3大在线分析工具:Enrichr、WebGestalt、gprofiler与R包clusterprofiler的比较
目前富集分析工具多种各样包括在线工具与R包等，富集到的结果以及分析的库也各不相同，昨天我在生信技能树介绍了：从基因名到GO注释一步到位，里面提到了其实有3个常见的网页工具也可以做到同样的分析，代码并没有任何神奇的地方，结果的解读才是重中之重！
但是网页工具我用起来毕竟还是有些丢脸，所以安排学徒比较了一下常用的3大在线分析工具:Enrichr、WebGestalt、gprofiler与R包clusterprofiler，有了这个笔记分享给大家。
Enrichr
网页版
enrichr由Ma'ayan Lab由2013年开发并维护，现引用量突破2500，而且很多都是CNS文章引用。现有来自164个库的332911个terms包括常规的GO KEGG以及Pfam，wikipathways等等。
除了进行人与老鼠的富集分析之外，还可以对其他物种例如飞虫等进行功能注释。
上传基因集
注意的是，Enrichr 只支持Entrez gene symbols 作为输入，所以其他格式需要进行ID转换，方法是Fisher exact test with the z-score of the deviation
上传数据后页面如下
分为表观修饰，通路，GO，疾病或药物中的表达以及杂项等，
用常见的GO举个例子，可分为BP CC MF三类
将3个库的所有富集到的terms下载下来做后续比较
R包版
install.packages("enrichR")
library(enrichR)
dbs <- listEnrichrDbs() ###列出164个库
dbs[1:4,1:4]
# geneCoverage genesPerTerm libraryName
# 1 13362 275 Genome_Browser_PWMs
# 2 27884 1284 TRANSFAC_and_JASPAR_PWMs
# 3 6002 77 Transcription_Factor_PPIs
# 4 47172 1370 ChEA_2013
# link
# 1 http://hgdownload.cse.ucsc.edu/goldenPath/hg18/database/
# 2 http://jaspar.genereg.net/html/DOWNLOAD/
# 3
# 4 http://amp.pharm.mssm.edu/lib/cheadownload.jsp
###从中选择你要富集的库
dbs$libraryName ###查看库名
dbs <- c("GO_Molecular_Function_2018", "GO_Cellular_Component_2018", "GO_Biological_Process_2018")###这里我选择GO库的3个process
library(clusterProfiler)
library(org.Hs.eg.db)
keytypes(org.Hs.eg.db)
symbol=read.csv("igno.txt",header = F)
#BiocManager::install("org.Hs.eg.db")
###id转换
df <- bitr(unique(symbol$V1), fromType = "ENSEMBL",
toType = c("SYMBOL","ENTREZID"),
OrgDb = org.Hs.eg.db)
enrichr<- enrichr(symbol, dbs)
###有点久，因为要联网
WebGestalt
WebGestalt同样是高引用率富集分析工具,现引用量超过 2,500(几版加起来)，支持3种算法进行富集：
Over-Representation Analysis (ORA)
Gene Set Enrichment Analysis (GSEA)
Network Topology-based Analysis (NTA)
截止2019年1月，现有321,251terms, 以及新添加了蛋白库 CORUM 与WikiPathway中肿瘤相关通路，重要的是有已经去除了GO的redundant terms
什么是redundant terms？
GO分为多级目录，也就是父目录下有很多分支子目录，nonredundant GO就是已经去除掉一 二级父目录了，富集结果更一目了然
网页版
运行后会得到一个汇总表：算法是BH，阈值为FDR=0.05，421个ID中有358个IDs被注释。总共有61506个 entrez gene IDs，只有16005个IDs用作这一次的功能聚类
第二个图告诉你，你的gene list 能比对到每个GO process中的三级目录的个数
最后就是常规的可视化，分为表格，条形图，火山图以及可以导入到cytoscape的DAG
最后通过右上角下载数据到本地
R包版
install.packages("WebGestaltR")
install.packages("gdtools")
library(WebGestaltR)
####ORA
head(listGeneSet())###列出所有的库的前6个
# name
# 1 geneontology_Biological_Process
# 2 geneontology_Biological_Process_noRedundant
# 3 geneontology_Cellular_Component
# 4 geneontology_Cellular_Component_noRedundant
# 5 geneontology_Molecular_Function
# 6 geneontology_Molecular_Function_noRedundant
enrichResult <- WebGestaltR(enrichMethod="ORA", organism="hsapiens", enrichDatabase=c("geneontology_Biological_Process_noRedundant","geneontology_Cellular_Component_noRedundant","geneontology_Molecular_Function_noRedundant")####GO 3个process
, interestGeneFile=df$SYMBOL,interestGeneType="genesymbol", isOutput=TRUE,
outputDirectory="./", projectName=NULL)
####也很慢，需要用外网
gprofiler
由爱沙尼亚的塔尔图大学开发，从2007到现在引用量800左右，个人觉得是网页配色最好看的一个。目前网页总共有4种功能，功能注释，ID转换，同源基因物种间转换以及snp id转换
网页版
我们主要用的就是g:GOSt这一个功能
g:GOSt
performs functional enrichment analysis, also known as over-representation analysis (ORA) or gene set enrichment analysis, on input gene list.
It maps genes to known functional information sources and detects statistically significantly enriched terms. We regularly retrieve data from Ensembl database and fungi, plants or metazoa specific versions of Ensembl Genomes, and parasite specific data from WormBase ParaSite.
In addition to Gene Ontology, we include pathways from KEGG Reactome andWikiPathways; miRNA targets from miRTarBase and regulatory motif matches from TRANSFAC; tissue specificity from Human Protein Atlas; protein complexes from CORUM and human disease phenotypes from Human Phenotype Ontology.
g:GOSt supports close to 500 organisms and accepts hundreds of identifier types.
可以看到主流的库基本囊括了。
run query后可以得到具体结果
R包版
install.packages("gprofiler2")
library(gprofiler2)
gostres <- gost(query = df$SYMBOL,
organism = "hsapiens", ordered_query = FALSE,
multi_query = FALSE, significant = TRUE, exclude_iea = FALSE,
measure_underrepresentation = FALSE, evcodes = FALSE,
user_threshold = 0.05, correction_method = "fdr",
domain_scope = "annotated", custom_bg = NULL,
numeric_ns = "", sources = GO, as_short_link = FALSE)
####这个也需要连接到外网
上面3个工具R包与网页版功能是一致的，但是在国内建议网页版，因为3个R包需要连接到外网，真的很慢~
cluster profiler
最后就是Y叔开发的R包cluster profiler，至今被引用率2518次，也可以用来做如GO、KEGG、DO（Disease Ontology analysis）、Reactome pathway analysis以及GSEA富集分析等，就是昨天jimmy老师介绍的：从基因名到GO注释一步到位，大家加油学习哦！
library(clusterProfiler)
library(org.Hs.eg.db)
keytypes(org.Hs.eg.db)
symbol=read.csv("igno.txt",header = F)
#BiocManager::install("org.Hs.eg.db")
df <- bitr(unique(symbol$V1), fromType = "ENSEMBL", ###输入只能说entrez id，所以需要iD转换
toType = c("SYMBOL","ENTREZID"),
OrgDb = org.Hs.eg.db)
go <- enrichGO(df$ENTREZID, OrgDb = "org.Hs.eg.db", ont="all",readable = T,qvalueCutoff = 0.05) ###默认阈值是pvalue=0,05，方法是"BH"，物种是人,这里需要将FDR设为0.05，因为这样4个阈值才统一
对4种工具筛选结果查看
###结果查看
###cluster profiler
cgo=go@result
ccc=cgo[cgo$ONTOLOGY=="CC",]
cbp=cgo[cgo$ONTOLOGY=="BP",]
cmf=cgo[cgo$ONTOLOGY=="MF",]
dim(ccc);dim(cbp);dim(cmf) ###分别为19 19 14个terms
# [1] 19 10
# [1] 19 10
# [1] 14 10
###enrich r
###enrichr
rbp=read.csv("mrna+lncrna/salmon/GO_Biological_Process_2018_table.txt",sep = "\t")
rcc=read.csv("GO_Cellular_Component_2018_table.txt",sep = "\t")
rmf=read.csv("GO_Molecular_Function_2018_table.txt",sep = "\t")
rbp=rbp[rbp$Adjusted.P.value<0.05,]
rcc=rcc[rcc$Adjusted.P.value<0.05,]
rmf=rmf[rmf$Adjusted.P.value<0.05,]
dim(rcc);dim(rbp);dim(rmf) ##可以看到如果根据fdr来筛选，cc与bp并没有显著terms，可能原因是其用的不是FDR or BH算法，而是fisher exact test
# [1] 0 9
# [1] 0 9
# [1] 4 9
###gprofiler
g=read.csv("gProfiler_hsapiens.csv")
gcc=g[g$source=="GO:CC",]
gbp=g[g$source=="GO:BP",]
gmf=g[g$source=="GO:MF",]
dim(gcc);dim(gbp);dim(gmf)
# [1] 20 10
# [1] 54 10
# [1] 8 10
### WebGestalt
####WebGestalt
wcc=read.csv("goslim_summary_wg_result1586186048_cc.txt",sep = "\t")
wbp=read.csv("goslim_summary_wg_result1586186048_bp.txt",sep = "\t")
wmf=read.csv("goslim_summary_wg_result1586173032_mf.txt",sep = "\t")
dim(wcc);dim(wbp);dim(wmf)
# [1] 21 3
# [1] 12 3
# [1] 17 3
####交集
library(venn)
cc=venn(list(ccc$ID,rcc$Term,gcc$term_id,wcc$V1),snames = c("cluster profiler","enrichr","gprofiler","WebGestalt"),zcolor = "style",sncs = 1,ellipse = T,box = F)
BP=venn(list(cbp$ID,rbp$Term,gbp$term_id,wbp$V1),snames = c("cluster profiler","enrichr","gprofiler","WebGestalt"),zcolor = "style",sncs = 1,ellipse = T,box = F)
mf=venn(list(cmf$ID,rmf$Term,gmf$term_id,wmf$V1),snames = c("cluster profiler","enrichr","gprofiler","WebGestalt"),zcolor = "style",sncs = 1,ellipse = T,box = F)
CC
BP
MF
感觉gprofiler与cluster profiler的结果比较相似，与其他2个工具分析结果相距甚远。（但是，总体来说，这些工具的一致性都好弱！！！）
upsetR
clcc=ifelse(ccu%in%ccc$ID,1,0)
rlcc=ifelse(ccu%in%rcc$Term,1,0)
glcc=ifelse(ccu%in%gcc$term_id,1,0)
wlcc=ifelse(ccu%in%wcc$V1,1,0)
cccc=data.frame(clcc,rlcc,glcc,wlcc)
rownames(cccc)=ccu
upset(cccc,nsets = 4)
CC
BP
MF
同样gprofiler与cluster profiler的结果比较相似，与其他2个工具分析结果相距甚远。（这里仅仅是把韦恩图替换成了upsetR的展现方式而已）
如果enrichR用p值=0.05筛选,一样是没有交集
CC
BP
MF
以上4种富集分析软件，gprofiler与cluster profiler结果比较相近，其他2种工具 enrichr可能是算法不一样，但是webgestalt算法是BH，FDR同样是0.05，不知道为什么结果相差甚远。接下来提取4个工具的GO库的gmt文件去做交集，看一看是不是本身的库文件本身就存在巨大差异。
3种原始库文件比较
因为webgestalt指向GO官网，没有找到2019.01.14的GO版本,暂时放到一边
myfun=function(x){
unlist(str_extract_all(x$V1,"GO:\\d+"))
}
###enrichr
rbp=read.csv("GO_Biological_Process_2018.txt",sep = "\n",header = F)
rbp=myfun(rbp)
rcc=read.csv("gocc.csv",sep = "\n",header = F)
rcc=myfun(rcc)
rmf=read.csv("GO_Molecular_Function_2018.txt",sep = "\n",header = F)
rmf=myfun(rmf)
length(rcc);length(rbp);length(rmf)
# [1] 446
# [1] 5103
# [1] 1151
###gprofiler
library(tidyr)
library(stringr)
gcc=read.csv("hsapiens.GO_CC.name.gmt",sep = "\n",header = F)
gcc=myfun(gcc)
gbp=read.csv("hsapiens.GO_BP.name.gmt",sep = "\n",header = F)
gbp=myfun(gbp)
gmf=read.csv("hsapiens.GO_MF.name.gmt",sep = "\n",header = F)
gmf=myfun(gmf)
length(gcc);length(gbp);length(gmf)
# [1] 2005
# [1] 16262
# [1] 4704
###cluster profiler
library(clusterProfiler)
cmf <- enrichGO(df$ENTREZID, OrgDb = "org.Hs.eg.db", ont="MF",readable = T,qvalueCutoff = 1,pvalueCutoff = 1)
cmf=names(cmf@geneSets)
ccc<- enrichGO(df$ENTREZID, OrgDb = "org.Hs.eg.db", ont="cc",readable = T,qvalueCutoff = 1,pvalueCutoff = 1)
ccc=names(ccc@geneSets)
cbp=enrichGO(df$ENTREZID, OrgDb = "org.Hs.eg.db", ont="BP",readable = T,qvalueCutoff = 1,pvalueCutoff = 1)
cbp=names(cbp@geneSets)
length(ccc);length(cbp);length(cmf)
# [1] 692
# [1] 4937
# [1] 872
CC
BP
MF
差异原因
从3个工具的数据库来看，gprofiler的 GOterms 是最多的，原因是什么呢？
1.不同工具的GO数据库更新时间不同
那最多的BP中一个为例子
bp=Reduce(setdiff,list(gbp,cbp,rbp))
bp[1]
#[1] "GO:0000019"
#获取该通路上所有基因
# ENSG00000020922
# ENSG00000076242
# ENSG00000104884
# ENSG00000113522
# ENSG00000132604
# ENSG00000180532
###放到cluster profiler去找
term=c("ENSG00000020922",
"ENSG00000076242",
"ENSG00000104884",
"ENSG00000113522",
"ENSG00000132604",
"ENSG00000180532")
###id转换
df_1 <- bitr(unique(term), fromType = "ENSEMBL",
toType = c("SYMBOL","ENTREZID"),
OrgDb = org.Hs.eg.db)
cbp_1=enrichGO(df_1$ENTREZID, OrgDb = "org.Hs.eg.db", ont="BP",readable = T,qvalueCutoff = 1,pvalueCutoff = 1)
cbp_1=cbp_1@result
cluster profiler 可以看到最多可以mapping到5/6，并没有完全mapping上
enrich r也是同样的结果，不能完全mapping上，更惨的是最多4/4
而webgestalt不仅完全mapping上，还多2个基因在这个term上
从各个数据库给出的GO库的更新时间来判断，enrichr的GO库是2018年的，cluster profiler与gprofiler没有给出，但是研发gprofiler的实验室在 2019年5月发表文章说更新了GO数据库，所以推定数据库是2019.05之前，而webgestalt官网给出的GO更新时间是2020.01，所以可以对4个工具的GO数据库做一个大致排序，
从新到旧：webgestalt > gprofiler >cluster profiler >enrichr
2. 冗余数据库
gprofiler 库并没有去除冗余的GO数据库，打开gprofiler的网页可以看到，例如比对同样的GO:0000019进行富集分析，能mapping到一个父目录GO0006259，打开之后会看到这个terms有上943个基因，这样反复多次父子级目录mapping，就会造成很多冗余的GO分集，这也是造成数量多的原因。而其他三个工具都是non-redudnant的GO数据，省去了很多时间，特别是webgestalt，既有去冗余也有所有GO。
综上所述，4个工具中GO分析最好用的还是webgestalt：
因为更新时间为最新（2020.01）
有冗余与非冗余2种GO库选择
国内最好还是使用网页版
如果想用R包就用Y叔的cluster profiler，其他3个工具R版真的慢，其他富集分析则需要具体情况具体考虑了。