别再手动整理KEGG基因集了！用R包KEGGREST和msigdbr一键搞定357条通路（附完整代码）

告别手动整理：用R包高效获取KEGG通路基因集的完整指南

在单细胞转录组分析中，基因集富集分析（GSEA）和基因集变异分析（GSVA）是揭示生物学通路活性的重要工具。然而，许多研究者常常在第一步——获取和整理KEGG通路基因列表时就陷入困境。手动收集不仅耗时耗力，还容易出错，特别是在处理数百条通路时。本文将介绍如何利用R包KEGGREST和msigdbr自动化这一过程，显著提升研究效率。

1. 为什么需要自动化获取KEGG基因集？

手动整理KEGG通路基因列表存在几个主要痛点：

• 耗时严重：KEGG数据库包含数百条通路，手动收集每条通路的基因成员极其繁琐
• 更新滞后：手动整理的列表难以实时同步KEGG数据库的最新更新
• 易出错：人工操作容易引入基因ID匹配错误或遗漏
• 格式转换复杂：将原始数据转换为GSVA/GSEA所需格式需要额外处理

相比之下，使用R包自动化获取具有明显优势：

# 比较手动与自动化方法的时间成本 time_cost <- data.frame( Method = c("Manual", "Automated"), Time_hours = c(8, 0.5), Error_rate = c("High", "Low"), Update_frequency = c("Rarely", "On-demand") )

表1：手动与自动化方法对比

2. 两种主流R包获取KEGG基因集的方法对比

目前主要有两种通过R获取KEGG基因集的方法，它们在通路数量和基因标识符类型上有所不同。

2.1 使用msigdbr包获取KEGG基因集

msigdbr包提供了MSigDB数据库的接口，包含多个预定义的基因集，其中就包含KEGG通路。

# 安装并加载msigdbr包 # install.packages("msigdbr") library(msigdbr) # 获取人类KEGG基因集 human_kegg <- msigdbr(species = "Homo sapiens", category = "C2", subcategory = "CP:KEGG")

msigdbr的特点：

• 获取速度快，数据已预处理
• 但只包含186条核心KEGG通路
• 基因标识符为Symbol格式
• 适合快速分析和初步探索

2.2 使用KEGGREST包获取完整KEGG基因集

KEGGREST提供了直接访问KEGG数据库的接口，能获取更全面的通路信息。

# 安装KEGGREST if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install("KEGGREST") library(KEGGREST)

获取完整流程：

1. 列出所有人类通路
2. 下载每条通路的基因列表
3. 转换为适合分析的格式

# 获取人类所有KEGG通路ID pathways <- keggList("pathway", "hsa") pathway_ids <- names(pathways) # 获取第一条通路的基因列表示例 keggGet(pathway_ids[1])[[1]]$GENE

KEGGREST的优势：

• 获取357条完整KEGG通路
• 数据直接来自KEGG，更新及时
• 可灵活选择不同基因ID类型
• 适合需要全面分析的研究

3. 从原始数据到分析就绪格式的完整流程

获取原始基因列表只是第一步，将其转换为分析工具所需的格式同样重要。

3.1 使用EnrichmentBrowser处理KEGGREST结果

EnrichmentBrowser包提供了便捷的函数来处理KEGG基因集：

library(EnrichmentBrowser) # 获取并保存基因集 hsa_kegg_genesets <- getGenesets(org = "hsa", db = "kegg", gene.id.type = "SYMBOL", cache = TRUE) # 查看前5条通路 length(hsa_kegg_genesets) # 应返回357 names(hsa_kegg_genesets)[1:5]

3.2 转换为GMT格式

GMT格式是GSEA等工具的标准输入格式，转换方法如下：

# 写入GMT文件 writeGMT(hsa_kegg_genesets, gmt.file = "kegg_hsa.gmt") # 保存为R对象 save(hsa_kegg_genesets, file = "hsa_kegg_genesets.RData")

3.3 ID类型转换技巧

不同数据库可能使用不同的基因标识符，统一ID类型至关重要。

# 将ENTREZID转换为SYMBOL library(org.Hs.eg.db) library(clusterProfiler) # 示例转换 gene_ids <- c("10", "100", "1000") bitr(gene_ids, fromType = "ENTREZID", toType = "SYMBOL", OrgDb = org.Hs.eg.db)

表2：常见基因ID类型及转换方法

4. 实战应用：将KEGG基因集用于GSVA分析

有了格式正确的基因集，就可以进行下游分析了。以下是GSVA分析的完整示例。

4.1 准备表达矩阵和基因集

library(GSVA) # 加载之前保存的基因集 load("hsa_kegg_genesets.RData") # 假设expr是准备好的表达矩阵 # 行名为基因Symbol，列名为样本ID # expr <- readRDS("expression_matrix.rds") # 运行GSVA gsva_results <- gsva(expr, hsa_kegg_genesets, method = "gsva", kcdf = "Gaussian", parallel.sz = 4)

4.2 参数优化建议

• kcdf选择：

  • RNA-seq数据：建议使用"Poisson"
  • 微阵列数据：建议使用"Gaussian"

• 并行计算：

  • 设置parallel.sz参数利用多核加速
  • 大数据集可考虑parallel.sz=detectCores()-1

• 基因集大小过滤：

  • min.sz：排除基因数过少的通路（建议≥5）
  • max.sz：排除基因数过多的通路（建议≤500）

4.3 结果解读与可视化

