1 什么是GSEA

基因集富集分析(Gene Set Enrichment Analysis, GSEA)是是一種計(jì)算方法，用于確定事先定義的一組基因是否在不同的樣品中差異表達(dá)。

GSEA首先將基因在樣品中的差異倍數(shù)值（logFC）由大到小排序，然后判斷來(lái)自功能注釋等預(yù)定義的基因集或自定義的基因集中的基因是富集在這個(gè)排序列表的頂部還是底部，如果在富集頂部，則該基因集是上調(diào)趨勢(shì)，反之，如果富集在底部，則是下調(diào)趨勢(shì)。

GSEA官網(wǎng)提供了詳細(xì)說(shuō)明，以及對(duì)應(yīng)軟件的下載地址。

2 GSEA特點(diǎn)

傳統(tǒng)的KEGG（通路富集分析）和GO（功能富集）分析時(shí)，針對(duì)總體的差異基因，不區(qū)分哪些差異基因是上調(diào)還是下調(diào)。所以會(huì)出現(xiàn)同一通路下富集到的的既有上調(diào)差異基因，也有下調(diào)差異基因，無(wú)法判斷這條通路表現(xiàn)形式究竟是怎樣。

而GSEA考慮了基因的表達(dá)水平，不需要明確指定差異基因閾值，檢驗(yàn)的是基因集而非單個(gè)基因的表達(dá)變化，算法會(huì)根據(jù)實(shí)際數(shù)據(jù)的整體趨勢(shì)進(jìn)行分析，以判斷這條通路的表達(dá)情況，激活或者抑制。

3 GSEA結(jié)果解讀

• 第1部分：Enrichment Score的折線圖

橫軸為排序后的基因，縱軸為對(duì)應(yīng)的ES, 在折線圖中出現(xiàn)的峰值就是這個(gè)基因集的富集分?jǐn)?shù)（Enrichment Score，ES）。ES是從排序后的表達(dá)基因集的第一個(gè)基因開(kāi)始，如果排序后表達(dá)基因列表中的基因出現(xiàn)在功能基因數(shù)據(jù)集中則加分，反之則減分。正值說(shuō)明在頂部富集，峰值左邊的基因?yàn)楹诵幕?，?fù)值則相反。

• 第2部分：基因位置圖

黑線代表排序后表達(dá)基因列表中的基因位于當(dāng)前分析的功能注釋基因集的位置，紅藍(lán)相間的熱圖是表達(dá)豐度排列，紅色越深的表示該位置的基因logFC越大 ，藍(lán)色越深表示logFC越小。如果研究的功能注釋基因集的成員顯著聚集在表達(dá)數(shù)據(jù)集的頂部或底部，則說(shuō)明功能基因數(shù)據(jù)集中的基因在數(shù)據(jù)集中高表達(dá)或低表達(dá)，若隨機(jī)分配，則說(shuō)明表達(dá)數(shù)據(jù)集與該通路無(wú)關(guān)。

• 第3部分：每個(gè)基因?qū)?yīng)的信噪比（Signal2noise）

以灰色面積圖展示?；疑幱暗拿娣e比，可以從整體上反映組間的Signal2noise的大小。

一般認(rèn)為校正后的富集分?jǐn)?shù)（NES）|NES|>1,p<0.05, qvalue(即FDR)<0.25的結(jié)果有意義。

4 GSEA實(shí)戰(zhàn)

#加載數(shù)據(jù) 數(shù)據(jù)鏈接：https://wwu.lanzouf.com/idmJB07bcefa
load(file = 'GSEA_test.Rdata')
colnames(deg)
head(deg)
#得到差異基因列表后取出 ，p值和logFC
nrDEG=deg[,c(2,1)]
colnames(nrDEG)=c('log2FoldChange','pvalue')
head(nrDEG)

       log2FoldChange        pvalue
LYZ         -1.812078 1.228136e-144
FCGR3A       2.617707 3.977801e-119
S100A9      -2.286734  2.481257e-95
S100A8      -2.610696  3.626489e-92
IFITM2       1.445772  7.942512e-87
LGALS2      -2.049431  1.275856e-75
#注：最后需要為文件如上：一列為基因名，一列為FC值，一列為p值

#加載Y叔的R包，把SYMBOL轉(zhuǎn)換為ENTREZID，后面可以直接做 KEGG 和 GO
library(org.Hs.eg.db)
library(clusterProfiler)
gene <- bitr(rownames(nrDEG),     #轉(zhuǎn)換的列是nrDEG的列名
             fromType = "SYMBOL",     #需要轉(zhuǎn)換ID類型
             toType =  "ENTREZID",    #轉(zhuǎn)換成的ID類型
             OrgDb = org.Hs.eg.db)    #對(duì)應(yīng)的物種，小鼠的是org.Mm.eg.db
