把MsigDB数据库的全部通路转为gsva分析要求的输入格式

首先了解一下超几何分布检验和GSEA富集分析的区别：

• 通常拿到了上下调差异基因列表，然后说的GO/KEGG数据库注释，指的是超几何分布检验。
• 但是如果我们并不是在差异分析结果里面的自定义阈值，定上下调差异基因列表，而是根据某个指标（比如logFC）把全部的基因排序，再去进行GO/KEGG数据库注释，一般来说就是GSEA分析啦。

无论是超几何分布检验和GSEA富集分析，都离不开生物学功能数据库，数据库不仅仅是GO/KEGG哦，目前最齐全的应该是属于 MSigDB（Molecular Signatures Database）数据库中定义了已知的基因集合：http://software.broadinstitute.org/gsea/msigdb 包括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: 免疫相关基因集合。

可以看到，GO/KEGG是最出名的，但不是唯一的，起码和kegg数据库并列的就有Reactome数据库哈。

也可以自定义基因集

而且有各种各样的参考文献基因列表，比如转录因子列表，关于转录因子列表我在生信菜鸟团公众号看到了有一个介绍：TCGA数据挖掘常见基因集合，首先是Cancer Manag Res. 2020的文章《Prognostic and Predictive Value of a 15 Transcription Factors (TFs) Panel for Hepatocellular Carcinoma》就提到了：

• The TF list was downloaded from the Human Transcription Factors website (http://humantfs.ccbr.utoronto.ca/.

然后是2021的文章《A Transcription Factor-Based Risk Model for Predicting the Prognosis of Prostate Cancer and Potential Therapeutic Drugs》提到一个出处：

• Atotal of 1665 transcription factors were obtained from the Animal TFDB database。（AnimalTFDB 3.0 ），链接：http://bioinfo.life.hust.edu.cn/AnimalTFDB/#!/

这两个数据库关于转录因子的收录，都是接近于2000个基因。

这些来源于参考文献基因列表往往是千奇百怪的格式，它们并不会遵循MSigDB的gmt文件标准（其实绝大部分人应该是都没有听说过这个标准），绝大部分都是Excel里面的列表格式。

要么是长表，如下所示：

pathway1  gene1 
pathway1  gene2 
pathway1  gene3
pathway2  gene4 
pathway2  gene2

要么是不整齐的宽表格，如下所示：

pathway1  gene1 gene2 gene3
pathway2  gene4 gene2

这些就需要读入到R里面进行整理，然后才能承接到后面的注释步骤。详见；基因集合的数据框，列表和对象形式

转化 MsigDB数据库的全部通路转为gsva分析要求的输入格式

如果是从 MsigDB数据库下载，通常是gmt文件格式， 可以读入。之前我们的学徒作业，都是以公众号推文的方式发布出来，希望更多人加入一起学习，前面两次的作业是：

• 写一个函数把基因集，写出成为gmt文件。详见： GSVA或者GSEA各种算法都是可以自定义基因集的
• 自行读取作者给出的Nanostring表达矩阵，做差异分析，然后跟作者的结果进行比较。 详见：Nanostring的表达矩阵分析也是大同小异
• 使用TCGA数据源，制作生存分析，看看显著与否！详见：并不是只有TCGA才有临床信息用来做生存分析

其中写一个函数把基因集，写出成为gmt文件，我看到学徒完成的作业惨不忍睹。

这个时候其实有一个取巧的办法，就是使用msigdbr这个包，比如msigdbr包提取 KEGG数据信息：

library(msigdbr)  #install.packages("msigdbr") 
## msigdbr包提取下载  KEGG数据集
KEGG_df_all <-  msigdbr(species = "Homo sapiens", # Homo sapiens or Mus musculus
                        category = "C2",
                        subcategory = "CP:KEGG") 
colnames(KEGG_df_all)
KEGG_df <- dplyr::select(KEGG_df_all,gs_exact_source,gs_exact_source,ensembl_gene)  
kegg_list <- split(KEGG_df$ensembl_gene, KEGG_df$gs_exact_source)  
kegg_list

不过呢， 值得提醒的是msigdbr包提取 KEGG数据受限于 MsigDB数据库本身，里面的kegg信息是过时的，所以仍然是建议使用kegg数据库官方来源哈。

gsva是针对表达量矩阵进行分析

这个时候表达量矩阵很容易获得，多个基因在多个样品的表达量行列式而已，但是这个msigdbr这个包里面的通路信息没办法直接被gsva函数使用，需要一点点转化，代码如下所示：

library(msigdbr)
all_gene_sets = msigdbr(species = "Homo sapiens",
                        category='H') #从Msigdb中获取信号通路，此处H是癌症的标志物

#?msigdbr #可改成相应的感兴趣的数据集：C1-C7;H
# https://www.jianshu.com/p/f750ddcc440d （MsigDB）

length(unique(table(all_gene_sets$gs_name))) #癌症Marker对应的H有41条通路
tail(table(all_gene_sets$gs_name))

# 制作gsva分析所需要的genelist
gs=split(all_gene_sets$gene_symbol,all_gene_sets$gs_name)
gs = lapply(gs, unique)
head(gs)
gsc <- GeneSetCollection(mapply(function(geneIds, keggId) {
  GeneSet(geneIds, geneIdType=EntrezIdentifier(),
          collectionType=KEGGCollection(keggId),
          setName=keggId)
}, gs, names(gs)))
head(gsc)
geneset <- gsc

这个时候如果你有表达量矩阵，就可以很容易跑gsva啦，比如如下所示的代码里面的X变量就是多个基因在多个样品的表达量行列式矩阵：

# 运行gsva
es.max <- gsva(X, geneset, 
               mx.diff=FALSE, verbose=FALSE, 
               parallel.sz=8)
# 保存gsva分析所需要的结果
write.csv(es.max,file='es.max.csv')
head(es.max)

gsea是针对单个样品的表达量排序

如果是多个基因在多个样品的表达量行列式矩阵，需要一个个样品独立跑gsea分析，可以是每个样品的全部的基因自己的表达量排序，也可以某次差异分析两个分组后全部的基因在这个差异分析的变化倍数排序 ，示例代码如下所示 。

geneList = DEG_DESeq2$log2FoldChange
names(geneList)= rownames(DEG_DESeq2)
geneList=sort(geneList,decreasing = T)
head(geneList)

all_gene_sets = msigdbr(species =  "Mus musculus", # Homo sapiens or Mus musculus
                        category = "H" ) 
egmt <- GSEA(geneList, TERM2GENE= all_gene_sets[,c('gs_name','gene_symbol')] , 
             minGSSize = 1,
             pvalueCutoff = 1,
             verbose=FALSE)
head(egmt)
egmt@result 
gsea_results_df <- egmt@result 
rownames(gsea_results_df)