Gene set enrichment analysis, 基因集富集分析

1、算法思想

對某基因Knock down和WT做差異表達分析之后尊搬，得到差異表達基因list，按照某個統(tǒng)計量土涝，比如fold change佛寿，從小到大排序，得到一個rank list但壮，記錄為L

假設(shè)某個通路的所有基因在L中隨機分布冀泻，假如我們能觀測到某通路的所有基因突然富集與L中的一端，計算其富集程度和統(tǒng)計顯著性茵肃，如果小于某個cutoff腔长，那么我們就可以拒絕原假設(shè)，認為該通路在L中富集验残，并且通過富集程度的打分捞附，如果為正，則該通路傾向于在上調(diào)的基因中富集，如果為負鸟召，則該通路傾向于在下調(diào)的基因中富集胆绊。

在GSEA對應(yīng)的軟件中，用normalized enrichment score(NES)作為富集程度的度量欧募，用p值和FDR作為統(tǒng)計顯著性的度量压状，

收集的基因集的數(shù)據(jù)庫叫做MSigDB；

特定的基因集合可以從GO跟继、KEGG种冬、Reactome、hallmark或MSigDB等基因集中獲取舔糖，其中MSigDB數(shù)據(jù)庫整合了上述所有基因集娱两。研究者也可自定義gene set (即新發(fā)現(xiàn)的基因集或其它感興趣的基因的集合)

image

H: hallmark gene sets （效應(yīng)）特征基因集合，共50組金吗；
C1: positional gene sets 位置基因集合十兢，根據(jù)染色體位置，共326個摇庙；
C2: curated gene sets：（專家）共識基因集合旱物，基于通路、文獻等（包括KEGG）卫袒；
C3: motif gene sets：模式基因集合宵呛，主要包括microRNA和轉(zhuǎn)錄因子靶基因兩部分；
C4: computational gene sets：計算基因集合玛臂，通過挖掘癌癥相關(guān)芯片數(shù)據(jù)定義的基因集合烤蜕；
C5: GO gene sets：Gene Ontology 基因本體論（包括BP/CC/MF）；
C6: oncogenic signatures：癌癥特征基因集合迹冤，大部分來源于NCBI GEO 未發(fā)表芯片數(shù)據(jù)讽营；
C7: immunologic signatures: 免疫相關(guān)基因集合。

2泡徙、與GO的區(qū)別

類似但又有所不同：GO分析更加依賴差異基因橱鹏，實則是對一部分基因的分析 (忽略差異不顯著的基因)，而GSEA是從全體基因的表達矩陣中找出具有協(xié)同差異 (concordant differences)的基因集堪藐，故能兼顧差異較小的基因

GO富集是定性的分析莉兰，GSEA考慮到了表達或其它度量水平的值的影響。GSEA分析不需要指定閾值（p值或FDR）來篩選差異基因礁竞，在沒有經(jīng)驗存在的情況下分析我們感興趣的基因集糖荒，而這個基因集不一定是顯著差異表達的基因。GSEA分析可以將那些GO/KEGG富集分信息中容易遺漏掉的差異表達不顯著卻有著重要生物學(xué)意義的基因包含在內(nèi)模捂。

另外捶朵，對于時間序列數(shù)據(jù)或樣品有定量屬性時蜘矢，GSEA的優(yōu)勢會更明顯，不需要每個分組分別進行富集综看，直接對整體進行處理品腹。可以類比于WGCNA分析红碑。

差異基因排序度量的選取

case-control樣本而言（每個樣本至少三個重復(fù)）

μ為均值舞吭， 為標準差

image

連續(xù)表型的樣本，如時間序列表達譜析珊，可以選下列統(tǒng)計量
Pearson 相關(guān)系數(shù)
Cosine
Manhattan measure
Euclidean measure

3羡鸥、 GSEA計算中的幾個關(guān)鍵概念

計算富集得分 (ES, enrichment score).

ES反應(yīng)基因集成員s在排序列表L的兩端富集的程度。計算方式是忠寻，從基因集L的第一個基因開始兄春，計算一個累計統(tǒng)計值。當遇到一個落在s里面的基因锡溯，則增加統(tǒng)計值。遇到一個不在s里面的基因哑姚，則降低統(tǒng)計值

每一步統(tǒng)計值增加或減少的幅度與基因的表達變化程度（更嚴格的是與基因和表型的關(guān)聯(lián)度祭饭，可能是fold-change，也可能是pearson corelation值）是相關(guān)的叙量，可以是線性相關(guān)倡蝙，也可以是指數(shù)相關(guān)

富集得分ES最后定義為最大的峰值。正值ES表示基因集在列表的頂部富集绞佩，負值ES表示基因集在列表的底部富集寺鸥。

評估富集得分(ES)的顯著性

通過基于表型而不改變基因之間關(guān)系的排列檢驗 (permutation test)計算觀察到的富集得分(ES)出現(xiàn)的可能性。若樣品量少品山，也可基于基因集做排列檢驗 (permutation test)胆建，計算p-value。

