一种用于代谢通路信息全基因组关联分析的通路关联分析工具

摘要

通过使用 Shiny 应用程序或 R 控制台运行通路关联分析工具（PAST），研究人员可以通过研究相关的代谢通路，更深入地理解全基因组关联研究（GWAS）结果的生物学意义。

摘要

最近，一种基于代谢通路分析解读全基因组关联研究（GWAS）数据的已有方法得到了新的实现并已发布。该通路关联分析工具（PAST）的开发旨在解决原有工具在用户友好性方面的不足以及分析速度缓慢的问题。这一新型用户友好型工具已发布于Bioconductor和Github平台。在测试中，PAST将原本需要二十四小时或更长时间的分析任务缩短至一小时以内完成。本文中，我们介绍了使用Shiny应用程序或R控制台运行PAST的操作方案。

引言

全基因组关联分析（GWAS）是研究复杂性状及其相关基因组区域的常用方法^1,2,3。在此类研究中，会检测数十万个单核苷酸多态性（SNP）标记与性状之间的关联，并评估这些关联的显著性。满足错误发现率（FDR）阈值（或其他类型的显著性阈值）的标记-性状关联将被保留用于后续分析，但真实的关联可能因此被过滤掉。对于复杂且多基因控制的性状，每个基因的效应可能很小（因而易被过滤），并且某些等位基因仅在特定条件下表达，而这些条件在研究中可能并不存在³。因此，尽管可能保留许多与性状相关的SNP，但每个SNP的效应可能非常微弱。过多缺失的SNP调用会导致对性状生物学意义及遗传结构的解释不完整且令人困惑。代谢通路分析可通过关注按生物学功能分组的基因的综合效应，帮助解决其中部分问题^4,5,6。

已有若干研究采用了本文所述方法的先前版本。先前版本曾被用于研究黄曲霉毒素积累⁷、玉米穗虫抗性⁸以及油脂生物合成⁹。尽管这些分析取得了成功，但分析过程复杂、耗时且繁琐，原因是分析工具由 R、Perl 和 Bash 多种语言编写，且流程未实现自动化。由于修改该方法以适应不同分析需要专业知识，目前开发出了一种新方法，可供其他研究人员共享使用。

通路关联分析工具（Pathway Association Study Tool, PAST）¹⁰ 旨在克服先前方法的不足，其设计减少了对编程语言知识的需求，并可在更短时间内完成分析。尽管该方法已在玉米中进行了测试，但PAST并不依赖于特定物种的假设。PAST可通过R控制台运行，也可作为Shiny应用程序使用，且预计不久后将在MaizeGDB上推出在线版本。

方案

1. 设置

1. 如果尚未安装 R，请先安装 R。
   注意：PAST 使用 R 编写，因此要求用户安装 R。在撰写本文时，直接从 Bioconductor 安装 PAST 需要 R4.0。对于 R3.6，可从 Bioconductor 安装较早版本的 PAST；而对于 R3.5 用户，可从 Github 安装 PAST。R 的安装说明可从以下链接下载：https://www.r-project.org/。
2. 安装最新版本的 RStudio Desktop 或更新现有 RStudio（可选）。
   注意：RStudio 是一个有助于使用 R 语言的开发环境，建议安装，特别是对于选择通过命令行而非 Shiny 图形用户界面运行 PAST 的用户。RStudio 及其安装说明可在以下链接获取：https://rstudio.com/products/rstudio/。
3. 按照 Bioconductor 上的说明，从 Bioconductor 安装 PAST¹¹。
   注意：通过 Bioconductor 安装将自动处理 PAST 所需依赖包的安装。此外，PAST 也可从 Github¹² 安装，但从 Github 安装不会自动安装依赖包。
4. 安装 PAST Shiny（可选）。从 Github 代码仓库的 Releases 页面下载文件 "app.R"：https://github.com/IGBB/PAST/releases/，并记住下载文件的存储位置。
   注意：用户可通过在 R 中直接调用 PAST 的方法来使用该工具，但对于不熟悉 R 的用户，可运行 PAST Shiny 应用程序，其提供引导式用户界面。PAST Shiny 是一个 R 脚本，位于 PAST Github 仓库的 shiny_app 分支中。PAST Shiny 在首次运行时会尝试自动安装其依赖包。
