一文掌握GSEA,超詳細教程!

生信寶典之前總結了一篇關於GSEA富集分析的推文——GSEA富集分析：從概念理解到界面實操，介紹了GSEA的定義、GSEA原理、GSEA分析、Leading-edge分析等，是全網最流行的原理+操作兼備教程，不太了解的朋友可以點擊閱讀先理解下概念 （為了完整性，下面也會摘錄一部分）。

GSEA案例解析

介紹GSEA分析之前，我們先看一篇Cell文章(https://sci-hub.tw/10.1016/j.cell.2016.11.033)的一個插圖 (SCI-HUB客戶端（文獻神器V4.0）——下載文獻如此簡單)。

以下是文章原文對圖的註解：GSEA analyses of genesets for cardiac (top) and endothelial/endocardial (bottom) development. NES, normalized enrichment score. FDR, false discovery rate. Positive and negative NES indicate higher and lower expression in iwt, respectively.

關於文章中使用的GSEA分析方法和參數，我們截取對應原文：Gene Set Enrichment Analysis was performed using the GSEA software (https://www.broadinstitute.org/gsea/) with permutation = geneset, metric = Diff_of_classes, metric = weighted, #permutation = 2500.

根據以上信息可知，上圖是研究者使用GSEA軟體所做的分析結果。文章通過GSEA分析，發現

與心臟發育有關的基因集 (影響心臟的收縮力、鈣離子調控和新陳代謝活力等)在iwt組 (GATA基因野生型)中普遍表達更高，而在G296S組 (GATA基因的一種突變體)中表達更低；

而對於參與內皮或內膜發育的基因集，在iwt組中表達更低，在G296S組中表達更高。

作者根據這個圖和其它證據推測iwt組的心臟發育更加完善，而G296S組更傾向於心臟內皮或內膜的發育，即GATA基因的這種突變可能導致心臟內皮或內膜的過度發育而導致心臟相關疾病的產生。

那麼GSEA分析是什麼？

參考GSEA官網主頁的描述：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 (e.g. phenotypes). 

在上述Cell文章中，作者更加關心參與心臟發育的基因集 (即a priori defined set of genes)與兩個狀態（突變體和野生型，狀態的度量方式是基因表達）的關係，因此利用GSEA對其進行分析後發現，參與心臟發育 (收縮力、鈣調控和新陳代謝)的基因集的表達模式更接近於iwt組的表型，而不是G296S組; 而參與心臟內皮或內膜發育的這些基因的表達模式更接近於G296S組的表型而不是iwt組的表型。

這就是GSEA分析所適用的主要場景之一。它能幫助生物學家在兩種不同的生物學狀態 (biological states)中，判斷某一組有特定意義的基因集合的表達模式更接近於其中哪一種。因此GSEA是一種非常常見且實用的分析方法，可以將數個基因組成的基因集與整個轉錄組、修飾組等做出簡單而清晰的關聯分析。

除了對特定gene set的分析，反過來GSEA也可以用於發現兩組樣本從表達或其它度量水平分別與哪些特定生物學意義的基因集有顯著關聯，或者發現哪些基因集的表達模式或其他模式更接近於表型A、哪些更接近於表型B。這些特定的基因集合可以從GO、KEGG、Reactome、hallmark或MSigDB等基因集中獲取，其中MSigDB資料庫整合了上述所有基因集。研究者也可自定義gene set (即新發現的基因集或其它感興趣的基因的集合)。

GSEA分析似乎與GO分析類似但又有所不同。GO分析更加依賴差異基因，實則是對一部分基因的分析 (忽略差異不顯著的基因)，而GSEA是從全體基因的表達矩陣中找出具有協同差異 (concordant differences)的基因集，故能兼顧差異較小的基因。因此二者的應用場景略有區別。另外GO富集是定性的分析，GSEA考慮到了表達或其它度量水平的值的影響。另外，對於時間序列數據或樣品有定量屬性時，GSEA的優勢會更明顯，不需要每個分組分別進行富集，直接對整體進行處理。可以類比於之前的WGCNA分析。

GSEA定義

Gene Set Enrichment Analysis (基因集富集分析)用來評估一個預先定義的基因集的基因在與表型相關度排序的基因表中的分布趨勢，從而判斷其對表型的貢獻。其輸入數據包含兩部分，一是已知功能的基因集 (可以是GO注釋、MsigDB的注釋或其它符合格式的基因集定義)，一是表達矩陣 (也可以是排序好的列表)，軟體會對基因根據其與表型的關聯度(可以理解為表達值的變化)從大到小排序，然後判斷基因集內每條注釋下的基因是否富集於表型相關度排序後基因表的上部或下部，從而判斷此基因集內基因的協同變化對表型變化的影響。

(The gene sets are defined based on prior biological knowledge, e.g., published information about biochemical pathways or coexpression in previous experiments. The goal of GSEA is to determine whether members of a gene set S tend to occur toward the top (or bottom) of the list L, in which case the gene set is correlated with the phenotypic class distinction.)

這與之前講述的GO富集分析不同。GO富集分析是先篩選差異基因，再判斷差異基因在哪些注釋的通路存在富集；這涉及到閾值的設定，存在一定主觀性並且只能用於表達變化較大的基因，即我們定義的顯著差異基因。而GSEA則不局限於差異基因，從基因集的富集角度出發，理論上更容易囊括細微但協調性的變化對生物通路的影響，尤其是差異倍數不太大的基因集。

GSEA原理

給定一個排序的基因表L和一個預先定義的基因集S (比如編碼某個代謝通路的產物的基因, 基因組上物理位置相近的基因，或同一GO注釋下的基因)，GSEA的目的是判斷S裡面的成員s在L裡面是隨機分布還是主要聚集在L的頂部或底部。這些基因排序的依據是其在不同表型狀態下的表達差異，若研究的基因集S的成員顯著聚集在L的頂部或底部，則說明此基因集成員對表型的差異有貢獻，也是我們關注的基因集。

GSEA計算中幾個關鍵概念：

計算富集得分 (ES, enrichment score). ES反應基因集成員s在排序列表L的兩端富集的程度。計算方式是，從基因集L的第一個基因開始，計算一個累計統計值。當遇到一個落在s裡面的基因，則增加統計值。遇到一個不在s裡面的基因，則降低統計值。

每一步統計值增加或減少的幅度與基因的表達變化程度（更嚴格的是與基因和表型的關聯度，可能是fold-change，也可能是pearson corelation值，後面有介紹幾種不同的計算方式）是相關的，可以是線性相關，也可以是指數相關 （具體見後面參數選擇）。富集得分ES最後定義為最大的峰值。正值ES表示基因集在列表的頂部富集，負值ES表示基因集在列表的底部富集。

評估富集得分(ES)的顯著性。通過基於表型而不改變基因之間關係的排列檢驗 (permutation test)計算觀察到的富集得分(ES)出現的可能性。若樣品量少，也可基於基因集做排列檢驗 (permutation test)，計算p-value。

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

Leading-edge subset，對富集得分貢獻最大的基因成員。

本文通過總結多人學習使用過程中遇到的問題進一步記錄軟體操作過程和結果解讀，力求講清每個需要注意的細節點。

從前文中我們了解到GSEA分析的目的是要判斷S集基因（基於先驗知識的基因注釋信息，某個關注的基因集合）中的基因是隨機分布還是聚集在排序好的L基因集的頂部或底部（這便是富集分析）。

與GO富集分析的差異在於GSEA分析不需要指定閾值（p值或FDR）來篩選差異基因，我們可以在沒有經驗存在的情況下分析我們感興趣的基因集，而這個基因集不一定是顯著差異表達的基因。GSEA分析可以將那些GO/KEGG富集分信息中容易遺漏掉的差異表達不顯著卻有著重要生物學意義的基因包含在內。

下面來看看軟體具體操作和結果解讀。

一、軟體安裝

軟體下載地址：http://software.broadinstitute.org/gsea/downloads.jsp

使用官方推薦的第一個軟體javaGSEA Desktop Application，根據分析數據的大小和電腦內存多少可以選擇下載不同內存版本的軟體。該軟體是基於java環境運行的，而且需要聯網。若會出現打不開的現象（小編就是就碰到了），要麼是沒有安裝java，要麼是java版本太低了，安裝或更新下java就能打開。也可能是網速太慢，或Java安全性問題，這時選擇官網提供的第二個軟體javaGSEA Java Jar file，同樣依賴java運行，但不需聯網，啟動快。

軟體啟動界面如下：

二、數據準備

所有矩陣的列以tab鍵分割，不同類型的數據格式和後綴要求見下表。

Data FileContentFormatSourceExpression datasetContains features (genes or probes), samples, and an expression value for each feature in each sample. Expression data can come from any source (Affymetrix, Stanford cDNA, and so on).res, gct, pcl, or txtYou create the file. 一般的基因表達矩陣整理下格式就可以。如果是其它類型數據或自己計算rank也可以，後面有更多示例。(如果後綴為txt格式，傳統的基因表達矩陣就可以，第一列為基因名字，名字與待分析的功能注釋數據集一致，同為GeneSymbol或EntrezID或其它自定義名字，第一行為標題行，含樣品信息。gct文件需要符合下面的格式要求。)Phenotype labelsContains phenotype labels and associates each sample with a phenotype.clsYou create the file or have GSEA create it for you. 一般是樣品分組信息或樣品屬性度量值或時間序列信息。Gene setsContains one or more gene sets. For each gene set, gives the gene set name and list of features (genes or probes) in that gene set.gmx or gmtYou use the files on the Broad ftp site, export gene sets from the Molecular Signature Database (MSigDb) or create your own gene sets file. 欲檢測是否富集的基因集列表。注意基因ID與表達矩陣基因ID一致。自己準備的基因集注意格式與官網提供的gmt格式一致。Chip annotationsLists each probe on a DNA chip and its matching HUGO gene symbol. Optional for the gene set enrichment analysis.ChipYou use the files on the Broad ftp site, download the files from the GSEA web site, or create your own chip file. 主要是為晶片探針設計的轉換文件。如果表達矩陣的基因名與注釋集基因名一致，不需要這個文件。

1. 表達數據集文件

GESA提供有Example Datasets，下載地址：http://software.broadinstitute.org/gsea/datasets.jsp。

在這裡可以下載表達矩陣Expression dataset(gct文件，常見txt格式也可以)和樣品分組信息Phenotype labels（cls文件）

數據示例中兩個gct文件都是表達矩陣，其中*hgu133a.gct文件第一列是探針名字，*collapsed.gct文件的第一列是gene symbol。

第一行：#1.2，表示版本號，自己準備文件時照抄就行；

第二行：兩個數分別表示gene NAME的數量和樣本數量（矩陣列數-2）；

矩陣：第一列是NAME；第二列Description，沒有的話可以全用na或任意字符串填充；後面的就是基因在不同樣本中標準化後的表達數據了 （部分統計量metrics for ranking genes計算需要log轉換後的數據，後面會有提及。其它情況是否為log轉換的數據都可用，GSEA關注的是差異，只要可比即可）。

2. 樣品分組信息

第一行：三個數分別表示：34個樣品，2個分組，最後一個數字1是固定的；

第二行：以#開始，tab鍵分割，分組信息（有幾個分組便寫幾個，多個分組在比較分析時，後面需要選擇待比較的任意2組）；

（樣品分組中NGT表示正常耐糖者，DMT表示糖尿病患者，自己使用時替換為自己的分組名字）

第三行：樣本對應的組名。樣本分組信息的第三行，同一組內的不同重複一定要命名為相同的名字，可以是分組的名字。例如相同處理的不同重複在自己試驗記錄裡一般是Treat6h_1、Treat6h_2、Treat6h_3，但是在這裡一定都要寫成一樣的值Treat6h。與表達矩陣的樣品列按位置一一對應，名字相同的代表樣品屬於同一組。如果是樣本分組信息，上圖中的0和1也可以對應的寫成NGT和DMT，更直觀。但是，如果想把分組信息作為連續表型值對待，這裡就只能提供數字。

3. 功能基因集文件（gene sets）

GSEA官網提供了8種基因分類資料庫，都是關於人類的數據，包括Marker基因，位置臨近基因，矯正過的基因集，調控motif基因集，GO注釋，癌基因，免疫基因，最新一次更新是在2018年7月，下載地址：http://software.broadinstitute.org/gsea/downloads.jsp#msigdb。

官網提供的gmt文件有兩種類型，*.symbols.gmt中基因以symbols號命名，*.entrez.gmt中基因以entrez id命名。注意根據表達矩陣的基因名字命名方式選擇合適的基因集。表達數據和通路數據能關聯在一起依賴的是基因名字相同，所以一定保證基因命名方式的統一。

 gmt格式是多列注釋文件，第一列是基因所屬基因集的名字，可以是通路名字，也可以是自己定義的任何名字。第二列，官方提供的格式是URL，可以是任意字符串。後面是基因集內基因的名字，有幾個寫幾列。列與列之間都是TAB分割。

Pathway_description Anystring Gene1 Gene2 Gene3
Pathway_description2 Anystring Gene4 Gene2 Gene3 Gene5

GSEA官網只提供了人類的數據，但是掌握了官網中基因表達矩陣和注釋文件的數據格式，就可以根據自己研究的物種，在公共資料庫下載對應物種的注釋數據，自己製作格式一致的功能基因集文件，這樣便就可以做各種物種的GSEA富集分析了。

4. 晶片注釋文件

如果分析的表達數據是晶片探針數據就需要用到晶片注釋文件（chip），用來做ID轉換，把探針名字轉換為基因名字。如果我們的表達數據文件中已經是基因名了就不再需要這個文件了。

三、分析參數設置和軟體運行

演示使用的數據來自GSEA官網：

表達矩陣：Diabetes_collapsed_symbols.gct

樣品分組信息：Diabetes.cls

基因功能分類數據選擇GO資料庫：c5.all.v6.2.symbols.gmt

因為表達矩陣與注釋中基因名字可以直接對應，第四個文件不需要

1. 數據導入

按照上圖步驟依次點擊Load data——Browse for file——在彈出文件框中找到待導入的文件，選中點擊打開即可；

若文件格式沒問題會彈出一個提示There were no error的框，證明文件上傳成功，並且會顯示在5所示的位置；若出錯，請仔細核對文件格式。

注意：1）本地文件存放路徑不要有中文、空格（用_代替空格）和其他特殊字符；2）所有用到的文件都需要通過上述方式先上傳至軟體；3）數據上傳錯誤後可以通過點擊工具欄file——clear recent file history進行清除。

2. 指定參數

點擊軟體左側Run GSEA，將跳出參數選擇欄。參數設置分為三個部分Require fields（必須設置的參數項）、Basic fields（基本參數設置欄）和Advanced fields（高級參數設置欄），後面兩欄的參數一般不做修改，使用默認的就行。後面兩部分參數設置，如果涉及到需要根據實驗數據做調整的地方，會在後面的分析中會提到。

1）Require fields

Expression dataset： 導入表達數據集文件，點擊後自動顯示上一步中從本地導入軟體內的文件，所以一定要確認上一步導入數據是否成功；

Gene sets database： 基因功能集資料庫，可以從本地導入（上一步）；

在聯網的情況下軟體也可以為自動下載GSEA官網中的gene sets文件；

Number of permutations： 置換檢驗的次數，數字越大結果越準確，但是太大會佔用太多內存，軟體默認檢驗1000次。

軟體分析時會得到一個基因富集的評分（ES），但是富集評分是否具有統計學意義，軟體就會採用隨機模擬的方法，根據指定參數隨機打亂1000次，得到1000個富集評分，然後判斷得到的ES是否在這1000個隨機產生的得分中有統計學意義。測試使用時建議填一個很小的數如10，先讓程序跑通。真正分析時再換為1000。

Phenotype labels： 選擇比較方式，如果文件只有2個組別的話就比較方便了，任意選一個就行，哪個在前在後全在自己怎麼解釋方便；如果數據有多組的話，GSEA會提供兩兩間比較的組合選項或者某一組與剩下所有組的比較。選擇好後，GSEA會在分析過程中根據組別信息自動到表達數據集文件中提取對應的數據作比較。

Collapse dataset to gene symbols： 如果表達數據集文件中NAME已經與gene sets database中名字一致，選擇FALSE，反之選擇TRUE。

Permutation type： 選擇置換類型，phenotype或者gene sets。

每組樣本數目大於7個時 ，建議選擇phenotype，否則選擇gene sets。

Chip platform： 表達數據集為晶片數據時才需要，目的是對ID進行注釋轉換，如果已經轉換好了就不需要了。應該也適用於其它需要轉換ID的情況，不過事先轉換最方便。

2）Basic fields

通常選擇默認參數即可，在此簡單介紹一下

Analysis name： 取名需要注意不能有空格，需要用_代替空格。如果做的分析多，最好選擇一個有意義的名字，比如shengxinbaodian （生信寶典全拼），方便查找。

Enrichment statistic： 基因集富集分析(PNAS)的最後一部分給出了GSEA中所用方法的數學描述，感興趣的可以查看一下論文。在此給出每種富集分析不同算法的參數情況：▪ classic: p=0 若基因存在，則ES值加1；若基因不存在，則ES值減1▪ weighted (default): p=1 若基因存在，則ES加rank值；若基因不存在，則ES減rank值▪ weighted_p2: p=2 基因存在，ES加rank值的平方，不存在則減rank值的平方▪ weighted_p1.5: p=1.5 基因存在，ES加rank值得1.5次方，不存在則減rank值得1.5次方

備註：如果想用其它加權，就自己計算rank值，使用preranked mode。

Metric for ranking genes： 基因排序的度量

下面提到的均值也可以是中位數。

如果表型是分組信息，GSEA在計算分組間的差異值時支持5種統計方式，分別是signal2noise、t-Test、ratio_of_class、 diff_of_class(log2轉換後的值計算倍數)和log2_ratio_of_class。

下面公式很清楚。

如果表型是連續數值信息（定量表型）: GSEA通過表型文件（cls）和表達數據集文件（gct），使用pearson相關性、Cosine、Manhattan 或Euclidean指標之一計算兩個配置文件之間的相關性。

（注意：若是分組表型文件想轉換為定量表型，cls文件中分類標籤應該指定為數字）

Gene list sorting mode： 對表達數據集中的基因進行排序，按照排序度量的真實值（默認）或者絕對值排序；

Gene list ordering mode： 使用此參數確定表達數據集中基因是按照降序（默認）或者升序排列；

Max size & Min size： 從功能基因集中篩選出不屬於表達數據集中的基因後，剩下基因總數在此範圍內則保留下來做後續的分析，否則將此基因集排除；一般太多或太少都沒有分析意義。

Save results in this folder： 在此可以選擇分析文件在本地電腦的存儲地址。

3）Advanced fields

Collapsing mode for probe sets => 1 gene：多個探針對應一個基因時的處理方式。

Normalized mode: 富集得分的標準化方式。

Randomization mode：只用於phenotype permutation。

Median for class metrics: 計算metrics ranking時用中值而不是平均值。

Number of markers：紅蝶圖中展示的Gene Marker數目。

Plot graphs for the top sets of each phenotype：繪製多少GSEA plot，默認top 20，其它不繪製。一般會把這個值調高。

Seed for permutation：隨機數種子，如果想讓每次結果一致，這裡需要設置同樣的一個整數。

以上參數都設置好後點擊參數設置欄下方的一個綠色按鈕Run，若軟體左下方GSEA reports處的狀態顯示Running的話則表示運行成功，此過程大概需要十分鐘左右，視數據大小而定。

四、結果解讀

數據分析完後的結果會保存到我們設置的路徑下，點開文件夾中的index.html就可以查看網頁版結果，更加方便。

結果報告分為多個子項目，其中最重要的是前面兩部分，基因富集結果就在這裡。從第三部分開始其實是軟體在分析數據的過程產生的中間文件， 也很重要，讀懂後可以加深對GSEA分析的認識，理解我們是如何從最初的基因表達矩陣得到最終的結果（即報告的前兩個項目）。建議先從Dataset details看起，然後再返回看第一部分的結果報告。

1. Enrichment in phenotype

以正常人組NGT的17個樣本數據為例解析最終結果。

報告首頁文字總結信息表示：

經過條件篩選後還剩下3953個GO條目，其中1697個GO條目在NGT組中富集；

有36個GO基因條目在FDR<25%的條件下顯著富集，這部分基因最有可能用於推進後續實驗；

在統計檢驗p<0.01, p<0.05的條件下分別有19和114個GO條目顯著富集；

結果有多種顯示方式：圖片快照（snapshot）、網頁（html）和表格（Excel）形式；

點擊Guide to可以查看官方幫助解讀結果的文檔。

1） 點擊enrichment results in html，在網頁查看富集結果，如下：

GS：基因集的名字，GO條目的名字

SIZE：GO條目中包含表達數據集文中的基因數目（經過條件篩選後的值）；

ES：富集評分；

NES：校正後的歸一化的ES值。

由於不同用戶輸入的基因資料庫文件中的基因集數目可能不同，富集評分的標準化考慮了基因集個數和大小。

其絕對值大於1為一條富集標準。

計算公式如下：

NOM p-val：即p-value，是對富集得分ES的統計學分析，用來表徵富集結果的可信度；

FDR q-val：即q-value，是多重假設檢驗校正之後的p-value，即對NES可能存在的假陽性結果的概率估計，因此FDR越小說明富集越顯著；

RANK AT MAX：當ES值最大時，對應基因所在排序好的基因列表中所處的位置；（註：GSEA採用p-value<5%，q-value<25%進行數據過濾）

LEADING EDGE：該處有3個統計值，tags=59%表示核心基因佔該基因集中基因總數的百分比；list=21%表示核心基因佔所有基因的百分比；signal=74%，將前兩項統計數據結合在一起計算出的富集信號強度，計算公式如下：

其中n是列表中的基因數目，nh是基因集中的基因數目

點擊Details跳轉至對應的詳情結果。只有前20個GO富集詳情可以查看，想要生成的結果報告可以查看更多的富集信息，可以通過在Advanced fields處設置參數Plot graphs for the top sets of each phenotype。

2）Details for gene set首先是一個選定GOset下的匯總信息表，每一部分意思在上面已做解釋，其中Upregulated in class表示該基因集在哪個組別中高表達，這個主要看富集分析後的leading edge分布位置。

接下來是富集分析的圖示，該圖示分為三部分，在圖中已做標記：

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

第二部分，用線條標記了基因集合中成員出現在基因排序列表中的位置，黑線代表排序基因表中的基因存在於當前分析的功能注釋基因集。leading edge subset 就是（0,0）到綠色曲線峰值ES出現對應的這部分基因。

第三部分是排序後所有基因rank值得分布，熱圖紅色部分對應的基因在NGT中高表達，藍色部分對應的基因在DMT中高表達，每個基因對應的信噪比（Signal2noise，前面選擇的排序值計算方式）以灰色面積圖顯展示。

在上圖中，我們一般關注ES值，峰出現在排序基因集的前端還是後端（ES值大於0在前端，小於0在後端）以及Leading edge subset（即對富集貢獻最大的部分，領頭亞集）；在ES圖中出現領頭亞集的形狀，表明這個功能基因集在某處理條件下具有更顯著的生物學意義；對於分析結果中，我們一般認為|NES|>1，NOM p-val<0.05，FDR q-val<0.25的通路是顯著富集的。最後還有一個該GO基因集下每個基因的詳細統計信息表，RANK IN GENE LIST表示在排序好的基因集中所處的位置；RANK METRIC SCORE是基因排序評分，我們這裡是Signal2noise；RUNNING ES是分析過程中動態的ES值；CORE ENRICHMENT是對ES值有主要貢獻的基因，即Leading edge subset，在表中以綠色標記。

2. Dataset details

晶片原始數據和去重後的數據；如果分析的時候沒有用到晶片數據或沒涉及到名字轉換則前後基因數目一樣。

3. Gene set details

我們分析提供的gmt文件中有多個GO條目，每個GO條目裡又有多個基因；GSEA分析軟體會在每個GO條目中搜索表達數據集gct文件中的基因，並判斷有多少個在GO條目中；若經過篩選後保留在GO條目中的基因在15-500（閉區間）時該GO條目才被保留下來進行後續的分析。

此結果顯示我們從5917個GO條目中淘汰了了1964個GO，剩下3953個GO條目用作後續分析。

點擊gene sets used and their sizes可以下載詳細Excel表。

Excel第一列是GO名稱，第二列是GO條目中包含的基因數目，第三列是篩選後每個GO中還有多少基因屬於表達數據集文件中的基因，不滿足參數（15-500）的條目被拋棄，顯示為Rejected不納入後續分析。

備註： 此處的篩選範圍15-500是可調參數，在軟體的參數basic fields處的Max size和Min size處更改。

4. Gene markers for the NGT versus DMT comparison

這部分展示的是我們提供的表達數據集文件中的基因在兩個組別中的表達情況。

輸入的文件中總共有15056個基因，其中有7993個基因在正常人（NGT）中表達更高，佔總基因數的53.1%；有7063個基因在糖尿病患者（DMT）中表達更高，佔總基因數的46.9%。後面一個面積百分比，稍後看圖的時候再做解釋。

點擊rank ordered gene list可以下載一個排序好的基因集Excel表，排序原則是根據Basic fields參數設置處的Metric for ranking genes決定的。我們選的是信噪比（signal2noise），顯示在表格中的最後一列。根據NGT_vs_DMT評分得到一個降序排列的基因集，之後便可以做基因的富集分析了。

GSEA基因富集分析的原理就是基於該排列好的基因集，從第一個基因開始判斷該基因是否存在於經過篩選的GO功能基因集中，如果存在則加分，反之減分。所以評分過程是一個動態的過程，最終我們會得到一個評分峰值，那就是GO功能富集的評分。加分規則通過Basic fields參數設置處的Enrichment statistic決定的。

接著有一個分析的結果的熱圖和gene list相關性的圖。

熱圖中展示了分別在兩組處理中高表達的前50個基因，總共100個基因的表達情況。

gene list相關性圖如下。橫坐標是已經排序好的基因，縱坐標是signal2noise的值。虛線左側的基因是在NGT中高表達，右側的基因在DMT中高表達。這部分結果報告中的面積比就是基於該圖計算的，可以看出面積百分比和基因數目百分比有一定的差異，面積百分比可以從整體上反映組間信噪比的大小。

Butterfly plot顯示了基因等級與排名指標評分之間的正相關（左側）和負相關性（右側）。左側藍色虛線和右側紅色虛線是真實的信噪比結果，其他顏色的線是軟體對數據做了隨機重排後的結果。默認情況下，圖形只顯示前100個基因，也就是排名第一和最後100個基因。可以使用運行GSEA頁面上Advanced fields處的Number of markers來更改顯示的基因數量。

5. Global statistics and plots

這部分包含兩個圖：1) p值與歸一化富集分數(NES)的對比圖，這提供了一種快速、直觀的方法來掌握有意義的富集基因集的數量。2) 通過基因集的富集分數統計圖，提供了一種快速、直觀的方法來掌握富集的基因集的數量。

