參考資料來源：
http://software.broadinstitute.org/gsea/doc/GSEAUserGuideFrame.html
https://cloud.tencent.com/developer/article/1426130
GSEA圖的理解：
含義概述：
豎線表示：某個GO過程（比如：脂肪代謝）的基因。
我測了10000個基因在treat組和control組中的表達量。那么，與脂肪代謝相關的基因，在這兩個組中是有差異表達的嗎？
通過此圖，可以知道這個信息。
豎線的基因集中在左端，說明這些基因在treat組中高表達。
得出結論：與脂肪代謝相關的基因，在treat組和control組中，是差異表達的。這些基因在treat組高表達。
若黑線不聚集在某端，而是散落地排列。說明：與脂肪代謝相關的基因，在treat組和control組的表達沒有差異。
含義詳述：
給定一個排序的基因表L和一個預先定義的基因集S (比如編碼某個代謝通路的產物的基因, 基因組上物理位置相近的基因，或同一GO注釋下的基因，即MSigDB數據庫中的基因)，GSEA的目的是判斷S里面的成員s在L里面是隨機分布還是主要聚集在L的頂部或底部。這些基因排序的依據是其在不同表型狀態下的表達差異，若研究的基因集S的成員顯著聚集在L的頂部或底部，則說明此基因集成員對表型的差異有貢獻，也是我們關注的基因集。
比如：根據數據庫，上圖的GO過程有一個先驗的gene集（即S），我也有一個自己的50000+個排序基因集（即L）中，此GO過程中有多少個基因落入我的基因L中呢？答：圖中間的黑色條帶表示GO過程基因落入排序基因集L中的基因。這些基因在ctrl中聚集還是sd中聚集呢？答：在ctrl中聚集。這些黑條帶基因是隨機分布在50000+的基因集中，還是顯著聚集在50000+基因集的頂部或底部？答：聚集在頂部。
1. 圖中下面灰色圖的含義：http://software.broadinstitute.org/gsea/doc/GSEAUserGuideFrame.html?_Interpreting_GSEA_Results。所有的灰色圖都是一樣的，即結果文件ranked_gene_list_treat_versus_control_1572940764277.xls中的值。
圖中的橫坐標表示你輸入的gene列表。
1）縱坐標含義：縱坐標表示“Ranked list metric”，GSEA官網對此的解釋是ranking metric，即signal-to-noise ratio。
2）計算方法：參考下面的Signal2noise的解釋。
3）ranking metric的含義：它表示基因與表型(phenotype)的關系。即：這個基因與treat相關，還是與control相關，相關的度量值是多少。即ranking metric的含義。
positive ranking metric表示與phenotype1相關，negative ranking metric表示與phenotype2相關。
表型是什么？比如，我們對給葯組(treat)和正常組(control)做RNA-seq，得到每個基因在給葯組和正常組中的表達值。treat和control就是表型。
主要概念 ：
Signal2noise：信噪比。信號強度/噪聲強度。值大於1，表示信號強於噪聲；值小於1，表示信號弱於噪聲。舉例：gene A在treat組和和control組的表達量的比值。上圖中，對該值取log，值大於0，表示在前表型中表達量高；值小於0，表示在后表型中表達量高。即：result文件中“ranked_gene_list_treat_versus_control_1572940764277.xls”的含義。
ES（enrichment score，富集得分）：它是一個累積值。表示GO過程的基因集S在排序基因集L中的富集程度。（之前，我總是想不過來這一點。因為我想的是：排序基因集L中有多少個基因落入GO過程的基因集中呢？其實，應該反着想：GO過程的基因有多少個落入排序基因集L中呢？是散落在兩種phenotype之間，還是富集在一種phenotype中呢？）正值ES表示該GO過程的基因集在前表型（phenotype1）中富集，負值ES表示基因集在后表型（phenotype2）中富集。
計算方法：首先，對基因集L按照Signa2noise排序（除了Signal2noise，GSEA也有其它的計算方法），得到排序基因集L。
然后，對L中的每個基因，根據這兩個條件：基因是否出現在GO基因集S中；該基因的Signal2noise值。計算得到L中每個基因的Running ES值。（這兩個條件體現了基因集L與phenotype的關聯度）
最后，對L中的每個基因，計算累積值。即得到上圖的第一個圖。
在上圖第一個圖中，ES的最高或最低值，作為該GO過程的ES值。
NES（Normalized Enrichment Score）：它用來比較各個GO過程之間的ES值。上面的ES值是一個GO過程內部的每個Gene的ES值。而，對於每個GO過程都有一個NES值。
FDR ：q-value。gene sets的錯誤概率。（理解：該GO過程包含這些基因的錯誤概率）。對NES值的估計。FDR=25%表示在4次驗證中，有3次是成功的。如果樣本數少，用gene_set作為permutation，那么，需要更樣的FDR，比如<5%.
gene set：某GO過程的基因集。
p-value：對每個gene set都有一個p-value。它表示某個gene-set的ES值的顯著性。
置換檢驗結合GSEA說明（參考”置換檢驗—結合GSEA解釋“的文章）：
隨機抽取次數（比如上文的1000），即是GSEA中提到的number of permutation。次數越多，得到的統計量的值越多，抽樣分布的值越多，抽樣分布的分布圖更接近真實值。看真實數據在此抽樣分布圖上的分布時，就更准確。
P-value：
由以上可得知P-value的含義。拒絕零假設的概率、上圖的末梢的面積、”GO過程在排序基因集L中富集“這一結論錯誤的概率。
通常，我們取P值的閥值為0.05。即：P值小於0.05的GO過程都認為在排序基因集中是富集的。
那么，既然錯誤率是0.05，也就是說：檢測100次，有5次是出錯的。如果檢測1萬次，有500次是出錯的。也就是說，當檢測次數很多時，出錯的情況就會多。結合GSEA：GSEA給出的GO過程（P值小於0.01）中，100個GO過程中有5個是出錯的。即：100個GO過程中，有5個GO過程在兩種phenotype中是不富集的。10000個GO過程中，有500個是不富集的。當GSEA給出富集的GO過程很多時，假陽性的GO過程就很多。
適用場景：
例子：cell，2016年發表的文章(https://sci-hub.tw/10.1016/j.cell.2016.11.033)，使用了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).
我理解的GSEA分析所適用的場景之一：對特定的基因集，在兩種phenotype中的富集程度。
舉例說明：有兩種類型的細胞，treated和control。想知道哪種功能的基因集在這兩種類型的細胞上有差異。比如：與髓鞘形成相關的一個功能，這個功能對應一個基因集（即：GSEA數據庫中的基因集），想知道與髓鞘形成相關的基因在treat組上表達，還是在control組上表達。要實現這個功能，需要知道這些數據集：treat和control組的與髓鞘形成相關基因的表達值（RNA-seq得到）；與髓鞘形成相關的基因都有哪些（GSEA MSigDB的GO數據庫）。
我理解的GSEA分析所適用的場景之二：對於兩種phenotype，哪些gene set在這兩種phenotype中有差異。
參考資料：
原文發布於微信公眾號 - 生信寶典（Bio_data）
原文發表時間：2019-04-26
工具下載：
下載網址：http://software.broadinstitute.org/gsea/downloads.jsp
命令：
執行命令：./gsea.sh ，打開了界面版的GSEA工具。
簡要步驟：
Load data，需要load：expression文件、phenotype文件、下載的MSigDB的文件。
設置Required fields的參數，點擊最下面的Run。
詳細步驟：
使用：
准備輸入文件：
1. Expression dataset：是RNA-seq(或chipseq)得到的數據。
通常有這么幾列：
例1：
NAME DESCRIPTION AML1 AML2 AML3 ALL1 ALL2 ALL3
TP53 na 681.3 638.0 665.0 240.0 587.0 737.0
例2：
NAME DESCRIPTION G C
TP53 na 681.3 638.0
注意：
1）必須用tab分隔；
2）第一行的名稱必須是NAME和DESCRIPTION；
3）文件前兩行必須是：
#1.2（固定的內容）
53796（gene個數，即行數） 2（sample個數，不是phenotype數。如例1中phenotype數是2，sample數是6） 參考網址：http://software.broadinstitute.org/cancer/software/gsea/wiki/index.php/Data_formats
2. Phenotype labels文件：共三行。
2 2 1 分別為：sample數，phenotype數，最后的1是固定的
# G C phenotype的名稱
G C phenotype的名稱
以上兩個文件是需要自己准備的，是自己的的數據。
3. Gene sets database選項：
可以用GSEA ftp網站上的，也可以從GSEA 下載MSigDB的文件，load到工具中使用。
4. Collapse/Remap to gene symbols選項：
我選擇的是“No_Collapse”。
5.Permutation type選項：
GSEA文檔中，給出的說明是 ：樣本數大於7，建議用“phenotype”選項；小於7個樣本的數據，建議用“gene_set” 的選項。
6. Number of permutations選項：
GSEA官網建議1000。剛開始運行時，建議10；運行成功后，再設置為1000。
術語說明：
GSEA文檔中的Phenotype、sample、gene sets等的含義。
Phenotype：
sample：
gene sets：
問題：
1.error：“Read timed out！”
這個錯誤，gsea官網沒有說明。自己查的解決方法，費勁...。
原因：Gene sets database中，我選的是website的dataset，比如：ftp/...。這樣，程序運行時，去GSEA的ftp服務器尋找並使用該文件。我懷疑可能網絡不好，導致讀取超時。
解決：下載gene set到本地。其實，就是把MSigDB的數據文件（.gmt文件）下載到本地。下載地址：http://software.broadinstitute.org/gsea/downloads.jsp。
注意：load本地的gene set文件時，只能通過“Load Data”功能執行。我在gene set database打開的界面找了好久都沒找到load文件的圖標，最后，靈機一動，才想到如何load。浪費時間%>_<%
2. error：After pruning, none of the gene sets passed size thresholds.
這個錯誤，gsea官網有說明。
原因：下載的c5.all.v7.0.symbols.gmt文件中的gene symbol全是大寫，而我提供的expression文件中的基因名是小寫。
解決：將expression文件的基因名改為大寫即可。
具體解決方法：awk '{print toupper($3)}’，使用shell腳本的toupper函數即可。
3.運行時間：
我輸入兩個sample，兩種phenotype，共53000+個gene的數據，permutation選擇默認值1000，選擇一個bp.gmtwenjian，運行時間長，大約十幾分鍾。
