使用人类差异表达基因列表进行下游通路富集分析和靶点优先排序

摘要

目前的工作描述了一种运行 Pathway2Targets 算法的协议，这是一种 R 脚本，它根据通过比较来自批量 RNA 测序实验的病例样本与对照样本生成的细胞内信号通路的概况来预测治疗靶点并确定治疗靶点的优先级。

摘要

该协议概述了一个多步骤计算管道，用于从 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 种针对 TBNC 的单独药物，其中 37 种正在测试¹⁷.在淋巴瘤研究中，该算法预测了 915 种药物，其中 461 种已获得 FDA 批准¹⁹。

目前工作的目的是描述一种计算协议，使更多的研究人员能够有效地使用最近开发的 Pathway2Targets 算法，他们可能会受益于访问有关在命令行上运行程序的更具描述性的指令（图 1）。Pathway2Targets 通过结合差异表达数据、基因-疾病关联、临床试验信息、公共靶点数据²²、通路信息和其他指标来预测给定病症的靶点。重要的是，该算法结合了一种独特且可定制的加权方案，允许用户确定他们喜欢在分析中强调的~20个与靶点相关的指标，例如疾病关联的数量，信号通路的数量，独特药物的数量，临床试验每个阶段的治疗方法数量等²³。作为该协议的示例用例，我们将重新分析现有的结直肠癌数据集²⁴。

方案

本研究中分析的批量RNA测序数据是从公开可用的数据库（NCBI基因表达综合和序列读取档案）^25,26中获取的。因此，原始数据收集者确保从知情和同意的人类受试者那里合乎道德和适当地收集这些样本。

1. 下载并安装 R 软件

1. 通过单击综合 R 存档网络 （CRAN） https://cran.r-project.org/mirrors.html 中的相应链接来安装 R（版本 4.0 或更高版本）并使用选项 0-Cloud，然后按照计算机作系统的相应说明进行作。这个过程通常需要 5-10 分钟。
2. 从 https://posit.co/download/rstudio-desktop/ 安装 R Studio（版本 2024 或更高版本），然后按照下载页上的说明进行作。安装 RStudio 通常需要 ~10 分钟。
   注意：安装 RStudio 是可选的（但强烈建议），因为它提供了一个集成的开发环境，可以简化代码执行、提供语法突出显示、促进包管理并帮助用户可视化输出，这对于不太熟悉 R 的用户特别有益。

2. 下载并安装相关工具的 R 脚本

1. 从 GitHub 存储库 https://github.com/bpickett/Pathway2Targets 下载以下必要的 R 脚本。此网址仅供参考。
2. 使用以下链接下载脚本：SPIA 1.0 版（通过单击 下载原始文件下载）：https://github.com/bpickett/Pathway2Targets/blob/main/SPIA_Code.Rmd（提交 ID：60fcd46）;Pathway2Targets版本 3.1（通过单击 下载原始文件下载）：https://github.com/bpickett/Pathway2Targets/blob/main/Pathway2Targets.R（提交 ID：8e4c7c8）

3. 下载相关工具的 R 库

1. 在运行 R 时（在 RStudio 中或在终端窗口中），输入以下命令以下载并安装运行软件所需的附加 R 库。
   1. 启动 RStudio 程序。默认情况下，控制台面板位于 RStudio 的左下角。单击控制台面板窗口中的任意位置，键入光标应出现在箭头“>”符号之后的底部。
   2. 在控制台区域复制并粘贴以下命令，然后按 Enter 键：
      install.packages(c("RCurl", "stringr", "jsonlite", "httr")). 
      成功安装后，将出现一条状态消息，显示“下载的二进制包位于 ......”。
   3. 将以下命令复制并粘贴到控制台面板区域，然后按 Enter 键。
      BiocManager::install(c("SummarizedExperiment", "EnrichmentBrowser", "biomaRt", "org.Hs.eg.db")).
      成功安装后，将出现类似的消息，指出“下载的二进制包位于 ......”。
      注意：第一组文库（步骤 3.1.1）由典型的 R 文库组成，而第二组文库（2.1.2）由 BioConductor 文库组成。因此，需要单独输入命令。这些库的相应版本应根据计算机上安装的 R 版本自动下载。这些命令中的每一个都应需要 ~5 分钟才能完成。

4. 文件处理

