傳統(tǒng)的富集分析的步驟是：

• 計算出差異表達基因列表
• 運行GO
• 運行KEGG
  當然槽畔，也可以運行其他富集分析哲银，比如DO，Wikipathways或者Reactome等等初肉，然后再可視化結(jié)果酷鸦，這里所有的操作都不能一次性運行。
  在R里面有一個list的數(shù)據(jù)格式，對于相同的需求臼隔，其實是可以后臺一步運行的嘹裂，只是之前沒有相應(yīng)的R包和代碼。

最近Y叔開發(fā)了一個gson 格式的包摔握，制定了一種新的背景基因集.gson （類似于.gmt）寄狼，并且嵌入到了clusterProfiler包里，通過gson_GO()氨淌、gson_KEGG()泊愧、gson_WP()函數(shù)可以在線爬取相應(yīng)的GSON背景基因集對象，然后可以進行后續(xù)的GSEA或者普通的富集分析盛正。

新開發(fā)的gson包中有一個gsonList()函數(shù)删咱，可以進行多組富集分析，以 list的形式進行組合分析蛮艰，然后通過enrichplot包中的autofacet()函數(shù)就可以將多個list的結(jié)果可視化腋腮。這個教程在一次搞定所有的富集分析這個公眾號里面。結(jié)果如下：

image.png

但是壤蚜，這個教程剛出的時候，我就想試試GSEA能不能行徊哑，結(jié)果發(fā)現(xiàn)當時的enrichplot包中的autofacet()函數(shù)并不支持GSEA結(jié)果袜刷，所以我當天就在Y叔公眾號和Github進行了提問，經(jīng)過一個月的等待莺丑，終于在前幾天收到了解決的郵件著蟹。還給我發(fā)了一個測試pdf結(jié)果，確實可以分面多個GSEA的結(jié)果了梢莽。

image.png

不過卵蛉，我還是發(fā)現(xiàn)一個問題潮峦，那就是GSEA本身就包含一個分面，根據(jù)NES值分為激活和抑制（這個功能可以通過ggplot2的facet_grid(~.sign)實現(xiàn)），那么想實現(xiàn)多個富集分析列表的激活和抑制卻又卡住了滤愕。

我發(fā)現(xiàn)：

使用了autofacet()，就不能使用facet_grid(~.sign)梧却；

使用了facet_grid(~.sign)挡毅，就不能使用autofacet()。

當然這個Bug可以通過圖片的拼接功能實現(xiàn)仑扑，但雙分面才是我要的結(jié)果览爵，這個bug我也已經(jīng)反映了，不知道什么時候能夠解決镇饮。

說完了前言蜓竹，現(xiàn)在就開始正式演示，使用gson包可以在線獲取或讀取本地的.gson文件，我們可以一次性分析BP和KEGG俱济，也可以一次性運行KEGG和Reactome這兩種信號通路司蔬，或者一次性運行KEGG和WikiPathways兩種信號通路。

獲得GSON對象

GSON庫

Y叔的工作網(wǎng)頁中制作了一些gson格式的基因集庫姨蝴，我們可以直接下載好并且使用俊啼，網(wǎng)址是https://yulab-smu.top/gson-files/

Y叔的庫

然而，這個基因集庫中還沒有收錄WikiPathways的結(jié)果左医，并且還需要訪問Github才能下載授帕。

所以我重新制作了一下這個庫，順便更新了一下時間浮梢，將結(jié)果導(dǎo)入到碼云Gitee上了跛十，地址是

http://swcyo.gitee.io/gson-file/

我修改的庫

你也可以這樣點擊直接下面的鏈接直接下載保存到本地文件夾：

Gene Ontology;ALL

Gene Ontology;BP

Gene Ontology;CC

Gene Ontology;MF

KEGG

MKEGG

reactome pathway

WikiPathways

在線GSON

• 我們使用最常用的BP生物學(xué)功能和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") 

如果你想獲得目前可獲得的所有結(jié)果，可以通過clusterProfiler的gson_GO()秕硝，gson_KEGG()芥映，gson_WP()和ReactomePA的gson_Reactome()去爬相應(yīng)背景基因集的最新數(shù)據(jù)，然后保存到本地使用远豺。代碼如下：

# 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()函數(shù)進行加載。

library(gson)
KEGG<-read.gson("KEGG_human.gson")
WP<-read.gson("WikiPathways_human.gson")

運行GSEA

我們使用DOSE包中g(shù)eneList數(shù)據(jù)列表進行GSEA演示

library(clusterProfiler)
library(gson)
library(enrichplot)
# 構(gòu)建gson列表
gsonlist <- gsonList(BP = BP, KEGG=KEGG)
data(geneList,package = 'DOSE')
# 一行運行多組GSEA
GSEA <- GSEA(geneList, 
                 gson = gsonlist,
                 eps=0, # 設(shè)置這個可以獲得更好的p值
                 pvalueCutoff = 0.9 # p值調(diào)大一點
             )

可視化

使用dotplot()對多組GSEA結(jié)果繪制點圖躯护，但是這個結(jié)果把BP和KEGG的兩種結(jié)果全部都包含了惊来，不能體現(xiàn)出具體的結(jié)果，見Figure 1所示棺滞。

dotplot(GSEA)

Figure 1: 多組GSEA結(jié)果同時顯示

自動分面

使用autofacet()可以很好的顯示分面裁蚁，默認是按行分面，見Figure 2所示继准。

dotplot(GSEA,showCategory = 8)+autofacet()

Figure 2: 按行對兩種結(jié)果進行分面顯示

我們也可以按列分面

dotplot(GSEA)+autofacet(by='col')

Figure 3: 按列對兩種結(jié)果進行分面顯示

Bug

我們知道GSEA是支持以NES為界枉证，分為激活和抑制兩個結(jié)果的，傳統(tǒng)方法是+facet_grid(~.sign)移必，但是autofacet()和+facet_grid(~.sign)目前只能顯示一種結(jié)果室谚。

可能是數(shù)據(jù)有誤，只跑出列激活的結(jié)果避凝，見Figure 4所示舞萄。

library(ggplot2)
dotplot(GSEA)+facet_grid(~.sign)

Figure 4: 顯示激活抑制通路

解決方法

目前可以通過拼圖的方法對兩個結(jié)果進行處理，見Figure 5所示管削。

library(patchwork)
p1<-dotplot(GSEA$`Gene Ontology;BP`)+facet_grid(~.sign)
p2<-dotplot(GSEA$KEGG)+facet_grid(~.sign)
p1/p2

Figure 5: 結(jié)果拼圖

通過gsonlist倒脓，加上Y叔開發(fā)的包包，以后做富集分析就越來越覺得含思，越來越快了崎弃。