MSigDB（Molecular Signatures Database）数据库里面的gmt文件超级多，是broad研究所为他们开发的gsea分析定义的文本文件规范，就是每一行都是一个通路（基因集合），每个行所代表的通路可以是不限制的列。但是第一列必须是通路的名字或者ID，第二个是通路的描述，第三列以及之后的全部的列都是基因名字或者ID即可。
不得不说，数据资源真心丰富啊：The 32880 gene sets in the Molecular Signatures Database (MSigDB) are divided into 9 major collections, and several sub-collections.
交流群的小伙伴神秘兮兮的给大家分享了他从文章附件supplementary pdf一个个抠出来的280 genes的 Splicing factor 基因列表 ，并且制作好了如下所示的 gmt文件 ，其实就是普通文本文件啦，编辑器打开可以看到就是一行的内容，如下所示：
Splicing_factors_geneset NA TARDBP SRRM1 PPP1R8 PPIE 后面的基因省略掉
写出gmt文件
假如你目前的基因列表在R里面，就可以自己写一个函数，比如 write.gmt 进行输出到gmt文件（本质上仍然是文本文件）：
library(msigdbr)
all_gene_sets = msigdbr(species = "Homo sapiens",
category='H')
length(unique(table(all_gene_sets$gs_name)))
tail(table(all_gene_sets$gs_name))
gcSample = split(all_gene_sets$gene_symbol,
all_gene_sets$gs_name)
names(gcSample)
file="Homo-H-examp.txt"
gs=gcSample
write.gmt <- function(gs,file){
sink(file)
lapply(names(gs), function(i){
cat( paste(c(i,'tmp',gs[[i]]),collapse='\t') )
cat('\n')
})
sink()
}
write.gmt(gs,file)
读取gmt文件
假如你拿到了gmt文件，很容易读取它并且去做分析，下面有两个不同包的函数：
file="Homo-H-examp.txt"
gmt <- GSEABase::getGmt( file ) #GeneSetCollection 对象
gmt
clusterProfiler::read.gmt( file ) # 普通数据框
跨越gmt文件
其实完全没必要写出gmt文件，然后再读取的，直接把自己的基因集合制作成为 GeneSetCollection 对象即可。
library(msigdbr)
all_gene_sets = msigdbr(species = "Homo sapiens",
category='H')
length(unique(table(all_gene_sets$gs_name)))
tail(table(all_gene_sets$gs_name))
gcSample = split(all_gene_sets$gene_symbol,
all_gene_sets$gs_name)
names(gcSample)
gs = lapply(gcSample, unique)
gsc <- GeneSetCollection(mapply(function(geneIds, keggId) {
GeneSet(geneIds, geneIdType=EntrezIdentifier(),
collectionType=KEGGCollection(keggId),
setName=keggId)
}, gs, names(gs)))
gsc
有了前面的 GeneSetCollection 对象 就可以很容易去做GSVA分析：
# 需要一个普通的表达量矩阵，X
es.max <- gsva(X, gsc,
mx.diff=FALSE, verbose=FALSE,
parallel.sz=8)
如果你使用clusterProfiler的read.gmt和GSEA
首先呢，clusterProfiler重新改写了gmt文件和做gsea的方法，所以代码稍微有一点点不同。
仍然是需要在msigdb数据库网页可以下载全部的基因集，我这里方便起见，仅仅是下载 h.all.v7.2.symbols.gmt文件：
### 对 MsigDB中的全部基因集 做GSEA分析。
# http://www./2105.html
# http://www./2102.html
# https://www./gsea/msigdb
# https://www./gsea/msigdb/collections.jsp
{
geneList= deg$avg_logFC
names(geneList)= toupper(rownames(deg))
geneList=sort(geneList,decreasing = T)
head(geneList)
library(ggplot2)
library(clusterProfiler)
library(org.Hs.eg.db)
#选择gmt文件（MigDB中的全部基因集）
gmtfile ='MSigDB/symbols/h.all.v7.2.symbols.gmt'
# 31120 个基因集
#GSEA分析
library(GSEABase) # BiocManager::install('GSEABase')
geneset <- clusterProfiler::read.gmt( gmtfile )
length(unique(geneset$term))
egmt <- GSEA(geneList, TERM2GENE=geneset,
minGSSize = 1,
pvalueCutoff = 0.99,
verbose=FALSE)
head(egmt)
egmt@result
gsea_results_df <- egmt@result
rownames(gsea_results_df)
write.csv(gsea_results_df,file = 'gsea_results_df.csv')
library(enrichplot)
gseaplot2(egmt,geneSetID = 'HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION',pvalue_table=T)
gseaplot2(egmt,geneSetID = 'HALLMARK_MTORC1_SIGNALING',pvalue_table=T)
}
代码都一脉相承，注意认真看文档，看数据结构即可。
关于 MSigDB
MSigDB（Molecular Signatures Database）数据库中定义了已知的基因集合：http://software./gsea/msigdb 需要注册才能下载。但是 msigdbr 包，它内置了MSigDB（Molecular Signatures Database）数据库的全部基因集合。安装方法非常简单：
install.packages("msigdbr")
包括H和C1-C7八个系列（Collection），每个系列分别是：
- H: hallmark gene sets （癌症）特征基因集合，共50组，最常用；
- C1: positional gene sets 位置基因集合，根据染色体位置，共326个，用的很少；
- C2: curated gene sets：（专家）校验基因集合，基于通路、文献等：
- C3: motif gene sets：模式基因集合，主要包括microRNA和转录因子靶基因两部分
- C4: computational gene sets：计算基因集合，通过挖掘癌症相关芯片数据定义的基因集合；
- C5: GO gene sets：Gene Ontology 基因本体论，包括BP（生物学过程biological process，细胞原件cellular component和分子功能molecular function三部分）
- C6: oncogenic signatures：癌症特征基因集合，大部分来源于NCBI GEO 发表芯片数据
- C7: immunologic signatures: 免疫相关基因集合。
因为我们这里针对的是大家熟知的pbmc3k数据集，所以这里选择C7: immunologic signatures: 免疫相关基因集合。如果你是肿瘤相关数据，建议选择H: hallmark gene sets （癌症）特征基因集合，共50组，最常用。