5. 通过以下三种方式之一启动应用程序，开始分析。
   1. PAST Shiny 与 RStudio 配合使用
      1. 在 RStudio 中，于 app.R 所在的文件夹中创建一个新项目。点击文件 | 新建项目，并选择该文件夹。
      2. 创建新项目后，打开之前下载的 app.R 文件。RStudio 会识别出 app.R 是一个 Shiny 应用，并在显示的源代码上方工具栏中生成一个 运行应用 按钮。点击 运行应用。RStudio 随即会启动一个窗口，显示 PAST Shiny 应用程序。
   2. PAST Shiny 与 R 控制台配合使用
      1. 启动 R，并运行以下代码以启动 PAST Shiny 应用程序：shiny::runApp('path/to/folder/with/shiny/app.R')。请将引号内的文本替换为 app.R 文件的下载路径，并保留引号。
   3. 不使用 R Shiny 的 PAST 使用方式
      1. 在 R 控制台中运行 library(PAST) 以加载 PAST。

2. 自定义 Shiny 分析（可选）

1. 将分析标题从“New Analysis”更改为更能反映所执行分析类型的名称，以便于跟踪多个分析（见图1）。

图1。 请点击此处查看此图的放大版本。

2. 修改核心数量和模式。将核心数设置为 1 至机器总核心数之间的任意数值，但需注意，为 PAST 分配更多资源可能会减慢机器上其他操作的运行速度。根据第 6 节中的描述设置模式。

3. 加载 GWAS 数据

注意：请验证 GWAS 数据是否为制表符分隔。确保关联文件包含以下列：性状、标记名称、位点或染色体、染色体上的位置、p 值以及标记的 R² 值。确保效应文件包含以下列：性状、标记名称、位点或染色体、染色体上的位置以及效应值。这些列的顺序无关紧要，因为用户在加载数据时可指定各列的名称。任何额外的列将被忽略。可使用 TASSEL¹³ 生成这些文件。

1. 使用 PAST Shiny 加载 GWAS 数据。
   1. 使用关联文件和效应文件选择框分别选择一个关联文件和一个效应文件。在文件选择框下方的关联列名称和效应列名称输入框中，修改列名称以匹配数据中的实际列名。

图 2. 请点击此处查看此图的放大版本。

