寫在前面的話：GSEA對(duì)于生物信息學(xué)技能的要求不高，以至于不需要太多的生信知識(shí)就能學(xué)會(huì)，這也成為了生物研究生的必備技能，如果你還不會(huì)的話，這篇教程將讓你100%學(xué)會(huì)，只要你從頭到尾跟著我的步驟走。
我按照網(wǎng)上的教程學(xué)習(xí)GSEA的時(shí)候繞了不少?gòu)澛?，參考了很多不同的教程，并且不斷地Debug，花了一周時(shí)間才跑出結(jié)果。之后我再做GSEA就如行云流水般順暢。這篇教程，保證看完包會(huì)，文末我甚至給出了可以直接拿來(lái)跑的文件，不耽誤大家時(shí)間，100%能重復(fù)出結(jié)果！
下面進(jìn)入正題：
GSEA（Gene Set Enrichment Analysis）,基因集富集分析，與GO分析類似，是一種基因差異表達(dá)分析。但是在原理上，與GO分析有很大區(qū)別。GO分析是針對(duì)差異基因，即看差異表達(dá)的基因在哪些GO分類或者pathway中富集，從而說明這些差異基因影響了哪些pathway。GSEA分析是針對(duì)所有基因的表達(dá)（不能用差異基因來(lái)跑GSEA,否則GSEA就失去了意義），這些基因的表達(dá)矩陣稱為dataset。dataset一般是經(jīng)過處理的RNA測(cè)序結(jié)果的基因表達(dá)矩陣，或者是從數(shù)據(jù)庫(kù)中下載得到的表達(dá)矩陣。這些表達(dá)矩陣必須手動(dòng)分組，比如分為control組與treat組。然后GSEA會(huì)自動(dòng)分析出control組與treat組之間每個(gè)基因的差異并按照差異大小給每個(gè)基因算分，然后GSEA會(huì)根據(jù)這些基因的分值以及他們?cè)谀承┗蚣╣eneset)中的富集程度來(lái)計(jì)算geneset富集分?jǐn)?shù)（EnrichmentScore)并繪圖（geneset由GSEA軟件自帶也可自行到數(shù)據(jù)庫(kù)下載（后面由鏈接）），這樣就得到了geneset的富集信息，從而推斷出control組與treat組究竟在哪些geneset所代表的通路中存在差異。以上為GSEA簡(jiǎn)單的原理，算法什么的我不太懂，也不需要去了解，只要會(huì)用GSEA進(jìn)行數(shù)據(jù)分析就可以了。
一. GSEA的前期準(zhǔn)備
1.電腦一臺(tái)，華為matebook，win10，8G內(nèi)存。
2.GSEA軟件：
http://software.broadinstitute.org/gsea/downloads.jsp
右邊選擇內(nèi)存大小，最好選小于電腦內(nèi)存的最大內(nèi)存，比如我電腦8G內(nèi)存，我就選擇4GB(for 64-bit Java only)。內(nèi)存大小決定了GSEA能處理的數(shù)據(jù)大小以及計(jì)算的精細(xì)程度。這里有兩個(gè)坑我曾經(jīng)栽過——請(qǐng)注意：第一是Java環(huán)境是必須的，請(qǐng)到j(luò)ava官網(wǎng)下載64位java8?（https://www.java.com/zh_CN/，官網(wǎng)只會(huì)推薦最適合你系統(tǒng)的java，沒有其他版本的下載選項(xiàng)，但這不一定是我們所需要的版本！）或者用我文末分享的網(wǎng)盤鏈接進(jìn)行下載，GSEA軟件不支持java9或更高版本。第二是：如果選擇2GB以上內(nèi)存的java需要64位系統(tǒng)，如果你的電腦是32位系統(tǒng)，請(qǐng)選擇1GB，win7-8系統(tǒng)一般會(huì)安裝64與32兩個(gè)系統(tǒng)，而win10不再有64與32位的區(qū)別，win10用戶可放心下載。點(diǎn)launch進(jìn)行下載，下載完成后可直接運(yùn)行，如無(wú)法運(yùn)行或沒有相應(yīng)請(qǐng)檢查java是否安裝正確。
下載完成后的圖標(biāo)為下圖所示，如果圖標(biāo)未顯示，請(qǐng)檢查java是否安裝正確。
這樣，軟件就下載完成了。
3. 實(shí)驗(yàn)數(shù)據(jù)（dataset）的準(zhǔn)備。
NCBI的GEO數(shù)據(jù)庫(kù)收錄的絕大多數(shù)已發(fā)表文章的測(cè)序數(shù)據(jù)，而且下載起來(lái)很方便。本次使用的是人肝臟單細(xì)胞測(cè)序的數(shù)據(jù)，GEO編號(hào)GSE115469。
由于該測(cè)序數(shù)據(jù)集太大，EXECEL無(wú)法打開，需要用R語(yǔ)言進(jìn)行整理，所以我將樣本數(shù)縮小為100個(gè)，基因縮小為5000個(gè)，重新構(gòu)建了一個(gè)表達(dá)矩陣，可以直接在EXCEL里進(jìn)行排版（注意啊，我是方便教大家用excel操作才這樣搞的，實(shí)際操作千萬(wàn)不能隨意縮小數(shù)據(jù)啊，如果有幾萬(wàn)基因幾十個(gè)樣本就用R處理）。這個(gè)矩陣我放在網(wǎng)盤里，文件名“human liver dataset”，文末有下載鏈接。打開之后如下圖所示，第一列是基因名，第一行是樣本名。這種基因名（ID)的專業(yè)術(shù)語(yǔ)叫做Gene symbol，我們平時(shí)接觸到的Actin等基因名都是Gene symbol，此外還有很多其他的基因名，比如Ensemble ID，不同的測(cè)序平臺(tái)對(duì)基因的編號(hào)也不同，這里不詳細(xì)討論Gene ID的轉(zhuǎn)換問題，本篇主要帶大家重復(fù)一遍GSEA的流程。下圖是normalize后的單細(xì)胞測(cè)序表達(dá)矩陣，可用作GSEA的dataset。事實(shí)上，如果不是單細(xì)胞的數(shù)據(jù)，普通轉(zhuǎn)錄組測(cè)序的數(shù)據(jù)樣本數(shù)一般只有10個(gè)以內(nèi)，基因數(shù)至少2萬(wàn)個(gè)，完全可以用EXCEL進(jìn)行處理。
對(duì)dataset的下一步處理是分組。我們這次就研究一下APOA2基因高表達(dá)的細(xì)胞與低表達(dá)的細(xì)胞有什么差異吧，首先我們將數(shù)據(jù)按照APOA2的表達(dá)水平排序。搜索到APOA2，在1454行
選擇除基因名外的所有列，然后點(diǎn)擊工具欄的“排序”
在新的對(duì)話框中，點(diǎn)擊選項(xiàng)，選擇按行排序，確定。
選擇主要關(guān)鍵字“更多行”
點(diǎn)擊1454行標(biāo)，點(diǎn)擊確定。
然后升序排列，點(diǎn)擊確定。
排序后的APOA2
我們?cè)谛闹邪亚?0列定義為APOA2low，把后10列定義為APOA2high，保存EXCEL。
下一步，將.csv格式的ECXEL轉(zhuǎn)換為.gct格式——可被GSEA軟件識(shí)別的格式。這里使用既可以繪制熱圖又可以轉(zhuǎn)換格式的在線軟件：Morpheus。
https://software.broadinstitute.org/morpheus/
界面如下
然后上傳上文處理好的表達(dá)矩陣。
點(diǎn)OK，會(huì)出現(xiàn)一個(gè)熱圖，無(wú)視這個(gè)熱圖，對(duì)我們沒有幫助，我們是利用這個(gè)軟件轉(zhuǎn)換格式而不是畫熱圖。然后點(diǎn)擊工具欄File-Save dataset
命名文件，切記選擇GCT1.2,?點(diǎn)OK，然后保存文件到本地
到這里我們對(duì)數(shù)據(jù)集(dataset)的處理就全部完成了。接下來(lái)進(jìn)行最后一步準(zhǔn)備工作。
4. 數(shù)據(jù)的Annotation文件的準(zhǔn)備。
打開excel，輸入
第一行20表示20個(gè)樣本，2表示分為2組，1為識(shí)別字符（必須為1）
第二行表示分組信息，分為APOA2low，APOAhigh兩組，注意#號(hào)不要丟。
第三行對(duì)應(yīng)于dataset中每一列的分組，前10列輸入APOA2low，后10列輸入APOA2high。
然后點(diǎn)擊保存，文件名結(jié)尾加上.cls，可被GSEA識(shí)別，儲(chǔ)存為制表符分隔的文本文件（.txt)
打開文件檢查一下，應(yīng)該是這樣的。
至此，GSEA前期準(zhǔn)備就完成了。
下一步運(yùn)行GSEA軟件。
二. GSEA軟件的運(yùn)行
上傳文件，點(diǎn)擊load data, 使用Method1依次上傳.gct與.cls結(jié)尾的txt文件。
成功會(huì)有以下提示，如果失敗則可能是文件格式不對(duì)或者數(shù)據(jù)的排列不對(duì)。
然后我們成功上傳了兩個(gè)文件，如下圖。
接下來(lái)我們?cè)O(shè)置GSEA output的儲(chǔ)存位置，file-preferences。根據(jù)自己習(xí)慣設(shè)置。
2.設(shè)置參數(shù)。
我們進(jìn)入run GSEA界面
第一行選擇上傳的dataset，第二行g(shù)enesets選擇c5.bp（biological process)。genesets詳細(xì)介紹見MsigDB:http://software.broadinstitute.org/gsea/msigdb.
下一個(gè)選項(xiàng)，permutations內(nèi)存夠大CPU夠強(qiáng)，就選1000，這里選100。越大越好，越大結(jié)果越精確。
phenotype二選一均可。
這里選false，因?yàn)槲覀兊幕蛎呀?jīng)是gene symbol了，不需要轉(zhuǎn)換。
然后選phenotype，最后一行不選，如果基因名是測(cè)序平臺(tái)的ID的話最后一行要選對(duì)應(yīng)的測(cè)序平臺(tái)并將上面的false改成ture
最后點(diǎn)擊最下面的run.
成功的提示：
三. 結(jié)果解讀
打開儲(chǔ)存的文件夾，可以看到很多GSEA結(jié)果圖
值得注意的是一個(gè)名為index的網(wǎng)頁(yè)文件
打開它我們看以看到有446個(gè)基因sets在APOA2low中上調(diào)，1008個(gè)在APOA2high中上調(diào)，
點(diǎn)擊enrichment results in?html，得到下圖的列表，GS為基因集的名字，SIZE代表該基因集下的基因總數(shù)，ES代表Enrichment score, NES代表歸一化后的Enrichment score, NOM p-val代表pvalue，表征富集結(jié)果的可信度，F(xiàn)DR q-val代表qvalue, 是多重假設(shè)檢驗(yàn)矯正后的p值，注意GSEA采用pvalue < 5%, qvalue < 25% 對(duì)結(jié)果進(jìn)行過濾。Rank at max表示ES分?jǐn)?shù)達(dá)到最大值時(shí)最后一個(gè)基因所處的位置。說的通俗一點(diǎn)，比如說在我們此次GSEA的5000個(gè)基因里，有3個(gè)基因?qū)τ趯?duì)于cell cycle通路的ES分?jǐn)?shù)有正貢獻(xiàn)，三個(gè)基因分別位于第1，2，3位（本文開頭說過，GSEA會(huì)根據(jù)基因表達(dá)差異大小對(duì)基因進(jìn)行排序打分），他們對(duì)于ES分?jǐn)?shù)的貢獻(xiàn)分別為0.1，0.01，0.001，ES分?jǐn)?shù)在第3位達(dá)到最大值0.111, 此時(shí)的Rank at max就為3,第三位往后的基因?qū)ell cycle通路的ES分?jǐn)?shù)是負(fù)貢獻(xiàn)。
GSEA的結(jié)果圖：
峰值的頂點(diǎn)為最大富集分?jǐn)?shù)，下面的黑色豎條表示基因的位置，在最大富集分?jǐn)?shù)對(duì)應(yīng)橫坐標(biāo)位置左側(cè)的黑線表示leading edge subset，意思是這些gene的排序靠前，對(duì)ES分?jǐn)?shù)有主要貢獻(xiàn)。總而言之，這張圖主要是為了展示GSEA結(jié)果，說明我們確實(shí)通過GSEA分析發(fā)現(xiàn)了low和high兩個(gè)表型在這條通路存在差異，這張圖上的其他參數(shù)不必細(xì)究。
補(bǔ)充一下：
1. GSEA軟件自帶的geneset只有人的，如果是鼠的測(cè)序數(shù)據(jù)則要用鼠的geneset，把它下載下來(lái)再導(dǎo)入GSEA，很遺憾，我沒有找到鼠的geneset?；蛘邔⑹蟮幕蜣D(zhuǎn)換為人的進(jìn)行GSEA分析，這種轉(zhuǎn)換方法我后續(xù)會(huì)進(jìn)行介紹。
2. GSEA軟件還有很多可供調(diào)整的參數(shù)，文中沒有提到，還有一些其他功能，如leading edge analysis，這些功能的介紹我放在了文末的云盤中。
最后祝大家都能跑出結(jié)果！
資料鏈接：
數(shù)據(jù)集（dataset)以及可直接用來(lái)跑的的.cls與.gct文件，不過我還是建議用dataset從頭到尾自己分析一遍：
鏈接：https://pan.baidu.com/s/1mTlXDXSUe9CRMZKu8P74uw
提取碼：oyl2
java64位：
鏈接：https://pan.baidu.com/s/1qXoT_ocDA1F8zWlxbIJpRg
提取碼：i1xz
參考資料與深度學(xué)習(xí)資料：
鏈接：https://pan.baidu.com/s/1jNc3Ho8lRoONuKyCDbgJHw
提取碼：lvds
