Reactome通路基因集本地化：从GMT文件到MySQL数据库的完整方案-开发者社区

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.db（Bioconductor项目）和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_PASSWORD=yourpassword -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. 你的基因列表包含了很多新发现的基因或非编码RNA，Reactome未收录。	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通路基因集，看似是一个数据准备的“脏活累活”，但它奠定了后续所有富集分析结果的可靠性与效率基础。花时间搭建好这个基础设施，之后的分析工作就会顺畅得多。从我自己的经验来看，在项目初期投入时间做好这类基础数据工程，远比在分析中途被不匹配的基因名、过时的通路或网络超时问题打断要划算得多。