生信寶典之前總結了一篇關於GSEA富集分析的推文——《GSEA富集分析 - 界面操作》，介紹了GSEA的定義、GSEA原理、GSEA分析、Leading-edge分析等，不太了解的朋友可以點擊閱讀先理解下概念 （下面摘錄一部分）。
GSEA案例解析
介紹GSEA分析之前，我們先看一篇Cell文章(https://sci-hub.tw/10.1016/j.cell.2016.11.033)的一個插圖。
以下是文章原文對圖的注解：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 File Content Format Source
Expression dataset Contains 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 txt You create the file. 一般的基因表達矩陣整理下格式就可以。如果是其它類型數據或自己計算rank也可以，后面有更多示例。(如果后綴為txt格式，傳統的基因表達矩陣就可以，第一列為基因名字，名字與待分析的功能注釋數據集一致，同為GeneSymbol或EntrezID或其它自定義名字，第一行為標題行，含樣品信息。gct文件需要符合下面的格式要求。)
Phenotype labels Contains phenotype labels and associates each sample with a phenotype. cls You create the file or have GSEA create it for you. 一般是樣品分組信息或樣品屬性度量值或時間序列信息。
Gene sets Contains 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 gmt You 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 annotations Lists each probe on a DNA chip and its matching HUGO gene symbol. Optional for the gene set enrichment analysis. Chip You 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
1
2
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的話則表示運行成功，此過程大概需要十分鍾左右，視數據大小而定。
Command： 顯示運行這個分析的命令行，以后就可以批量運行類似分析了。
四、結果解讀
數據分析完后的結果會保存到我們設置的路徑下，點開文件夾中的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到多個富集分析得分的最大值。
GSEA分析自定義功能注釋集
軟件安裝和數據格式准備不清晰的請見前文。
下面以擬南芥的轉錄組數據為例介紹GSEA的使用。待分析的數據是擬南芥的兩組實驗 (mut_vs_wt)，每組各三個生物重復，故一共6個樣本。各個樣本有唯一的名稱：mut_1、mut_2、mut_3、wt_1、wt_2、wt_3。
准備數據
在分析之前需要准備三個文件：樣本分組文件 (Phenotype labels, cls格式)、基因表達文件 (Expression dataset)和基因集數據庫文件(gene sets database)
Phenotype labels file
本質上是一個文本文件，需按照以下格式書寫：
其中每兩個字符串之間均用空格隔開
第一行6 表示共有6個樣本；2 表示是兩個分組；1 是固定寫法，不要改。
第二行的#后面沒有空格，mut和wt為2個分組的名字，可以自行更改為自己試驗設定的分組名字，名字內部不只允許有字母、數字、下划線。
第三行的書寫方式需要格外注意，這里寫的是6個字符串(代表6個樣本)，樣本重復必須賦予相同的名字，以便GSEA判斷樣品所屬的組，並且其排列的相對順序必須與稍后要准備的基因表達文件(或基因表達矩陣)中對應的樣本書寫順序相同。
如果你使用windows來書寫這個文件，建議使用NotePadd++或UltraEdit等。不建議使用Excel或記事本(空格、換行及Tab鍵等格式不易把握)。寫完之后注意需要將文件后綴名保存為.cls。
如果你使用bash環境的vi或vim來寫這個文件(使用windows的請忽略此處)，那么保存為.cls文件后使用cat -A查看這個文件，輸出到屏幕上的應是 （$代表換行符，不需要自己書寫）：
6 2 1$
#mut wt$
mut mut mut wt wt wt$
1
2
3
基因表達文件(Expression dataset)
該文件建議使用Excel文件生成(測序公司一般都會提供這個文件)。只是需要按照下面幾項稍作修改：
第一行是固定寫法：#1.2
第二行25000表示此表格有25000個基因，注意25000不是這個表的總行數，而是從第四行開始數的基因數，並且無重復的基因名稱。6表是有6個樣本(與樣本設置文件中的6必須對應)。
第三行NAME和DESCRIPTION為固定寫法。DESCRIPTION列下面可以寫na或任意字符串。
第四行及以后的行中，基因的表達值最好是沒有取過對數的。比如TPM、RPKM/FPKM或DESeq2標准化后的值均可。
將此文件另存為文本文件(制表符分隔)(*.txt)，然后再將后綴名稱改為.gct即可
若在bash環境中生成該文件(使用windows的請忽略此處)，使用cat -A查看這個文件，輸出到屏幕上的應是：
#1.2$
25000^I6$
NAME^IDESCRIPTION^Imut_1^Imut_2^Imut_3^Iwt_1^Iwt_2^Iwt_3$
AT3G45264^Ina^I296472.8029^I298638.6905^I280473.9967^I385382.3187^I355358.3136^I377348.3648$
AT5G38382^Ina^I164223.8048^I169543.267^I164432.3066^I124346.784^I120357.6296^I123873.0918$
1
2
3
4
5
符號^I表示此處是一個Tab鍵。生成文件后將后綴名稱改為.gct
基因集數據庫文件(gene sets database)
如果研究的物種是人，在MSigDB數據庫中下載感興趣的基因集即可。地址：http://software.broadinstitute.org/gsea/downloads.jsp。
其中含有GO、KEGG、Reactome、hallmark、BioCarta及MSigDB匯總的基因集，其中GO基因集又分為mf、bp和cc，且基因名稱包含Entrez IDs和gene symbols等。因此在選擇的時候需要仔細甄別一下再下載，且注意保證注釋集中的基因名稱與表達矩陣中基因名稱類型一致。
基因集數據庫文件的格式：
第一列是某一個生物過程、細胞組分、分子功能或信號通路等等，是一個基因集(gene set)的名稱
第二列是該基因集的編號或其它描述信息均可
第三列及以后的每列是從屬於該基因集的基因，有的是一個，有的是兩個或多個。你也可以輸入自己感興趣的基因，只需滿足以上格式即可。
列與列之間以Tab鍵隔開，基因與基因之間也以Tab鍵隔開。
文件后綴名寫為.gmt
本文分析的對象是擬南芥，MSigDB並無提供，因此需要自己去GO數據庫或擬南芥信息資源網站(The Arabidopsis Information Resource, TAIR)上下載。再將其加工修改為一個gmt文件。以bash環境為例：
wget https://www.arabidopsis.org/download_files/GO_and_PO_Annotations/Gene_Ontology_Annotations/ATH_GO_GOSLIM.txt
# 提取biological process分類
# 注意按照自己的文件篩選合適的列，修改下面的$8,$1,$5,$6.
awk 'BEGIN{OFS=FS="\t"}{if($8=="P") print $1,$5,$6}' ATH_GO_GOSLIM.txt | sort -u | awk 'BEGIN{OFS=FS="\t"}{anno=$2"\t"$3; a[anno]=a[anno]==""?$1:a[anno]"\t"$1}END{for(i in a) print i,a[i]}' > ATH_GO_GOSLIM_LocusName_bp_useme.gmt
1
2
3
4
5
這個ATH_GO_GOSLIM_LocusName_bp_useme.gmt就是一個新建的基因集數據庫文件。注意其中的基因名稱需要與已經做好的基因表達文件中的基因名稱對應。
使用javaGSEA軟件開始分析
Java安裝或更新后即可雙擊gsea_4096m.jnlp啟動GSEA,如下圖：
點擊1處的Load data后，點擊2處的Browse for files；將已經准備好的樣本分組文件(.cls)、基因表達文件(.gct)和基因集數據庫文件(.gmt)一並上傳。注意：如果分析的是人的基因數據，且調用GSEA數據庫中的gmt文件，軟件會自動下載，此處可以無需上傳本地的gmt文件。
上傳后點擊3處的Run GSEA，彈出下圖。在4處的Expression dataset選擇上傳的基因表達文件(.gct文件)。在gene sets database右側點擊復選框，選擇Gene Matrix(local gmx/gmt)，選擇已上傳的基因集數據庫文件(.gmt)。然后回來，在Phenotype labels右側點擊復選框，選擇已上傳的樣本分組文件(.cls)並選擇一個phenotype。
圖6
在Collapse dataset to gene symbols處選擇false，這是因為在基因表達文件(.gct)和基因集數據庫文件(.gmt)中已經有對應好的基因名稱，不需要再做轉換。
在Permutation type處選擇gene_set。
在6處即可點擊Run來運行GSEA。
7處會顯示運行的狀態running，當出現success后點點擊 (等同於點擊Show result folder並打開文件夾中的index.html)，即可彈出分析報告，如下圖。
點擊圖中Detailed enrichment results in html format可查看排序好的基因集富集結果的表格(依照NES值排序)，並可點擊Details..查看詳情。
點擊Snapshot即可看到GSEA分析的結果圖 （后面解釋）。
在Snapshot頁面點擊圖片即可進入對應基因集的詳細頁面，其中的信息如GSEA Results Summary、Enrichment plot和GSEA details(可看到該基因集的基因數目)等。
圖中1說明描繪的是對核糖體生物合成(ribosome biogenesis)基因集的GSEA分析結果。
首先GSEA將所有基因按照與表型的相關度排序，獲得一個排序后的基因列表，排序值以熱圖形式展示 (圖中的4)。
然后在排序好的基因列表從前向后挨個查看每個基因是否出現在此圖對應的基因集中，若初選用豎線表示（圖中的3(共119個基因)），並且增加一個統計值( a running-sum statistic)，反之則減少一個統計值，最終每個位置的基因各自獲得一個累計值（圖中的綠色曲線）。
每個基因對應的累計值就叫做富集得分 (Enrichment score, ES) （圖8中的2），而這個基因集的富集得分 (ES)則定義為遍歷基因列表時遇到的離零的最大偏差，即峰值 (此圖中為0.597)。峰值為正值表示基因集富集在列表的頂部(mut)，負值表示富集在底部(wt)。
標准化富集評分(NES)是檢驗基因集富集結果的主要統計指標。GSEA富集分析時，不同的用戶所輸入的基因集數據庫文件中的基因集數目可能不同。富集得分的標准化考慮了基因集個數和大小的差異。因此NES可以用來比較不同基因集之間的分析結果。NES計算方法如下:
因此，GSEA對基因集富集分析的排序 (如Snapshot頁面或Detailed enrichment results in html format頁面)是依據NES的值，而不是ES或p-val等值。這也就是為什么文章開頭的那篇cell文章需要在兩個基因集富集分析結果中(圖1)明確標明NES值。
NOM p-value，即Nominal p value，是對富集得分(ES)的統計學分析，但未根據基因集的大小或多重假設檢驗來校正，因此在比較不同基因集的作用上有限。
FDR q-value，即假陽性率(False discovery rate)。是對標准化富集得分(NES)可能存在的假陽性結果的概率估計。因此FDR越小說明富集越可靠。FDR比NOM p-value更有參考意義 (圖1也只標了FDR)。
參考文獻
[1] Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES, Mesirov JP. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005 Oct 25;102(43):15545-50. Epub 2005 Sep 30.
[2] Ang YS, Rivas RN, Ribeiro AJS, Srivas R, Rivera J, Stone NR, Pratt K, Mohamed TMA, Fu JD4, Spencer CI, Tippens ND, Li M, Narasimha A, Radzinsky E, Moon-Grady AJ, Yu H, Pruitt BL, Snyder MP, Srivastava D. Disease Model of GATA4 Mutation Reveals Transcription Factor Cooperativity in Human Cardiogenesis. Cell. 2016 Dec 15;167(7):1734-1749.e22. doi: 10.1016/j.cell.2016.11.033.
————————————————
版權聲明：本文為CSDN博主「生信寶典」的原創文章，遵循CC 4.0 BY-SA版權協議，轉載請附上原文出處鏈接及本聲明。
原文鏈接：https://blog.csdn.net/qazplm12_3/article/details/83474140