如何實現GSEA-基因富集分析?

基因富集分析(Gene Set Enrichment Analysis,GSEA)是一種針對全基因組表達譜晶元數據的分析方法,將基因與預定義的基因集進行比較。即綜合現有的對基因的定位、性質、功能、生物學意義等信息基礎,構建一個分子標籤資料庫,在此資料庫中將已知基因按照染色體位置、已建立基因集、模序、腫瘤相關基因集和GO基因集等多個功能基因集進行分組與歸類。通過分析基因表達譜數據,了解它們在特定的功能基因集中的表達狀況,以及這種表達狀況是否存在某種統計學顯著性。

統計過程:1.計算富集分數。2.估計富集分數的顯著程度。3.校正多重假設檢驗。

載入超時,點擊重試

流程圖

工具:

GSEA軟體下載:software.broadinstitute.org 要下載到Java,這個是在Java基礎上運行的軟體,根據你的數據大小,選擇不同內存的版本,2G內存開始的GSEA版本需要的是64位的Java 1.8版。

(操作教程:GSEA | Desktop Tutorial)

軟體界面

數據準備:主要準備一個表達矩陣和一個分組說明的cls文件,軟體界面如上圖,操作簡單,按照步驟Load data and run就行了,比較需要注意的是準備表達矩陣,如果選取的是GEO的公共數據集,就要將數據集進行預處理(採用R/bioconductor Affy和affyPLM程序包對數據集原始CEL文件進行質量控制後,使用Affy程序包中rma演算法對該數據集進行進行預處理。),因為GSEA只支持特定的格式,所以要剔除不必要的信息,將癌組織和對應的癌旁組織的數據分別提取出來分別作為兩組的表達矩陣(gct文件)以及分組文件(cls文件)(此步驟可以手動excel整理也可以找個代碼模板用R來操作)

data preparing:

1.如果是自己已經排序好了的基因,可以直接拿來做GSEA分析的見: GSEAPreranked Page in the GSEA User Guide.

2.如果是affymetrix的表達矩陣,不需要提前進行Present/Marginal/Absent Calls. 來過濾掉一些表達探針,GSEA需要各種情況的表達數據。

3.如果是gct and pcl 的表達矩陣,缺失值空著就好了。但是如果缺失值太多了,這樣在計算signal-to-noise的時候,不同group的樣本數就不一致了,mean和SD都會變好,最好是避免這樣的情況,可以考慮進行插值,或者過濾掉這樣的探針。

我是表達矩陣

我是分組文件

txt文檔格式會不一樣,GSEA有給出模板,照著修改就OK,如果格式有誤或數據有問題GSEA會報錯的。(格式參考說明書:Data formats - GeneSetEnrichmentAnalysisWiki)

load data

設置參數

成功導入數據後,點擊RUN GSEA,這時候要指定幾個參數的選擇,就是你要用哪些標籤資料庫來進行分析,以及如何分組等。

1. Expression dataset:輸入的表達矩陣

2. Gene sets database:分析的資料庫

3. Number of permutations:置換檢驗的次數

4. Phenotype labels:選擇比較組,如果你輸入的文件就只有2個組別的話,這個就很方便選一個就行了;如果你輸入的有三個組別及以上的話,則這裡就要跟你的需要選擇兩個組別的比較組,而且GSEA也會根據你的組別信息去表達矩陣中提取相對應的數據。

5. Collapse dataset to gene symbols: 如果你已經ID轉化為HUGO gene symbol,那麼這裡選FALSE,否則選擇TRUE。

6. Permutation type:選擇置換的類型,是random phenotype還是random gene sets,一般每組樣本數目大於7個時,建議選擇phenotype,否則選擇gene sets。

Chip platform:選擇晶元類型,是對ID進行注釋,即ID轉化,選擇ID對應的chip文件即可,如果已自行轉化了ID的話,則空著就行(那麼Collapse dataset to gene symbols應選擇否)

提交之後,如果運行失敗會出error提示,成功的話直接進入success的界面。

結果的解讀:

431/899表示在WT這一分組中,一共有899個功能基因集,其中421個上升

99個基因集的FDE小於25%

118個基因的名義P值小於1%

118個基因的名義P值小於5%

點擊snapshot可以看富集結果,就是下圖Enrichment plot

點擊enrichment result in html 可以查看所有的富集分析結果,進去之後可以點開查看每個Enrichment plot的參數。

點擊enrichment result in excel就可以直接下載附帶結果的excel。

SIZE:表示基因集里的基因數量

ES(enrichment score):富集分數

NES(normalized enrichment score):表示校正後的富集分數

NOM p-val (nominal p value ): 名義P值

FDR q-val(false discovery rate):錯誤發現率

FWER p-val:用bonferonni校正後的P值

RANK AT AMX:ES值對應的通路基因排名

Leading-edge subset:對富集貢獻最大的基因成員,即領頭亞集,用於定義Leading-edge subset的參數有:Tags,List,Signal。

Enrichment plot

當Enrichment plot過多的時候,可以整理成如下的表格形式展現在文章中。

在這些enrichment plot中,我們最關注的四個指標為ES值、NES、NOM p-val、FDR。

綠色曲線就是gene set裡面對應的每個基因的enrichment score值(ES),開始時為零,從左到右每遇到一個基因就計算出一個ES值,連成一條綠線。當ES值大於0時,表示某一功能基因富集在排序序列的前端,若為小於0時,則某一功能基因富集在排序序列的後端,ES值越高說明這些基因在通路中有富集,非散在分布。中間條形碼似的黑線是gene set裡面的基因在背景基因里的位置,每條豎線代表該通路下的基因,從左到右按照表達水平排序。Leading-edge subset(對富集貢獻最大的基因成員,即領頭亞集);在ES圖中出現領頭亞集的形狀,表明這個功能基因集在某處理條件下具有更顯著的生物學意義;對於結果的分析,通常認為|NES|>1,NOM p-val<0.05,FDR q-val<0.25的通路下的基因集合是有意義的;NES的絕對值越大,FDR值就越小,說明分析的結果可信度越高。NOM p-val是針對某一功能基因集得到的ES值的統計顯著性,P值越小,說明基因的富集性越好,但P值很小時,FDR值也可能很大,這說明和其他功能基因子相比較,它的富集並不是很顯著,原因可能是數據樣本量較少、雜交信號微弱或者是選擇的功能基因子集並未很好得反映樣本的物理學意義。

ES score的演算法

基因富集的熱圖

熱圖用5種顏色來表示基因表達水平的高低水平

蝴蝶圖顯示的是基因順序和排序度量得分之間的正相關和負相關的關係。

參考文獻:

1.From the Cover: Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles

2.GSEA (GSEA小組官網)

3.software.broadinstitute.org (說明書)

4.GSEA學習筆記

5.GSEA富集分析 - 界面操作

6. 基因探針富集分析(GSEA)翻譯+心得 (作者為為)

7.3.GSEA-基因富集分析

推薦閱讀:

TAG:生物統計學 | 基因 | 基因組 |