1. 项目概述为什么需要批量下载Reactome通路基因集如果你正在做生物信息学分析特别是涉及到通路富集分析、基因功能注释或者多组学数据整合那么Reactome这个数据库对你来说一定不陌生。它是一个高质量、手工注释的通路知识库涵盖了从代谢、信号传导到疾病相关的各种生物学过程。很多时候我们手头有一份差异表达基因列表或者是从单细胞测序中鉴定出的细胞类型特异性基因下一步就是想看看这些基因在哪些通路上富集。这时候你可能会打开网页在Reactome官网一个通路一个通路地查基因列表或者用一些R包比如ReactomePA在线获取。但你想过没有如果网络不稳定或者你需要一次性分析成百上千个基因集甚至想本地构建一个通路基因库供后续反复查询这种“现用现查”的方式就非常低效了。这就是“下载所有通路基因集”这个需求的由来。它不是一个简单的数据搬运而是一个提升分析效率、确保分析可重复性、并深入理解数据背后生物学意义的基础性工作。通过本地化存储Reactome的全套通路-基因对应关系你可以实现离线分析摆脱网络依赖在封闭环境或计算集群上也能快速进行富集分析。批量处理一次性完成大量基因列表的富集分析效率远超交互式查询。自定义分析基于本地的基因集文件你可以轻松地与其他数据库如GO、KEGG的基因集进行合并、比较或自定义筛选构建更符合你研究背景的注释集。版本控制固定使用某一特定版本的Reactome数据确保多年后你的分析结果依然可以精确复现。简单来说这就像把一本厚厚的、随时可能更新的“通路百科全书”下载到自己的书房里随时翻阅、批注而不必每次都跑去图书馆。2. 核心思路与方案选型要实现这个目标核心思路很清晰从Reactome官方数据源获取结构化的通路与基因对应关系然后整理成方便程序尤其是R语言读取的格式。这里有几个关键的技术选型点直接决定了后续操作的复杂度和结果的可用性。2.1 数据源选择官网下载 vs. 生物信息学R包Reactome的数据主要通过两种方式公开官方FTP/下载页面Reactome官网提供了多种格式的数据下载包括MySQL数据库dump文件、BioPAX格式、SBML格式等。这是最原始、最完整的数据源。生物信息学R包例如reactome.dbBioconductor项目和ReactomePA。这些包通常内置了某个时间点的通路注释数据或者提供了API来在线查询。为什么我推荐从官方数据源下载虽然R包用起来方便一句library(reactome.db)就行但它有几个局限数据版本滞后Bioconductor的发布周期导致包内数据可能不是最新的。数据内容受限R包通常只提供核心的基因-通路映射关系而官方数据源包含更丰富的元信息如通路层级结构、反应细节、参与分子类型蛋白质、复合物、小分子等。灵活性差你无法定制化地提取信息比如只提取人类Homo sapiens的通路或者同时提取小鼠的直系同源基因信息。因此为了获得最全面、最新且可定制的数据直接从Reactome官网下载是最佳选择。本项目将聚焦于处理官网提供的MySQL数据库dump文件和GMT格式文件这两种最实用的数据源。2.2 输出格式设计GMT文件是金标准下载和解析的最终目的是生成一个可用的基因集文件。在生物信息学领域GMT (Gene Matrix Transposed) 格式是基因集富集分析的事实标准。它被GSEA、clusterProfiler等主流工具广泛支持。一个GMT文件是纯文本格式每一行代表一个基因集通路其结构如下通路标识符(Tab)通路描述(Tab)基因1(Tab)基因2(Tab)基因3(Tab)...例如一条人类细胞周期通路的记录可能看起来像R-HSA-1640170 Cell Cycle CDC20 CDK1 CCNB1 PLK1 ...我们的核心任务就是将Reactome中成千上万个通路以及每个通路对应的基因列表通常使用NCBI Gene ID或Ensembl Gene ID整理成这样的GMT文件。有了这个文件后续使用clusterProfiler的enricher()函数或GSEA软件进行分析就变得轻而易举。2.3 技术路线规划基于以上选择我规划了两条主要的技术路线适合不同技术背景的从业者路线一基于预处理的GMT文件推荐给大多数R用户这是最快捷的路径。Reactome官网直接提供了按物种划分的GMT格式文件。我们只需要下载、解压并可能进行简单的ID转换如从Ensembl ID转成常用的Symbol即可投入使用。这条路线省时省力适合快速启动分析。路线二基于MySQL数据库Dump文件适合需要深度定制或最新数据的用户这条路线更底层也更强大。我们需要下载Reactome的MySQL数据库备份文件在本地或服务器上恢复成一个临时数据库然后执行SQL查询精确地提取我们需要的通路-基因关系并生成GMT文件。这个过程虽然步骤较多但让你能完全掌控数据的提取逻辑例如可以关联查询通路的层级结构、筛选特定类型的参与者等。接下来我将详细拆解这两条路线的实操步骤。3. 实操详解路线一基于官网GMT文件这条路线是“开箱即用”的典范特别适合希望快速获得结果的分析者。3.1 获取数据文件访问Reactome下载页面打开浏览器访问 Reactome 的官方下载页面通常路径为https://reactome.org/download-data。定位GMT文件在页面中寻找名为 “Gene sets” 或 “Pathway GMT files” 的板块。Reactome 为多个物种提供了GMT文件。选择物种和标识符类型你会看到类似以下的文件列表ReactomePathways.gmt.zip(可能包含所有物种)ReactomePathways_Homo_sapiens_Ensembl.gmt.zip(人类Ensembl Gene ID)ReactomePathways_Mus_musculus_Ensembl.gmt.zip(小鼠Ensembl Gene ID)ReactomePathways_Homo_sapiens_Symbol.gmt.zip(人类基因Symbol)选择建议对于大多数基于RNA-seq差异表达基因的分析基因Symbol是最直观、最常用的标识符。因此直接下载ReactomePathways_Homo_sapiens_Symbol.gmt.zip通常是最佳选择。如果你需要与其他使用Ensembl ID的数据集整合则选择对应的Ensembl版本。3.2 在R中加载与初步探索下载的ZIP文件解压后得到一个.gmt文件。在R中我们可以使用clusterProfiler包提供的read.gmt()函数轻松读取。# 安装并加载必要包 if (!requireNamespace(clusterProfiler, quietly TRUE)) BiocManager::install(clusterProfiler) library(clusterProfiler) # 读取GMT文件 gmt_file - ~/Downloads/ReactomePathways_Homo_sapiens_Symbol.gmt reactome_gmt - read.gmt(gmt_file) # 查看数据结构 head(reactome_gmt) # 输出通常有两列ont通路ID/名称和 gene基因Symbol # ont gene # 1 R-HSA-109581 ABL1 # 2 R-HSA-109581 ABL2 # 3 R-HSA-109581 CRK # 4 R-HSA-109581 CRKL # 5 R-HSA-109581 FYN # 6 R-HSA-109581 GRB2 # 查看共有多少条通路唯一的ont项 pathway_count - length(unique(reactome_gmt$ont)) cat(sprintf(总共下载了 %d 条Reactome通路基因集。\n, pathway_count)) # 查看某个特定通路的基因 library(dplyr) apoptosis_genes - reactome_gmt %% filter(ont R-HSA-109581) %% # 以“Apoptosis”通路为例 pull(gene) head(apoptosis_genes)注意read.gmt()读入的是长格式数据框。对于clusterProfiler::enricher()函数这本身就是理想的输入格式。如果你需要将其转换为列表格式通路名称为名基因向量为值可以这样做reactome_list - split(reactome_gmt$gene, reactome_gmt$ont)3.3 关键步骤基因标识符检查与转换这是最容易出错的一步。Reactome提供的Symbol文件使用的是官方基因名来自HUGO Gene Nomenclature Committee, HGNC。而你的表达矩阵或差异基因列表中的基因名可能存在以下问题别名如TP53与P53。过时符号一些旧的基因名可能已被更新。大小写不一致RNA-seq流程输出的可能是全大写而GMT文件中可能是首字母大写。强烈建议进行基因名匹配检查# 假设你的基因列表是 deg_genes deg_genes - c(TP53, MYC, BRCA1, SOME_GENE) # 检查有多少基因能在Reactome基因集中找到 matched_genes - deg_genes[deg_genes %in% reactome_gmt$gene] unmatched_genes - setdiff(deg_genes, reactome_gmt$gene) cat(sprintf(匹配到的基因数: %d\n, length(matched_genes))) cat(sprintf(未匹配的基因数: %d\n, length(unmatched_genes))) cat(未匹配的基因:, paste(unmatched_genes, collapse, ), \n)如果未匹配的基因较多你需要进行基因标识符转换。可以使用biomaRt或AnnotationDbi包配合org.Hs.eg.db等物种注释包。# 使用 org.Hs.eg.db 进行Symbol的同步和转换 if (!requireNamespace(org.Hs.eg.db, quietly TRUE)) BiocManager::install(org.Hs.eg.db) library(org.Hs.eg.db) # 将你的基因Symbol转换为最新的官方Symbol # 首先获取所有已知的别名和官方Symbol的映射 # 这里假设你的基因名是“ALIAS”类型尝试映射到“SYMBOL” map - select(org.Hs.eg.db, keys unmatched_genes, columns c(SYMBOL, ALIAS), keytype ALIAS) # 查看映射结果 print(map) # 根据映射结果更新你的 deg_genes 列表实操心得在实际项目中我通常会准备一个“基因名清洗”的预处理脚本。这个脚本利用org.Hs.eg.db将所有输入基因名尽可能统一为官方Symbol并记录下无法映射的基因。这能极大提高后续富集分析的基因检出率。4. 实操详解路线二基于MySQL Dump文件当预生成的GMT文件不满足需求例如需要特定版本、需要包含非编码基因、需要提取通路层级信息时这条路线的价值就凸显出来了。4.1 环境准备与数据下载安装MySQL或MariaDB在你的本地机器或服务器上安装一个MySQL数据库服务。对于快速测试使用Docker容器是极佳的选择因为它不会污染本地环境用完即删。# 拉取MySQL官方镜像并运行一个临时容器 docker run --name reactome_db -e MYSQL_ROOT_PASSWORDyourpassword -d mysql:8.0 # 进入容器内的bash环境 docker exec -it reactome_db bash在容器内你可以使用mysql客户端。下载数据库Dump文件回到Reactome下载页面找到 “MySQL database” 部分下载最新版本的数据库dump文件通常是一个巨大的.sql.gz文件如reactome_yourversion.sql.gz。4.2 数据库恢复与探索将数据文件复制到容器内如果使用Docker# 在宿主机执行将下载的sql.gz文件复制到容器中 docker cp ~/Downloads/reactome_yourversion.sql.gz reactome_db:/tmp/解压并恢复数据库# 在容器内执行 cd /tmp gunzip reactome_yourversion.sql.gz mysql -uroot -pyourpassword reactome_yourversion.sql这个过程可能会持续几分钟到十几分钟取决于文件大小和机器性能。连接数据库并探索关键表mysql -uroot -pyourpassword reactome进入MySQL命令行后查看有哪些表SHOW TABLES;对我们最重要的几张表是Pathway存储通路的基本信息ID、名称、物种等。DatabaseObject所有对象的基类信息。ReferenceGeneProduct关联基因产物蛋白质和外部数据库如NCBI Gene, Ensembl的ID。Pathway_2_hasEvent/ReactionLikeEvent_2_hasEvent描述通路和反应/子通路之间的层级关系。ReactionLikeEvent_2_input/ReactionLikeEvent_2_output/ReactionLikeEvent_2_requiredInputComponent描述反应的事件参与者。4.3 编写SQL查询提取通路-基因关系这是最核心的一步。我们需要编写一个SQL查询将通路与其对应的基因这里以NCBI Gene ID为例关联起来。思路是通路(Pathway) - 包含的反应/事件(ReactionLikeEvent) - 事件的参与者(PhysicalEntity) - 对应的基因产物(ReferenceGeneProduct) - NCBI Gene ID。以下是一个经过简化的查询示例用于提取人类speciesId 48887通路的基因SELECT DISTINCT p.stId AS pathway_id, p.displayName AS pathway_name, rgp.identifier AS ncbi_gene_id FROM Pathway p -- 连接通路与其包含的直接事件反应或子通路 JOIN Pathway_2_hasEvent phe ON p.DB_ID phe.DB_ID JOIN ReactionLikeEvent rle ON phe.hasEvent rle.DB_ID -- 连接反应与其输入的物理实体参与者 JOIN ReactionLikeEvent_2_input rlei ON rle.DB_ID rlei.DB_ID JOIN PhysicalEntity pe ON rlei.input pe.DB_ID -- 连接物理实体到其对应的基因产物如果是蛋白质 JOIN PhysicalEntity_2_hasComponent phc ON pe.DB_ID phc.DB_ID JOIN ReferenceGeneProduct rgp ON phc.hasComponent rgp.DB_ID -- 筛选物种为人类且标识符类型为NCBI Gene的记录 WHERE p.speciesId 48887 AND rgp.identifier REGEXP ^[0-9]$ -- 简单筛选数字型的NCBI Gene ID AND rgp.identifierDatabaseName NCBI Gene -- 按通路和基因排序 ORDER BY p.stId, rgp.identifier;请注意实际的Reactome数据模型非常复杂上述查询是一个高度简化的版本可能无法涵盖所有情况例如复合物、多聚体、非蛋白质分子。一个更健壮的方案可能需要递归查询来处理通路的层级嵌套并联合查询_2_output和_2_requiredInputComponent等表。对于生产环境建议仔细研究Reactome的数据模型文档。4.4 在R中连接数据库并生成GMT文件你可以直接在MySQL命令行中将查询结果导出为CSV也可以在R中使用DBI和RMySQL/RMariaDB包来连接数据库并处理数据。# 安装必要的R包 install.packages(c(DBI, RMariaDB)) library(DBI) library(RMariaDB) # 创建数据库连接假设MySQL运行在本地 con - dbConnect(MariaDB(), host localhost, port 3306, user root, password yourpassword, dbname reactome) # 执行我们编写好的复杂查询这里用简化版示例 query - SELECT ... -- 将上面完整的SQL语句粘贴在这里 pathway_gene_df - dbGetQuery(con, query) # 断开连接 dbDisconnect(con) # 查看获取的数据 head(pathway_gene_df) # pathway_id pathway_name ncbi_gene_id # 1 R-HSA-109581 Intrinsic Pathway for Apoptosis 25 # 2 R-HSA-109581 Intrinsic Pathway for Apoptosis 100 # ... ... # 将数据转换为GMT格式所需的列表 library(dplyr) library(tidyr) # 首先将NCBI Gene ID转换为基因Symbol如果需要 # 这里需要用到 org.Hs.eg.db 进行ID转换 if (!requireNamespace(org.Hs.eg.db, quietly TRUE)) BiocManager::install(org.Hs.eg.db) library(org.Hs.eg.db) # 映射NCBI Gene ID 到 Symbol gene_ids - unique(pathway_gene_df$ncbi_gene_id) id_map - select(org.Hs.eg.db, keys as.character(gene_ids), columns SYMBOL, keytype ENTREZID) # 注意可能会有一些ENTREZID映射不到SYMBOL或者一对多的情况需要处理 # 合并映射关系 pathway_gene_df - pathway_gene_df %% left_join(id_map, by c(ncbi_gene_id ENTREZID)) %% filter(!is.na(SYMBOL)) # 过滤掉没有Symbol的基因 # 聚合生成GMT列表格式 gmt_list - pathway_gene_df %% group_by(pathway_id, pathway_name) %% summarise(genes paste(unique(SYMBOL), collapse \t), .groups drop) %% mutate(gmt_line paste(pathway_id, pathway_name, genes, sep \t)) # 写入GMT文件 writeLines(gmt_list$gmt_line, Reactome_Custom_Human_Symbol.gmt)重要提示直接从数据库提取并生成GMT文件是一个复杂的过程上述代码示例是一个概念性框架。你需要根据实际的Reactome数据库表结构和你的具体需求如是否包含复合物成员、是否包含小分子来调整SQL查询和后续的R处理逻辑。这可能涉及多次迭代和验证。5. 质量验证与常见问题排查无论采用哪种路线获得GMT文件后都必须进行质量验证这是保证后续分析可靠性的关键一步。5.1 基础完整性检查文件格式检查用文本编辑器或head/tail命令查看GMT文件开头和结尾几行确保是标准的制表符分隔格式且每行以通路ID开头。通路数量核对与Reactome官网公布的该物种通路总数进行粗略比对。例如人类通路大约有2500多条。如果数量差异巨大如少了一个数量级说明提取过程可能有误。基因数量与唯一性统计总共有多少唯一的基因Symbol。对于人类Reactome覆盖的基因通常在10000-15000个左右。如果只有几千个可能丢失了大量基因。5.2 逻辑一致性检查选取已知通路验证挑选几个你熟悉的经典通路如“Cell Cycle (R-HSA-1640170)”、“Apoptosis (R-HSA-109581)”检查你的GMT文件中这些通路下的基因列表是否合理。可以去Reactome官网搜索该通路对比其“Proteins”列表。检查基因重复同一个基因出现在同一个通路下应该是唯一的。检查并确保每个通路内的基因列表没有重复项。检查空通路确保没有基因列表为空的通路行。5.3 常见问题与解决方案速查表问题现象可能原因排查与解决方案GMT文件行数远少于预期1. SQL查询条件过严丢失了数据。2. 从官网下载的GMT文件本身就不全可能性低。3. ID转换过程丢失了大量基因。1.简化查询先写一个最简单的查询只关联Pathway和ReferenceGeneProduct的核心表确保能取出大量数据再逐步添加关联和筛选条件。2.检查ID映射在R中检查org.Hs.eg.db映射前后基因数量的变化处理一对多或映射失败的情况。富集分析时很多基因匹配不上1. 基因标识符不统一如Symbol vs. Ensembl ID。2. GMT文件使用的基因名版本过旧。3. 你的基因列表包含了很多新发现的基因或非编码RNAReactome未收录。1.统一ID确保你的差异基因列表与GMT文件使用同一种ID系统。优先使用官方基因Symbol。2.更新GMT下载最新版本的Reactome GMT文件。3.接受局限性了解Reactome主要关注蛋白质编码基因和已知通路未匹配是正常的。从数据库提取时SQL查询报错或极慢1. 表连接错误导致笛卡尔积数据量爆炸。2. 缺少关键索引。3. 查询逻辑过于复杂。1.检查JOIN条件确保每个JOIN都有明确的ON条件。2.使用EXPLAIN在MySQL中使用EXPLAIN命令分析查询执行计划看是否用上了索引。3.分步查询将复杂查询拆分成几个临时表分步执行最后再整合。生成的GMT文件通路名称乱码或包含特殊字符数据库中的文本字段包含非ASCII字符如希腊字母α、β或制表符、换行符。1.导出时指定编码在SQL查询中使用CONVERT(column USING utf8mb4)。2.在R中清洗使用stringr::str_remove_all或gsub移除控制字符。gmt_line - gsub(“[\r\n\t]”, “ “, gmt_line)注意保留作为分隔符的\t。5.4 高级技巧构建带层级信息的通路集有时我们不仅需要通路和基因的对应关系还希望知道通路之间的上下级关系例如“细胞周期”通路下包含“有丝分裂”子通路。这在进行层级富集分析或可视化时非常有用。你可以通过查询Pathway_2_hasEvent表来获取这种关系。一个简单的递归查询或是在R中用循环处理可以构建出通路的父子关系树。然后你可以选择扁平化将所有下层子通路的基因也合并到上层父通路中。分层存储将层级关系单独存储为一个数据框在富集分析后用enrichplot或ggplot2绘制树状图。-- 示例获取人类通路的直接父子关系非递归 SELECT parent.stId AS parent_id, parent.displayName AS parent_name, child.stId AS child_id, child.displayName AS child_name FROM Pathway parent JOIN Pathway_2_hasEvent phe ON parent.DB_ID phe.DB_ID JOIN Pathway child ON phe.hasEvent child.DB_ID WHERE parent.speciesId 48887 AND child.speciesId 48887;在R中处理这种层级数据data.tree或tidyr的嵌套功能会很有帮助。6. 集成到分析流程与自动化脚本手动操作一次是可行的但最好的实践是将这个过程脚本化、自动化。这样当Reactome数据库更新时你可以一键重新生成最新的基因集。我建议创建一个R脚本例如update_reactome_gmt.R将路线一或路线二的核心步骤封装进去。脚本可以包含以下功能自动下载使用download.file()从Reactome官网获取最新的GMT或SQL dump文件。自动处理执行数据提取、清洗、ID转换和GMT文件生成。版本记录在生成的GMT文件头或单独的文件中记录Reactome数据版本号、生成日期和脚本版本。完整性检查内置上一节提到的质量验证步骤如果检查不通过则报错并停止。对于路线二你甚至可以编写一个Dockerfile构建一个包含MySQL和自动恢复、查询、导出功能的镜像实现完全自动化的流水线。最终将这个脚本或流水线集成到你的项目分析流程最开始的部分。确保每个使用这些通路基因集的分析报告都明确标注了所使用的Reactome数据版本和基因集生成方法这极大地提升了研究的可重复性。本地化Reactome通路基因集看似是一个数据准备的“脏活累活”但它奠定了后续所有富集分析结果的可靠性与效率基础。花时间搭建好这个基础设施之后的分析工作就会顺畅得多。从我自己的经验来看在项目初期投入时间做好这类基础数据工程远比在分析中途被不匹配的基因名、过时的通路或网络超时问题打断要划算得多。