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