GSEA縮寫
GSEA的全稱是Gene Set Enrichment Analysis，中文翻譯就是基因集富集分析，最早提出GSEA的文章是下面的這篇文獻：
基因集富集分析是分析基因表達信息的一種方法，富集是指將基因按照先驗知識，也就是基因組注釋信息進行分類（注：我個人理解先驗知識就是先于經(jīng)驗的知識，例如我們從書本上知道地球是圓的，這就是先驗知識，我自己沒有親自去測量，也無法觀察到，就知道地球是圓的?；蚪M注釋就是對堿基序列的注釋，給出一段堿基序列，根據(jù)前人研究的結(jié)果，我知識它是什么基因，發(fā)揮了什么功能，以上就是對先驗知識的理解）。最早提出GSEA的文章是下面的這篇文獻：
Subramanian, A., et al. (2005). "Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles." Proceedings of the National Academy of Sciences 102(43): 15545-15550.
常規(guī)分析的局限
使用芯片或測序技術(shù)已經(jīng)成了生物研究的常規(guī)手段，獲取某個樣本的整個表達譜已經(jīng)不再困難，但是，如何從這些表達譜中發(fā)掘有用的信息，以及對生物學(xué)現(xiàn)象進行解釋則是一個挑戰(zhàn)。在一個典型的研究中，例如我們研究某個抗腫瘤藥物對于兩種腫瘤細(xì)胞（一個對藥物敏感，一個不敏感）的效應(yīng)，然后檢驗這兩類腫瘤細(xì)胞的表達譜，可以獲得數(shù)千個基因的變化信息。我們根據(jù)這兩類細(xì)胞之間基因表達差異的情況，例如根據(jù)它們之間的差異倍數(shù)，對它們進行排序，得到一個列表L，就可以從中獲得一些信息。常規(guī)的做法就關(guān)注那些表達量差異極大的基因，也就是列表L頂部或底部的基因，用于發(fā)現(xiàn)這兩類細(xì)胞不同變化的一些線索，不過這種方法有一些局限。
- 經(jīng)過多重假設(shè)檢驗校正，或許沒有一個基因滿足統(tǒng)計學(xué)顯著性的閾值，因為相關(guān)的生物學(xué)差異相對于實驗誤差來說，并不明顯，也就是不好區(qū)分出基因的差異與實驗誤差。
- 也許會有很多差異基因，但是這些差異基因無法與生物學(xué)內(nèi)容整合起來，也就是說用這些差異基因不太好解釋生物學(xué)現(xiàn)象。
- 單基因分析方法也許會錯過許多對通路有重要影響的基因。細(xì)胞過程經(jīng)常會協(xié)同影響一批基因。例如，編碼一個代謝通路的一批基因的表達量都升高20%，就會造成非常大的影響，這種影響有可能比一個基因的表達量增加20倍都有效。
- 當(dāng)兩個不同的團隊對同一個生物學(xué)體系都進行研究時，他們得到的差異基因或許就不太一致。
GSEA就是為了解決上述的問題而提出來的。此方法可以在基因集（也就是多個基因）水平上來處理表達譜數(shù)據(jù)。這里所說的基因集（gene sets）是指基于當(dāng)前積累的知識，例如關(guān)于生物通路知識或以前得到的共表達數(shù)據(jù)來定義的一組基因。GSEA主要研究一個基因集S（這里的S指的是我們定義的基因集）出現(xiàn)在目標(biāo)基因列表L（也就是實驗得到的按表達倍數(shù)排序的基因列表L）的頂部或底部是否與實驗中的不同表型有關(guān)。
因此，我們就可以知道，GSEA分析有三個特點：①分析的基因集合而不是單個基因；②將基因與預(yù)定義的基因集進行比較；③富集分析。
與GO（Gene Ontology）和KEGG pathway分析相比，GSEA分析的主要優(yōu)勢在于：
一般的差異分析（GO和Pathway）往往側(cè)重于比較兩組間的基因表達差異，集中關(guān)注少數(shù)幾個顯著上調(diào)或下調(diào)的基因，這容易遺漏部分差異表達不顯著卻有重要生物學(xué)意義的基因，忽略一些基因的生物特性、基因調(diào)控網(wǎng)絡(luò)之間的關(guān)系及基因功能和意義等有價值的信息。而GSEA不需要指定明確的差異基因閾值，算法會根據(jù)實際數(shù)據(jù)的整體趨勢， 為研究者們提供了一種合理地解決目前芯片分析瓶頸問題的方法,即使在沒有先驗經(jīng)驗存在的情況下也能在表達譜整體層次上對數(shù)條基因進行分析,從而從數(shù)理統(tǒng)計上把表達譜芯片數(shù)據(jù)與生物學(xué)意義很好地銜接起來,使得研究者們能夠更輕松、更合理地解讀芯片或測序結(jié)果。
GSEA的原理
首先是對每個樣本里面的基因表達值在樣本內(nèi)部進行排序，這里可以理解為，實驗組與對照組的差異基因進行排序。但是，這個差異如何量化，有多種方法，可以是Signal2noise，或者是T檢驗的值，或者是fold change的值，或者是logFC等。而GSEA默認(rèn)的就是signal-to-noise metric來對基因進行排序。如果要想計算Signal2Noise ，每個group必須要有3個及以上的samples。
現(xiàn)在我們定義一下各種參數(shù)。
其中我們把自己測出來的差異基因的排序列表稱為目標(biāo)基因列表L，把根據(jù)先驗知識預(yù)先定義的基因集稱為功能基因集S（這個基因集可以是與免疫有關(guān)的，可以是與某個代謝通路有關(guān)的，也可以是與miRNA有關(guān)的，具體的看GSEA官網(wǎng)），把這個基因集中的成員稱為s（前者是大寫，后者是小寫）。GSEA的運行原理就是判斷功能基因集S里面的成員s在目標(biāo)基因列表L里面是隨機分布的，還是主要聚集在目標(biāo)基因列表L的頂部或底部。如果我們研究的功能基因集S的成員顯著聚集在目標(biāo)基因列表L的頂部或底部，那么就說明此基因集成員對實驗干預(yù)造成的結(jié)果有作用，就是我們要關(guān)注的基因集。
以下是GSEA計算中的幾個概念。
- Enrichment Score，即ES，中文翻譯為富集得分。它反應(yīng)的是基因集成員
s在目標(biāo)基因列表L端富集的程度，計算方法是，從目標(biāo)基因列表L的第一個基因開始，計算一個累計統(tǒng)計值。當(dāng)遇到一個落在功能基因集S里面的基因，則增加統(tǒng)計值。遇到一個不在功能基因集S里面的基因，則降低統(tǒng)計值。每一步統(tǒng)計值增加或減少的幅度與基因的表達變化程度是相關(guān)的。富集得分ES是從沒有遇到時候開始計算，直到最大值。正值ES表示基因集在列表的頂部富集，負(fù)值ES表示基因集在列表的底部富集，它的算法如下所示： - ES score的算法如下圖所示：
其中目標(biāo)基因列表L表示的就是差異基因已經(jīng)排序后的序列，S指的是功能基因集S，我們看左邊的Ranked List這個部分，它右邊是FC，即表達倍數(shù)（明顯能看出來是從大小到排列的），紅色表示的是Hit，也就是說這個基因在功能基因集S里面，如果在，就加分，也就是說在右邊的表格里的Contribution to running sum for ES加上相應(yīng)的分?jǐn)?shù)。如果是藍色，表示這個基因不在，就減去相應(yīng)的分?jǐn)?shù)。
評估富集得分(ES)的顯著性。顯著性是通過置換檢驗 (permutation test)的方法來進行檢驗的。具體過程是，我們轉(zhuǎn)換不同分?jǐn)?shù)下的數(shù)據(jù)，并且再一次計算ES值，使之形成一個新的ES分布（假分布），從經(jīng)驗上來講，交換之后，ES的p值相對于新的ES值（統(tǒng)計分布）來說若是顯著的變化，則有理由說明此基因集是有一定的生物學(xué)意義的。若樣品量少，也可基于基因集做置換檢驗 (permutation test)，計算p-value。
多重假設(shè)檢驗校正。當(dāng)評估了所有基因的數(shù)據(jù)之后，我們要使用多重假設(shè)檢驗來評估它們的顯著性。首先把每一個基因的ES值根據(jù)基因集的大小進行標(biāo)準(zhǔn)化，得到均一化后的ES值，即NES，全稱為Normalized Enrichment Score，之后針對NES計算假陽性發(fā)現(xiàn)率（FDR），并以此劃出假陽性部分對應(yīng)每一個NES值。FDR是評估一個NES表達值中所發(fā)現(xiàn)的假陽性可能性大??；它是由NES的觀測值和零分布時比較得出的。
Leading-edge subset，對富集得分貢獻最大的基因成員。
GSEA分析所需要的文件
GCT文件
GCT格式是由分割符分開的一個數(shù)據(jù)文件，它記錄的是基因表達的數(shù)據(jù)集，也就是我們要進行GSEA計算的原始數(shù)據(jù)，它的格式如下所示：
從左到右，順時針看到下，GCT文件包含以下這些信息：
-
Always "#1.2"：這個單元格里一直標(biāo)注的是#1.2，不用管它。 -
# of samples：這個是樣本的數(shù)目，在上面的表格中，這個數(shù)目是130。 -
Third column onwards are sample names. These must be UNIQUE，它是樣本的名稱，樣本的名稱必須保證是唯一的，不能出現(xiàn)重復(fù)。 -
Data stars on line：第4行，這是每個基因具體的數(shù)值。 -
Each column contains expression values from 1 sample. Missing values are allowed(leave empty)：基因表達的數(shù)值可以是空。 -
Column 2:Row descriptions.：第2列，對基因的描述，可以為空（填上na）。 -
Column 1:Rowidentifiers：這一列是基因的標(biāo)識，通常是探針集的id（probe set ids）或克隆號，此列數(shù)據(jù)必須唯一，不能出現(xiàn)重復(fù)。 -
A2單元格：The # of rows：探針的數(shù)目。
CLS文件
CLS的全稱是Categorical (e.g tumor vs normal) class file format，它表示的是分組的信息，下面是一個CLS文件的內(nèi)容：
以上面的圖片為例說明，第1行有3個數(shù)字，分別為38， 2， 1。其中38是樣本的數(shù)目，2是分組的數(shù)目，1不用管。第2行，分別是ALL和AML，這表示的是分組的信息，即有2個組，分別為ALL和AML。第4行，分別是0和1，0表示ALL，1表示AML，它們的數(shù)目表示各自的樣本有多少個。再看一個案例：
在這個案例中，我們可以看到一共有9個樣本，分了3組；這3組分別為KRAS_MUT，WT，MYC_MUT，第3行是具體的分組數(shù)目。
常用的定義基因集的文件主要有GMT與GMX。
GMT文件
GMT文件的全稱是Gene Matrix Transposed file format，中文是“基因矩陣轉(zhuǎn)換文件格式”，縮寫為GMT。GMT文件是由制表符分割的文件，它用于描述基因集，可以理解為高通量測序的注釋文件。下圖是GMT文件的一個案例，每行代表一個基因集：
在上面的圖片中，我們可以看到，第1行表示的一個是基因集，第1列表示的是是基因集的名稱，這一列不能出現(xiàn)重復(fù)的基因集。第2列是對基因集的一個簡單描述，可以選擇不填（寫上na），第1列的長度可以不等，如果要在Excel表格中編輯GMT文件，一定要注意Excel的自動糾正功能（例如，Excel有可能把SEP8基因自動糾正為8-Sep）。
GMX文件
GMX文件是另外一種描述基因集的文件格式，下圖是GMX的一個案例，它與GMT格式正好相反，它的每列代表一個基因，如下所示：
RNK文件
RNK文件的全稱是Ranked list file format，這個文件是一個有序的基因列表，第1列是基因名稱（A1單元格忽略），第2列是數(shù)值（可以是表達倍數(shù)差異，也可以使用類似于t檢驗進行的排序），如下所示：
GSEA使用案例
這個案例使用的是GESA官網(wǎng)的數(shù)據(jù)，數(shù)據(jù)是p53突變型與野生型的肝細(xì)胞系的表達譜，這里使用的是芯片的表達數(shù)據(jù)，芯片現(xiàn)在已經(jīng)不怎么用了，這里只是講一下GSEA官網(wǎng)這個軟件的使用。
第1案例
第一步，去官網(wǎng)地址下載GSEA軟件，如下所示：
下載后，打開，界面如下所示：
GSEA的官網(wǎng)提供了一些數(shù)據(jù)的案例，我們可以拿來練習(xí)一下，地址點擊這里，這些數(shù)據(jù)如下所示：
這里以第2行的p53為例說明一下，這一行中有三個數(shù)據(jù)文件，分別為P53-hug95av2.gct，P53_collapsed.gct和P53.cls，都把它們下載下來。其中P53-hug95av2.gct，P53_collapsed.gct這兩個文件都是基因表達文件，其中，帶有collapsed.gct的文件表示這個表達譜文件中的標(biāo)記符（通常是一串?dāng)?shù)字）已經(jīng)被基因名（symbol）替換掉了，另外一個文件是原始文件，如下所示：
使用Excel打開P53-hug95av2.gct這個數(shù)據(jù)文件（其實也可以使用Notepad++或PilotEdit Lite這類文本編輯器來打開gct文件），這個文件是分級以及表達量的數(shù)據(jù)，如下所示：
現(xiàn)在用Excel打開P53.cls這個文件，這個文件是注釋信息，如下所示：
從注釋信息中我們可以得知，有50個樣本，分為兩組，分別為2和1表示這兩組，2表示MUT組，1表示W(wǎng)T組。將上述的數(shù)據(jù)導(dǎo)入到GSEA中，過程如下所示：
打開Load data，如下所示：
此時右側(cè)出現(xiàn)新的頁畫，打開Browse for file，找到原始數(shù)據(jù)文件夾，如下所示：
導(dǎo)入要分析的文件，分別為cls，gct文件（gmt文件后面可以直接在軟件中下載），如下所示：
此時，單擊左側(cè)的Run GSEA，如下所示：
此時出現(xiàn)一些參數(shù)的選擇頁面，如下所示：
在第一行的Expression dataset中選擇表達的數(shù)據(jù)文件，即P53_collapsed_symbols.gct；
在第二行的Gene sets database中選擇基因集文件，即gmt文件，單擊右側(cè)的省略號，此時出現(xiàn)基因集文件的選擇框，如下所示：
根據(jù)官網(wǎng)的的信息，選擇c2.all.v6.2.symbols.gmt這個基因集，如下所示：
第三行，Number of permutations，此處設(shè)為10（我開始設(shè)了1000，沒跑出來）；
第四行，Phenotype labels，選擇MUT_versus_WT或WT_versus_MUT，順序不同會導(dǎo)致最終生成的ES圖不同，我開始的時候，設(shè)為MUT_versus_WT，發(fā)現(xiàn)ES得分的圖是倒的，不好看，現(xiàn)在就設(shè)為WT_versus_MUT；
第五行，Collapse dataset to gene symbols，此處填false，這個參數(shù)用于說明是否要把探針的編號轉(zhuǎn)換為基因ID（gene symbols），由于P53_collapsed.gct這個文件本身就已經(jīng)是使用基因ID了，因此不需要轉(zhuǎn)換，設(shè)為false；
第六行，Permutation type：填入phenotype。
第七行，Clip platfrom，芯片類型，前面Collapse dataset to gene symbols已經(jīng)是false了，因此這里不用選，再看第二個頁面，也就是Basic fields，填入要分析的樣本名稱，保存結(jié)果的文件夾，如下所示：
填完之后，單擊最下面的Run，如下所示：
[圖片上傳失敗...(image-895646-1551524317172)]
此時，GSEA開始運行，如下所示：
運行結(jié)束后，系統(tǒng)有提示，如下所示：
現(xiàn)在進入保存結(jié)果的文件夾，如下所示：
其中，我們選擇index這個html文件，就可以看到分析整體結(jié)果，如下所示：
打開，如下所示：
GSEA結(jié)果解讀
[圖片上傳失敗...(image-d98572-1551524317172)]
第1部分
現(xiàn)在分別查看這些結(jié)果的信息，先看第一部分：
在這一部分中，從上到下，可以看到：
431/899表示在WT這一分組中，一共有899個的功能基因集，其中421個上升；
99個基因集的FDR小于25%；
118個基因集的名義p值小于1%；
118個基因集的名義p值小于5%；
-
Snapshot of enrichment results中單擊
Snapshot可以看富集的結(jié)果，默認(rèn)情況下，只顯示前20個，如下所示：image現(xiàn)在單擊第1張圖，就可以查看富集分析的結(jié)果，如下所示：
image
再回到主頁面，如下所示：
單擊enrichment results in html就可以查看所有的富集分析結(jié)果，如下所示：
現(xiàn)在把標(biāo)題欄放大，如下所示：
其中，GS follow link to MSgnDB表示具體的功能基因集S的名稱；
GS DETAILS可以查看具體的某個功能基因集S的富集分析結(jié)果；
SIZE表示基因集里的基因數(shù)；
ES表示富集分?jǐn)?shù)的ES值；
NES表示校正后的ES值；
NOM p=val是名義P值；
FDR q-val用FDR法校正后的P值，即Qvalue；
FWER p-val用FWER法（Bonferonni校正）校正后的P值
RANK AT AMX表示ES值達到最大進度，對應(yīng)的通路基因的排名；
LEADING EDGE，顯示了用于定義Leading edge subset的一些參數(shù)，分別有Tags，List，Signal。
Tags指的是，在ES最大值這個點之前，功能基因集S中的成員s在目標(biāo)基因列表中的概率（此時ES為正值），或者說是在ES最大值這個點之后，功能基因集S中的成員s在目標(biāo)基因列表中的概率（此時ES為負(fù)值）。這個指標(biāo)看的是，對ES有貢獻的基因數(shù)目占功能基因集S中總數(shù)目的百分比。
List指的是，在ES最大值這個點之前，在排序的目標(biāo)基因列表L中的基因數(shù)目（此時ES值為正）占L總基因數(shù)目的百分比，或者說是在ES最大值這個點之后，在排序的目標(biāo)基因列表L中的基因數(shù)目（此時ES值為正）占L總基因數(shù)目的百分比，它反映的是，目標(biāo)基因列表L中多少個基因參與了ES的計算。
-
Signal指的是富集信號的強度，它由Tags和List構(gòu)成，如下所示：
其中，N是目標(biāo)基因集L中的基因數(shù)目，Nh指的是功能基因集S中的基因數(shù)目。如果功能基因集中的所有基因都位于目標(biāo)基因集L中的前列，那么信號強度就是最大，或者說是100%，如果功能基因集S中的基因都分散在整個目標(biāo)基因列表L中，那么信號強度就下降到了0%（這一處不太理解，Gene%代表的是什么文中沒有提，可能是List，原文如下所示：
where N is the number of genes in the list and Nh is the number of genes in the gene set.If the gene set is entirely within the first Nh positions in the list, then the signal strength is maximal or 100%. If the gene set is spread throughout the list, then the signal strength decreases towards 0%.
現(xiàn)在單擊一下GS DETAILS這一欄中的Details..，就會看到某個基因集的詳細(xì)結(jié)果，如下所示：
拉到中間區(qū)域，就會看到這樣的結(jié)果，即；
現(xiàn)在解釋一下這個列表：
PROBE，基因探針名；
DESCRIPTION，對這個基因的描述，如果沒有，則是na
GENE SYMBOL，基因名
GENE_TITLE，來源于芯片的注釋文件，如果不是芯片，就沒有，空白；
RANK IN GENE LIST，基因集S的成員在目標(biāo)基因列表L中的位置；
RANK METRIC SCORE，基因集S的成員在目標(biāo)基因列表L位于中的打分情況，貌似最高位于的分?jǐn)?shù)最高，然后遞減；
-
RUNNING ES，運行的富集分?jǐn)?shù)，也就是在排序基因列表L中某個點上的富集分?jǐn)?shù)，可以看到最高點的分?jǐn)?shù)是0.7459，這個可能就是ES的得分，如下所示：
image CORE ENRICHMENT，這一欄中，有的標(biāo)記為Yes，背景是淺綠色，這些標(biāo)記為Yes的基因，則是功能基因集S中的Leading-edge subset。它們對ES的貢獻最大。
如果在GSEA軟件中運行Leading Edge Analysis，就會出現(xiàn)下面的結(jié)果，如下所示：
其中，先看第1張圖，如下所示：
這是一張熱圖，它顯示了leading edge subset中的基因情況，在熱圖中，基因的表達值按5種顏色進行區(qū)分，分別為red, pink, light blue, dark blue，它們代表基要因的表達水平為high, moderate,low, lowest。
再看第2張圖，即右上，這是Set-to-Set圖，如下所示：
這些顏色塊顯示了subset之間的重疊情況，如果顏色塊越深（仔細(xì)看，有些顏色非常淺），表明subsets之間的重疊性越高，一個小方格中的A與B的密度對應(yīng)了X/Y的比例，其中X是A中l(wèi)eading edge genes的數(shù)目，Y表示，A與B的并集。顏色越深，表明，A與B含有的相同的leading edge gens越多，如果方格是白色的，表明A與B不含有相同的leading edge genes。
再看第3張圖，即左下的圖，如下所示：
這張圖的最底部是基因的名稱，縱坐標(biāo)表示，含有這個基因的功能基因集S的數(shù)目。
再看第4張圖，也就是右下圖，如下所示：
這個圖是直方圖，其中橫坐標(biāo)Jacquard表示，每對leading edge subset的交集除以并集（可以看到，它就是個小數(shù)，畢竟交集的數(shù)目要小于并集的數(shù)目）。Number of Occurrences表示在一個指定的bin內(nèi)，配對的leading edge subset數(shù)目，在這個案例中，大多數(shù)的subset配對的Jaccard結(jié)果都位于0到0.02之間，可以說，大部分的配對的Subsets都沒有重合的部分。
下面是原文：
The last plot is a histogram, where the Jacquard is the intersection divided by the union for a pair of leading edge subsets. Number of Occurrences is the number of leading edge subset pairs in a particular bin. In this example, most subset pairs have no overlap (Jacquard = 0).
GSEA的主要觀察指標(biāo)
GSEA主要看四個指標(biāo)，分別為ES，NES，名義p值，F(xiàn)DR值，現(xiàn)在分別解釋一下。
- ES，即Enrichment Score，富集分?jǐn)?shù)，看一張典型的GSEA分析圖，如下所示：
- ES，全稱為Enrichment Score，即富集分?jǐn)?shù)，一開始都是0，就是綠色長蛇狀的粗線，這個綠色的曲線就是ES的值計算過程，從左到右每遇到一個基因，就計算出一個ES值，連成線，最高峰為富集得分（ES）。
- 中間的長得像條形碼的圖形，圖片中長得像條形碼的東西就是功能基因集S在目標(biāo)基因列表L中所在的位置。
- 紅色直線是綠色曲線的頂點，它左側(cè)的數(shù)據(jù)是Leading edge subs，et，也就是對富集得分貢獻最大的基因成員，若富集得分為正值，則則表示這基因集位于L的頂部，如果富集得分為負(fù)值，則位于L的底部。
- 當(dāng)用功能基因集S從上到下，遍歷排序好的目標(biāo)基因列表L時，此時最下面的灰色區(qū)域就是不同基因的排序結(jié)果，它與分組情況相關(guān)，排序的結(jié)果從正值到負(fù)值進行排列，正值是與第1個分組有關(guān)（例如WT組），負(fù)值與第2個分組有關(guān)（例如MUT組），它是按從高到低進行排序的基因次數(shù)，排序用的是Signal2noise。
-
NES，即Normalized Enrichment Score。
這是校正后的ES數(shù)值。由于ES的計算原理是，一個
目標(biāo)基因列表L中的基因是否在一個功能基因集S中出現(xiàn)來計算的，但是各個功能基因集S中所包括的基因數(shù)不同，而且不同的功能基因集S與要分析的目標(biāo)基因列表L的相關(guān)性也不同，置換次數(shù)也不同，因此在比較數(shù)據(jù)集在不同功能基因集中的富集程度時，需要對ES進行標(biāo)準(zhǔn)化處理。
通過校正原始ES值，NES可以解釋功能基因集S大小的差異，以及功能基因集S和目標(biāo)基因列表L之間的相關(guān)性，因此NES可以用于比較基因集之間的結(jié)果，GSEA定義的NES公式如下所示：
NES基于所有數(shù)據(jù)集置換(permutation)的基因集富集分?jǐn)?shù)進行計算的，因此通過改變置換方法，以及置換數(shù)目，或者是基因表達數(shù)據(jù)量的大小，都能影響NES。例如，現(xiàn)在我們看2個分析：
第一，你分析一個表達數(shù)據(jù)集，GSEA會產(chǎn)生一個有序列表，并對這個列表進行分析；
第二，你使用GSEA對上面這個有序列表重排序，并分析。
如果你使用相同的參數(shù)設(shè)置（這里的相同是指相同的功能基因集S，相同的目標(biāo)基因列表L），那么你得到的的ES就是相同的；但是如果你在兩次分析時，設(shè)定了不同的置換次數(shù)（permuation）的話，那么校正后的富集分?jǐn)?shù)，即NES則不同。這兩次分析能反映出用于置換的數(shù)據(jù)集的很大不同（表達數(shù)據(jù)集與排序的基因列表），如下所示（第1次的置換次數(shù)是149，第2次是10）：
- FDR，全稱為false discovery rate，即錯誤發(fā)現(xiàn)率。
FDR描述的是一個估計的可能性，即當(dāng)一個功能基因集的NES值確定后，判斷其中可能包含的錯誤的陽性發(fā)現(xiàn)率，例如FDR=25%意味著，對此NES值的確定，4次中可能有1次是錯的，在GSEA的結(jié)果報告中，高亮顯示了FDR值小于25%的富集功能基因集，因為從這些富集功能基因集中最可能產(chǎn)生有意義的科學(xué)假設(shè)，方便進一步深入研究，在多數(shù)情況下，選擇FDR為25%來作為閾值，因為通常用于分析的芯片表達數(shù)據(jù)之間，大部分都缺乏一致性，而且每次分析的功能基因集數(shù)據(jù)不多，但是當(dāng)分析的芯片數(shù)據(jù)集較小，分析時選擇的是探針之間的隨機組合，而不是表型間的隨機組合，P值采用的嚴(yán)格度不高時，應(yīng)該選擇更加嚴(yán)格的FDR閾值，例如FDR=5%。一般而言，NES的絕對值越大，F(xiàn)DR的值就越小，說明寶座度越高的功能基因集，分析的結(jié)果可信度就越高。
- 名義P值，全稱為Nominal P Value。
名義p值描述提針對某一個功能基因集S得到的富集得分的統(tǒng)計顯著性，p值越小，越能說明基因的富集性越好。但在GSEA分析結(jié)果的四個參數(shù)中，只有FDR值進行了功能基因集S大小和多重假設(shè)檢驗的校正，而p值沒有，因此，當(dāng)分析結(jié)果中出現(xiàn)一個高度富集的功能基因集S，具有很小的名義p值，卻有很大的FDR值時，往往意味著其實和其它功能基因子相比較，它的富集并不是很顯著，原因可能是用于分析的芯片數(shù)據(jù)樣本量較少，或者是雜交信號微弱，又或是選擇的功能基因子集并沒有很好地反映樣本的物理學(xué)意義。在GSEA的報告中，當(dāng)p=0時，這就意味著實際的p值小于1/隨機組合的次數(shù)，例如，當(dāng)選擇的隨機組合數(shù)為100時，報告中的p=0實際上是說p值小于0.01，所以提高隨機數(shù)，就能得到更精確的p值，通常情況下，運行GSEA的隨機組合次數(shù)都會選擇1000，但不要超過1000（個人電腦一般選10），否則系統(tǒng)內(nèi)存會被耗盡。
通常情況下，我們一般認(rèn)為|NES|>1，NOM p-val<0.05，F(xiàn)DR q-val<0.25的通路下的基因集合是有意義的。
現(xiàn)在再回到主頁面，如下所示：
Detailed enrichment results in excel format (tab delimited text)這一行是使用Excel表格呈現(xiàn)的富集分析結(jié)果；
最后一行，即Guide to interpret results，則是對結(jié)果的說明。
第2部分
再看結(jié)果的第二部分，即如下所示：
這是對MUT組的富集分析結(jié)果，內(nèi)容與WT組的差不多。
第3部分
現(xiàn)在看第3部分，如下所示：
這里的信息提示，一共分析了10100個基因；
第4部分
看第4部分，如下所示：
這是是功能基因集S的一些信息，第一行提示，根據(jù)min=15,max=500的標(biāo)準(zhǔn)來篩選，1329個基因集中，不能用的有430個，能用的有899個；
第二行信息說明，能用的基因集是899個（WT組有421個，MUT組有478個）；
第三行說明，基因集的內(nèi)容，也就是不同的功能基因集中有多少個基因。
第5部分
再看第5部分，如下所示：
從這一部分中，我們可以看到這些信息：
第1行：一共有10100個基因；
第2行：WT組的基因是5707個，占56.5%，它所點相關(guān)關(guān)系面積的59.8%（這個59.8%我的理解就是在ES富集分析圖中，最下面的那個灰色區(qū)域）；
第3行：MUT組的基因是4393個，占43.5%，它所點相關(guān)關(guān)系面積的40.2%；
第4行：目標(biāo)基因列表的排序，是個Excel表格，打開，如下所示：
一共是10100個，最后一部分如下所示：
第5行：熱圖和相關(guān)關(guān)系的示意圖（只顯示了前50個基因），其中，熱圖用5種顏色來表示基因表達水平的高低，即red, pink,light blue, dark blue分別表示high, moderate, low, lowest，如下所示：
（熱圖內(nèi)容：Heat Map of the top 50 features for each phenotype in P53_collapsed_symbols.P53.cls#WT_versus_MUT ）
（圖片內(nèi)容：Fig 2: Ranked Gene List Correlation Profile Ranked list correlations for P53_collapsed_symbols.P53.cls#WT_versus_MUT ）
此外，在GSEA的使用的手冊上發(fā)現(xiàn)，還有一種蝴蝶圖(Butterfly plot)不過，在我的這個分析中，貌似沒有找到。蝴蝶圖顯示的是基因順序和排序度量得分（ranking metric score）之間的正相關(guān)和負(fù)相關(guān)的關(guān)系。默認(rèn)情況下，蝴蝶圖只顯示排名前100的基因，這也就是說，它只顯示目標(biāo)基因列表的前100個和最后100個基因。不過，可以在Run GSEA Page這個頁面中設(shè)置Number of markers參數(shù)來設(shè)置顯示的基因數(shù)。富集圖下面的部分顯示了觀察到的基因排列和ranking metric score得分之間的相關(guān)關(guān)系，蝴蝶圖顯示了觀察到的相關(guān)關(guān)系，以及置換(1%,5%,50%)正相關(guān)，負(fù)暗相關(guān)的關(guān)系。蝴蝶圖提供了一些哪些功能基因集的置換能夠改變基因排序和ranking metric score的信息，如下所示（這張圖是從使用手冊中截取的，不清楚）：
第6部分
再看第6部分，如下所示：
第1行表示，p值與NES的圖，如下所示：
第2行顯示了全局的ES值，如下所示：
第7-8部分
剩下的內(nèi)容如下所示：
Other則是這次分析的一些參數(shù)，Comments作為隨機種子的時間戳。
參考資料
6.Quick Tour of the GSEA Java Desktop Application
10.Quick Tour of the GSEA Java Desktop Application
11.馮春瓊, 鄒亞光, 周其趙,等. GSEA在全基因組表達譜芯片數(shù)據(jù)分析中的應(yīng)用[J]. 現(xiàn)代生物醫(yī)學(xué)進展, 2009, 9(13):2553-2557.
