GSEA的分析汇总
学习GSEA
GSEA的统计学原理试讲
GSEA
GSEA这个java软件使用非常方便，只需要根据要求做好GCT/CLS格式的input文件就好了。我以前也写过用法教程:
原理
但说到统计学原理，就有点麻烦了，我试着用自己的思路阐释一下：
假设芯片或者其它测量方法测到了2万个基因，那么这两万个基因在case和control组的差异度量(六种差异度量，默认是signal 2 noise，GSEA官网有提供公式，也可以选择大家熟悉的foldchange)肯定不一样。
那么根据它们的差异度量，就可以对它们进行排序，并且Z-score标准化，在下图的最底端展示的就是
image那么图中间，就是我们每个gene set里面的基因在所有的2万个排序好基因的位置，如果gene set里面的基因集中在2万个基因的前面部分，就是在case里面富集，如果集中在后面部分，就是在control里面富集着。
而最上面的那个ES score的算法，大概如下：
image仔细看，其实还是能看明白的，每个基因在每个gene set里面的ES score取决于这个基因是否属于该gene set，还有就是它的差异度量，上图的差异度量就是FC（foldchange）,对每个gene set来说，所有的基因的ES score都要一个个加起来，叫做running ES score，在加的过程中，什么时候ES score达到了最大值，就是这个gene set最终的ES score！
算法解读我参考的PPT，反正我是看懂了，但不一定能讲清楚：
软件还有大把的参数可以调整，自行摸索！
用GSEA来做基因集富集分析
how to use GSEA?
这个有点类似于pathway（GO,KEGG等）的富集分析，区别在于gene set（校验好的基于文献的数据库）的概念更广泛一点
how to download GSEA ?
what's the input for the GSEA?
说明书上写的输入数据是：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.
实际上没那么复杂，一个表达矩阵即可！然后做一个分组说明的cls文件即可。
表达矩阵我这里下载GSE1009数据集做测试吧！
cls的样本说明文件，就随便搞一搞吧，下面这个是例子：
6 2 1
# good bad
good good good bad bad bad
文件如下，六个样本，根据探针来的表达数据，分组前后各三个一组。
imagestart to run the GSEA !
首先载入数据
image确定无误，就开始运行，运行需要设置一定的参数！
imagewhat's the output ?
输出的数据非常多，对你选择的gene set数据集里面的每个set都会分析看看是否符合富集的标准，富集就出来一个报告。
点击success就能进入报告主页，里面的链接可以进入任意一个分报告。
最大的特色是提供了大量的数据集：
参考文献
有些文献是基于GSEA的：
批量运行GSEA，命令行版本
之前用过有界面的那种，那样非常方便，只需要做好数据即可，但是如果有非常多的数据，每次都要点击文件，点击下一步，也很烦，不过，，它既然是java软件，就可以以命令行的形式来玩转它！
一、程序安装
二、输入数据
image三、运行命令
hgu95av2的芯片数据，只有一万多探针，所以很快就可以出结果
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
四、输出数据
image自己看官网去理解这些结果咯！
需要下载的数据如下：
都是gmt格式的文件！
| CP: CANONICAL PATHWAYS() | Gene sets from the pathway databases. Usually, these gene sets are canonical representations of a biological process compiled by domain experts. | Download GMT Files |
|---|---|---|
| CP:BIOCARTA: BIOCARTA GENE SETS() | Gene sets derived from the BioCarta pathway database (). | Download GMT Files |
| CP:KEGG: KEGG GENE SETS() | Gene sets derived from the KEGG pathway database (). | Download GMT Files |
| CP:REACTOME: REACTOME GENE SETS() | Gene sets derived from the Reactome pathway database (). | Download GMT Files |
然后做出表达数据gct文件和cls表型文件~
然后就可以直接运行了
制作自己的gene set文件给gsea软件
gene set
image我这里提供一个2列的文件，直接转换成gmt的R代码！
image首先在R里面赋值一个变量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)
这个变量kegg2GeneID_list是一个list，因为是entrez gene ID，需要转换成symbol，我就不多说了，转换后的数据，就是kegg2symbol_list 。
最后对 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()
}