ixxmu
commented
11 months ago 任意物种的go和kegg数据库注释（视频号直播互动答疑） by 生信技能树
本文分享内容如下：
六大安装R包的方法 解决此次clusterProfiler包安装的方法 物种注释的三层境界 举个小绵羊注释的例子
一、六大安装R包的方法
1 从CRAN安装
# 1 directly download from CRAN
# mirror for cran：https://mirrors.ustc.edu.cn/CRAN/; https://mirrors.tuna.tsinghua.edu.cn/CRAN/
options(repos="https://mirrors.tuna.tsinghua.edu.cn/CRAN/")
install.packages("ggplot2", repos="https://mirrors.ustc.edu.cn/CRAN/")
install.packages("tidyverse")
library(ggplot2)
2 从Bioconductor安装
# 2 downlaod packages from bioconductor
# mirror for bioc: https://mirrors.ustc.edu.cn/bioc/; http://mirrors.ustc.edu.cn/bioc/
# options(BioC_mirror="http://mirrors.ustc.edu.cn/bioc/")
if (!requireNamespace("BiocManager")) install.packages("BiocManager")
BiocManager::available() # find packages
BiocManager::install() # update packages
BiocManager::install(c("ape", "ggtree")) # install packages
3 从Github安装
# 3 download packages from github
install.packages("devtools", repos="https://mirrors.ustc.edu.cn/CRAN/")
install.packages("githubinstall")
library(devtools)
library(githubinstall)
# There is an install_github function to install R packages hosted on GitHub in the **devtools** package. But it requests developer’s name.
install_github("DeveloperName/PackageName")
eg:
devtools::install_github("twitter/AnomalyDetection")
devtools::install_github('dviraran/SingleR')
# The **githubinstall** package provides a function githubinstall. It does not need developer’s name.
githubinstall("PackageName")
eg: githubinstall("AnomalyDetection")
4 从网页安装
# 4 download from URL
packageurl <- "http://cran.r-project.org/src/contrib/Archive/ggplot2/ggplot2_0.9.1.tar.gz"
packageurl <- "http://cran.r-project.org/src/contrib/Archive/gridExtra/gridExtra_0.9.1.tar.gz"
install.packages(packageurl, repos=NULL, type="source")
packageurl <- "http://www.bioconductor.org/packages/2.11/bioc/src/contrib/ggbio_1.6.6.tar.gz"
packageurl <- "http://cran.r-project.org/src/contrib/Archive/ggplot2/ggplot2_1.0.1.tar.gz"
install.packages(packageurl, repos=NULL, type="source")
5 用函数从本地安装
# 5 download from local
download.file("http://bioconductor.org/packages/release/bioc/src/contrib/BiocInstaller_1.20.1.tar.gz","BiocInstaller_1.20.1.tar.gz")
install.packages("BiocInstaller_1.20.1.tar.gz", repos = NULL, type="source")
## Install the package from local
devtools::install_local("C:/Users/whbio/Desktop/temp/SingleR-master.zip")
library(SingleR)
Notice：下载压缩包到本地，指定压缩包目录，然后用R函数安装
6 用R从本地安装
# 6 Download from R
# sudo R CMD INSTALL package.tar.gz
Notice：下载压缩包到本地，然后用R软件安装
二、解决此次clusterProfiler包安装的方法
Error如下
> library(clusterProfiler)
Error: package or namespace load failed for ‘clusterProfiler’:
'namespace:HPO.db'没有出口‘HPOGENE’这个对象
Try1：从Bioconductor下载HPO压缩到到本地，FAILED
在确认一系列依赖包顺利安装并下载压缩包后，从本地安装此包。
# dependenct packages: AnnotationDbi, AnnotationHub, BiocFileCache, DBI, methods, utils
dependency_packages <- c('AnnotationDbi', 'AnnotationHub', 'BiocFileCache', 'DBI', 'methods', 'utils')
inst = lapply(dependency_packages, library, character.only = TRUE)
# install from local
install.packages("D:/2023_softwares/Installation_softwares/R/R-4.3.0/library/HPO.db_0.99.2.tar.gz",
repos = NULL, type="source")
Try2: 下载包的旧版本到本地，再次安装，SUCCESS
包的下载地址为：https://github.com/huerqiang/HPO.db_MPO.db
# https://github.com/huerqiang/HPO.db_MPO.db
install.packages('D:/2023_softwares/Installation_softwares/R/R-4.3.0/library/HPO.db_MPO.db-main/HPO.db_0.99.0.tar.gz')
install.packages('D:/2023_softwares/Installation_softwares/R/R-4.3.0/library/HPO.db_MPO.db-main/MPO.db_0.99.0.tar.gz')
library(clusterProfiler)
三、物种注释的三层境界
模仿张敬信老师写循环的三层境界，注释也有三层境界。
第一层：经典模式生物
如果是常见人鼠物种，直接指定物种名称即可，so easy！
# 直接指定物种名
go_organism <- "org.Hs.eg.db" # org.Hs.eg.db / org.Mm.eg.db / org.Gg.eg.db
kegg_organism <- 'hsa' # hsa / mmu / gga
第二层：非模式生物，但OrgDb库可检索
网页查找OrgDb库
1、Bioconductor提供19种物种OrgDb包 （图1）
https://bioconductor.org/packages/release/BiocViews.html#___OrgDb
2、GOC（Gene Ontology Consortium）提供39种不同模型生物的GAF格式的注释信息 （图2）
http://current.geneontology.org/products/pages/downloads.html
3、Bioconductor的Annotationhub数据库提供了6126种物种 （图3）
https://annotationhub.bioconductor.org/species直接搜索自己所研究的物种的拉丁学名，即可找到对应的sqlite文件，下载到本地即可。如何使用sqlite格式文件呢？
4、KEGG物种注释库提供真核生物968种，细菌8073种，古菌421种，共9462种物种信息 （图4）
https://www.genome.jp/kegg/catalog/org_list.html
除此之外，还有下面的查找方法：
NCBI的基因ftp数据库：ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/
包含gene2go,gene2ensembl,gene2pubmed等，每天或每周更新，包含所有物种的信息。
图1
图2
图3
图4
R函数查找OrgDb库
1、GO注释库
# 1 用R包AnnotationHub检索
library(AnnotationHub)
ah = AnnotationHub()
# 2 非模式物种，以毛果杨为例
query(ah, "Populus_trichocarpa")
# title
# AH10434 | hom.Populus_trichocarpa.inp8.sqlite
# AH66282 | org.Populus_trichocarpa.eg.sqlite
# 3 下载对应的文件org.xx.eg.db.sqlite
populus <- ah[['AH66282']]
# > populus <- hub[['AH66282']]
# downloading 1 resources
# retrieving 1 resource
# 4 保存变量
saveDb(populus, file = "populus.orgdb")
# loadDb(file = "populus.orgdb")
# 5 查看变量
columns(populus)
length(keys(populus)) # 37273
head(keys(populus, keytype="REFSEQ"))
# 6 使用变量
# mydf <- bitr(x$SYMBOL, fromType="SYMBOL", toType="ENTREZID", populus)
# ====================================================================
2、KEGG注释库
# KEGG物种名检索
library(clusterProfiler)
search_kegg_organism('ece', by='kegg_code')
ecoli <- search_kegg_organism('Escherichia coli', by='scientific_name')
dim(ecoli)
Taxonomy 分类数据库是NCBI公共序列数据库中所有生物的策划分类和命名法，包含地球上大概10%的物种。可以看到不同的分类下的分布情况，总体包含有789855条物种信息。
目前可以直接检索解决的物种中GO注释有6184种（19+39+6126），KEGG注释有9462种。
如果你的物种着实罕见，欢迎分享物种名，一起找一找。
第三层：非模式生物，自建OrgDb库
如果物种实在稀缺，用R包AnnotationForge自建OrgDb库。
参见参考链接[3]和[4]，尤其是[4]。因为写得太好，所以不画蛇添足了。
四、举个小绵羊注释的例子
可爱的小绵羊名字为Ovis aries（sheep）
1、R函数查找OrgDb库
library(AnnotationHub)
hub = AnnotationHub::AnnotationHub()
一个大大的Error出现，ok，尝试网页搜索，直接下载。
在Annotationhub数据库中输入Ovis aries，下载org.Ovis_aries.eg.sqlite
或者shell下执行axel -n 20 https://bioconductorhubs.blob.core.windows.net/annotationhub/ncbi/uniprot/3.15/org.Ovis_aries.eg.sqlite。
如果网络顺利后，可以尝试R中如下安装
sheep <- AnnotationHub::AnnotationHub()[["AH101101"]]
saveDb(sheep, file = "sheep.orgdb")
loadDb(file = "sheep.orgdb")
KEGG注释物种输入Ovis aries，得到物种名为oas
2、基因ID转换
library(clusterProfiler)
library(AnnotationHub)
library(AnnotationDbi)
library(ggplot2)
library(ggpubr)
library(ggsci)
# Read genelist
setwd(r"(E:\Year_2023\Annual_small_work\Work_November\2023-11-07_test)")
load('BigData.Rdata')
head(target_gene)
## [1] "XP_011988802.3" "XP_042099641.1" "XP_042099639.1" "XP_042099640.1" "XP_011988809.3" "XP_011988797.3"
# Read OrgDb
sheep <- loadDb(file='./OrgDb/org.Ovis_aries.eg.sqlite')
# 指定物种名 Ovis aries (sheep)test
go_organism <- sheep # org.Hs.eg.db / org.Mm.eg.db / org.Gg.eg.db
kegg_organism <- 'oas' # hsa / mmu / gga
columns(go_organism)
## [1] "ACCNUM" "ALIAS" "CHR" "ENSEMBL" "ENTREZID" "EVIDENCE" "EVIDENCEALL" "GENENAME" "GID" "GO"
## [11] "GOALL" "ONTOLOGY" "ONTOLOGYALL" "PMID" "REFSEQ" "SYMBOL" "UNIGENE"
head(keys(go_organism, keytype="REFSEQ"))
## [1] "NM_001009189.1" "NM_001009191.1" "NM_001009192.2" "NM_001009194.1" "NM_001009195.1" "NM_001009196.1"
# 转基因ID
eg <- bitr(target_gene, fromType="REFSEQ", toType="ENTREZID", OrgDb=go_organism)
head(eg)
## REFSEQ ENTREZID
## 1 XP_011988802.3 101105639
## 2 XP_042099641.1 101105639
## 3 XP_042099639.1 101105639
## 4 XP_042099640.1 101105639
## 5 XP_011988809.3 101105639
## 6 XP_011988797.3 101105639
genelist <- eg$ENTREZID
# Notice: fromType类型
3、GO富集分析
ego <- enrichGO(genelist,
OrgDb = go_organism,
keyType = "GID",
ont = "ALL",
qvalueCutoff = 0.05,
pvalueCutoff = 0.05)
dim(ego@result)
## [1] 120 10
ego@result[1:2, ]
## ONTOLOGY ID Description GeneRatio BgRatio pvalue p.adjust qvalue
## GO:0006091 BP GO:0006091 generation of precursor metabolites and energy 13/51 120/4831 1.529984e-10 1.040389e-07 8.181391e-08
## GO:0005975 BP GO:0005975 carbohydrate metabolic process 12/51 149/4831 2.681084e-08 9.115685e-06 7.168372e-06
## geneID Count
## GO:0006091 101115777/101115434/443005/101113198/443086/101109163/101108778/443091/101120193/442998/101111501/100270719/100216440 13
## GO:0005975 101118831/101108630/101115777/443005/101113198/443086/101109163/443091/443089/442998/101111501/100216440 12
4、KEGG富集分析
kk_df <- enrichKEGG(genelist,
organism = kegg_organism,
keyType = 'kegg', # one of "kegg", 'ncbi-geneid', 'ncbi-proteinid' and 'uniprot'
pvalueCutoff = 0.05,pAdjustMethod = 'BH',
minGSSize = 10,maxGSSize = 500,
qvalueCutoff = 0.2,use_internal_data = FALSE)
dim(kk_df@result)
## [1] 168 11
kk_df@result[1:2, ]
## category subcategory ID Description GeneRatio BgRatio pvalue p.adjust qvalue
## oas05012 Human Diseases Neurodegenerative disease oas05012 Parkinson disease 25/118 332/9556 1.887259e-13 3.170594e-11 2.026320e-11
## oas00010 Metabolism Carbohydrate metabolism oas00010 Glycolysis / Gluconeogenesis 11/118 72/9556 1.018531e-09 8.555661e-08 5.467904e-08
## geneID
## oas05012 101118037/101115434/101116188/101105620/101118985/100270717/105614945/101110188/101115230/443016/443296/100135694/443297/101122127/101108778/101112274/101113030/121820456/443230/105614854/101116338/101115216/101121036/101116512/100270719
## oas00010 443005/101113460/101113198/443086/114108635/114110581/101108127/443089/101108565/443489/100216440
## Count
## oas05012 25
## oas00010 11
5、保存运行结果
res.go = DOSE::setReadable(ego, OrgDb = go_organism, keyType = 'ENTREZID') # 按需替换ENTREZID为SYMBOL
res.kegg = DOSE::setReadable(kk_df, OrgDb = go_organism, keyType = 'ENTREZID') # 按需替换ENTREZID为SYMBOL
writexl::write_xlsx(list('GO注释结果' = res.go@result, 'KEGG注释结果' = res.kegg@result), 'GO-KEGG注释结果.xlsx')
6、简单可视化
dotplot(ego)
dotplot(kk_df)ego@result %>%
dplyr::group_by(ONTOLOGY) %>%
# 为便于图形展示，提取部分子集
dplyr::slice_head(n = 5) %>%
ggplot(aes(x = reorder(Description, Count), y = Count)) +
geom_bar(stat = "identity", aes(fill = ONTOLOGY), alpha = 1) +
facet_grid(ONTOLOGY ~ . , scale="free", space="free_x", margins = F) +
labs(y = "Number of Genes", x = "Term")+
coord_flip() + theme_classic(base_size = 16) +
scale_y_continuous(expand = c(0,0)) +
scale_fill_simpsons() +
geom_text(aes(label=Count), hjust=1.2, color="white")
查找的资料越多，越发现自己的无知，推荐一个哈佛大学转录组的课程吧，见参考链接[6]。
只是简单地写个总结，没曾想被狗粮喂饱(狗粮见参考链接[4]的主页)。呼朋引伴后，被投喂了更多狗粮。Ok，为人类高质量爱情喝彩！
Reference
[1] 生信分析：常用文件数据库地址（GO,uniref）https://www.jianshu.com/p/7c108b851fa3
[2] 用AnnotationHub获取非模式物种注释信息https://www.bioinfo-scrounger.com/archives/512/
[3] 富集分析：（四）clusterProfiler：不同物种的GO+KEGG富集分析https://zhuanlan.zhihu.com/p/536082841
[4] 208-详细回顾非模式物种注释构建过程https://www.jieandze1314.com/post/cnposts/208/
[5] Biomedical Knowledge Mining using GOSemSim and clusterProfilerhttps://yulab-smu.top/biomedical-knowledge-mining-book/index.html
[6] Differential gene expression workshop using Salmon countshttps://hbctraining.github.io/DGE_workshop_salmon_online/
视频号直播互动答疑
本文代码是可以复制粘贴后在自己的项目中使用的，如果大家测试过程遇到了困难，欢迎明晚七点参与直播互动哈
文末友情宣传
强烈建议你推荐给身边的博士后以及年轻生物学PI，多一点数据认知，让他们的科研上一个台阶：
生物信息学马拉松授课（买一得五） ，你的生物信息学入门课 时隔5年，我们的生信技能树VIP学徒继续招生啦 144线程640Gb内存服务器共享一年仍然是仅需800 千呼万唤始出来的独享生物信息学云服务器 生信技能树知识整理实习生又又又开放申请啦（不招了，谢谢）
https://mp.weixin.qq.com/s/av5TWUyqWMpM8xpbFp2NlQ