GSVA结果是一个矩阵，行是通路，列是样本，值是通路活性分数。

# 简单可视化 library(pheatmap) pheatmap(gsva_results, show_rownames = FALSE, clustering_method = "complete")

对于差异分析，可以结合limma包：

library(limma) # 假设group是分组因子 design <- model.matrix(~0 + group) colnames(design) <- levels(group) fit <- lmFit(gsva_results, design) contrast.matrix <- makeContrasts(treatment-control, levels=design) fit2 <- contrasts.fit(fit, contrast.matrix) fit2 <- eBayes(fit2) # 获取差异通路 topTable(fit2, adjust="BH")

5. 常见问题与解决方案

在实际应用中，研究者常遇到以下几类问题：

5.1 物种匹配问题

• 问题表现：使用错误物种代码导致获取不到数据
• 解决方案：
  • 使用keggList("organism")查看所有可用物种
  • 人类代码为"hsa"，小鼠为"mmu"

# 查看所有KEGG支持的物种 kegg_organisms <- keggList("organism") head(kegg_organisms)

5.2 基因ID不一致

• 问题表现：基因集与表达矩阵基因ID类型不匹配
• 解决方案：
  • 统一使用SYMBOL或ENSEMBL等常用ID
  • 使用clusterProfiler进行ID转换

5.3 通路数量差异

• 问题表现：不同方法获取的通路数量不同
• 原因分析：
  • msigdbr提供186条核心通路
  • KEGGREST提供357条完整通路
• 选择建议：
  • 初步分析可使用msigdbr
  • 全面分析推荐KEGGREST

5.4 网络连接问题

• 问题表现：KEGGREST访问超时或失败
• 解决方案：
  • 设置超时时间：setOption(timeout=300)
  • 使用cache=TRUE缓存结果
  • 考虑非高峰时段访问

# 设置超时时间 options(timeout=300) # 使用缓存 getGenesets(org="hsa", db="kegg", cache=TRUE)

6. 进阶技巧与性能优化

对于大规模数据集或频繁分析的需求，以下技巧可以进一步提升效率。

6.1 本地缓存策略

频繁从KEGG数据库下载会增加时间成本和服务器负担，建议建立本地缓存。

# 设置本地缓存目录 cache_dir <- "kegg_cache" if(!dir.exists(cache_dir)) dir.create(cache_dir) # 自定义缓存函数 get_cached_kegg <- function(pathway_id, cache_dir) { cache_file <- file.path(cache_dir, paste0(pathway_id, ".rds")) if(file.exists(cache_file)) { return(readRDS(cache_file)) } else { result <- keggGet(pathway_id) saveRDS(result, cache_file) return(result) } }

6.2 并行获取通路信息

使用parallel包可以显著加速多条通路的获取过程。

library(parallel) # 设置并行集群 cl <- makeCluster(detectCores() - 1) clusterExport(cl, c("keggGet")) # 并行获取多条通路 pathway_ids <- names(keggList("pathway", "hsa"))[1:50] # 前50条通路 pathway_data <- parLapply(cl, pathway_ids, function(id) keggGet(id)) stopCluster(cl)

6.3 自定义基因集过滤

根据研究需求，可以针对性地筛选通路。

# 示例：筛选与癌症相关的通路 cancer_terms <- c("pathway in cancer", "signaling", "metabolism") cancer_pathways <- hsa_kegg_genesets[ grepl(paste(cancer_terms, collapse="|"), names(hsa_kegg_genesets), ignore.case=TRUE) ] length(cancer_pathways) # 查看筛选后的数量

6.4 与其他数据库整合

将KEGG通路与其他数据库（如GO）结合，可以获得更全面的生物学见解。

# 获取GO基因集 go_genesets <- getGenesets(org="hsa", db="go", onto="BP") # 合并KEGG和GO combined_genesets <- c(hsa_kegg_genesets, go_genesets) # 确保名称唯一 names(combined_genesets) <- make.unique(names(combined_genesets))

7. 从分析到发表：完整工作流建议

为确保研究的可重复性和高效性，建议遵循以下工作流程：

1. 数据获取阶段：

   • 使用KEGGREST获取最新完整通路
   • 保存原始数据和GMT文件
   • 记录获取日期和参数

2. 预处理阶段：

   • 统一基因ID类型
   • 过滤低质量通路
   • 保存处理后的R对象

3. 分析阶段：

   • 使用GSVA/GSEA进行富集分析
   • 保存中间结果和脚本
   • 记录所有参数设置

4. 可视化阶段：

   • 选择适当的可视化方法
   • 确保图形清晰可读
   • 保存高分辨率图片

5. 报告阶段：

   • 整理完整分析代码
   • 包含版本信息（R、包版本）
   • 提供示例数据

# 示例工作流文档 workflow_doc <- list( date = Sys.Date(), r_version = R.version.string, packages = installed.packages()[,c("Package", "Version")], parameters = list( organism = "hsa", gene_id_type = "SYMBOL", min_genes = 5, max_genes = 500 ), input_files = "expression_matrix.rds", output_files = c("kegg_genesets.RData", "gsva_results.rds") ) saveRDS(workflow_doc, "analysis_workflow.rds")

在实际项目中，这套自动化流程相比手动整理可以节省约90%的时间，同时显著降低错误率。一位长期从事单细胞分析的同行在采用这种方法后，单在数据准备阶段每周就能节省10-15小时，使得可以将更多精力投入到结果解读和生物学发现上。