說(shuō)到富集,，富集是將基因根據(jù)一些先驗(yàn)的知識(shí)（也就是常見(jiàn)的注釋?zhuān)┻M(jìn)行分類(lèi)的過(guò)程。我們一般會(huì)想到最常見(jiàn)的是GO/KEGG富集，其思路是先篩選差異基因,，然后確定這些差異基因的GO/KEGG注釋?zhuān)缓笸ㄟ^(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

GSEA的使用

這里不詳細(xì)說(shuō)算法，具體可看GSEA的文章,，因?yàn)槲乙彩且恢虢?。?！?/p>

下載地址http://software./gsea/downloads.jsp,，PS.會(huì)驗(yàn)證下你的郵箱,，先注冊(cè)下

第一次使用的話(huà)，而且數(shù)據(jù)不大的話(huà),，建議使用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)后的界面如下：

數(shù)據(jù)準(zhǔn)備

因?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:///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./cancer/software/gsea/wiki/index.php/Data_formats

如果網(wǎng)絡(luò)不佳的話(huà),，接下來(lái)最好將Gene sets file（也就是GSEA軟件上需要輸入的Gene sets database），作者將gene sets都儲(chǔ)存在Signature Database (MSigDb)中,，去官網(wǎng)下載即可http://software./gsea/downloads.jsp,，比如下載個(gè)c2.cp.kegg.v6.1.symbols.gmt文件作為Gene sets database

如果數(shù)據(jù)是芯片數(shù)據(jù)或者需要GSEA的chip文件做ID轉(zhuǎn)化的話(huà)，則也可以先將chip文件下載下來(lái),，F(xiàn)TP地址：ftp://ftp./pub/gsea/annotations

軟件使用

因?yàn)槭莣indows桌面式軟件,，使用就比較簡(jiǎn)單了,。首先將BRCA_fpkm.txt和group.cls文件導(dǎo)入，如果網(wǎng)速不好的話(huà)（比如我自己）,，那也把上述下載下來(lái)的c2.cp.kegg.v6.1.symbols.gmt文件也導(dǎo)入

接下來(lái)點(diǎn)擊RUN GSEA,，就是幾個(gè)指定參數(shù)的選擇了，如下圖所示:

1. Expression dataset：輸入的表達(dá)矩陣
2. Gene sets database：gene sets文件,，這里是c2.cp.kegg.v6.1.symbols.gmt文件,；PS.這個(gè)文件也可以自己創(chuàng)建哦
3. Number of permutations：置換檢驗(yàn)的次數(shù)
4. Phenotype labels：選擇比較組，如果你輸入的文件就只有2個(gè)組別的話(huà),，這個(gè)就很方便選一個(gè)就行了,；如果你輸入的有三個(gè)組別及以上的話(huà)，則這里就要跟你的需要選擇兩個(gè)組別的比較組,，而且GSEA也會(huì)根據(jù)你的組別信息去表達(dá)矩陣中提取相對(duì)應(yīng)的數(shù)據(jù)
5. Collapse dataset to gene symbols： 如果你已經(jīng)ID轉(zhuǎn)化為HUGO gene symbol,，那么這里選FALSE，否則選擇TRUE
6. Permutation type：選擇置換的類(lèi)型,，是random phenotype還是random gene sets,，一般每組樣本數(shù)目大于7個(gè)時(shí)，建議選擇phenotype,，否則選擇gene sets（這句話(huà)一直沒(méi)在官網(wǎng)上找到。,。,。似乎在文章里？）
7. Chip platform：早些時(shí)候主要都是芯片數(shù)據(jù),，那時(shí)必須要將芯片id轉(zhuǎn)化為HUGO gene symbol,，所以這個(gè)參數(shù)一直叫做chip（我猜的。,。,。），其實(shí)就是對(duì)ID進(jìn)行注釋?zhuān)碔D轉(zhuǎn)化,，選擇你ID對(duì)應(yīng)的chip文件即可,，如果已自行轉(zhuǎn)化了ID的話(huà)，則這里空著就行（記得Collapse dataset to gene symbols選擇否）

除了Required field參數(shù)外,，下面還有Basic fields和Advanced fields,，具體參見(jiàn)官網(wǎng)吧（注：或者鼠標(biāo)懸浮在對(duì)應(yīng)參數(shù)名稱(chēng)上，有簡(jiǎn)單的參數(shù)介紹哦）

最后點(diǎn)擊RUN,，等待左下角的Running變成Success,，然后點(diǎn)擊Success即可查看完整的結(jié)果，也可以點(diǎn)擊Show results folder,，GSEA將所有結(jié)果都放在一個(gè)文件夾中了?。,。?/p>

分析結(jié)果

來(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ù)專(zhuān)題-GSEA
GSEA富集分析 – 界面操作
http://software./gsea/doc/desktop_tutorial.jsp
http://software./gsea/doc/GSEAUserGuideFrame.html
GSEA使用（初級(jí)）

本文出自于http://www.轉(zhuǎn)載請(qǐng)注明出處