#讓基因名、ENTREZID、logFC對(duì)應(yīng)起來(lái)
gene$logFC <- nrDEG$log2FoldChange[match(gene$SYMBOL,rownames(nrDEG))]
head(gene)
geneList=gene$logFC
names(geneList)=gene$ENTREZID 
#按照l(shuí)ogFC的值來(lái)排序geneList
geneList=sort(geneList,decreasing = T)
head(geneList)

以上即完成了數(shù)據(jù)的準(zhǔn)備工作，開(kāi)始進(jìn)行GSEA分析

GSEA-KEGG

library(clusterProfiler)
#clusterProfiler包的妙用
kk_gse <- gseKEGG(geneList     = geneList,
                  organism     = 'hsa',
                  nPerm        = 1000,
                  minGSSize    = 10,
                  pvalueCutoff = 0.9,
                  verbose      = FALSE)
#提取GSEA-KEGG結(jié)果
tmp=kk_gse@result
kk=DOSE::setReadable(kk_gse, OrgDb='org.Hs.eg.db',keyType='ENTREZID')
tmp=kk@result
pro='test_gsea'
write.csv(kk@result,paste0(pro,'_kegg.gsea.csv'))

#按照enrichment score從高到低排序，便于查看富集通路
  kk_gse=kk
  sortkk<-kk_gse[order(kk_gse$enrichmentScore, decreasing = T),]
  head(sortkk)
  dim(sortkk)
  #write.table(sortkk,"gsea_output2.txt",sep = "\\t",quote = F,col.names = T,row.names = F)
  #可以根據(jù)自己想要的通路畫(huà)出需要的圖
  library(enrichplot)
  gseaplot(kk_gse, "hsa04510")
  gseaplot2(kk_gse, "hsa04510", color = "firebrick",
            rel_heights=c(1, .2, .6))#改變更多參數(shù)，為了美觀
  
  #同時(shí)展示多個(gè)pathways的結(jié)果。
  #畫(huà)出排名前四的通路
  gseaplot2(kk_gse, row.names(sortkk)[1:4])
  
  #上圖用的是ES排名前4個(gè)畫(huà)圖，也可以用你自己感興趣的通路畫(huà)圖
  # 自己在剛才保存的txt文件里挑選就行。
  paths <- c("hsa04510", "hsa04512", "hsa04974")
  paths <- row.names(sortkk)[5:8]
  paths
  gseaplot2(kk_gse, paths)
  
  #這里的GSEA分析其實(shí)由三個(gè)圖構(gòu)成，GSEA分析的runningES折線圖
  # 還有下面基因的位置圖，還有所謂的蝴蝶圖。如果不想同時(shí)展示，還可以通過(guò)subplots改變。
  gseaplot2(kk_gse, paths, subplots=1)#只要第一個(gè)圖
  gseaplot2(kk_gse, paths, subplots=1:2)#只要第一和第二個(gè)圖
  gseaplot2(kk_gse, paths, subplots=c(1,3))#只要第一和第三個(gè)圖
  
  #如果想把p值標(biāo)上去，也是可以的。
  gseaplot2(kk_gse, paths, color = colorspace::rainbow_hcl(4),
            pvalue_table = TRUE)
  
  #最后的總結(jié)代碼
  gseaplot2(kk_gse,#數(shù)據(jù)
            row.names(sortkk)[1:5],#畫(huà)哪一列的信號(hào)通路
            title = "Prion disease",#標(biāo)題
            base_size = 15,#字體大小
            color = "green",#線條的顏色
            pvalue_table = TRUE,#加不加p值
            ES_geom="line")#是用線，還是用d點(diǎn)

###當(dāng)然，這里不知道具體需要什么通路，就全部都畫(huà)出來(lái)
# 這里找不到顯著下調(diào)的通路，可以選擇調(diào)整閾值，或者其它。
down_kegg<-kk_gse[kk_gse$pvalue<0.1 & kk_gse$enrichmentScore < -0.4,];down_kegg$group=-1
up_kegg<-kk_gse[kk_gse$pvalue<0.1 & kk_gse$enrichmentScore > 0.4,];up_kegg$group=1
dat=rbind(up_kegg,down_kegg)
colnames(dat)
dat$pvalue = -log10(dat$pvalue)
dat$pvalue=dat$pvalue*dat$group 
dat=dat[order(dat$pvalue,decreasing = F),]

