书名：Microarray Data Analysis
编辑：Giuseppe Agapito
年份：2022
出版社：HUMAN PRESS
ISBN：ISBN 978-1-0716-1839-4

Pathway Enrichment Analysis of Microarray Data

作者

• Chiara Pastrello, Yun Niu, and Igor Jurisica

摘要

• 微阵列分析通常产生一个差异基因列表，这些基因需要进行注释以：

  • 链接到所研究的表型
  • 帮助规划验证实验
  • 解释结果

• 通路富集分析经常用于上述目的，其中通路是人工创建的分子活动和过程模型。

• 虽然有不同类型的通路富集分析可用，本协议重点介绍最常见的类型——过度表示分析。

• 数据库多样性：

  • 许多数据库收集不同的通路集
  • 为相同的通路策划不同的基因集
  • 因此，选择最合适的通路来源进行富集分析非常重要

• 综合通路分析工具：

  • 本协议将使用pathDIP
  • pathDIP整合了22个主要通路数据库
  • 支持全面的富集分析

• 可视化步骤：

  • 描述使用GSOAP进行富集通路可视化的步骤

关键词

• Pathway enrichment analysis, Pathway orphans, Pathway consolidation, Database coverage, Database overlap, Enrichment visualization

引言

• 微阵列分析的目标是识别与研究条件相关的基因。
• 完成预处理和归一化后，常见步骤是识别不同条件下差异显著的基因（如疾病与正常、处理与未处理等）。
• 即使使用严格的统计过滤，研究人员通常也会得到一个较长的基因列表，需要对其进行优先级排序以进行后续验证，并结合表型进行解释。

基因注释的多样性

• 基因有多种注释类型，尤其在人类和模式生物中，这些注释广泛且丰富。

• 注释类型包括：

  • 组织注释（例如，TissueEnrich [1]）
  • 疾病注释（例如，Disease Ontology [2]）
  • 突变注释（例如，COSMIC [3]）
  • 基因本体论（Gene Ontology [4]）
  • 通路注释（最常用）

通路数据库的发展

• 通路是研究人员创建的复杂生物事件的简化表示，用于总结和存储当前知识。
• 1990年代中期，首两个通路数据库EcoCyc [5]和KEGG [6]被创建。
• 机器可读形式收集通路数据的重要性显现。
• 根据 pathguide.org 数据，目前有超过两百个数据库用于策划、收集、存储和/或可视化特定或全面的通路集合。
• 然而，研究人员通常只使用少数几个通路数据库进行注释，选择不总是基于“最佳”数据库，因为“最佳”在不同场景和需求下定义不同。

通路数据库的异质性与重叠

• 通路不是静态实体，缺乏“真实值”作为参考。
• 策划自文献的通路数据库之间预期会有一定的重叠。
• 某些数据库专注于特定反应和信号级联，如代谢通路（Biocyc [8]）或信号传导通路（SignaLink [9]），因此不与其他类型数据库重叠。
• 一般通路数据库涵盖广泛类别，预期有较高的重叠率，但实际最高重叠率仅约50% [7]。
• 数据库在包含的通路数量、通路名称和每个通路包含的基因上非常异质，导致选择合适的通路数据库具有挑战性。

[图片上传失败...(image-9b0764-1727822732738)]

<figcaption>image</figcaption>

基因覆盖率的差异

• 通路数据库在基因覆盖率上存在差异：

  • Reactome [10]注释约一万个人类基因
  • KEGG注释约七千个
  • WikiPathways [11]注释约六千个

• Reactome和KEGG之间的最高基因重叠率为六千多个基因。

• 注释缺失可能导致通路富集分析中的偏差和相关通路排名的差异 [12]。

• 例如，WNT信号通路在不同数据库中的注释数量和基因内容差异显著，导致富集分析结果差异较大。

[图片上传失败...(image-5adce1-1727822732738)]

<figcaption>image</figcaption>

解决方案：通路数据库的整合

• 整合不同的数据库可以扩展注释基因数量和涵盖的通路类型。
• 例如，pathDIP [13]集成了22个通路数据库，涵盖17个物种，注释人类超过13,000个基因到5,380个通路。
• 虽然整合有助于解决覆盖率问题，但也增加了返回的富集通路数量，导致多重比较调整后显著性p值急剧下降。
• 许多“相似”通路可能被返回，难以确定其具体来源和内容相似性。

