? ? ? ?上一期中給大家介紹了擬時(shí)序分析的意義及具體的分析過(guò)程，本期繼續(xù)給大家?guī)?lái)轉(zhuǎn)錄組個(gè)性化分析——GSEA。廢話不多說(shuō)，干貨直接奉上！
1 GSEA基本概念
? ? ? ?基因集富集分析（Gene Set Enrichment Analysis，GSEA）：用一個(gè)預(yù)先定義的基因集中的基因來(lái)評(píng)估在與表型相關(guān)度排序的基因表中的分布趨勢(shì)，從而判斷其對(duì)表型的貢獻(xiàn)。
2 GSEA原理
(1)?背景基因排序：將全部基因按照某種指標(biāo)（差異分析p值，表型相關(guān)性，表達(dá)量等）進(jìn)行排序，比如log2FC排序。
(2)?目標(biāo)基因富集：將某個(gè)特定類型的基因在排序表中標(biāo)出，目標(biāo)基因可以是某個(gè)通路或GO terms的基因等。
(3)?計(jì)算富集分?jǐn)?shù)：使用加權(quán)法，計(jì)算ES值變化。對(duì)位于中部（與性狀相關(guān)性低）的部分采用較小的權(quán)值，所以越集中在兩端，與表型的相關(guān)性越高。ES曲線最大值為富集分?jǐn)?shù)（Enrichment Score）。
(4) Permutation test：對(duì)基因集的ES值進(jìn)行顯著性檢驗(yàn)及多重假設(shè)檢驗(yàn)，從而計(jì)算出顯著富集的基因集。
3 GSEA分析的作用
? ? ? ?GSEA和常規(guī)的GO、KEGG的差異在于，GSEA使用的是基因集，傳統(tǒng)的富集分析不需要考慮基因表達(dá)量的變化趨勢(shì)，其算法的核心只關(guān)注這些差異基因的分布是否和隨機(jī)抽樣得到的分布一致，即使后期在可視化時(shí)，我們?cè)谕穲D上用不同顏色標(biāo)記了上下調(diào)的基因，但是由于沒(méi)有采用有效的統(tǒng)計(jì)學(xué)手段去分析這條通路下所有差異基因的總體變化趨勢(shì)，這使得傳統(tǒng)的富集分析結(jié)果無(wú)法回答如下問(wèn)題：一個(gè)富集到的通路下，既有上調(diào)差異基因，也有下調(diào)差異基因，那么這條通路總體的表現(xiàn)形式究竟是怎樣呢，是被抑制還是激活？或者更直觀點(diǎn)說(shuō)，這條通路下的基因表達(dá)水平在實(shí)驗(yàn)處理后是上升了呢，還是下降了呢？
? ? ? ?想要知道進(jìn)行了差異分析的兩組有什么功能和通路的差別，手上有大部分的功能分子以及對(duì)應(yīng)的值，這個(gè)值可以是 logFC。可以用這個(gè) logFC作為分子的排序，從而來(lái)評(píng)估在預(yù)先定義的基因集中是否顯著富集。
4 GSEA分析代碼
library(DO.db)
require(DOSE)
library(clusterProfiler)
library(AnnotationHub)
library(readr)
library("genefilter")
library(pheatmap)
library(tidyverse)
library(DESeq2)
library(ggplot2)
library(export)
library(enrichplot)
library(Rgraphviz)
#構(gòu)造圖片輸出函數(shù)，need input filename width height
#函數(shù)依賴export包
out_img <- function(filename,pic_width=5,pic_height=7){
?graph2png(file=filename,width=pic_width,height=pic_height)
?graph2ppt(file=filename,width=pic_width,height=pic_height)
?graph2tif(file=filename,width=pic_width,height=pic_height)
}
#all_entrez.csv的第一列是entrezid,第二列是FoldChange的值。
#gse需要單獨(dú)做數(shù)據(jù)格式
d <- read.csv("all_entrez.csv")
geneList <- d[,2]
names(geneList)=as.factor(d[,1])
geneList <- sort(geneList,decreasing = TRUE)
#gseGO進(jìn)行GSEA分析
#參考連接https://yulab-smu.github.io/clusterProfiler-book/chapter12.html
###gseBP <- gseGO(geneList=geneList,ont="BP",OrgDb=maize,keyType = 'ENTREZID',nPerm = 50000,minGSSize = 100,maxGSSize = 6000,pvalueCutoff = 0.05,verbose = FALSE)
############# GSEA CC 模式 start
ego3 <- gseGO(geneList = geneList,OrgDb = maize,ont = "CC",nPerm = 1000,minGSSize = 100,maxGSSize = 1000,pvalueCutoff = 0.05,verbose = FALSE)
write.csv(ego3,file = "GESA-GO_CC.csv")
#ridgeline plot for expression distribution of GSEA result
ridgeplot(ego3)
out_img(filename = "ridgeplot_CC",pic_width = 12,pic_height = 12)
#只顯示值最高的一組的信息
#gseaplot(ego3,geneSetID = 1,by="runningScore",title=ego3$Description[1])
#gseaplot(ego3,geneSetID = 1,by="preranked",title=ego3$Description[1])
#gseaplot(ego3,geneSetID = 1,title=ego3$Description[1])
#顯示前4組信息
gseaplot2(ego3,geneSetID = 1:4, ES_geom = "dot",pvalue_table = TRUE)
out_img(filename = "gseaplot_CC",pic_width=12,pic_height = 10)
#gsearank(ego3,1,title=ego3[1,"Description"])
############GSEA CC 模式end
############# GSEA BP 模式 start
ego2 <- gseGO(geneList = geneList,OrgDb = maize,ont = "BP",pvalueCutoff = 0.05,verbose = FALSE)
write.csv(ego2,file = "GESA-GO_BP.csv")
#ridgeline plot for expression distribution of GSEA result
ridgeplot(ego2)
out_img(filename = "ridgeplot_BP",pic_width = 12,pic_height = 12)
#只顯示值最高的一組的信息
#gseaplot(ego3,geneSetID = 1,by="runningScore",title=ego3$Description[1])
#gseaplot(ego3,geneSetID = 1,by="preranked",title=ego3$Description[1])
#gseaplot(ego3,geneSetID = 1,title=ego3$Description[1])
#顯示前4組信息
gseaplot2(ego2,geneSetID = 1:4, ES_geom = "dot",pvalue_table = TRUE)
out_img(filename = "gseaplot_BP",pic_width=12,pic_height = 10)
#gsearank(ego3,1,title=ego3[1,"Description"])
############GSEA BP 模式end
############# GSEA MF 模式 start
ego4 <- gseGO(geneList = geneList,OrgDb = maize,ont = "MF",pvalueCutoff = 0.05,verbose = FALSE)
write.csv(ego4,file = "GESA-GO_MF.csv")
#ridgeline plot for expression distribution of GSEA result
ridgeplot(ego4)
out_img(filename = "ridgeplot_MF",pic_width = 12,pic_height = 12)
#只顯示值最高的一組的信息
#gseaplot(ego3,geneSetID = 1,by="runningScore",title=ego3$Description[1])
#gseaplot(ego3,geneSetID = 1,by="preranked",title=ego3$Description[1])
#gseaplot(ego3,geneSetID = 1,title=ego3$Description[1])
#顯示前4組信息
gseaplot2(ego4,geneSetID = 1:4, ES_geom = "dot",pvalue_table = TRUE)
out_img(filename = "gseaplot_MF",pic_width=12,pic_height = 10)
#gsearank(ego3,1,title=ego3[1,"Description"])
############GSEA MF 模式end
#gsaKEGG基因富集分析
kk2 <- gseKEGG(geneList = geneList,organism = 'zma',pvalueCutoff = 0.05,verbose = FALSE)
write.csv(kk2,file="GSEA_KEGG.csv")
gseaplot2(kk2,geneSetID = 1:4,ES_geom = "dot",pvalue_table = TRUE)
out_img(filename="GSEA_KEGG",pic_width = 12,pic_height = 10)
ridgeplot(kk2)out_img(filename="ridgeplot_GSEA_KEGG",pic_width = 12,pic_height = 12)
5 結(jié)果解讀
? ? ? ?典型結(jié)果圖由上、中、下三個(gè)部分組成：
? ? ? ?上：為富集評(píng)分的情況，如果 NES 為正，則峰出現(xiàn)在左側(cè)（頭部富集）（高表達(dá)組富集）基因集中核心分子主要集中在左側(cè)高表達(dá)組中；如果 NES 為負(fù)，則尾部會(huì)出現(xiàn)谷（尾部富集）（低表達(dá)組富集），基因集中核心分子主要集中在右側(cè)低表達(dá)組中。
? ? ? ?中：每一根豎線代表基因集中一個(gè)分子，上傳數(shù)據(jù)的分子根據(jù)給定的值進(jìn)行排序，排序后單獨(dú)提取當(dāng)前基因集中的定義的分子，分子的位置情況即為中間部分的所示。
? ? ? ?下：把上傳數(shù)據(jù)分子給定的值進(jìn)行歸一化后的值進(jìn)行可視化。下部分的結(jié)果可以不用怎么關(guān)注。
? ? ? ?從該圖中可以看出，這個(gè)基因集是在MUT這一組高表達(dá)的，在折線圖中有個(gè)峰值，該峰值就是這個(gè)基因集的Enrichemnt score，峰值之前的基因就是該基因集下的核心基因。基因集下的所有基因在這個(gè)排序列表的頂部富集，我們可以說(shuō)，從總體上看，該基因集是上調(diào)趨勢(shì)，反之，如果在底部富集，則是下調(diào)趨勢(shì)。
