說到富集，富集是將基因根據(jù)一些先驗(yàn)的知識（也就是常見的注釋）進(jìn)行分類的過程。我們一般會想到最常見的是GO/KEGG富集，其思路是先篩選差異基因，然后確定這些差異基因的GO/KEGG注釋，然后通過超幾何分布計算出哪些通路富集到了，通常會選擇一個閾值來卡一下，比如p值和FDR等。因此這會涉及到人為的閾值選擇，具有一定的主觀性，而且只能用于差異較大的基因，所以結(jié)果可能有一定的局限性。

根據(jù)上述情況，有了GSEA（Gene Set Enrichment Analysis），其思路是發(fā)表于2005年的Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles，主要是要有兩個概念：預(yù)先定義的基因集S（基于先驗(yàn)知識的基因注釋信息）和待測基因集L（一般是表達(dá)矩陣）；然后GSEA目的就是為了判斷S基因集中的基因是隨機(jī)分布于L（排序后的數(shù)據(jù)集），還是聚集分布在L的頂部或者底部（這也就是富集）。如果待測基因集中的某些基因顯著富集在L的頂部或者底部，這說明這些基因的表達(dá)（因?yàn)槠涫歉鶕?jù)表達(dá)譜數(shù)據(jù)）對你定義的分組（預(yù)先分組）的差異有顯著影響（一致性），從而找到我們關(guān)注的基因集；在富集分析的理論中，GSEA可以認(rèn)為是第二代，即Functional Class Scoring (FCS) Approaches

GSEA的使用

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

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

第一次使用的話，而且數(shù)據(jù)不大的話，建議使用javaGSEA
Desktop Application，Java的圖形化界面，使用較為友好點(diǎn)，根據(jù)你的數(shù)據(jù)大小，選擇不同內(nèi)存的版本[1-8G不等]，PS.記得看看系統(tǒng)的Java版本，2G內(nèi)存開始的GSEA版本需要的是64位的Java 1.8版

打開后的界面如下：

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

因?yàn)镚SEA分析一般只作用于人物種的，所以我準(zhǔn)備以TCGA的BRCA的mRNA數(shù)據(jù)作為測試數(shù)據(jù)，正好也試下UCSC xena 瀏覽器才是最簡單的TCGA數(shù)據(jù)下載途徑這個方法來下載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的文件（整理好的），正好可以拿來用，因?yàn)殡m然GSEA有EnsemblID轉(zhuǎn)化的chip文件，但是感覺有些數(shù)據(jù)有點(diǎn)問題（可能是由于Ensembl的版本一直在更新的緣故），HUGO gene symbol最好

然后用R處理下，將癌組織和對應(yīng)的癌旁組織的數(shù)據(jù)分別提取出來分別作為兩組的表達(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個癌組織樣本和對應(yīng)的113個癌旁樣本的表達(dá)譜數(shù)據(jù)，并且將Ensembl ID均轉(zhuǎn)化為了Gene symbol（避免之后用GSEA時，再用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照抄；第二行井號注釋說明分組信息；第三行為每個樣本對應(yīng)的組名，空格分隔

上述文件的詳細(xì)格式可參照網(wǎng)站：http://software./cancer/software/gsea/wiki/index.php/Data_formats

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

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

軟件使用

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

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

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

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

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

分析結(jié)果

來看下文章里最常見的GSEA的結(jié)果圖片，如下圖所示：

從圖上，我們一般關(guān)注ES值，峰出現(xiàn)在前端還是后端（ES值大于0在前端，小于0在后端）以及Leading-edge subset（即對富集貢獻(xiàn)最大的部分,領(lǐng)頭亞集）；在ES圖中出現(xiàn)領(lǐng)頭亞集的形狀，表明這個功能基因集在某處理?xiàng)l件下具有更顯著的生物學(xué)意義；對于分析結(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)中的解釋說明啦

多于多個樣本的批處理，GSEA也有服務(wù)器版本，通過命令行即可操作，適合批處理操作；其還提供了R腳本可供使用（但官網(wǎng)上說似乎并一定可行，需要自己調(diào)整？），反正我也正準(zhǔn)備都試試看。。。

參考資料：

功能數(shù)據(jù)庫專題-GSEA
GSEA富集分析 – 界面操作
http://software./gsea/doc/desktop_tutorial.jsp
http://software./gsea/doc/GSEAUserGuideFrame.html
GSEA使用（初級）

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