单细胞分析避坑指南：KEGGREST与msigdbr获取KEGG基因列表的深度对比与实战解析
在单细胞转录组分析中，基因集富集分析（GSEA）和基因集变异分析（GSVA）是揭示细胞状态和功能的重要工具。而获取高质量的KEGG通路基因列表，往往是整个分析流程的第一个关键步骤。面对R语言生态中多个可选的工具包，初学者常陷入选择困境：msigdbr提供的186条通路是否够用？KEGGREST返回的357条通路又该如何处理？本文将带您深入解剖这两个主流工具的差异，并提供可直接落地的解决方案。
1. 工具选择的核心考量维度
1.1 数据来源与更新机制
msigdbr的数据来源于MSigDB数据库的C2子集，其KEGG通路是通过定期抓取KEGG官网信息构建的。而KEGGREST则是直接通过KEGG官方API实时获取数据，这意味着：
- 时效性差异：KEGGREST能获取最新添加的通路（如2023年新增的COVID-19相关通路），而msigdbr可能存在数月延迟
- 数据完整性：KEGGREST包含357条人类通路，msigdbr仅186条，缺失的通路主要分布在：
- 药物代谢（如hsa00983）
- 罕见疾病（如hsa04925）
- 新兴信号通路（如hsa04932）
提示：如果研究涉及最新发现的生物过程，建议优先使用KEGGREST
1.2 物种覆盖范围对比
虽然两者都支持多物种分析，但实现方式截然不同：
| 特性 | msigdbr | KEGGREST |
|---|---|---|
| 内置物种数 | 17个常用模式生物 | 全物种支持（通过KEGG代码） |
| 非模式生物支持 | 有限 | 完整（如山羊：hsa→bta） |
| 基因ID类型 | Symbol/Entrez | 需自行转换 |
| 安装复杂度 | 简单（CRAN） | 需Bioconductor依赖 |
# msigdbr物种列表示例 library(msigdbr) unique(msigdbr_species()$species_name) # KEGGREST物种查询 KEGGREST::keggList("organism")[1:5,]1.3 输出格式兼容性
与下游分析工具的适配性直接影响工作效率：
msigdbr天然输出为data.frame，与clusterProfiler无缝衔接：
genesets <- msigdbr(species = "Homo sapiens", category = "C2", subcategory = "KEGG") head(genesets[, c("gs_name", "gene_symbol")])KEGGREST需要额外处理才能用于GSVA：
library(GSEABase) kegg_genesets <- getGmt("kegg_hsa.gmt")
2. 实战操作全流程解析
2.1 使用msigdbr的标准化流程
对于快速分析场景，msigdbr提供开箱即用的体验：
安装与加载：
install.packages("msigdbr") library(msigdbr)获取基因集：
human_kegg <- msigdbr( species = "Homo sapiens", category = "C2", subcategory = "CP:KEGG" ) %>% dplyr::select(gs_name, gene_symbol)转换为GSVA所需格式：
kegg_list <- split(human_kegg$gene_symbol, human_kegg$gs_name)
2.2 KEGGREST的完整工作流
需要更全面数据时的解决方案：
library(KEGGREST) library(EnrichmentBrowser) # 获取所有人类通路ID pathways <- keggList("pathway", "hsa") pathway_ids <- names(pathways) # 逐条获取基因列表 kegg_genes <- lapply(pathway_ids, function(id) { pathway <- keggGet(id) genes <- pathway[[1]]$GENE if(!is.null(genes)) { gene_symbols <- genes[seq(1, length(genes), 2)] return(gene_symbols) } }) names(kegg_genes) <- gsub("path:", "", pathway_ids) kegg_genes <- kegg_genes[!sapply(kegg_genes, is.null)]注意：KEGGREST的API有访问频率限制（约每秒1次请求），建议添加Sys.sleep(1)避免被封禁
2.3 混合使用策略
结合两者优势的进阶方案：
# 先用msigdbr获取基础集合 base_kegg <- msigdbr(species = "hsa", category = "C2", subcategory = "KEGG") # 识别缺失的通路 all_kegg <- names(keggList("pathway", "hsa")) missing <- setdiff( gsub("path:", "", all_kegg), unique(base_kegg$gs_name) ) # 针对性补充缺失通路 supp_kegg <- lapply(missing, function(p) { pathway <- keggGet(paste0("hsa", p)) genes <- pathway[[1]]$GENE[seq(1, length(pathway[[1]]$GENE), 2)] data.frame(gs_name = p, gene_symbol = genes) })3. 与单细胞分析管道的整合
3.1 Seurat兼容性适配
将基因列表转换为单细胞分析友好格式：
# 创建GMT文件 writeGMT <- function(gene_sets, file) { conn <- file(file, "w") lapply(names(gene_sets), function(n) { writeLines( paste(n, "NA", paste(gene_sets[[n]], collapse = "\t")), conn ) }) close(conn) } # 在Seurat中应用 library(Seurat) seurat_obj <- AddModuleScore( object = seurat_obj, features = kegg_list[1:50], name = "KEGG_" )3.2 GSVA分析最佳实践
不同工具获取的基因集在GSVA中的表现差异：
| 参数 | msigdbr基因集 | KEGGREST基因集 |
|---|---|---|
| min.sz | 建议≥5 | 建议≥3 |
| max.sz | 可设200 | 可能需要调整到500 |
| 计算速度 | 较快（约10分钟） | 较慢（约25分钟） |
| 内存占用 | 中等 | 较高 |
# 基准测试代码 library(microbenchmark) microbenchmark( msigdbr = gsva(expr, msigdbr_sets, method="gsva"), keggrest = gsva(expr, keggrest_sets, method="gsva"), times = 3 )4. 常见问题排查指南
4.1 基因符号匹配失败
当出现"None of the genes in the gene sets are present in the expression matrix"警告时：
检查基因命名风格：
# 转换基因符号 library(HGNChelper) new_symbols <- checkGeneSymbols(rownames(expr))处理过时符号：
library(org.Hs.eg.db) mapIds( org.Hs.eg.db, keys = rownames(expr), column = "SYMBOL", keytype = "ENSEMBL" )
4.2 通路覆盖不全的补救措施
当关键通路缺失时，可交叉验证其他来源：
# 从Reactome获取补充 library(reactome.db) reactome_pathways <- as.list(reactomePATHID2EXTID)4.3 大规模分析的性能优化
处理大型单细胞数据集时：
# 并行化设置 library(BiocParallel) register(SnowParam(workers = 8)) # 内存优化 gsva_res <- gsva( expr, kegg_genes, mx.diff=TRUE, parallel.sz=8, BPPARAM=SnowParam(workers=8) )在实际项目中，我发现对心肌细胞单细胞数据的分析中，使用KEGGREST补充的"Cardiac muscle contraction"通路（hsa04260）能捕获更多功能差异。而将msigdbr生成的GMT文件与GSEA软件配合使用时，需要注意将基因符号全部转为大写字母以避免匹配失败。