多重假設(shè)檢驗校正
首先對每個基因子集s計算得到的ES根據(jù)基因集的大小進行標準化得到Normalized Enrichment Score (NES)肘交。隨后針對NES計算假陽性率笆载。（計算NES也有另外一種方法，是計算出的ES除以排列檢驗得到的所有ES的平均值）

Leading-edge subset涯呻，對富集得分貢獻最大的基因成員凉驻，領(lǐng)頭亞集

該處有3個統(tǒng)計值，tags=59%表示核心基因占該基因集中基因總數(shù)的百分比复罐；list=21%表示核心基因占所有基因的百分比涝登；signal=74%，將前兩項統(tǒng)計數(shù)據(jù)結(jié)合在一起計算出的富集信號強度效诅，計算公式如下：

image

其中n是列表中的基因數(shù)目胀滚，nh是基因集中的基因數(shù)目

4趟济、結(jié)果解釋

富集分析的圖示，該圖示分為三部分蛛淋，在圖中已做標記：

image

第一部分是Enrichment score折線圖：顯示了當分析沿著排名列表按排序計算時咙好，ES值在計算到每個位置時的展示。最高峰處的得分 (垂直距離0.0最遠)便是基因集的ES值褐荷。

第二部分勾效，用線條標記了基因集合中成員出現(xiàn)在基因排序列表中的位置，黑線代表排序基因表中的基因存在于當前分析的功能注釋基因集叛甫。leading edge subset 就是（0,0）到綠色曲線峰值ES出現(xiàn)對應(yīng)的這部分基因层宫。在ES圖中出現(xiàn)領(lǐng)頭亞集的形狀，表明這個功能基因集在某處理條件下具有更顯著的生物學(xué)意義

第三部分是排序后所有基因rank值得分布其监，熱圖紅色部分對應(yīng)的基因在NGT中高表達萌腿，藍色部分對應(yīng)的基因在DMT中高表達，每個基因?qū)?yīng)的信噪比（Signal2noise抖苦，前面選擇的排序值計算方式）以灰色面積圖顯展示毁菱。

image

一般認為|NES|>1，NOM p-val<0.05锌历，F(xiàn)DR q-val<0.25的通路是顯著富集的贮庞。

5、操作

南方醫(yī)科大的余光創(chuàng)教授（人稱Y叔）開發(fā)出了特別好用的進行各種富集分析的R包clusterProfiler究西，這個包里也集成了GSEA分析的功能窗慎，有兩種算法可選，一種是Y叔自己寫的DOSE包里的卤材，根據(jù)GSEA原始文獻的算法實現(xiàn)的函數(shù) 遮斥，一種是俄國人開發(fā)的fgsea包實現(xiàn)的快速GSEA算法，默認為fgsea算法扇丛。

clusterProfiler

gsea.raw<-GSEA(geneList = gene_List.entrez,nPerm=1000,pvalueCutoff = 1,TERM2GENE=gene_set)
print(head(gsea.raw@result)) 

gsea.go<- gseGO(geneList = gene_List.entrez,ont="BP",OrgDb = org.Hs.eg.db,keyType = "ENTREZID",pvalueCutoff = 1,nPerm = 2)
 print(head(gsea.go@result))

gsea.kegg<- gseKEGG(geneList = gene_List.entrez,organism = "hsa",keyType = "kegg",pvalueCutoff = 1,nPerm = 2)
 print(head(gsea.kegg@result)) 

可視化

 ridgeplot(Go_Reactomeresult, 10) # 輸出前十個結(jié)果术吗，查看terms對應(yīng)基因的上下調(diào)方向

image

barplot() #圖的呈現(xiàn)效果和常規(guī)GO分析一樣

 dotplot() #圖的呈現(xiàn)效果和常規(guī)GO分析一樣

cowplot::plot_grid(p1, p2, p3) # use plot_grid() to combine multiple plot in one figure.

cnetplot() # 獲取不同富集terms間的共同基因，加上表達數(shù)據(jù)也可以看上或下調(diào)方向

image

heatplot() # 獲取不同富集terms間的共同基因晕拆，加上表達數(shù)據(jù)也可以看上或下調(diào)方向

image

emapplot() # 不同terms之間的關(guān)系

image

compareCluster()

emapplot()

image

upsetplot() # 不同富集terms之間的共同基因

image

browseKEGG(kk, "hsa04934")

image

也可以繪制特定通路

pathview(gene.data  = geneList,
                     pathway.id = "hsa04750", #上述結(jié)果中的hsa04750通路
                     species    = "hsa",
                     limit      = list(gene=max(abs(geneList)), cpd=1))

寫在后面

圖可以各種海量得畫藐翎，終歸是術(shù)，但要想給這些圖一個有趣的靈魂实幕，還需要設(shè)計出好的故事吝镣，修煉道的水準才能真正提高文章的質(zhì)量～