基本概念
Gene Set Enrichment Analysis (GSEA)是一種的計算方法，用于評估在一個預設訂的基因集中，兩個生物學數(shù)據(jù)集間，是否存在統(tǒng)計學上顯著且一致的差異。
Gene Set Enrichment Analysis (GSEA) is a computational method that determines whether an a priori defined set of genes shows statistically
significant, concordant differences between two biological states.
GSEA富集分析和傳統(tǒng)富集分析的比較
傳統(tǒng)富集分析特點：
傳統(tǒng)的富集分析的主要分析過程包括：采取固定閾值的篩選差異基因（p值，差異倍數(shù)），然后對其進行功能/通路富集分析（GO/KEGG），簡單來說就是一種歸類的方法。這種方法存在一個問題，一個通路的各個基因的表達量差異很大，生物調控也是遞進關系，上游基因的微小變化可能就會導致下游基因的明顯變化，如果直接一刀切用定閾值篩選，會損失很多基因信息，導致后續(xù)的富集不顯著（基因數(shù)量太少做這種富集分析沒有生物學意義）。
GSEA富集分析的特點
GSEA的分析流程主要包括：計算指標、排序、標注基因集、計算基因集得分、置換檢驗。從流程就可以看出，GSEA富集分析不需要用固定閾值去過濾基因，是基于所有基因分析的一種方法，規(guī)避了傳統(tǒng)富集分析方法的短板。
因此，簡單來說GSEA分析適用于，傳統(tǒng)分析方法篩選后樣本過少的數(shù)據(jù)集。
工具
GSEA軟件下載：https://www.gsea-msigdb.org/gsea/index.jsp 直接選擇你操作系統(tǒng)對應的版本，下載安裝即可。
數(shù)據(jù)準備
基因表達量表（.gct）
橫坐標為樣本，縱坐標為基因，陳列基因表達量信息；第一行#1.2為信息注釋，不需要改動；第二行兩個數(shù)字分別為：基因數(shù)量和樣本數(shù)量（如實填寫）；.gct文件，直接把表格用excel保存為制表符分隔文件（.txt），然后直接改后綴為gct即可。
樣本分類關系表（.cls）
第一行三個數(shù)字分別為：樣品總數(shù)，分組數(shù)，恒定數(shù)字1；第二行為樣品分組名稱；第三行為樣品分組信息。需要注意：第二行和第三行需嚴格和.gct文件對應。
基因集合（.gmt）
第一列為基因集編號，第二列為基因集名稱，第三列以后為集合中所包含的所有基因名稱，注意基因名稱需要與.gct文件中的名稱一致。
基因集文件的準備相對比較復雜，我用到的是KEGG的通路基因集，也可以用GO的的基因集，甚至可以自己創(chuàng)造一個，只需要滿足上圖所示的文件格式。下面介紹一下如何從KEGG網(wǎng)站上下載某物種（小鼠）的基因集信息。
首先，在KEGG數(shù)據(jù)庫中下載小鼠的KEGG通路數(shù)據(jù)庫文件
大家可以看到，這個.keg文件和我們需要的.gmt文件還有不小的差異，那么，我們需要按如下方法，先在Linux下解析keg文件，然后用R制作.gmt文件：
#解析信息
perl -alne '{if(/^C/){/PATH:mmu(\d+)/;$kegg=$1}else{print "$kegg\t$F[1]" if /^D/ and $kegg;}}' mmu00001.keg > mmu00001_kegg2gene.txt
grep '^C' mmu00001.keg | grep "mmu" | sed 's/^C //g' | sed 's/\([0-9]*\)\s/\1\t/' > mmu00001_kegg2description.txt
#轉換ID同時制作.gmt文件
#安裝轉換ID需要用到的各個R包
#安裝data.table
install.packages("data.table")
#安裝org.Xx.eg.db系列包中小鼠資源包
BiocManager::install("org.Mm.eg.db")
#安裝clusterProfiler包，安裝時需要“科學上網(wǎng)”
BiocManager::install("clusterProfiler")
#調用各個R包
library(data.table)
library(org.Mm.eg.db)
library(clusterProfiler)
#讀入數(shù)據(jù)
path2gene_file<-"mmu00001_kegg2gene.txt"
tmp=read.table(path2gene_file,sep="\t",colClasses=c('character'))
dt <- tmp
dt$geneid <- rownames(dt)
# transform id
map_dt <- bitr(dt$V2, fromType = "ENTREZID",toType = c( "SYMBOL"),OrgDb = org.Mm.eg.db)
dt_merge <- merge(map_dt,dt, by.y = "V2", by.x = "ENTREZID")
# 可以把轉換ID后的表輸出
#write.table(dt_merge, file = "example42.txt", sep = "\t", col.names = T, row.names = T, quote = F)
# first column is kegg ID, second column is entrez ID
GeneID2kegg_list<- tapply(tmp[,1],as.factor(tmp[,2]),function(x) x)
kegg2symbol_list<- tapply(dt_merge$SYMBOL,as.factor(dt_merge$V1),function(x) x)
#輸出數(shù)據(jù)，但描述列為NA
write.gmt <- function(geneSet=kegg2symbol_list,gmt_file='kegg2symbol.gmt'){
sink( gmt_file )
for (i in 1:length(geneSet)){
cat(names(geneSet)[i])
cat('\tNA\t')
cat(paste(geneSet[[i]],collapse = '\t'))
cat('\n')
}
sink()
}
write.gmt()
GSEA分析
導入數(shù)據(jù)
打開GSEA軟件，點擊左側菜單欄的Load data選項，將數(shù)據(jù)按如下圖所示方法導入。
設置參數(shù)
按照下圖標記的批注設置合適的參數(shù)。
運行GSEA
設置好參數(shù)后點擊Run運行即可。
GSEA結果展示
GSEA的所有結果都保存在結果文件夾中，完整報告為index.html文件，雙擊打開即可查看。
172/230表示在這一分組中，一共有230個功能基因集，其中172個上升
97個基因集的FDE小于25%
76個基因的名義P值小于1%
87個基因的名義P值小于5%
點擊snapshot可以看富集結果
點擊enrichment result in html 可以查看所有的富集分析結果，進去之后可以查看每個Enrichment plot的參數(shù)
點擊enrichment result in excel就可以直接下載附帶結果的excel
我們重點關注enrichment result in html中的內容，點開后會出現(xiàn)如下圖所示的表格：
圖中：
SIZE：表示基因集里的基因數(shù)量
ES（enrichment score）：富集分數(shù)（重點）
NES（normalized enrichment score）:表示校正后的富集分數(shù)（重點）
NOM p-val (nominal p value ): 名義P值（重點）
FDR q-val（false discovery rate）：錯誤發(fā)現(xiàn)率（重點）
FWER p-val：用bonferonni校正后的P值
RANK AT AMX：ES值對應的通路基因排名
Leading-edge subset：對富集貢獻最大的基因成員,即領頭亞集，用于定義Leading-edge subset的參數(shù)有：Tags，List，Signal。
點擊每個通路的Details即可看到如下圖所示的典型GSEA數(shù)據(jù)展示圖：
綠色曲線就是gene set里面對應的每個基因的enrichment score值（ES），開始時為零，從左到右每遇到一個基因就計算出一個ES值，連成一條綠線。當ES值大于0時，表示某一功能基因富集在排序序列的前端，若為小于0時，則某一功能基因富集在排序序列的后端，ES值越高說明這些基因在通路中有富集，非散在分布。
中間條形碼似的黑線是gene set里面的基因在背景基因里的位置，每條豎線代表該通路下的基因。
蝴蝶圖：按相關系數(shù)對所有基因排序，大于0正相關，小于0位負相關，越靠近兩級相關性越弱。
對于結果的分析，通常認為|NES|>1，NOM p-val<0.05，F(xiàn)DR q-val<0.25的通路下的基因集合是有意義的；NES的絕對值越大，F(xiàn)DR值就越小，說明分析的結果可信度越高。NOM p-val是針對某一功能基因集得到的ES值的統(tǒng)計顯著性，P值越小，說明基因的富集性越好，但P值很小時，F(xiàn)DR值也可能很大，這說明和其他功能基因子相比較，它的富集并不是很顯著，原因可能是數(shù)據(jù)樣本量較少、雜交信號微弱或者是選擇的功能基因子集并未很好得反映樣本的物理學意義。