解决方案：通路整合与合并

• 通路整合不仅包括数据库的整合，还包括相似通路的合并，以减少通路数量并扩展通路特定覆盖。
• 每个原始数据库中，相同通路可能有不同名称，相似名称可能指不同通路。
• 需要通过结合通路内容和名称相似性来分组功能相似的通路（例如，Pathcards [14]）。
• 最佳策略可能包括文本和内容相似性的结合。

通路孤立基因的预测

• 即使经过整合，经过策划的数据库最多覆盖基因组的三分之二，许多基因可能无注释，成为“通路孤立基因”。
• 预测通路注释可以弥补这一空白。
• 在pathDIP中，使用策划的通路进行预测，增加了人类超过18,000个基因的覆盖率，以及17个物种中超过120,000个基因的覆盖率。

计算需求

• 本协议执行的第一组分析将使用基于网络的软件。
• 用户需要互联网连接和网页浏览器。
• 对于可视化，我们将使用 R 中的 GSOAP。

获取差异表达基因

1. 为了从微阵列数据中获取差异表达的基因，我们将使用 GEO (https://www.ncbi.nlm.nih.gov/geo)，这是一个收集了数千个针对许多不同疾病的微阵列数据集的数据库。
2. 拥有自己微阵列数据和差异表达基因列表的研究人员可以跳过此步骤。

执行通路富集分析

1. 为了执行通路富集分析，我们将使用 pathDIP 4.0 (http://ophid.utoronto.ca/pathDIP)。
2. 有兴趣的读者可参考 http://ophid.utoronto.ca/pathDIP/API.jsp，将此类分析集成到使用 Java、R 或 Python 的生物信息学工作流程中。

在 R 中使用富集分析结果

• 富集分析的结果将用于 R。

• 本协议中的示例在 R 4.0.3 (https://www.r-project.org) 上运行。

• 需要安装 GSOAP [15] 包。要安装它，请在 R 中运行以下命令：

  require(devtools)install_github("tomastokar/gsoap", dependencies=TRUE)
  

方法

3.1 从微阵列数据获取差异基因

1. 连接到 GEO 并搜索疾病：

   • 研究人员在 GEO 网站（https://www.ncbi.nlm.nih.gov/geo）的搜索框中输入所选疾病。在我们的案例中，输入“骨关节炎”（“Osteoarthritis”），点击搜索，选择“GEO DataSets 数据库中共有 1292 个“骨关节炎”结果。点击数字将打开搜索结果页面。

2. 筛选搜索结果：

   • 根据研究兴趣筛选结果。在左侧，选择“Series”作为条目类型（见注1）和“Expression profiling by array”作为研究类型；在右侧，选择“人类”（“Homo sapiens”）作为主要生物体。对于本协议，选择数据集“GSE55584”。

3. 使用 GEO2R 进行差异基因分析：

   • 在数据集页面，点击“Analyze with GEO2R”。这将打开用于获取差异基因的页面。

   • 设置分析中使用的组：

     • 点击“define groups”，在自由文本字段中输入“OA”并按回车，然后输入“RA”并按回车（见注2）。

     • 选择适当的样本并将其链接到对应组：

       • 选择所有“临床状态”为“骨关节炎”的样本，点击“OA”。
       • 选择所有“临床状态”为“类风湿性关节炎”（“Rheumatoid arthritis”）的样本，点击“RA”。

     • 定义组将自动显示每个标签的样本数量（本例中为 6 和 10）。

   • 点击页面底部的“Analyze”按钮，GEO2R 将运行差异基因表达分析。

4. 下载和筛选差异基因：

   • 分析结果显示在页面底部的表格中。
   • 选择“Download full table”以获取平台上所有基因的结果（包含在数据 S2 中）。
   • 筛选最显著的基因，选择调整后的 p 值 < 0.01，折叠变化大于 1 或小于 -1，并具有基因符号注释的基因（见注3）。

3.2 通路富集分析

1. 在 pathDIP 中进行富集分析：

   • 将前一步骤中获得的 56 个基因符号列表粘贴到 pathDIP 的搜索框中。列表可以由基因符号、Entrez 基因 ID 或 UniProt ID 组成。
   • 保持默认选项或使用页面上的筛选器。注意通路来源，通路数据库根据存储的通路类型进行颜色编码（例如，信号传导、代谢、多种）（见注4）。
   • 在本分析中，选择“文献策划集”（literature curated set）。用户也可以为注释不足的基因选择“扩展通路关联”（Extended pathway association），并选择是否仅使用实验蛋白质相互作用或实验和高置信度预测蛋白质相互作用生成的集，以及预测的置信度水平。
   • 选择是下载还是显示结果。

2. 查看和下载富集分析结果：

   • pathDIP 将使用所选数据库执行过度表示分析（见注5），并返回三组结果：

     1. 注释基因表：列出在 pathDIP 中具有注释的基因，帮助用户确认所用集和数据库的适用性（见表2）。
     2. 富集通路表：包含富集通路及其 p 值、调整后的 p 值，以及指向原始数据库的链接。
     3. 基因-通路注释表。

   • 如果选择下载结果，三个表格将包含在一个 txt 文件中。富集通路的列表通常用于生物学解释（见注6）并作为出版物中的表格。

   • 页面还允许用户获取通路/基因矩阵，并执行词汇富集。词汇富集结果以每个有意义的通路词的 p 值和调整后的 p 值的表格形式展示，并提供用于词云软件（例如 Wordle—https://www.wordle.net）的词列表。数据 S4 显示了本示例的词汇富集结果。

[图片上传失败...(image-a60ad8-1727822732738)]

<figcaption>image</figcaption>

3.3 GSOAP 绘图

[图片上传失败...(image-1d96b8-1727822732738)]

<figcaption>image</figcaption>

1. 准备数据文件：

   • 从 pathDIP 获得的通路列表将在 R 中使用。

   • 为简化操作，将下载的文件分成两部分：

     • 富集通路集：包含“Pathway Source”、“Pathway Name”、“p-value”、“q-value (FDR: BH-method)”、“q-value (Bonferroni)”、“External Source ID”标题的部分，命名为“pathways_to_be_plotted.txt”。
     • 基因-通路注释：包含“UniProt”、“Gene Symbol”、“Entrez Gene”、“Pathway Source”、“Pathway Name”标题的部分，命名为“genes_annotation.txt”。避免使用 MS Excel 处理此部分（见注7）。

   • 在 R 中读取文件：

     pathways_to_plot = read.delim("pathways_to_be_plotted.txt",stringsAsFactors = FALSE, header = TRUE)genes_annotations = read.delim("genes_annotation.txt",stringsAsFactors = FALSE, header = TRUE)
     

2. 格式化数据以符合 GSOAP 要求：

   • GSOAP 需要一个至少包含“Pathway”、“p-value”或“adjusted p-value”以及每个通路中识别出的基因符号以斜杠分隔的文件。

   • 操作步骤：

     1. 聚合注释文件，使基因按通路分组并以斜杠分隔：

        genes_annotations = aggregate(Gene.Symbol ~ Pathway.Source +Pathway.Name, genes_annotations, paste, collapse = "/")
        

     2. 合并两个文件：

        pathways_to_plot = merge(x = pathways_to_plot, y = genes_annotations, by = c("Pathway.Source", "Pathway.Name"), all.x = TRUE)
        

     3. 筛选 p 值经过 Bonferroni 校正 < 0.01 的通路，并设置数据框的行名为通路名称：

        pathways_to_plot = pathways_to_plot[which(pathways_to_plot$q.value..Bonferroni. < 0.01),]rownames(pathways_to_plot) = make.names(pathways_to_plot[,2], unique = TRUE)
        

3. 使用 GSOAP 计算和绘图：

   • GSOAP 将使用 Jaccard 距离计算通路之间的基因重叠，并使用通路显著性（计算为 -log10(provided p-value)）来计算通路之间的接近度。显著性还用于对通路进行排序，因此在图中只显示最显著的标签（本例中为 5 个）。

   • 使用以下命令进行计算、排序和绘图：

     layout_path = gsoap_layout(pathways_to_plot, 'Gene.Symbol', 'q.value..FDR..BH.method.')layout_path = layout_path[order(layout_path$significance, decreasing = TRUE),]gsoap_plot(layout_path, as.color = 'cluster', as.alpha = 'significance', which.label = 1:5)
     

   • 图2显示了本示例的结果。

注释

1. 注1：Entry type 中的 “Series” 指的是 GEO 数据集中的系列类型，用于描述一组相关的实验数据。
2. 注2：定义组时，确保正确标注每个样本所属的组，以确保差异分析的准确性。
3. 注3：Gene Symbol 注释确保基因列表中的每个基因都有对应的标准符号，便于后续分析和解读。
4. 注4：通路数据库的颜色编码有助于快速识别不同类型的通路，提升分析效率。
5. 注5：过度表示分析（Over-Representation Analysis, ORA）是一种常用的通路富集方法，通过检测目标基因集中某些通路中的基因是否显著多于预期来识别富集通路。
6. 注6：富集通路的生物学解释有助于理解差异基因在生物过程中的潜在作用和机制。
7. 注7：避免使用 MS Excel 处理基因-通路注释文件，以防止数据格式和内容的意外更改。
8. 注8：设置数据框的行名为通路名称有助于后续的可视化和数据操作，确保每个通路都有唯一标识。

笔记

1. 注1：GEO 系列是由原始提交者提供的数据集。当 GEO 员工策划该数据时，它就成为了一个策划过的 GEO 数据集。每个系列和数据集都包括在特定平台（本例中为 Illumina HumanHT-12 V3.0 表达珠芯片）上运行的一组样本（其 ID 以 GSM 开头）。
2. 注2：设置组时使用的顺序很重要，因为 GEO2R 对所选组不敏感。分析将始终独立于组的名称执行“第一组输入” vs “第二组输入”。
3. 注3：在微阵列中，一个探针集应设计为靶向一个基因，但当一个基因家族包含保守区域时，同一个探针集可能会靶向不同的相关基因。另一方面，按设计，一个基因被多个探针集靶向，这些探针集对基因的不同区域对齐。这可能导致不同的探针集靶向基因的不同转录本，可能具有不同的表达水平 [14]。在某些情况下，探针集不映射到已知基因，需要从基因注释分析等基于基因的分析中排除。尤其重要的是，在出版物中跟踪并提供生成感兴趣基因列表的探针集。
4. 注4：必须注意用户正在处理的数据类型。与通路数据库类似，不同的微阵列具有不同的基因覆盖率。一些旨在靶向整个基因组（例如 Affymetrix Human Genome U133 Plus 2.0 Array），而另一些则专门针对一组基因（例如 CustomArray Signaling Pathway platform version 3 是针对关注信号传导通路的约 2000 个基因的阵列 [16]）。在后一种情况下，使用仅限于信号传导特定数据库进行富集分析会更合适。同样，取决于实验设计，某个数据库可能比其他数据库更合适。例如，专注于代谢的研究将受益于使用代谢特定数据库，而研究与疾病相关的基因则可能受益于使用通用数据库。
5. 注5：通路富集分析可以使用不同的技术进行，这些技术可以分为三大类：过度表示分析（ORA）、基因集富集分析（GSEA）和基于拓扑的通路富集分析（TPEA），如 [17] 中所述并进行比较。简而言之，ORA 是最常用的方法，因为它简单；它只需要一个感兴趣基因列表，计算列表中属于某个通路的基因比例（通常使用超几何检验或 Fisher 检验）。GSEA 考虑通路内的所有表达变化，以提供该通路受到影响的概率。提供的基因列表需要进行排序，以及排序值（例如，折叠变化）。TPEA 涵盖使用通路内基因之间的连接来计算通路富集和激活状态的方法。这些方法特别容易受到数据库选择的影响，因为拓扑层增加了异质性 [12]。
6. 注6：获得的通路数量可以是可管理的，如本例，或可能过多，如 [18]。在后一种情况下，将单个通路映射到更高层次的通路本体（例如，Pathway Ontology [19] 或 Reactome）可以帮助减少复杂性并突出显示感兴趣的过程。
7. 注7：在第一个表格中，因为我们仅处理通路名称，任何工具/软件都可以使用，包括 MS Excel。然而，由于第二个表格包含基因-通路对注释，MS Excel 不应被使用，因为它会默认将特定的基因符号转换为日期 [20]。
8. 注8：设置行名的常用命令是 rownames。然而，在来自不同数据库的通路名称的情况下，可能有少数不同的通路具有相同的名称，这会导致 R 中出现错误。因此，我们使用 make.names 命令，它将为行赋予与通路相同的名称，除非有两个（或更多）通路具有相同的名称（如本例中的 ‘Chemokine signaling’）。在后一种情况下，具有相同名称的行名将带有索引以进行区分（在本例中，一个行名为 ‘Chemokine.signaling’，另一个为 ‘Chemokine.signaling.1’）。