基因集合的数据框，列表和对象形式

通常拿到了上下调差异基因列表，然后说的GO/KEGG数据库注释，指的是超几何分布检验。但是如果我们并不是在差异分析结果里面的自定义阈值，定上下调差异基因列表，而是根据某个指标（比如logFC）把全部的基因排序，再去进行GO/KEGG数据库注释，一般来说就是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里面进行整理，然后才能承接到后面的注释步骤。如下所示就是长短不一的Excel，读取就考验大家的代码能力了：

数据框

这个大概是基因集合最容易看人看懂的形式了，

library(msigdbr)
all_gene_sets = msigdbr(species =  "Homo sapiens", # Homo sapiens or Mus musculus
                        category = "H" ) 

length(unique(table(all_gene_sets$gs_name)))

可以看到是长型数据框哦，因为数据框不能是不整齐的，所以没办法是宽的，每个基因集合里面的基因个数不一样，大概率都是不整齐的。

这种数据框格式的基因列表适合于 clusterProfiler::GSEA( 函数：

names(all_gene_sets)
egmt <- GSEA(geneList, TERM2GENE= all_gene_sets[,c('gs_name','gene_symbol')] , 
             minGSSize = 1,
             pvalueCutoff = 0.99,
             verbose=FALSE)
head(egmt)
egmt@result 
gsea_results_df <- egmt@result  

列表

在单细胞的标记基因可视化我们经常会遇到：

Myo=c("Krt17", "Krt14", "Krt5", "Acta2", "Myl9", "Mylk", "Myh11")
Lum=c("Krt19", "Krt18", "Krt8")
Hs=c("Prlr", "Cited1", "Pgr", "Prom1", "Esr1")  
AV=c("Mfge8", "Trf", "Csn3", "Wfdc18", "Elf5", "Ltf")
Lp=c("Kit", "Aldh1a3", "Cd14")
Fib=c("Col1a1", "Col1a2", "Col3a1", "Fn1")
genes_to_check =list(
  Myo=Myo,
  Lum=Lum,
  Hs=Hs, 
  AV=AV,
  Lp=Lp,  
  Fib=Fib 
  
) 

一般来说做gsea或者gsva这样的打分都不支持这样的列表格式，不过Seurat自带的AddModuleScore打分函数支持：

library(Seurat)
sce = AddModuleScore(sce,genes_to_check ,name = names(glist))

这样的列表如果想转换成为前面的数据框也很容易：

 TERM2GENE = do.call(rbind, lapply(names(genes_to_check), function(x){
  data.frame(gs_name=x,gene_symbol=glist[[x]])
}))

对象（遵循MSigDB的gmt文件标准）

前面的数据框或者列表，要弄成对象就比较麻烦了，需要做一些转换：

library(GSVA) # BiocManager::install('GSVA')
library(GSEABase)

gs=split(all_gene_sets$gene_symbol,all_gene_sets$gs_name)
gs = lapply(gs, unique)

geneset <- GeneSetCollection(mapply(function(geneIds, keggId) {
  GeneSet(geneIds, geneIdType=EntrezIdentifier(),
          collectionType=KEGGCollection(keggId),
          setName=keggId)
}, gs, names(gs))) 
# 这个 gsva 函数可以根据前面的  geneset对象，对任意表达量矩阵进行分析
es.max <- gsva(as.matrix( sce@assays$RNA@counts), geneset, 
               mx.diff=FALSE, verbose=FALSE, 
               parallel.sz=8) 

这个函数非常重要，GeneSetCollection，自己去理解吧。

写在文末

我在《生信技能树》，《生信菜鸟团》，《单细胞天地》的大量推文教程里面共享的代码都是复制粘贴即可使用的， 有任何疑问欢迎留言讨论，也可以发邮件给我，详细描述你遇到的困难的前因后果给我，我的邮箱地址是 jmzeng1314@163.com

如果你确实觉得我的教程对你的科研课题有帮助，让你茅塞顿开，或者说你的课题大量使用我的技能，烦请日后在发表自己的成果的时候，加上一个简短的致谢，如下所示：

We thank Dr.Jianming Zeng(University of Macau), and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes.

十年后我环游世界各地的高校以及科研院所（当然包括中国大陆）的时候，如果有这样的情谊，我会优先见你。