传统的富集分析的步骤是：
- 计算出差异表达基因列表
- 运行GO
- 运行KEGG
当然，也可以运行其他富集分析，比如DO，Wikipathways或者Reactome等等，然后再可视化结果，这里所有的操作都不能一次性运行。
在R里面有一个list的数据格式，对于相同的需求，其实是可以后台一步运行的，只是之前没有相应的R包和代码。
最近Y叔开发了一个gson 格式的包，制定了一种新的背景基因集.gson （类似于.gmt），并且嵌入到了clusterProfiler包里，通过gson_GO()、gson_KEGG()、gson_WP()函数可以在线爬取相应的GSON背景基因集对象，然后可以进行后续的GSEA或者普通的富集分析。
新开发的gson包中有一个gsonList()函数，可以进行多组富集分析，以 list的形式进行组合分析，然后通过enrichplot包中的autofacet()函数就可以将多个list的结果可视化。这个教程在一次搞定所有的富集分析这个公众号里面。结果如下：
image.png
但是，这个教程刚出的时候，我就想试试GSEA能不能行，结果发现当时的enrichplot包中的autofacet()函数并不支持GSEA结果，所以我当天就在Y叔公众号和Github进行了提问，经过一个月的等待，终于在前几天收到了解决的邮件。还给我发了一个测试pdf结果，确实可以分面多个GSEA的结果了。
image.png
不过，我还是发现一个问题，那就是GSEA本身就包含一个分面，根据NES值分为激活和抑制（这个功能可以通过ggplot2的facet_grid(~.sign)实现），那么想实现多个富集分析列表的激活和抑制却又卡住了。
我发现：
使用了
autofacet()，就不能使用facet_grid(~.sign)；使用了
facet_grid(~.sign)，就不能使用autofacet()。
当然这个Bug可以通过图片的拼接功能实现，但双分面才是我要的结果，这个bug我也已经反映了，不知道什么时候能够解决。
说完了前言，现在就开始正式演示，使用gson包可以在线获取或读取本地的.gson文件，我们可以一次性分析BP和KEGG，也可以一次性运行KEGG和Reactome这两种信号通路，或者一次性运行KEGG和WikiPathways两种信号通路。
获得GSON对象
GSON库
Y叔的工作网页中制作了一些gson格式的基因集库，我们可以直接下载好并且使用，网址是https://yulab-smu.top/gson-files/
Y叔的库
然而，这个基因集库中还没有收录WikiPathways的结果，并且还需要访问Github才能下载。
所以我重新制作了一下这个库，顺便更新了一下时间，将结果导入到码云Gitee上了，地址是
http://swcyo.gitee.io/gson-file/
我修改的库
你也可以这样点击直接下面的链接直接下载保存到本地文件夹：
在线GSON
- 我们使用最常用的BP生物学功能和KEGG信号通路进行演示
library(clusterProfiler)
BP <- gson_GO(OrgDb = org.Hs.eg.db, keytype = 'ENTREZID', ont = "BP")
KEGG <- gson_KEGG(species = "hsa", KEGG_Type="KEGG", keyType="kegg")
如果你想获得目前可获得的所有结果，可以通过clusterProfiler的gson_GO()，gson_KEGG()，gson_WP()和ReactomePA的gson_Reactome()去爬相应背景基因集的最新数据，然后保存到本地使用。代码如下：
# GO
library(clusterProfiler)
library(org.Hs.eg.db)
library(gson)
gson_BP_human <- gson_GO(OrgDb = org.Hs.eg.db, keytype = 'ENTREZID', ont = "BP")
gson_MF_human <- gson_GO(OrgDb = org.Hs.eg.db, keytype = 'ENTREZID', ont = "MF")
gson_CC_human <- gson_GO(OrgDb = org.Hs.eg.db, keytype = 'ENTREZID', ont = "CC")
gson_ALL_human <- gson_GO(OrgDb = org.Hs.eg.db, keytype = 'ENTREZID', ont = "ALL")
write.gson(gson_BP_human, file = "GO_BP_human.gson")
write.gson(gson_MF_human, file = "GO_MF_human.gson")
write.gson(gson_CC_human, file = "GO_CC_human.gson")
write.gson(gson_ALL_human, file = "GO_ALL_human.gson")
# KEGG
KEGG_human <- gson_KEGG(species = "hsa", KEGG_Type="KEGG", keyType="kegg")
MKEGG_human <- gson_KEGG(species = "hsa", KEGG_Type="MKEGG", keyType="kegg")
write.gson(KEGG_human, file = "KEGG_human.gson")
write.gson(MKEGG_human, file = "MKEGG_human.gson")
# WikiPathways
WikiPathways_human <- gson_WP("Homo sapiens")
write.gson(WikiPathways_human, file = "WikiPathways_human.gson")
# Reactome
library(ReactomePA)
Reactome_human <- gson_Reactome(organism = "human")
write.gson(Reactome_human, file = "Reactome_human.gson")
加载本地GSON对象
如果通过上述代码把gson文件下载到了本地文件夹，我们可以通过gson的read.gson()函数进行加载。
library(gson)
KEGG<-read.gson("KEGG_human.gson")
WP<-read.gson("WikiPathways_human.gson")
运行GSEA
我们使用DOSE包中geneList数据列表进行GSEA演示
library(clusterProfiler)
library(gson)
library(enrichplot)
# 构建gson列表
gsonlist <- gsonList(BP = BP, KEGG=KEGG)
data(geneList,package = 'DOSE')
# 一行运行多组GSEA
GSEA <- GSEA(geneList,
gson = gsonlist,
eps=0, # 设置这个可以获得更好的p值
pvalueCutoff = 0.9 # p值调大一点
)
可视化
使用dotplot()对多组GSEA结果绘制点图，但是这个结果把BP和KEGG的两种结果全部都包含了，不能体现出具体的结果，见Figure 1所示。
dotplot(GSEA)
Figure 1: 多组GSEA结果同时显示
自动分面
使用autofacet()可以很好的显示分面，默认是按行分面，见Figure 2所示。
dotplot(GSEA,showCategory = 8)+autofacet()
Figure 2: 按行对两种结果进行分面显示
我们也可以按列分面
dotplot(GSEA)+autofacet(by='col')
Figure 3: 按列对两种结果进行分面显示
Bug
我们知道GSEA是支持以NES为界，分为激活和抑制两个结果的，传统方法是+facet_grid(~.sign)，但是autofacet()和+facet_grid(~.sign)目前只能显示一种结果。
可能是数据有误，只跑出列激活的结果，见Figure 4所示。
library(ggplot2)
dotplot(GSEA)+facet_grid(~.sign)
Figure 4: 显示激活抑制通路
解决方法
目前可以通过拼图的方法对两个结果进行处理，见Figure 5所示。
library(patchwork)
p1<-dotplot(GSEA$`Gene Ontology;BP`)+facet_grid(~.sign)
p2<-dotplot(GSEA$KEGG)+facet_grid(~.sign)
p1/p2
Figure 5: 结果拼图
通过gsonlist，加上Y叔开发的包包，以后做富集分析就越来越觉得，越来越快了。
网友评论