說(shuō)到富集，富集是將基因根據(jù)一些先驗(yàn)的知識(shí)（也就是常見(jiàn)的注釋）進(jìn)行分類的過(guò)程。我們一般會(huì)想到最常見(jiàn)的是GO/KEGG富集，其思路是先篩選差異基因，然后確定這些差異基因的GO/KEGG注釋，然后通過(guò)超幾何分布計(jì)算出哪些通路富集到了，通常會(huì)選擇一個(gè)閾值來(lái)卡一下，比如p值和FDR等。因此這會(huì)涉及到人為的閾值選擇，具有一定的主觀性，而且只能用于差異較大的基因，所以結(jié)果可能有一定的局限性。
根據(jù)上述情況，有了GSEA（Gene Set Enrichment Analysis），其思路是發(fā)表于2005年的Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles，主要是要有兩個(gè)概念：預(yù)先定義的基因集S（基于先驗(yàn)知識(shí)的基因注釋信息）和待測(cè)基因集L（一般是表達(dá)矩陣）；然后GSEA目的就是為了判斷S基因集中的基因是隨機(jī)分布于L（排序后的數(shù)據(jù)集），還是聚集分布在L的頂部或者底部（這也就是富集）。如果待測(cè)基因集中的某些基因顯著富集在L的頂部或者底部，這說(shuō)明這些基因的表達(dá)（因?yàn)槠涫歉鶕?jù)表達(dá)譜數(shù)據(jù)）對(duì)你定義的分組（預(yù)先分組）的差異有顯著影響（一致性），從而找到我們關(guān)注的基因集；在富集分析的理論中，GSEA可以認(rèn)為是第二代，即Functional Class Scoring (FCS) Approaches
這里不詳細(xì)說(shuō)算法，具體可看GSEA的文章，因?yàn)槲乙彩且恢虢狻！！?/p>
下載地址http://software.broadinstitute.org/gsea/downloads.jsp，PS.會(huì)驗(yàn)證下你的郵箱，先注冊(cè)下
第一次使用的話，而且數(shù)據(jù)不大的話，建議使用javaGSEA
Desktop Application，Java的圖形化界面，使用較為友好點(diǎn)，根據(jù)你的數(shù)據(jù)大小，選擇不同內(nèi)存的版本[1-8G不等]，PS.記得看看系統(tǒng)的Java版本，2G內(nèi)存開(kāi)始的GSEA版本需要的是64位的Java 1.8版
打開(kāi)后的界面如下：
因?yàn)镚SEA分析一般只作用于人物種的，所以我準(zhǔn)備以TCGA的BRCA的mRNA數(shù)據(jù)作為測(cè)試數(shù)據(jù)，正好也試下UCSC xena 瀏覽器才是最簡(jiǎn)單的TCGA數(shù)據(jù)下載途徑這個(gè)方法來(lái)下載TCGA數(shù)據(jù)（數(shù)據(jù)還蠻新的，2017年的）
數(shù)據(jù)下載地址：https://xenabrowser.net/datapages/?dataset=TCGA-BRCA%2FXena_Matrices%2FTCGA-BRCA.htseq_fpkm.tsv&host=https%3A%2F%2Fgdc.xenahubs.net
其還提供了ID/Gene Mapping的文件（整理好的），正好可以拿來(lái)用，因?yàn)殡m然GSEA有EnsemblID轉(zhuǎn)化的chip文件，但是感覺(jué)有些數(shù)據(jù)有點(diǎn)問(wèn)題（可能是由于Ensembl的版本一直在更新的緣故），HUGO gene symbol最好
然后用R處理下，將癌組織和對(duì)應(yīng)的癌旁組織的數(shù)據(jù)分別提取出來(lái)分別作為兩組的表達(dá)矩陣（gct文件）以及或者分組文件（cls文件）
library(dplyr)library(stringr)data <- read.table(file = "TCGA-BRCA.htseq_fpkm.tsv", sep = "\t", header = T, stringsAsFactors = F, check.names = F)data_df <- data[, -1]normal_df <- data_df[, str_sub(names(data_df), 14, 15) %in% 11:19]length(names(normal_df))tumor_df <- data_df[, grepl(paste(str_sub(names(normal_df), 1, 12), collapse = "|"), names(data_df)) & !names(data_df) %in% names(normal_df)]length(names(tumor_df))idmapping <- read.table(file = "gencode.v22.annotation.gene.probeMap", sep = "\t", header = T, stringsAsFactors = F)geneid <- data.frame(id = data$Ensembl_ID, stringsAsFactors = F)geneid2symbol <- left_join(geneid, idmapping, by = "id")all_df <- cbind(Name = geneid2symbol$gene, DESCRIPTION = "na", tumor_df, normal_df)group <- c(rep("Tumor", 118), rep("Normal", 113))group <- paste(group, collapse = " ")group <- c(paste(c(118+113, 2, 1), collapse = " "), "# Tumor Normal", group)write.table(file = "BRCA_fpkm.txt", all_df, sep = "\t", col.names = T, row.names = F, quote = F)write.table(file = "group.cls", group, col.names = F, row.names = F, quote = F)從上述代碼，我獲得118個(gè)癌組織樣本和對(duì)應(yīng)的113個(gè)癌旁樣本的表達(dá)譜數(shù)據(jù)，并且將Ensembl ID均轉(zhuǎn)化為了Gene symbol（避免之后用GSEA時(shí)，再用chip做ID轉(zhuǎn)化）；然后可以直接將txt文件作為輸入，也可以將BRCA_fpkm.txt文件做一些處理變成GSEA標(biāo)準(zhǔn)的gct文件，如下圖所示：分別加上前2行內(nèi)容，第一行照抄就行，第二行則是geneid數(shù)目和樣本數(shù)目
接著是Phenotype labels文件（上述代碼直接出了），即cls文件，格式如下圖所示：第一行231代表樣本數(shù)目，2代表分2組，空格間隔，1照抄；第二行井號(hào)注釋說(shuō)明分組信息；第三行為每個(gè)樣本對(duì)應(yīng)的組名，空格分隔
上述文件的詳細(xì)格式可參照網(wǎng)站：http://software.broadinstitute.org/cancer/software/gsea/wiki/index.php/Data_formats
如果網(wǎng)絡(luò)不佳的話，接下來(lái)最好將Gene sets file（也就是GSEA軟件上需要輸入的Gene sets database），作者將gene sets都儲(chǔ)存在Signature Database (MSigDb)中，去官網(wǎng)下載即可http://software.broadinstitute.org/gsea/downloads.jsp，比如下載個(gè)c2.cp.kegg.v6.1.symbols.gmt文件作為Gene sets database
如果數(shù)據(jù)是芯片數(shù)據(jù)或者需要GSEA的chip文件做ID轉(zhuǎn)化的話，則也可以先將chip文件下載下來(lái)，F(xiàn)TP地址：ftp://ftp.broadinstitute.org/pub/gsea/annotations
因?yàn)槭莣indows桌面式軟件，使用就比較簡(jiǎn)單了。首先將BRCA_fpkm.txt和group.cls文件導(dǎo)入，如果網(wǎng)速不好的話（比如我自己），那也把上述下載下來(lái)的c2.cp.kegg.v6.1.symbols.gmt文件也導(dǎo)入
接下來(lái)點(diǎn)擊RUN GSEA，就是幾個(gè)指定參數(shù)的選擇了，如下圖所示:
- Expression dataset：輸入的表達(dá)矩陣
- Gene sets database：gene sets文件，這里是
c2.cp.kegg.v6.1.symbols.gmt文件；PS.這個(gè)文件也可以自己創(chuàng)建哦- Number of permutations：置換檢驗(yàn)的次數(shù)
- Phenotype labels：選擇比較組，如果你輸入的文件就只有2個(gè)組別的話，這個(gè)就很方便選一個(gè)就行了；如果你輸入的有三個(gè)組別及以上的話，則這里就要跟你的需要選擇兩個(gè)組別的比較組，而且GSEA也會(huì)根據(jù)你的組別信息去表達(dá)矩陣中提取相對(duì)應(yīng)的數(shù)據(jù)
- Collapse dataset to gene symbols： 如果你已經(jīng)ID轉(zhuǎn)化為HUGO gene symbol，那么這里選FALSE，否則選擇TRUE
- Permutation type：選擇置換的類型，是random phenotype還是random gene sets，一般每組樣本數(shù)目大于7個(gè)時(shí)，建議選擇phenotype，否則選擇gene sets（這句話一直沒(méi)在官網(wǎng)上找到。。。似乎在文章里？）
- Chip platform：早些時(shí)候主要都是芯片數(shù)據(jù)，那時(shí)必須要將芯片id轉(zhuǎn)化為HUGO gene symbol，所以這個(gè)參數(shù)一直叫做chip（我猜的。。。），其實(shí)就是對(duì)ID進(jìn)行注釋，即ID轉(zhuǎn)化，選擇你ID對(duì)應(yīng)的chip文件即可，如果已自行轉(zhuǎn)化了ID的話，則這里空著就行（記得Collapse dataset to gene symbols選擇否）
除了Required field參數(shù)外，下面還有Basic fields和Advanced fields，具體參見(jiàn)官網(wǎng)吧（注：或者鼠標(biāo)懸浮在對(duì)應(yīng)參數(shù)名稱上，有簡(jiǎn)單的參數(shù)介紹哦）
最后點(diǎn)擊RUN，等待左下角的Running變成Success，然后點(diǎn)擊Success即可查看完整的結(jié)果，也可以點(diǎn)擊Show results folder，GSEA將所有結(jié)果都放在一個(gè)文件夾中了！??！
來(lái)看下文章里最常見(jiàn)的GSEA的結(jié)果圖片，如下圖所示：
從圖上，我們一般關(guān)注ES值，峰出現(xiàn)在前端還是后端（ES值大于0在前端，小于0在后端）以及Leading-edge subset（即對(duì)富集貢獻(xiàn)最大的部分,領(lǐng)頭亞集）；在ES圖中出現(xiàn)領(lǐng)頭亞集的形狀，表明這個(gè)功能基因集在某處理?xiàng)l件下具有更顯著的生物學(xué)意義；對(duì)于分析結(jié)果中，我們一般認(rèn)為|NES|>1，NOM p-val<0.05，F(xiàn)DR q-val<0.25的通路下的基因集合是有意義的
除了上述的結(jié)果外，GSEA還提供了Running the Leading Edge Analysis等操作，也可以看看
GSEA的結(jié)果解讀我也不是太熟悉，還是得多看看文獻(xiàn)中的解釋說(shuō)明啦
多于多個(gè)樣本的批處理，GSEA也有服務(wù)器版本，通過(guò)命令行即可操作，適合批處理操作；其還提供了R腳本可供使用（但官網(wǎng)上說(shuō)似乎并一定可行，需要自己調(diào)整？），反正我也正準(zhǔn)備都試試看。。。
參考資料：
功能數(shù)據(jù)庫(kù)專題-GSEA
GSEA富集分析 – 界面操作
http://software.broadinstitute.org/gsea/doc/desktop_tutorial.jsp
http://software.broadinstitute.org/gsea/doc/GSEAUserGuideFrame.html
GSEA使用（初級(jí)）
本文出自于http://www.bioinfo-scrounger.com轉(zhuǎn)載請(qǐng)注明出處
Bioinformatics for Proteomics Data2018年2月26日在“Proteomics”中
甲基化芯片入門學(xué)習(xí)-基礎(chǔ)知識(shí)（一）2018年1月9日在“Microarray”中
淺談蛋白組的差異蛋白分析2018年4月12日在“Proteomics”中
聯(lián)系客服
