通常拿到了上下调差异基因列表，然后说的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里面的列表格式。
要么是长表，如下所示：
代码语言：javascript复制pathway1 gene1
pathway1 gene2
pathway1 gene3
pathway2 gene4
pathway2 gene2
要么是不整齐的宽表格，如下所示：
代码语言：javascript复制pathway1 gene1 gene2 gene3
pathway2 gene4 gene2
这些就需要读入到R里面进行整理，然后才能承接到后面的注释步骤。如下所示就是长短不一的Excel，读取就考验大家的代码能力了：
数据框
这个大概是基因集合最容易看人看懂的形式了，
代码语言：javascript复制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( 函数：
代码语言：javascript复制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
列表
在单细胞的标记基因可视化我们经常会遇到：
代码语言：javascript复制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打分函数支持：
代码语言：javascript复制library(Seurat)
sce = AddModuleScore(sce,genes_to_check ,name = names(glist))
这样的列表如果想转换成为前面的数据框也很容易：
代码语言：javascript复制 TERM2GENE = do.call(rbind, lapply(names(genes_to_check), function(x){
data.frame(gs_name=x,gene_symbol=glist[[x]])
}))
对象（遵循MSigDB的gmt文件标准）
前面的数据框或者列表，要弄成对象就比较麻烦了，需要做一些转换：
代码语言：javascript复制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.
十年后我环游世界各地的高校以及科研院所（当然包括中国大陆）的时候，如果有这样的情谊，我会优先见你。