理解了上面各個部分的結果後，再回過頭看這張GSEA分析原理圖就簡單了。

Cytoscape富集網路可視化

在GSEA軟體的左側提供了Enrichment Map Visualization的功能，點擊後GSEA軟體會自動調用Cytoscape，建議等待Cytoscape啟動後再進行接下來的操作，且保證在分析過程中Cytoscape是處於開啟狀態。

選擇一個GSEA分析結果，點擊Load GSEA Results，其他項為默認值就行，點擊Build Enrichment Map以展示基因富集結果的網絡圖。（備註：GSEA分析結果用的是和上面演示數據不同的文件，可自行更改）

運行成功之後會彈出下面的提示框，結果直接展現在了Cytoscape中，如下圖所示：

Graphpad作圖比較多個ES

GSEA富集分析可視化結果是給每個功能基因集富集情況單獨出一張圖，有的時候我們想要比較基因集在兩個不同的GO中的富集情況，利用GSEA軟體分析得到的Excel結果表，提取有用的數據結果，在graphpad裡進行加工再出圖，可以達到我們想要的結果！

效果圖如下：

《Graphpad，經典繪圖工具初學初探》一文中介紹了graphpad入門的基礎知識，基本操作可以單擊回看。最近使用graphpad發現其多圖排版功能十分強大，不僅可以實現多個圖形排版還能實現圖層疊加。上面這個圖的作圖思路也就是把該圖拆分為兩部分，Enrichment score和基因位置分布條帶圖。

在GSEA分析結果文件夾裡隨便找一個感興趣的GO條目分析結果Excel表，作圖需要提取的信息即圖中標黃的部分，RANK IN GENE LIST和RUNNING ES。

加工一下已有數據，添加一列high取值都為0.1，設置高度，黃色部分的數據就是用來繪製基因位置分布條帶圖的；綠色部分用來繪製動態的ES評分曲線。

打開graphpad之後，我們在XY類圖下選擇Enter and plot a single Y value for each point，將兩部分數據分開粘貼到軟體不同數據表格中（如下圖左側所示），下圖中間展示兩個圖選擇的不同繪圖方式，調整參數後最終得到右側的結果。

在左側目錄樹處點擊layout創建一個圖形排版界面，將Graphs下的圖形複製粘貼到layout1下，拖拉移動位置很快就能將兩部分圖對齊。

之後用同樣地方式畫另外一個富集結果，粘貼到layout1中便得到最開始展示的圖。

注意：設置X軸的範圍是1到總排序基因數，Y軸是0到多個富集分析得分的最大值。

往期精品(點擊圖片直達文字對應教程)

後臺回復「生信寶典福利第一波」或點擊閱讀原文獲取教程合集