本研究介绍了一种运行Pathway2Targets算法的实验方案，该算法是一个R语言脚本，可根据批量RNA测序实验中病例组与对照组样本间细胞内信号通路谱的比较结果，预测并筛选潜在的治疗靶点。
需要JoVE订阅才能观看此内容。 请登录或开始免费试用
本文内容
摘要
摘要
本方案概述了一个多步骤的计算流程，用于从RNA测序数据中识别潜在的治疗靶点，包括相关软件的安装、环境配置验证，以及使用edgeR进行差异表达分析。随后，我们展示了如何利用信号通路影响分析（SPIA）算法来预测具有统计学意义的通路。为确保结果的可靠性，我们重点关注显著性通路（p < 0.05），以减少假阳性结果。与传统的基因集不同，这些通路反映了蛋白质-蛋白质相互作用网络，能够为细胞周期、免疫应答和代谢等细胞过程提供机制性见解。随后，利用Pathway2Targets算法对这些通路进行分析，该算法通过应用程序编程接口（API）与OpenTargets.org数据库对接。该算法采用一种新颖的加权方法，对已知药物靶点在所识别通路中的重要性进行评分，并实时提供分析进度。运行时间取决于通路的复杂性和靶点密度。输出结果包含两个排序文件：第一个文件列出了预测的药物靶点及其加权评分，第二个文件则包含相关治疗药物的多种详细信息。该流程整体上有助于在疾病特异性基因表达谱背景下，对可成药靶点和治疗方案进行优先排序。
引言
大规模RNA测序能够比较病例细胞群体与对照细胞群体中数千个基因的表达水平。实验通常设计为至少包含三个重复样本，理想情况下为生物学重复，但技术重复亦可满足要求。该设计可考虑生物学变异性，并降低离群样本的影响。对这些表达模式的分析可深入了解目标疾病对正常细胞过程的影响，并可能有助于预测相关治疗药物。
批量RNA测序数据的预处理通常包括：对测序读段进行质量控制（重复序列、测序接头、GC含量等）、读段修剪与接头去除、读段比对/定量1,2,3，以及差异表达分析4,5,6。幸运的是，已有多种分析流程实现了自动化，以减少这些步骤中所需的手工操作7,8,9。预处理完成后，常见的下游分析通常包括基因本体功能富集分析、信号通路富集分析以及剪接变异分析。这些下游分析能够对差异表达结果进行更高层次的归纳，相较于单纯的基因列表，更有利于结果的解释。
已开发出多种工具，旨在将现有治疗药物重新用于具有明确定义的疾病类型 或亚型。这些工具通过针对目标疾病的多组学数据类型进行算法训练来实现药物重定位。然而，此类旨在提高特定疾病中特异性和敏感性的努力，往往导致这些工具在更广泛的应用场景中表现欠佳10,11。另一类工具则具有更广泛的适用性，可用于将基因表达谱与已有的基因表达特征谱12,13，或与当前治疗药物的量化效应进行匹配14,15。然而，这些适用范围更广的工具通常在多种疾病中表现出较低的特异性和敏感性，和/或其训练所用数据已过时。
相比之下，Pathway2Targets 算法此前已被用于预测B细胞淋巴瘤、牙周炎、雌激素受体阳性乳腺癌、三阴性乳腺癌以及基孔肯雅病毒的潜在治疗靶点16,17,18,19,20,21。这些研究结果表明，该工具能够预测出稳健且具有生物学意义的靶点。值得注意的是，Pathway2Targets 预测出三阴性乳腺癌的392个潜在药物靶点，其中60个已进入临床试验；同时还预测出针对三阴性乳腺癌的828种药物，其中37种已进入临床试验阶段17。在淋巴瘤的研究中，该算法共预测出915种药物，其中461种已获美国食品药品监督管理局（FDA）批准19。
本研究的目的是描述一种计算协议，旨在为更多研究人员提供更详尽的命令行程序操作指导，使其能够有效使用近期开发的Pathway2Targets算法 （图1）。Pathway2Targets通过整合差异表达数据、基因-疾病关联、临床试验信息、公共靶点数据22、通路信息及其他指标，预测特定条件下的潜在靶点。重要的是，该算法采用一种独特且可自定义的加权方案，允许用户确定其分析中希望重点强调的约20个与靶点相关的指标，例如疾病关联数量、信号通路数量、独特药物数量、各阶段临床试验中治疗药物的数量等23。作为本协议的一个应用示例，我们将重新分析一个已有的结直肠癌数据集24。
访问受限。请登录或开始试用以查看此内容。
方案
本研究中分析的批量RNA测序数据来自公开可用的数据库（NCBI基因表达综合数据库和序列读取档案）25,26。因此，原始数据收集者确保了这些样本来自知情并同意的受试人类个体，且采集过程符合伦理规范。
1. 下载并安装 R 软件
- 安装 R（版本 4.0 或更高版本），方法是从综合 R 存档网络（CRAN）https://cran.r-project.org/mirrors.html 选择合适的链接并点击 选择选项 0-Cloud，然后根据计算机的操作系统遵循相应的说明。此过程通常需要 5-10 分钟。
- 从 https://posit.co/download/rstudio-desktop/ 安装 R Studio（版本 2024 或更高版本），然后按照下载页面上的说明进行操作。安装 RStudio 通常需要约 10 分钟。
注意：安装 RStudio 是可选操作（但强烈推荐），因为它提供了一个集成开发环境，可简化代码执行、提供语法高亮、便于包管理，并帮助用户可视化输出结果，这对于不熟悉 R 的用户尤为有益。
2. 下载并安装相关工具的 R 脚本
- 从 GitHub 仓库 https://github.com/bpickett/Pathway2Targets 下载以下必要的 R 脚本。此 URL 仅作参考。
- 使用以下链接下载脚本：SPIA 版本 1.0（点击 Download Raw File 下载）：https://github.com/bpickett/Pathway2Targets/blob/main/SPIA_Code.Rmd（提交 ID：60fcd46）；Pathway2Targets 版本 3.1（点击 Download Raw File 下载）：https://github.com/bpickett/Pathway2Targets/blob/main/Pathway2Targets.R（提交 ID：8e4c7c8）
3. 下载相关工具的 R 软件包
- 在运行 R（无论是在 RStudio 还是终端窗口中）时，输入以下命令以下载并安装运行该软件所需的额外 R 库。
- 启动 RStudio 程序。默认情况下，控制台面板位于 RStudio 的左下角。单击控制台面板窗口中的任意位置，输入光标应在箭头 ">" 符号后出现在底部。
- 将以下命令复制并粘贴到控制台区域，然后按 Enter 键：
install.packages(c("RCurl", "stringr", "jsonlite", "httr")).
安装成功后，将出现一条状态消息，内容为 "The downloaded binary packages are in ...."。 - 将以下命令复制并粘贴到控制台面板区域，然后按 Enter 键。
BiocManager::install(c("SummarizedExperiment", "EnrichmentBrowser", "biomaRt", "org.Hs.eg.db")).
安装成功后，将出现类似的消息，内容为 "The downloaded binary packages are in ...."。
注意：第一组库（步骤 3.1.1）为常规 R 库，而第二组（2.1.2）为 BioConductor 库。因此，这些命令需分别输入。这些库的适当版本将根据计算机上安装的 R 版本自动下载。每个命令的执行时间约为 5 分钟。
4. 文件处理
- 从本地计算机下载先前由 ARMOR 软件（或类似软件）生成的 edgeR 差异表达输出文件（RDS 格式）。该文件的名称通常为 edgeR_dge.rds。
- 手动审阅 edgeR（或类似的差异表达分析）结果，以开始对结果进行具有生物学意义的解释。为此，至少应根据校正后的 p 值 < 0.05 进行筛选，并可进一步根据 log2 倍数变化值的绝对值 > 1.5 进行筛选。若要试用该软件，可在 Zenodo 上获取示例 edgeR_dge.rds 文件，地址如下：https://doi.org/10.5281/zenodo.15186609
注意：审阅经上述筛选后保留的基因列表，有助于初步解释与病例样本（相对于对照样本）相关表型背后的分子机制。需要认识到，以无偏倚的方式解读基因列表极为困难，原因在于人们能够快速回忆的基因符号数量相对较少。因此，信号通路分析是一种有效的方法，可基于这些筛选出的基因在细胞内相互作用和/或通讯的方式对其进行归纳总结。 - 批量 RNA-seq 数据的预处理可能需要数小时至数天的计算时间，具体取决于所分析数据集的大小。请将此 .rds 文件保存在计算机的“下载”文件夹中。请注意，该 .rds 文件格式不可直接被人阅读。
5. 运行 SPIA 通路富集算法
- 若使用 R 运行，请从 GitHub 或 补充代码文件 1 中获取 R 脚本。假设 edgeR_dge.rds 文件位于 Downloads 文件夹中，输入以下命令：
Rscript --vanilla SPIA_Code.Rmd ~/Downloads/edgeR_dge.rds- 如果 edgeR_dge.rds 文件位于其他文件夹（或目录）中，请将上述命令替换为以下内容：
Rscript --vanilla SPIA_Code.Rmd <edgeR_dge.rds 文件的路径>
- 如果 edgeR_dge.rds 文件位于其他文件夹（或目录）中，请将上述命令替换为以下内容：
- 若使用 RStudio 运行，请从 GitHub 或 补充代码文件 2 中获取 R 脚本。
注意：在 R 语言中，在代码行前添加井号 # 符号可临时禁用该行代码。这些脚本最初设计用于命令行环境而非 RStudio。启用或禁用特定代码行是重新配置输入文件设置的最简便方法。- 在 RStudio 中，通过点击“文件”菜单中的 打开文件 选项，然后选择脚本名称，打开 SPIA_Code.Rmd 脚本。默认情况下，RStudio 的代码窗口位于左上方面板。
- 选中文件中的所有代码行，然后点击 运行 按钮（或 运行所选行 按钮），该按钮位于代码窗口的右上方。成功运行后，将在下载目录中生成一个类似以下名称的文件：
edgeR_dge.rds-TreatmentTumor-TreatmentNativeTissue_2025-04-23_
10-56-45.12767_SPIA_Results.csv
该文件将包含具有统计学显著性的信号通路。 - 通过电子表格程序手动查看包含显著性结果的文件。该文件内容有助于总结由差异表达基因显著代表的细胞内信号级联反应。
注意：显著通路的计算完成时间可能从约 30 分钟到数小时不等，具体取决于所分析数据集中信号的强度。程序运行期间，控制台窗口将持续更新实时进度信息。频繁更新的消息表明程序正在正常运行。有关此步骤更详细的说明，请参见 GitHub 仓库：https://github.com/bpickett/Pathway2Targets/tree/main - 若使用示例输入文件，此步骤的输出文件将位于 Downloads 文件夹中，文件名为：
edgeR_dge.rds-TreatmentTumor-TreatmentNativeTissue_"时间戳"_SPIA_Results.csv
此命名方式反映了输入文件名、执行的处理过程及输出结果，有助于在处理多个文件时避免混淆。 - 根据以下说明调整其他参数。
- 本库中 SPIA 算法的默认参数为 1,000 次置换。可将置换次数增加至 2,000 次以提高结果的可信度。通过修改脚本第 84 和 85 行中的 perm = 2000 为所需置换次数来调整置换数。其他参数的调整方法如下所述。
- 通过删除第 84 和 85 行中的 padj.method = 'BH' 来调整 p 值校正方法。这将不对 p 值进行校正，可能导致假阳性结果增加。
6. 在 SPIA 输出结果上运行 Pathway2Targets 目标优先级排序算法
- 若使用 R 运行，请从 GitHub 或 补充代码文件 3 中获取 R 脚本。使用以下命令调用该算法：
Rscript --vanilla Pathway2Targets.R - 若使用 R Studio 运行，请从 GitHub 或 补充代码文件 4 中获取 R 脚本。通过点击“文件”菜单中的 打开文件 选项，然后选择脚本名称，将 Pathway2Targets.R 脚本在 R Studio 中打开。
- 在 RStudio 代码窗口（左上方面板）中，将第 22 行的文件名替换为 SPIA 结果文件的名称，例如（来自示例数据）：
infile <- "edgeR_dge.rds-TreatmentTumor-TreatmentNativeTissue_2025-04-23_
10-56-45.12767_SPIA_Results.csv" - 选中文件中的所有代码行，然后点击代码窗口上方偏右侧的 运行 按钮。实时进度状态信息将持续显示在右下面板中。成功运行后，将在下载目录中生成一个文件，名称类似于：
"edgeR_dge.rds-TreatmentTumor-TreatmentNativeTissue_2025-04-23_。文件的命名方式反映了输入、处理过程和输出
10-56-45.12767_SPIA_Results.csv-RankedTargets.tsv"
edgeR_dge.rds-TreatmentTumor-TreatmentNativeTissue_2025-04-23
_10-56-45.12767_SPIA_Results.csv
即为此类输出文件的输入。
注意：此步骤可能需要一小时至数小时不等，具体时间取决于具有显著 p 值的信号通路数量、这些显著通路中包含的基因产物数量，以及已知为药物靶点的基因产物数量。
- 在 RStudio 代码窗口（左上方面板）中，将第 22 行的文件名替换为 SPIA 结果文件的名称，例如（来自示例数据）：
- Pathway2Targets 算法的某些参数可进行调整。具体而言，可调整脚本中第 31–38 行的乘数数值，以自定义各指标的加权方案。有关此步骤更详细的说明，请参见对应的 GitHub 仓库：https://github.com/bpickett/Pathway2Targets/tree/main
注意：作为参考，在示例中，从包含数百个独立靶点的 132 条通路中识别潜在药物，大约需要 2 小时完成。基于此指标，可合理估算总计算时间。
7. 打开优先目标和治疗方案的结果文件
- 将生成包含优先级排序的靶点及其相关指标的文件。对于示例输入文件，请使用“下载”文件夹中的输出文件，文件名为：
edgeR_dge.rds-TreatmentTumor-TreatmentNativeTissue_"timestamp"_SPIA_Results.csv-RankedTargets.tsv- 默认情况下，该文件将根据自定义加权指标对靶点按降序排列。请手动审阅输出文件，确保结果具有生物学意义，且所列靶点与正在评估的表型相符。
- 还将生成包含优先级排序的治疗方案及其相关指标的文件。对于示例输入文件，请使用“下载”文件夹中的输出文件，文件名为：
edgeR_dge.rds-TreatmentTumor-TreatmentNativeTissue_"timestamp"_SPIA_Results.csv-Treatments.tsv- 同样，默认情况下，该输出文件将根据加权指标对各个靶点（见步骤 7.1.1）对应的治疗方案按降序排列。请结合对基础生物系统的充分背景知识手动审阅该文件，以判断是否有必要开展进一步实验。由于市场上许多治疗方案可能影响多个靶点，因此多个治疗方案具有相同加权指标的情况是预期之中的。
访问受限。请登录或开始试用以查看此内容。
结果
方案中第1-3步所述的设置对于后续执行SPIA和Pathway2Targets算法是必要的。每一步结束时，系统将生成一条消息，以确认软件已成功安装。第4步包括下载现有的差异表达结果集，这些结果可以是提供的示例文件、其他现有文件，或预处理自定义的RNA测序数据集。第4步的主要要求是工作流程使用edgeR作为差异表达分析算法，并将结果存储为rds文件中的SingleCellExperiment对象。
使用SPIA算法预测具有统计学意义的细胞内信号通路（补充代码文件1 和 补充代码文件2；约130行代码）在步骤5中进行描述。运行后，该算法通过实时日志显示当前正在计算的通路（图2），表明算法正在正常运行。尽管许多SPIA分析可在约1小时内完成，但实际所需时间取决于差异表达基因（DEG）的数量、自举重复次数及其他因素。表1展示了SPIA输出文件（补充表1）的部分内容，列出了相关p值<0.05的...
访问受限。请登录或开始试用以查看此内容。
讨论
本方案的第1-3步专门涉及安装底层R软件、脚本和依赖项，以成功运行下游软件。所需R软件包的详细列表见（补充表4）。方案的第4步涉及获取一个R数据文件（.rds格式），该文件包含差异表达分析的结果输出。常用的软件包括edgeR6、DESeq24和limma5。我们构建的工作流程兼容以rds格式输出的edgeR结果，所用参数如下：使用每百万序列的对数2计数进行标准化（log = TRUE），在标准化中考虑转录本长度偏移量（offset = dge0$offset），设定保留基因的过滤条件（min.count = 10），设定在所有样本中保留基因的过滤条件（min.total.count = 15），离散估计方法（trend.method = "locfit"），以及模型（glmQLFit）。目前已有多种计算工作流程可用于预处理RNA测序读段或读段计数文件，并生成edgeR输出，包括Galaxy平台27以...
访问受限。请登录或开始试用以查看此内容。
披露
BEP 持有 Pythia Biosciences 的股权。本研究未获得外部资金支持。
致谢
我们感谢杨百翰大学研究计算办公室在访问校园高性能计算环境期间提供的专业知识和支持。
访问受限。请登录或开始试用以查看此内容。
材料
| 姓名 | 公司 | 目录编号 | 评论 |
|---|---|---|---|
| Pathway2Targets R 脚本 | Brigham Young University（Pickett 实验室） | 版本 3.1 | https://github.com/bpickett/Pathway2Targets/blob/main/Pathway2Targets.R |
| R 软件 | 综合 R 存档网络（CRAN） | 版本：4.4.3 | https://cran.r-project.org |
| R Studio Desktop 软件 | posit | 版本：2024.12.1+563 | https://posit.co/download/rstudio-desktop/ |
| SPIA R 脚本 | Brigham Young University（Pickett 实验室） | 版本：3.1 | https://github.com/bpickett/Pathway2Targets/blob/main/SPIA_Code.Rmd |
参考文献
- Dobin, A., et al. Ultrafast universal RNA-seq aligner. Bioinformatics. 29 (1), 15-21 (2013).
- Kim, D., Paggi, J. M., Park, C., Bennett, C., Salzberg, S. L. Graph-based genome alignment and genotyping with hisat2 and hisat-genotype. Nat Biotechnol. 37 (8), 907-915 (2019).
- Patro, R., Duggal, G., Love, M. I., Irizarry, R. A., Kingsford, C. Salmon provides fast and bias-aware quantification of transcript expression. Nat Methods. 14 (4), 417-419 (2017).
- Love, M. I., Huber, W., Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15 (12), 550(2014).
- Ritchie, M. E., et al. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 43 (7), e47(2015).
- Robinson, M. D., McCarthy, D. J., Smyth, G. K. Edger: A bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 26 (1), 139-140 (2010).
- Orjuela, S., Huang, R., Hembach, K. M., Robinson, M. D., Soneson, C. Armor: An automated reproducible modular workflow for preprocessing and differential analysis of RNA-seq data. G3 (Bethesda). 9 (7), 2089-2096 (2019).
- Zhang, X., Jonassen, I. Rasflow: An RNA-seq analysis workflow with Snakemake. BMC Bioinformatics. 21 (1), 110(2020).
- Bhardwaj, V., et al. Snakepipes: Facilitating flexible, scalable, and integrative epigenomic analysis. Bioinformatics. 35 (22), 4757-4759 (2019).
- Chen, Y., Xu, R. Drug repurposing for glioblastoma based on molecular subtypes. J Biomed Inform. 64, 131-138 (2016).
- Xu, Y., Kong, J., Hu, P. Computational drug repurposing for Alzheimer's disease using risk genes from GWAS and single-cell RNA sequencing studies. Front Pharmacol. 12, 617537(2021).
- Subramanian, A., et al. A next-generation connectivity map: L1000 platform and the first 1,000,000 profiles. Cell. 171 (6), 1437-1452.e17 (2017).
- Wang, Z., Lachmann, A., Keenan, A. B., Ma'ayan, A. L1000fwd: Fireworks visualization of drug-induced transcriptomic signatures. Bioinformatics. 34 (12), 2150-2152 (2018).
- Chan, J., Wang, X., Turner, J. A., Baldwin, N. E., Gu, J. Breaking the paradigm: Dr insight empowers signature-free, enhanced drug repurposing. Bioinformatics. 35 (16), 2818-2826 (2019).
- Keenan, A. B., et al. The library of integrated network-based cellular signatures NIH program: System-level cataloging of human cells response to perturbations. Cell Syst. 6 (1), 13-24 (2018).
- Jackson, M., et al. Transcriptomic insights into gas6-induced placental dysfunction: Gene targets for preeclampsia therapy. Cells. 14 (4), 278(2025).
- Rapier-Sharman, N., et al. Secondary transcriptomic analysis of triple-negative breast cancer reveals reliable universal and subtype-specific mechanistic markers. Cancers (Basel). 16 (19), 3379(2024).
- Sutherland, L., Lang, J., Gonzalez-Juarbe, N., Pickett, B. E. Secondary analysis of human bulk RNA-seq dataset suggests potential mechanisms for letrozole resistance in estrogen-positive (ER+) breast cancer. Curr Issues Mol Biol. 46 (7), 7114-7133 (2024).
- Rapier-Sharman, N., Clancy, J., Pickett, B. E. Joint secondary transcriptomic analysis of non-Hodgkin's B-cell lymphomas predicts reliance on pathways associated with the extracellular matrix and robust diagnostic biomarkers. J Bioinform Syst Biol. 5 (4), 119-135 (2022).
- Moreno, C., Bybee, E., Tellez Freitas, C. M., Pickett, B. E., Weber, K. S. Meta-analysis of two human RNA-seq datasets to determine periodontitis diagnostic biomarkers and drug target candidates. Int J Mol Sci. 23 (10), (2022).
- Gray, M., et al. Chikungunya virus time course infection of human macrophages reveals intracellular signaling pathways relevant to repurposed therapeutics. PeerJ. 10, e13090(2022).
- Ochoa, D., et al. The next-generation open targets platform: Reimagined, redesigned, rebuilt. Nucleic Acids Res. 51 (D1), D1353-D1359 (2023).
- Dobbs Spendlove, M., et al. Pathway2targets: An open-source pathway-based approach to repurpose therapeutic drugs and prioritize human targets. PeerJ. 11, e16088(2023).
- Li, Q. L., et al. Genome-wide profiling in colorectal cancer identifies phf19 and tbc1d16 as oncogenic super enhancers. Nat Commun. 12 (1), 6407(2021).
- Clough, E., et al. Ncbi geo: Archive for gene expression and epigenomics data sets: 23-year update. Nucleic Acids Res. 52 (D1), D138-D144 (2024).
- Katz, K., et al. The sequence read archive: A decade more of explosive growth. Nucleic Acids Res. 50 (D1), D387-D390 (2022).
- Galaxy, C. The galaxy platform for accessible, reproducible and collaborative biomedical analyses: 2022 update. Nucleic Acids Res. 50 (W1), W345-W351 (2022).
- Tarca, A. L., et al. A novel signaling pathway impact analysis. Bioinformatics. 25 (1), 75-82 (2009).
- Kanehisa, M., Furumichi, M., Tanabe, M., Sato, Y., Morishima, K. Kegg: New perspectives on genomes, pathways, diseases and drugs. Nucleic Acids Res. 45 (D1), D353-D361 (2017).
- Gillespie, M., et al. The Reactome Pathway Knowledgebase 2022. Nucleic Acids Res. 50 (D1), D687-D692 (2022).
- Li, Z., et al. Construction and function analysis of the lncRNA-miRNA-mRNA competing endogenous RNA network in autoimmune hepatitis. BMC Med Genomics. 15 (1), 270(2022).
- Subramanian, A., et al. Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 102 (43), 15545-15550 (2005).
访问受限。请登录或开始试用以查看此内容。
重印与许可
