熟悉GSEA软件的都知道，它只需要GCT,CLS和GMT文件，其中GMT文件，GSEA的作者已经给出了一大堆！就是记录broad的Molecular Signatures Database (MSigDB) 已经收到了18026个geneset，MSigDB的确是多，但未必全，其实里面还有很多重复。而且有不少几乎没有意义的gene set。
了解GMT格式
那我想做自己的gene set来用gsea软件做分析，就需要自己制造gmt格式的数据。因为即使下载了MSigDB的gene set，本质上就是gmt格式的数据而已，如下图：
4
说明的非常清楚，每一行都是一个独立的gene set ， 第一列是该gene set的名字，第二列打酱油，后面所有的列都是该gene set 包含着的基因名字。每个gene set包含的基因个数不一样。
自己制作gene set文件
我们首先要拿到自己感兴趣的gene set里面的gene list，最好是以hugo规定的标准symbol。比如我感兴趣的是 ：http://www.cta.lncc.br/modelo.php 这里我提供一个2列的文件直接转换成gmt的R代码！
文件来自于：下载最新版的KEGG信息，并且解析好，如下：
img
在R里面赋值一个变量path2gene_file就是读取如上图所示的kegg2gene.txt文件
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()
}
得到的GMT格式的gene set文件如下：
5
网友评论