1. 将先前计算的由 ARMOR 软件（或类似软件）生成的 edgeR 生成的差分表达式输出文件（RDS 格式）下载到本地计算机。此文件的名称通常为 edgeR_dge.rds。
2. 手动查看 edgeR（或类似的差异表达）的结果，以开始对结果进行生物学相关解释。为此，至少按校正后的 p 值 < 0.05 进行过滤，并可能按对数₂ 倍变化值的绝对值 > 1.5 进行过滤。如果测试该软件，可以在Zenodo上找到示例edgeR_dge.rds文件：https://doi.org/10.5281/zenodo.15186609
   注意：回顾这些过滤结果中残留的基因可以开始解释与病例样本相关的表型的潜在分子机制（与对照样本相比）。重要的是要认识到，由于可以快速调用的基因符号数量相对较少，因此以公正的方式解释基因列表的能力是极其困难的。因此，信号通路分析是一种有用的方法，可以根据过滤后的基因在细胞中相互作用和/或相互交流的方式来总结它们。
3. 预处理批量 RNA-seq 数据可能需要几小时到几天的计算时间，具体取决于所分析数据集的大小。将此 .rds 文件保存在计算机上的“下载”文件夹中。请注意，此 .rds 文件类型不是人类可读的。

5. 运行SPIA通路富集算法

1. 如果使用 R 运行，请使用 GitHub 或 补充编码文件 1 中的 R 脚本。输入以下命令，假设 edgeR_dge.rds 文件位于 Downloads 文件夹中
   Rscript --vanilla SPIA_Code.Rmd ~/Downloads/edgeR_dge.rds
   1. 如果 edgeR_dge.rds 文件位于其他文件夹（或目录）中，则将此命令替换为以下命令
      Rscript --vanilla SPIA_Code.Rmd
2. 如果使用 RStudio 运行，请使用 GitHub 中的 R 脚本或 补充编码文件 2。
   注意：在语言 R 中，在一行代码之前添加哈希标签 # 符号会暂时禁用它。这些脚本最初设计为在命令行环境中运行，而不是在 RStudio 中运行。启用或禁用某些代码行是重新配置输入文件设置的最简单方法。
   1. 在 R Studio 中打开 SPIA_Code.Rmd 脚本，方法是单击“文件”菜单中的“ 打开文件 ”选项，然后选择脚本名称。在 R Studio 代码窗口中默认情况下，它位于左上角的面板。
   2. 选择文件中的所有代码行，然后单击位于代码窗口上方和右侧的“ 运行 ”按钮（或“ 运行所选行）”按钮。成功运行将生成一个名为
      edgeR_dge.rds-TreatmentTumor-TreatmentNativeTissue_2025-04-23_
10-56-45.12767_SPIA_Results.csv
      在下载目录中。该文件将包含具有统计学意义的信号通路
   3. 通过将文件作为电子表格打开，手动查看具有统计显着性结果的文件。该文件的内容应有助于总结由差异表达基因显着代表的潜在细胞内信号级联反应。
      注意：完成有效通路的计算可能需要~30分钟到几个小时，具体取决于正在分析的数据集中信号的强度。当程序运行时，实时进度消息将在控制台窗口中不断更新。经常更新的消息表明该程序正在成功运行。有关此步骤中发生的情况的更详细说明，请参阅 GitHub 存储库：https://github.com/bpickett/Pathway2Targets/tree/main
   4. 如果使用示例输入文件，则此步骤的输出文件将在“下载”文件夹中找到，并将命名为
      edgeR_dge.rds-TreatmentTumor-TreatmentNativeTissue_"timestamp"_SPIA_Results.csv.
      此命名样式反映了输入文件名、正在执行的过程和输出。这可以避免在处理多个文件时识别文件时出现混淆。
   5. 调整其他参数，如下所述。
      1. 该库中 SPIA 算法的默认参数为 1,000 个排列。将此值增加到 2,000 个排列，以提高对结果的置信度。通过将 perm = 2000 更改为所需的排列数量来调整此脚本第 84 行和第 85 行中的排列数量。调整其他参数，如下所述。
      2. 通过删除 padj.method = 'BH' 来调整第 84 行和第 85 行中的 p 值校正方法。这不会纠正 p 值，这可能会使假阳性结果更有可能。

6. 在 SPIA 输出上运行 Pathway2Targets 目标优先级算法

1. 如果运行 R，请使用 GitHub 中的 R 脚本或 补充编码文件 3。使用以下命令调用此算法
   Rscript --vanilla Pathway2Targets.R
2. 如果使用 R Studio 运行，请使用 GitHub 中的 R 脚本或 补充编码文件 4。在 R Studio 中打开 Pathway2Targets.R 脚本，方法是单击“文件”菜单中的“ 打开文件 ”选项，然后选择脚本名称。
   1. 在 RStudio 代码窗口（左上角面板）中，将第 22 行的文件名替换为 SPIA 结果文件名，例如（来自示例数据）
      infile <- "edgeR_dge.rds-TreatmentTumor-TreatmentNativeTissue_2025-04-23_
