GSEA的分析匯總
學(xué)習(xí)GSEA
GSEA的統(tǒng)計(jì)學(xué)原理試講
GSEA
GSEA這個(gè)java軟件使用非常方便，只需要根據(jù)要求做好GCT/CLS格式的input文件就好了彼硫。我以前也寫(xiě)過(guò)用法教程:
原理
但說(shuō)到統(tǒng)計(jì)學(xué)原理，就有點(diǎn)麻煩了乌助，我試著用自己的思路闡釋一下：
假設(shè)芯片或者其它測(cè)量方法測(cè)到了2萬(wàn)個(gè)基因溜在，那么這兩萬(wàn)個(gè)基因在case和control組的差異度量(六種差異度量陌知，默認(rèn)是signal 2 noise他托，GSEA官網(wǎng)有提供公式，也可以選擇大家熟悉的foldchange)肯定不一樣仆葡。
那么根據(jù)它們的差異度量赏参，就可以對(duì)它們進(jìn)行排序，并且Z-score標(biāo)準(zhǔn)化沿盅，在下圖的最底端展示的就是
那么圖中間把篓，就是我們每個(gè)gene set里面的基因在所有的2萬(wàn)個(gè)排序好基因的位置，如果gene set里面的基因集中在2萬(wàn)個(gè)基因的前面部分腰涧，就是在case里面富集韧掩，如果集中在后面部分，就是在control里面富集著窖铡。
而最上面的那個(gè)ES score的算法疗锐，大概如下：
仔細(xì)看坊谁，其實(shí)還是能看明白的，每個(gè)基因在每個(gè)gene set里面的ES score取決于這個(gè)基因是否屬于該gene set滑臊，還有就是它的差異度量口芍，上圖的差異度量就是FC（foldchange）,對(duì)每個(gè)gene set來(lái)說(shuō)，所有的基因的ES score都要一個(gè)個(gè)加起來(lái)雇卷，叫做running ES score鬓椭，在加的過(guò)程中，什么時(shí)候ES score達(dá)到了最大值关划，就是這個(gè)gene set最終的ES score小染！
算法解讀我參考的PPT，反正我是看懂了祭玉，但不一定能講清楚：
http://bioinformatics.mdanderson.org/MicroarrayCourse/Lectures09/gsea1_bw.pdf
https://bioinformatics.cancer.gov/sites/default/files/course_material/GSEA_Theory.pptx
http://compbio.ucdenver.edu/Hunter_lab/Phang/downloads/files/GSEA.ppt
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1239896/
http://www.baderlab.org/CancerStemCellProject/VeroniqueVoisin/AdditionalResources/GSEA
軟件還有大把的參數(shù)可以調(diào)整氧映，自行摸索！
用GSEA來(lái)做基因集富集分析
how to use GSEA?
這個(gè)有點(diǎn)類(lèi)似于pathway（GO,KEGG等）的富集分析脱货，區(qū)別在于gene set（校驗(yàn)好的基于文獻(xiàn)的數(shù)據(jù)庫(kù)）的概念更廣泛一點(diǎn)
how to download GSEA ?
http://software.broadinstitute.org/gsea/downloads.jsp
教程：http://software.broadinstitute.org/gsea/doc/desktop_tutorial.jsp 岛都，需要自己安裝好java環(huán)境！
what's the input for the GSEA?
說(shuō)明書(shū)上寫(xiě)的輸入數(shù)據(jù)是：GSEA supported data files are simply tab delimited ASCII text files, which have special file extensions that identify them. For example, expression data usually has the extension *.gct, phenotypes *.cls, gene sets *.gmt, and chip annotations *.chip. Click the More on file formats help button to view detailed descriptions of all the data file formats.
并且提供了測(cè)試數(shù)據(jù)：http://software.broadinstitute.org/gsea/datasets.jsp
實(shí)際上沒(méi)那么復(fù)雜振峻，一個(gè)表達(dá)矩陣即可臼疫！然后做一個(gè)分組說(shuō)明的cls文件即可。
主要是自己看說(shuō)明書(shū)扣孟，做出要求的數(shù)據(jù)格式：http://www.broadinstitute.org/cancer/software/gsea/wiki/index.php/Data_formats
表達(dá)矩陣我這里下載GSE1009數(shù)據(jù)集做測(cè)試吧烫堤！
- http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=gse1009
- ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE1nnn/GSE1009/matrix/GSE1009_series_matrix.txt.gz
cls的樣本說(shuō)明文件，就隨便搞一搞吧凤价，下面這個(gè)是例子：
6 2 1
# good bad
good good good bad bad bad
文件如下鸽斟，六個(gè)樣本，根據(jù)探針來(lái)的表達(dá)數(shù)據(jù)利诺，分組前后各三個(gè)一組富蓄。
start to run the GSEA !
首先載入數(shù)據(jù)
確定無(wú)誤，就開(kāi)始運(yùn)行慢逾，運(yùn)行需要設(shè)置一定的參數(shù)立倍！
what's the output ?
輸出的數(shù)據(jù)非常多，對(duì)你選擇的gene set數(shù)據(jù)集里面的每個(gè)set都會(huì)分析看看是否符合富集的標(biāo)準(zhǔn)侣滩，富集就出來(lái)一個(gè)報(bào)告口注。
點(diǎn)擊success就能進(jìn)入報(bào)告主頁(yè)，里面的鏈接可以進(jìn)入任意一個(gè)分報(bào)告君珠。
最大的特色是提供了大量的數(shù)據(jù)集：
You can browse the MSigDB from the Molecular Signatures Database page of the GSEA web site or the Browse MSigDB page of the GSEA application. To browse the MSigDB from the GSEA application:
還自己建立了wiki說(shuō)明主頁(yè)：http://www.broadinstitute.org/cancer/software/gsea/wiki/index.php/Main_Page
參考文獻(xiàn)
有些文獻(xiàn)是基于GSEA的：
http://www.ncbi.nlm.nih.gov/pubmed/16199517
http://stke.sciencemag.org/highwire/filestream/4681053/field_highwire_adjunct_files/1/2001966_Slides.zip
http://www.ingentaconnect.com/content/ben/cbio/2007/00000002/00000002/art00003
http://www.nature.com/articles/ng0704-663a
http://bioinformatics.oxfordjournals.org/content/23/23/3251.short
http://link.springer.com/article/10.1007/s00335-011-9359-x
批量運(yùn)行GSEA寝志，命令行版本
之前用過(guò)有界面的那種，那樣非常方便，只需要做好數(shù)據(jù)即可材部，但是如果有非常多的數(shù)據(jù)悠菜，每次都要點(diǎn)擊文件，點(diǎn)擊下一步败富，也很煩悔醋，不過(guò)，兽叮，它既然是java軟件芬骄，就可以以命令行的形式來(lái)玩轉(zhuǎn)它！
一鹦聪、程序安裝
直接在官網(wǎng)下載java版本軟件即可：http://software.broadinstitute.org/gsea/downloads.jsp
二账阻、輸入數(shù)據(jù)
需要下載gmt文件，自己制作gct和cls文件泽本，或者直接下載測(cè)試文件p53
http://www.broadinstitute.org/cancer/software/gsea/wiki/index.php/Data_formats
三淘太、運(yùn)行命令
hgu95av2的芯片數(shù)據(jù)，只有一萬(wàn)多探針规丽，所以很快就可以出結(jié)果
java -cp gsea2-2.2.1.jar -Xmx1024m xtools.gsea.Gsea -gmx c2.cp.kegg.v5.0.symbols.gmt \
-res P53_hgu95av2.gct -cls P53.cls -chip chip/HG_U95Av2.chip -out result -rpt_label p53_example
這個(gè)參數(shù)其實(shí)非常多蒲牧，見(jiàn)文件http://software.broadinstitute.org/gsea/doc/linkedFiles/GSEAParameters.txt
四、輸出數(shù)據(jù)
自己看官網(wǎng)去理解這些結(jié)果咯赌莺！
需要下載的數(shù)據(jù)如下：
http://software.broadinstitute.org/gsea/downloads.jsp
都是gmt格式的文件冰抢！
| CP: CANONICAL PATHWAYS(BROWSE 1330 GENE SETS) | Gene sets from the pathway databases. Usually, these gene sets are canonical representations of a biological process compiled by domain experts. details | Download GMT Filesoriginal identifiersgene symbolsentrez genes ids |
|---|---|---|
| CP:BIOCARTA: BIOCARTA GENE SETS(BROWSE 217 GENE SETS) | Gene sets derived from the BioCarta pathway database (http://www.biocarta.com/genes/index.asp). | Download GMT Filesoriginal identifiersgene symbolsentrez genes ids |
| CP:KEGG: KEGG GENE SETS(BROWSE 186 GENE SETS) | Gene sets derived from the KEGG pathway database (http://www.genome.jp/kegg/pathway.html). | Download GMT Filesoriginal identifiersgene symbolsentrez genes ids |
| CP:REACTOME: REACTOME GENE SETS(BROWSE 674 GENE SETS) | Gene sets derived from the Reactome pathway database (http://www.reactome.org/). | Download GMT Filesoriginal identifiersgene symbolsentrez genes ids |
然后做出表達(dá)數(shù)據(jù)gct文件和cls表型文件~
見(jiàn)：http://www.broadinstitute.org/cancer/software/gsea/wiki/index.php/Data_formats
然后就可以直接運(yùn)行了
如果是芯片數(shù)據(jù)，第一列是芯片探針艘狭，那么還需要下載chip數(shù)據(jù)：ftp://ftp.broadinstitute.org/pub/gsea/annotations
制作自己的gene set文件給gsea軟件
gene set
熟悉GSEA軟件的都知道挎扰，它只需要GCT,CLS和GMT文件，其中GMT文件巢音，GSEA的作者已經(jīng)給出了一大堆遵倦！就是記錄broad的Molecular Signatures Database (MSigDB) 已經(jīng)收到了18026個(gè)geneset， 但是我奇怪的是里面竟然沒(méi)有包括cancer testis的gene set官撼，MSigDB的確是多梧躺，但未必全，其實(shí)里面還有很多重復(fù)歧寺。而且有不少幾乎沒(méi)有意義的gene set燥狰。那我想做自己的gene set來(lái)用gsea軟件做分析棘脐，就需要自己制造gmt格式的數(shù)據(jù)斜筐。因?yàn)榧词瓜螺d了MSigDB的gene set，本質(zhì)上就是gmt格式的數(shù)據(jù)而已：http://software.broadinstitute.org/cancer/software/gsea/wiki/index.php/Data_formats#GMT:Gene_Matrix_Transposed_file_format.28.2A.gmt.29
我們首先要拿到自己感興趣的gene set里面的gene list蛀缝，最好是以hugo規(guī)定的標(biāo)準(zhǔn)symbol顷链。
比如我感興趣的是 ：http://www.cta.lncc.br/modelo.php
我這里提供一個(gè)2列的文件，直接轉(zhuǎn)換成gmt的R代碼屈梁！
文件來(lái)自于：下載最新版的KEGG信息嗤练，并且解析好榛了，如下：
首先在R里面賦值一個(gè)變量path2gene_file就是圖中的kegg2gene.txt文件，讀到R里面去
tmp=read.table(path2gene_file,sep="\t",colClasses=c('character'))
# tmp=toTable(org.Hs.egPATH)
# first column is kegg ID, second column is entrez ID
GeneID2kegg_list<<- tapply(tmp[,1],as.factor(tmp[,2]),function(x) x)
kegg2GeneID_list<<- tapply(tmp[,2],as.factor(tmp[,1]),function(x) x)
這個(gè)變量kegg2GeneID_list是一個(gè)list煞抬，因?yàn)槭莈ntrez gene ID霜大，需要轉(zhuǎn)換成symbol，我就不多說(shuō)了革答，轉(zhuǎn)換后的數(shù)據(jù)战坤，就是kegg2symbol_list 。
最后對(duì) kegg2symbol_list 輸出成gmt文件：
write.gmt <- function(geneSet=kegg2symbol_list,gmt_file='kegg2symbol.gmt'){
sink( gmt_file )
for (i in 1:length(geneSet)){
cat(names(geneSet)[i])
cat('\tNA\t')
cat(paste(geneSet[[i]],collapse = '\t'))
cat('\n')
}
sink()
}