2. 在 R 控制台中使用 PAST 加载 GWAS 数据。
   1. 修改并运行以下代码：
      gwas_data = load_GWAS_data("path/to/association_file.tsv", "path/to/effects_file.tsv", association_columns = c("Trait", "Marker", "Locus", "Site", "p", "marker_R2"), effects_columns = c("Trait", "Marker", "Locus", "Site", "Effect")
3. 注意：请将路径更改为 GWAS 文件的实际存储位置。association_columns 和 effects_columns 所提供的值为默认值。如果列名与默认值不一致，则需明确指定列名；否则可省略这些参数。

4. 加载连锁不平衡（LD）数据

注意：请验证连锁不平衡（LD）数据是否为制表符分隔，且包含以下类型的数据：位点（Locus）、位置1（Position1）、位点1（Site1）、位置2（Position2）、位点2（Site2）、位置1与位置2之间的碱基对距离，以及 R² 值。

1. 使用 PAST Shiny 加载 LD 数据。
   1. 选择包含 LD 数据的文件。如有必要，在文件选择框下方的 LD 列名称 输入框中修改列名称，使其与 LD 数据中的列名称一致。

图3。 请点击此处查看此图的放大版本。

2. 在 R 控制台中使用 PAST 加载 LD 数据。
   1. 修改并运行以下代码以加载 LD 数据：
      LD = load_LD("path/to/LD.tsv", LD_columns = c("Locus1", "Position1", "Site1", "Position2", "Site2", "Dist_bp", "R.2")
      注意：请将路径更改为 LD 文件的实际位置。LD_columns 提供的值为默认值。如果列名与这些默认值不匹配，请指定正确的列名；否则可省略此参数。

5. 将SNP分配给基因

注意：下载或以其他方式获取 GFF 格式的注释文件。这些注释通常可在特定生物体的在线数据库中找到。需谨慎对待质量较低的注释，因为注释数据的质量将直接影响通路分析的结果质量。请确认这些注释文件的第一列（染色体）与关联分析、效应值及 LD 数据中的位点/染色体命名格式一致。例如，如果 GWAS 和 LD 数据文件中将第一条染色体记为“1”，则注释文件不应将其记为“chr1”。

1. 使用 PAST Shiny 将 SNP 分配至基因。
   注意：有关确定合适 R² 截断值的更多信息，请参见 Tang 等人⁶ 的文章中名为“通路分析的 SNP 至基因算法”部分。
   1. 选择包含 GFF 注释的文件。根据所研究物种的特点，考虑最合适的窗口大小和 R² 截断值，若默认值不适用于上传的数据，则进行相应修改。
      注意：PAST 中的默认值主要适用于玉米。PAST Shiny 分析开始时设置的核数（步骤 2.2）将在此步骤中使用。

图4. 请点击此处查看此图的放大版本。

2. 在 R 控制台中使用 PAST 将 SNP 分配给基因。
   1. 修改并运行以下代码，以将 SNP 分配给基因：
      genes = assign_SNPs_to_genes(gwas_data, LD, "path/to/annotations.gff", c("gene"), 1000, 0.8, 2)
      注意：在此示例代码中，提供了若干默认参数：1000 表示在 SNP 周围搜索基因的窗口大小；0.8 为 R² 的阈值；2 为用于并行处理的计算核心数。注释文件的路径也应更改为注释文件的实际存储位置。

6. 发现显著性通路

注意：请验证通路文件是否包含以下以制表符分隔的数据，每条通路中的每个基因对应一行：通路 ID——例如“PWY-6475-1”的标识符；通路描述——对通路功能的详细说明，例如“反式番茄红素生物合成”；基因——通路中的一个基因，其名称应与注释中提供的名称一致。通路信息通常可在特定生物的在线数据库中找到，例如 MaizeGDB。第二个用户自定义选项是模式。“增加”指表型测量值越高越有利的情况，例如产量；而“减少”指测量值越低越有利的性状，例如虫害评分。通路的显著性通过先前描述的方法进行检验^4,6,14。

1. 使用 PAST Shiny 发现显著通路。
   1. 选择包含通路数据的文件，并确保在分析选项中已选择相应的模式。如有必要，可更改保留在分析中的通路所必须包含的基因数量，以及用于构建零分布以检验效应显著性的置换次数。

图 5．请点击此处查看此图的放大版本。

注意：此步骤使用在 PAST Shiny 分析开始时设置的核心数量和模式（步骤 2.2）。默认的基因数量目前设定为 5 个基因，因此已知基因数少于该值的通路将被剔除。用户可将此值降低至 4 或 3 以包含较短的通路，但这样做会增加假阳性结果的风险。提高该值可增强分析的效力，但会导致更多通路被排除在分析之外。改变所用置换次数的数值会相应地提高或降低检验的效力。

2. 在 R 控制台中使用 PAST 发现显著通路。
   1. 修改并运行以下代码以发现显著通路：
      rugplots_data <- find_pathway_significance(基因, "path/to/pathways.tsv", 5, "increasing", 1000, 2)
      ​注意：在此示例代码中，提供了若干建议的默认参数。其中，5 表示通路中必须包含的最少基因数，以便在分析中保留该通路；"increasing" 表示所测性状呈递增趋势（建议用户无论性状如何，均分别运行递增和递减方向；但两种方向的数据解释将有所不同）；1000 表示为确定零分布而进行抽样的次数；2 表示用于并行处理的计算核心数。请将路径修改为通路文件的实际存储位置。

7. 查看 Rug图

1. 使用 PAST Shiny 查看 Rugplots。
   1. 上传并设置所有输入后，点击 开始分析。此时将显示一个进度条，指示分析流程中最后完成的步骤。分析完成后，PAST Shiny 将切换至“结果”标签页。左侧列（标注为“通路”）将显示结果表格，右侧列（标注为“图谱”）将显示 Rugplots。
   2. 使用滑块调节筛选参数。当筛选水平达到满意时，点击左下角的 下载结果 按钮，将所有图像和表格分别下载为一个以分析标题命名的 ZIP 文件。该 ZIP 文件包含筛选后的表格、未筛选的表格，以及筛选后表格中每个通路对应的单张图像。

图 6。 请点击此处查看此图的放大版本。

图7。 请点击此处查看此图的放大版本。

2. 在 R 控制台中使用 PAST 查看 Rugplots
   1. 修改并运行以下代码以保存结果：
      plot_pathways(rugplots_data, "pvalue", 0.02, "increasing", "output_folder")
      注意：在此示例代码中提供了若干建议的默认参数。pvalue 提供可用于在用户选择显著性阈值后过滤不显著通路的数据；0.02 是用于过滤的默认值；increasing 表示所测量性状的增加趋势（建议用户无论性状如何，均分别运行 increasing 和 decreasing 两种情况；然而，两者的数据解释将有所不同）；output_folder 是图像和表格将被写入的文件夹（该文件夹必须在运行函数前已存在）。过滤后的结果表格、未过滤的结果以及每个过滤后通路的单独图像将被写入此文件夹。

结果

如果运行 PAST 软件工具后未生成结果，请检查所有输入文件的格式是否正确。使用 PAST 软件包中提供的示例数据（基于玉米籽粒颜色的全基因组关联分析）成功运行的结果如图8所示。该表格及生成的图像可通过“下载结果”按钮进行下载。下载后的图像示例如图2¹⁰所示。不正确的设置可能导致结果在生物学上无意义，但判断结果是否正确应由研究人员自行负责，需仔细核对所选设置的有效性，并综合考虑有关目标性状的所有已知证据。

图9¹⁰展示了利用包含288个自交系的玉米群体进行全基因组关联分析（GWAS）后，基于籽粒颜色表型数据所生成的通路分析结果的地毯图。该示例较为简单，表型仅分为“白色”或“黄色”两类，选择此性状的原因在于，已知类胡萝卜素黄色色素的合成通路明确，且应主导该表型的大部分变异。因此，我们预期会观察到反式番茄红素生物合成通路（即类胡萝卜素合成通路）与籽粒颜色显著相关，结果确实如此。图中顶部列出了通路的编号（Pathway ID）和名称。图的横轴表示分析中包含的所有基因，按其对性状影响效应从大到小从左至右排列。然而，仅在图中顶部以刻痕线形式标出属于反式番茄红素生物合成通路的基因（其位置对应于所有分析基因中的效应排序）。该通路共包含7个基因。纵轴表示累积富集评分（running enrichment score, ES）。每个基因的ES值按其效应大小顺序累加至总分中，并根据分析的基因总数进行调整。因此，随着沿横轴向右移动，评分不断变化，当效应较大的基因被纳入时，总分趋于上升；但当新增基因带来的效应增量小于因基因数量增加而进行的调整值时，总分开始下降。累积ES曲线的峰值由一条垂直虚线标出；该峰值即为整个通路的ES值，程序据此判断该通路是否被选中并以地毯图形式呈现。

图8：PAST Shiny 完成运行结果。 请点击此处查看该图的放大版本。

图9：PAST 完整运行后生成的（或从 Shiny 下载的）通路图。 本图引自 Thrash 等人¹⁰。请点击此处查看此图的放大版本。

讨论

PAST 的主要目标之一是将全基因组关联研究（GWAS）数据的代谢通路分析推广至更广泛的用户群体，尤其是非人类和非动物生物的研究者。现有的其他方法通常是基于命令行的程序，且主要集中于人类或动物。在开发 PAST 时，易用性是首要目标，这体现在选择开发 Shiny 应用程序以及采用 R 和 Bioconductor 来发布该应用程序。用户无需学习如何编译程序即可使用 PAST。

与大多数分析软件一样，PAST 的分析结果质量取决于输入数据的质量；如果输入数据存在错误或格式不正确，PAST 将无法运行或产生无意义的结果。确保全基因组关联研究（GWAS）数据、连锁不平衡（LD）数据、注释信息和通路文件的格式正确，是获得准确输出结果的关键。PAST 仅能分析双等位基因标记，并且每组输入数据只能分析一个性状。此外，由低质量基因分型或错误、不精确的表型鉴定所产生的 GWAS 数据，也不太可能产生清晰或可重复的结果。PAST 可辅助对 GWAS 结果进行生物学解释，但如果环境变异、实验误差或群体结构未得到适当校正，则难以澄清混乱的数据集。

用户可以选择更改分析的某些参数，既可以在 Shiny 应用程序中进行，也可以通过将这些参数传递给 R 控制台中的 PAST 函数来实现。这些参数可能会影响 PAST 报告的结果，因此用户在修改默认设置时应格外谨慎。由于连锁不平衡（LD）由用户自行测定，通常使用与全基因组关联分析（GWAS）相同的标记数据集，因此 LD 的测量结果具有群体特异性。对于所有研究，尤其是针对玉米以外的物种（特别是自花授粉、多倍体或高度异质性的物种），可能需要调整默认参数。

材料

本文使用的材料清单

姓名	公司	目录编号	评论
计算机	NA	NA	任何具有 8GB 内存的计算机均应足够
R	R Project	NA	需要 R 4.0 或更高版本才能从 Bioconductor 3.11 安装