library(ggplot2)
library(enrichplot)
gesa_res=kk@result

###畫(huà)出每張kegg通路
lapply(1:nrow(down_kegg), function(i){ 
  gseaplot2(kk,down_kegg$ID[i],
            title=down_kegg$Description[i],pvalue_table = T)
  ggsave(paste0(pro,'_down_kegg_',
                gsub('/','-',down_kegg$Description[i])
                ,'.pdf'))
})
lapply(1:nrow(up_kegg), function(i){ 
  gseaplot2(kk,up_kegg$ID[i],
            title=up_kegg$Description[i],pvalue_table = T)
  ggsave(paste0(pro,'_up_kegg_',
                gsub('/','-',up_kegg$Description[i]),
                '.pdf'))
})

GSEA-GO

GO和KEGG分析流程基本相同，除了函數(shù)名和變量名的變化

ego <- gseGO(geneList     = geneList,
             OrgDb        = org.Hs.eg.db,
             ont          = "ALL",
             nPerm        = 1000,   ## 排列數(shù)
             minGSSize    = 100,
             maxGSSize    = 500,
             pvalueCutoff = 0.9,
             verbose      = FALSE)  ## 不輸出結(jié)果


go=DOSE::setReadable(ego, OrgDb='org.Hs.eg.db',keyType='ENTREZID')
go=ego@result
pro='test_gsea'

go_gse=go
sortgo<-go_gse[order(go_gse$enrichmentScore, decreasing = T),]
head(sortgo)
dim(sortgo)
#write.table(sortkk,"gsea_output2.txt",sep = "\\t",quote = F,col.names = T,row.names = F)
#可以根據(jù)自己想要的通路畫(huà)出需要的圖
library(enrichplot)
gseaplot2(go_gse,#數(shù)據(jù)
          row.names(sortgo)[1:5],#畫(huà)那一列的信號(hào)通路
          title = "Prion disease",#標(biāo)題
          base_size = 15,#字體大小
          color = "green",#線條的顏色
          pvalue_table = TRUE,#加不加p值
          ES_geom="line")#是用線，還是用d點(diǎn)
write.csv(go,file = 'gse_go.csv') 

其他可視化方法

氣泡圖

dotplot(kk,split=".sign")+facet_wrap(~.sign,scales="free")

條形圖

down_kegg<-kk_gse[kk_gse$pvalue<0.1 & kk_gse$enrichmentScore < -0.4,];down_kegg$group=-1
up_kegg<-kk_gse[kk_gse$pvalue<0.1 & kk_gse$enrichmentScore > 0.4,];up_kegg$group=1
dat=rbind(up_kegg,down_kegg)
colnames(dat)
dat$pvalue = -log10(dat$pvalue)
dat$pvalue=dat$pvalue*dat$group 
dat=dat[order(dat$pvalue,decreasing = F),]
library(ggplot2)


g_kegg<- ggplot(dat, aes(x=reorder(Description,order(pvalue, decreasing = F)), y=pvalue, fill=group)) + 
  geom_bar(stat="identity") + 
  scale_fill_gradient(low="blue",high="red",guide = FALSE) + 
  scale_x_discrete(name ="Pathway names") +
  scale_y_continuous(name ="log10P-value") +
  coord_flip() + theme_bw(base_size = 15)+
  theme(plot.title = element_text(hjust = 0.5),  axis.text.y = element_text(size = 15))+
  ggtitle("Pathway Enrichment") 
g_kegg

網(wǎng)絡(luò)圖

library(ggplot2)
library(enrichplot)

cnetplot(kk,showCategory= 5, foldChange= geneList, colorEdge="TRUE")
#colorEdge不同的term展示不同的顏色，如果希望標(biāo)記節(jié)點(diǎn)的子集，可以使用node_label參數(shù)，它支持4種可能的選擇(即“category”、“gene”、“all”和“none”).

ego2<-pairwise_termsim(ego)
emapplot(ego2, showCategory= 10, color="pvalue", cex_category=1, layout="kk")
#cex_category調(diào)節(jié)節(jié)點(diǎn)大小，showCategory展示條目個(gè)數(shù)

參考來(lái)源

#section 3已更新#「生信技能樹(shù)」單細(xì)胞公開(kāi)課2021_嗶哩嗶哩_bilibili

史上最全GSEA可視化教程，今天讓你徹底搞懂GSEA！

詳解：基因集富集分析GSEA

enrichplot||基因富集結(jié)果可視化解決方案

致謝
I thank Dr.Jianming Zeng(University of Macau), and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes.

THE END