10-56-45.12767_SPIA_Results.csv"
   2. 选择文件中的所有代码行，然后单击位于代码窗口上方和右侧的“ 运行” 按钮。实时进度状态消息将持续显示在右下角面板中。成功运行将生成一个名为 （类似于） 的文件
      "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 值的信号通路数量、这些重要途径中的基因产物数量以及作为已知药物靶标的基因产物数量。
3. 可以调整 Pathway2Targets 算法的一些参数。具体来说，调整乘数值（脚本的第 31-38 行）以自定义每个指标的权重方案。有关此步骤中发生的情况的更详细说明，请参阅相应的 GitHub 存储库：https://github.com/bpickett/Pathway2Targets/tree/main
   注意：作为参考，在本例中，从由数百个单独靶标组成的 132 条途径中识别潜在药物大约需要 2 小时才能完成。根据这个指标，估计总计算时间是合理的。

7. 打开优先靶点和治疗药物的结果文件

1. 将生成包含优先目标及其指标的文件。对于示例输入文件，请使用 Downloads 文件夹中的输出文件，该文件名为：
   edgeR_dge.rds-TreatmentTumor-TreatmentNativeTissue_“timestamp”_SPIA_Results.csv-RankedTargets.tsv
   1. 默认情况下，此文件将根据自定义权重指标降序排序目标。手动查看输出文件，以确保结果具有生物学相关性，并且靶标对于正在评估的表型是合乎逻辑的。
2. 还将生成包含优先治疗及其指标的文件。对于示例输入文件，请使用 Downloads 文件夹中的输出文件，该文件名为：
   edgeR_dge.rds-TreatmentTumor-TreatmentNativeTissue_“时间戳”_SPIA_Results.csv-Treatments.tsv
   1. 同样，默认情况下，此输出文件也将与各种靶点的治疗药物（在步骤 7.1.1 中）进行排序，这些治疗药物根据加权指标降序排序。使用对底层生物系统的足够背景知识手动审查文件，以确定额外的实验是否合理。预计多种疗法可能具有相同的加权指标，因为许多靶点将受到市场上不止一种疗法的影响。

结果

协议步骤 1-3 中描述的设置对于启用 SPIA 和 Pathway2Targets 算法的后续执行是必要的。在每个步骤结束时，都会生成一条消息以确认软件安装成功。第 4 步包括下载一组现有的差异表达结果，其中可能包括提供的示例文件、不同的现有文件或预处理自定义 RNA 测序数据集。步骤 4 的主要要求是工作流使用 edgeR 作为差分表达式算法，并将结果存储为 rds 文件中的 SingleCellExperiment 对象。

步骤 5 描述了使用 SPIA 算法（补充 编码 文件 1 和补充 编码 文件 2;~130 行代码）预测具有统计学意义的细胞内信号通路。运行后，该算法通过正在计算的途径的实时日志来显示其进度（图 2

讨论

协议的步骤 1-3 专门涉及安装基础 R 软件、脚本和依赖项，以成功运行下游软件。提供了所需 R 库的详细列表（补充表 4）。该协议的第 4 步涉及检索 R 数据（.rds 格式）文件，其中包含差异表达分析的输出。此步骤的常用软件包括 edgeR⁶、DESeq2⁴ 和 limma⁵。我们使用以下参数构建了此工作流程，以与rds格式的edgeR输出兼容：使用每百万个log2计数进行归一化（log = TRUE），归一化时要考虑的转录长度偏移（offset = dge0$offset），过滤要保留的基因（min.count = 10），过滤要在所有样本中保留的基因（min.total.count = 15）， 色散估计方法 （trend.method = “locfit”） 和模型 （glmQLFit）。有各种计算工作流程可以预处理RNA测序读取或读取计数文件并生成edgeR输出，包括Galaxy平...

材料

本文使用的材料清单

姓名	公司	目录编号	评论
Pathway2Targets R 脚本	杨百翰大学（皮克特实验室）	版本 3.1	https://github.com/bpickett/Pathway2Targets/blob/main/Pathway2Targets.R
R 软件	综合R档案网络（CRAN）	版本：4.4.3	https://cran.r-project.org
R Studio 桌面软件	假设	版本：2024.12.1+563	https://posit.co/download/rstudio-desktop/
SPIA R 脚本	杨百翰大学（皮克特实验室）	版本：3.1	https://github.com/bpickett/Pathway2Targets/blob/main/SPIA_Code.Rmd