目前富集分析工具多种各样包括在线工具与R包等，富集到的结果以及分析的库也各不相同，昨天我在生信技能树介绍了： 从基因名到GO注释一步到位 ，里面提到了其实有3个常见的网页工具也可以做到同样的分析，代码并没有任何神奇的地方， 结果的解读才是重中之重！ 但是网页工具我用起来毕竟还是有些丢脸，所以 安排学徒 比较了一下常用的3大在线分析工具: Enrichr、WebGestalt、gprofiler 与 R包clusterprofiler，有了这个笔记分享给大家。

PART 01

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/
# 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)
###有点久，因为要联网

PART 02

WebGestalt

WebGestalt同样是高引用率富集分析工具,现引用量超过 2,500(几版加起来)，支持3种算法进行富集：

1. Over-Representation Analysis (ORA)
2. Gene Set Enrichment Analysis (GSEA)
3. 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)
####也很慢，需要用外网

PART 03

gprofiler

• 由爱沙尼亚的塔尔图大学开发， 从2007到现在引用量800左右，个人觉得是网页配色最好看的一个 。目前网页总共有4种功能，功能注释，ID转换，同源基因物种间转换以及snp id转换

网页版

我们主要用的就是g:GOSt这一个功能

g:GOSt

1. performs functional enrichment analysis, also known as over-representation analysis (ORA) or gene set enrichment analysis, on input gene list.
2. 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.
3. 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.
4. 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包需要连接到外网，真的很慢~

PART 04

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个阈值才统一

PART 05

对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文件去做交集，看一看是不是本身的库文件本身就存在巨大差异。

PART 06

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

PART 07

差异原因

从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"),