利用 CorExplorer 网络平台分析肿瘤基因表达因子

摘要

我们介绍 CorExplorer 网络门户，这是一个用于探索由机器学习算法 CorEx（相关性解释）发现的肿瘤 RNA 测序因子的资源平台，并展示如何将这些因子与生存数据、数据库注释、蛋白质-蛋白质相互作用以及彼此之间的关系进行分析，从而深入理解肿瘤生物学机制及治疗干预策略。

摘要

差异基因表达分析是理解疾病状态的重要技术。机器学习算法CorEx在分析肿瘤RNA-seq数据中基因群组的差异表达方面已显示出应用价值，可能有助于推动精准肿瘤学的发展。然而，CorEx会产生大量因子，这些因子在分析时可能难以与现有知识建立联系。为了促进此类关联，我们构建了一个名为CorExplorer的网站，使用户能够交互式地探索数据，并回答与其分析相关的常见问题。我们使用四种肿瘤类型的RNA-seq基因表达数据训练了CorEx，分别为卵巢癌、肺癌、黑色素瘤和结直肠癌。随后，我们将相应的生存数据、蛋白质-蛋白质相互作用、基因本体（GO）和京都基因与基因组百科全书（KEGG）通路富集结果以及热图整合到网站中，以关联因子图谱的可视化结果。本文通过示例方案说明如何利用该数据库结合这些外部数据，深入理解所学习到的肿瘤因子的生物学意义。

引言

自十余年前问世以来，RNA-seq 已成为测量基因表达的普遍工具¹。这得益于其能够快速且低成本地对样本的整个转录组进行从头分析。然而，RNA-seq 肿瘤数据反映的是内在复杂且常被采样不足的生物学背景，而数据本身又具有高维度和噪声大的特点。这为提取可靠信号带来了重大挑战。CorEx 算法利用多变量互信息，在此类情况下识别细微模式^2,3。该技术此前已被用于分析癌症基因组图谱（The Cancer Genome Atlas, TCGA）中的卵巢肿瘤 RNA-seq 样本，在此背景下，其表现明显优于更常用的分析方法⁴。

尽管RNA测序（RNA-seq）在包括肿瘤学在内的各类研究领域中已得到极为广泛的应用，但这些努力尚未推动其在临床干预中的广泛应用⁵。造成这一现状的部分原因在于，目前尚缺乏针对这些特定问题而设计的用户友好型算法和软件。为弥合这一差距，我们开发了CorExplorer网络平台，使来自不同背景的研究人员能够利用CorEx机器学习算法所发现的结果，研究肿瘤RNA-seq样本中的基因表达因子。CorExplorer平台支持对多种不同类型肿瘤（包括肺癌、结肠癌、黑色素瘤和卵巢癌）的因子进行交互式可视化和查询^6,7,8,9,10，旨在帮助研究人员梳理数据中的相关性，识别潜在通路，从而为治疗目的对患者进行分层。

我们预计 CorExplorer 门户网站可能对多种类型的用户具有帮助。该门户在设计时考虑了用户需求，旨在帮助用户理解公共数据库中驱动肿瘤基因表达差异的广泛因素，并可能将单个基因表达谱置于具有相似特征的肿瘤背景下进行分析。除了本文概述的代表性实验方案外，CorExplorer 的研究还可作为提出假设的起点，用于进一步验证；也可用于比较和对照 CorExplorer 以外数据集上的 CorEx 分析结果，以及将单个肿瘤中一个或少数几个基因的病理表达特征与可能受到协同影响的更大基因群组相联系。最后，对于刚进入该领域的研究人员而言，它还可作为将机器学习应用于 RNA-seq 分析的用户友好型入门工具。

方案

1. 探索包含目的基因的因素

1. 打开网页浏览器，访问 http://corex.isi.edu，进入 CorExplorer 主页。
2. 在右侧的 快速链接 下方，点击 Ovarian (TCGA-OV) 旁边的 + 展开按钮，查看基于 TCGA 卵巢癌数据训练得到的 CorEx 因子图的摘要（如 图 1 所示）。也可点击其他项目进行比较。
3. 完成因子图查看后，点击 Lung (TCGA-LUAD)，进入肺癌 RNA-seq 的 CorExplorer 页面。
   1. 使用 CorExplorer 的“因子图”窗口，探索目标基因的 CorEx 因子图。
      1. 将鼠标指针移至因子图显示窗口上，使用鼠标滚轮或触控板缩放因子图，以查看图中细节，例如每个因子中最重要的基因以及不同层级节点之间的连接。或者，点击并拖动以移动视图区域或任意节点。
      2. 要查找目标基因（此处以 BRCA1 为例），点击因子图窗口顶部的 基因 下拉菜单，输入“BRCA1”以在下拉列表中选择该基因，然后按 Return 键，视图将自动缩放至因子 26，即与 BRCA1 相关性最强的因子。
      3. 将鼠标重新移至图显示区域，滚动以缩小视图，查看二级节点 L2_8 及其与因子 26 相邻的关联因子。请注意，仅显示权重高于 最小连接权重 滑块所设阈值的基因。
      4. 要查看与该因子关联的所有基因，点击 L1_26 节点，并在弹出窗口中选择 加载更多基因。当显示“Done”（完成）时，关闭弹出窗口。
      5. 返回因子图窗口上方的标题区域，抓取并拖动 最小连接权重 调节器。当滑块向下移动至 0.05 时，因子 L1_26 中的其他基因（包括 BRCA2）将按权重顺序显示。可选操作：通过抓取并拖动节点以优化布局。
   2. 通过生存分析窗口查询，确定患者按因子分层对生存结果的影响。
      1. 在生存分析窗口中，取消勾选 按 p 值排序，然后在 单因子 下拉菜单中选择因子 26，以显示因子 26 的生存曲线。
      2. 向下滚动生存图，以在横轴上显示处于风险中的患者数量。
   3. 通过注释窗口查询，寻找与生物学功能的关联。
      1. 在注释窗口中，若要按因子编号而非错误发现率（FDR）对 因子 下拉菜单进行排序，请取消勾选 FDR 排序。
      2. 滚动并点击注释窗口的下拉菜单，选择因子 26，以显示该因子的功能富集注释。
      3. 向下滚动注释列表，直到看到 DNA 修复，点击该项，即可在图中立即看到相关基因以黄色高亮显示。参见 图 2 中间面板。
      4. 请注意，当选择不同的 GO 术语（例如“DNA 损伤响应中的内在凋亡信号通路”）时，因子会根据是否富集了具有所选注释的基因而出现或消失。
   4. 通过添加具有不同功能的窗口，进一步探索因子。
      1. 从顶部菜单栏中，选择 添加窗口 下拉菜单中的 PPI，然后点击 添加 按钮，在显示区域中添加一个蛋白质-蛋白质相互作用网络（PPI）窗口。在 PPI 图窗口中，选择因子“Layer1: 26”以显示蛋白质-蛋白质相互作用。注意连接的密度。
      2. 从顶部菜单栏中，选择 添加窗口 下拉菜单中的 热图 而非 PPI，然后点击 添加 按钮，在显示区域中添加一个热图窗口。在热图窗口中，选择因子“Layer1: 26”以显示基因表达模式。
      3. 抓取并重新定位热图窗口，使生存分析窗口也可见。观察热图顶部的橙色/蓝色/灰色条带，其对应于生存图中的患者风险分层。结果如 图 2 底部所示。

2. 使用基因权重、生存和注释数据过滤和解释 CorEx 因子

1. 使用生存数据和聚类质量筛选感兴趣的因子。
   1. 在页面顶部的“Dataset”下拉菜单中，选择 TCGA_OVCA，进入TCGA卵巢癌RNA-seq的CorExplorer页面。
   2. 页面加载完成后，从生存分析窗口中注意到，不同分层之间生存差异最大的因子是114。
   3. 在因子图谱窗口顶部的Factor下拉菜单中，选择“Layer1: 114”。
   4. 用鼠标拖动链接权重滑块至0.5。注意因子114中基因数量较多（共1609个），但没有任何基因的权重>0.35，表明其聚类相对较弱。
   5. 接着，在生存分析窗口中展开因子列表，并在下拉菜单中选择生存差异次优的因子39，以查看其对应的生存曲线。
   6. 在注释窗口中点击选择因子39，系统将显示其显著相关的GO和KEGG注释。
2. 为了更深入理解因子39中基因的生物学功能，可按以下步骤利用邻近注释信息对因子进行解读。
   1. 在因子图谱窗口顶部的因子下拉菜单中选择“Layer1: 39”。然后将鼠标移至因子图谱窗口并缩小视图，以完整显示包含6个因子的L2_14聚类：14、32、39、42、52和82（如图3所示）。
   2. 为理解与L2_14节点相连各因子的相对重要性，首先查看每个L2_14因子的生存差异。在生存分析窗口中取消勾选Sort by p-val，然后依次点击各个因子编号。通过此操作可注意到，仅有因子14、32和39显示出与生存相关的显著性。
   3. 现在从顶部菜单栏再次选择Add Window下拉菜单中的PPI选项，点击Add，在显示区域添加一个PPI图谱窗口。在PPI图谱窗口中选择因子“Layer1: 52”，以显示显著的蛋白质-蛋白质相互作用。此时窗口布局的一个示例如图3所示。
   4. 点击PPI窗口底部的View at StringDB链接，跳转至StringDB在线数据库。在首个页面点击Continue，然后像之前一样选择网络图下方的Analysis标签页，以获取PPI网络基因的在线GO分析结果。其中最显著的细胞组分为“MHC II类蛋白复合物”。
   5. 返回CorExplorer标签页及PPI窗口，从因子下拉菜单中选择因子32，再次点击View at StringDB链接跳转至StringDB分析页面。此时最显著的细胞组分为“MHC I类蛋白复合物”，与上一步中因子52的II类结果形成对比！
   6. 最后，回到PPI窗口，从顶部的因子下拉菜单中选择“Layer1: 39”。
   7. 点击View at StringDB链接跳转至StringDB分析页面。在首个页面点击Continue，然后选择网络图下方的Analysis标签页，获取PPI网络基因的在线GO分析结果。可观察到最显著的分子功能为“CXCR3趋化因子受体结合”。

3. 利用生存数据和数据库注释寻找有前景的治疗组合

1. 从Dataset下拉菜单中选择TCGA_SKCM，切换至TCGA黑色素瘤CorExplorer界面。
2. 注意，具有最大生存差异的因子是因子171。通过滚动查看因子171的注释，可注意到“免疫应答”和“细胞因子介导的信号通路”位于注释前列（与之前卵巢癌顶部因子的情况一致）。
3. 为寻找一个互补因子，请查看与生存显著相关的顶部因子及其对应的顶级注释术语。为此，请点击顶部菜单栏中的Dataset overview链接，打开一个包含数据集处理详情的独立标签页，其中还列出了按生存差异p值排序的顶部因子摘要。请注意，第一个非免疫相关因子为因子88。
4. 返回至TCGA_SKCM浏览器标签页。
5. 在生存、注释和图形窗口中选择因子88。其前几项GO术语与“rRNA加工”和“线粒体组织”相关，证实该因子与免疫相关因子明显不同。
6. 在生存窗口的成对因子下拉菜单中，选择“88_171”，以观察同时具有中等表达水平的因子171和因子88的患者群体其生存情况的改善。注释与生存比较结果如图4所示。

4. 使用搜索页面查找不同肿瘤类型间基因表达变异的共性与差异

1. 点击 CorExplorer 返回首页
2. 点击 搜索 在顶部菜单栏上点击，进入可搜索 CorExplorer 网站所有数据集的页面。
3. 在 基因 搜索框中输入“FLT1”（VEGFR1）并点击 返回 或按压 搜索FLT1 在以下因素中权重相对较高：OVCA - 76、LUAD - 162、SKCM - 195 和 SKCM - 184，以及 COAD - 112 和 COAD - 74。
4. 或者，在所有数据集中搜索相关的 GO 术语。您可以在“GO 搜索”框中输入“angiogenesis”并按下回车进行尝试 返回 或按压 搜索所有FLT1因子（除SKCM-195外）均被列为在“血管生成”基因中具有统计学富集意义；因子195实际上也具有该注释，但低于默认的10-8阈值。本步骤及前一步骤的搜索结果如图所示。 图5.
5. 作为进一步的示例，在GO搜索框中首先输入“epidermal growth factor receptor”。仅有LUAD对该术语富集，这是肺癌中一个众所周知的分层因子。接着，在搜索框中输入“mesenchymal”。该术语在OVCA的基因表达组中富集，而在OVCA中，它是一个被广泛研究的分层因子。

结果

在肺癌数据集中搜索基因“BRCA1”发现，该基因与CorEx因子26关联最强（图2）。该因子的GO术语富集程度极高，其中DNA修复的FDR仅为1 × 10^-19。该选择还突出了第二层聚类L2_8，其下包含六个密切相关的子因子。在GO术语注释或因子图的GO富集下拉菜单中选择“DNA修复”，可高亮显示各因子中相关的基因，其中因子26所关联的基因数量远超其他因子，符合预期¹¹。蛋白质-蛋白质相互作用网络显示出强连接性，进一步支持因子26中基因功能的高度关联性。相应的生存曲线图提示其可能与患者生存相关，但需在更大规模的数据集中进一步验证。

从生存分析入手，可以剖析特定基因表达组别与生存改善之间的关联原因。例如，卵巢癌生存率的首要影响因素为第39因子，该因子显著富集了与免疫系统相关的基因（图3）。另有五个与同一二级节点相关的因子也被识别为免疫相关，但它们对生存的影响差异显著，其中因子39的影响最强，而因子52的影响最弱。为某一因子添加蛋白质-蛋白质相互作用窗口，可显示其直接的相互作用网络，并支持链接至StringDB¹²网站，以查询该PPI网络基因的多种富集情况。依次对L2_14的每个因子进行此操作后发现，StringDB对PPI网络基因的富集分析提示了以下可能的生存关联解释：因子32包含编码主要组织相容性复合体（MHC）I类蛋白复合物的基因，该复合物可被细胞毒性T淋巴细胞识别；因子39对应于细胞因子信号传导和CXCR3受体结合，与CD8+ T淋巴细胞相关。这两个因子在相应基因表达水平较高的患者中均表现出显著的生存优势。细胞毒性CD8+ T淋巴细胞主要负责抗肿瘤免疫。相反，因子52由编码MHC II类复合体蛋白的基因组成，主要被CD4+辅助性T细胞识别，而非直接被细胞毒性T淋巴细胞识别。其余L2_14因子反映的是广义的免疫系统激活，未区分两种淋巴细胞亚群。生存关联特异性地体现在细胞毒性T淋巴细胞对MHC I类细胞抗原的识别，这一结果与我们对肿瘤免疫（包括黑色素瘤等其他癌症）的普遍认知一致^13,14。

该网络门户支持发现具有互补功能的因子对，这些因子对可能提示有效的肿瘤特异性联合疗法。可通过扫描数据集概览，寻找与生存相关但具有不同基因本体（GO）富集特征的因子。以黑色素瘤（TCGA_SKCM；图4）为例，排名最高的生存相关因子171与免疫相关，而列表中靠后的因子88则在与线粒体组织相关的基因中表现出富集。事实上，已有研究提出将后者作为黑色素瘤的潜在靶点¹⁵。在CorExplorer页面中添加生存时间窗口，可比较使用因子对与单独使用每个因子进行分层的效果，结果显示，同时具备两组有利基因表达模式的患者群体，其生存趋势优于仅具备单一因子有利模式的患者。然而，最高风险分层的生存情况并未明显改善，提示对于部分患者，单独使用免疫治疗可能是最佳选择。

通过在数据集中搜索基因或基因本体（GO）术语，可以观察到肿瘤之间的共性与差异图5）。例如，FLT1（又称VEGFR1）是一个被广泛研究的促血管生成标志物^16,17当将其输入搜索栏时，所有肿瘤均含有FLT1发挥主要作用的因素。相反，当在搜索页面输入GO术语“angiogenesis”（血管生成）时，6个FLT1组中有5个显示出该富集结果。除SKCM-195外，所有FLT1因素均被列为在“血管生成”基因上具有统计学意义的富集。第六个因素实际上也具有该注释，但低于默认的10-8阈值。当在替代的富集分析工具（例如基因集富集分析，Gene Set Enrichment Analysis, GSEA）中利用因子列表内的权重进行计算时¹⁸第六个因子也被发现显著富集了“血管生成”相关基因。

检查热图以确保基因表达模式具有足够的质量来支持生物学解释非常重要。显示强烈且清晰变异的热图可能表现出因子基因从低到高协调表达的模式，或者更复杂的模式，其中一些基因的低表达与另一些基因的高表达相关（图6）。高质量聚类的一个关键标志是存在多个基因，其表达随因子评分呈平滑变化。因子热图中的样本按因子评分排序，因此应从左到右呈现平滑的梯度变化。然而，这种情况至少可能在两种不同情况下无法实现。最常见的是，相关性可能极度嘈杂（图5C），从而质疑关于生存和/或生物学功能推论的稳健性和实用性。此外，仅在极少数样本中出现的模式可能不符合CorEx算法所假设的三种表达状态模型，导致对样本的分类产生误导（图5D右侧）。

图 1：CorExplorer 首页。在快速链接下的卵巢癌旁边点击+后，将显示因子图的详细信息。CorEx 层次模型由底层的输入变量（本例中为基因表达）和上层推断出的潜在因子组成。请点击此处查看此图的放大版本。

图 2：使用基因名称引导探索。 本图展示了一系列截图，说明如何探索与 BRCA1 高度相关的 CorEx 肺癌因子。首先，在因子图的 基因 下拉框中选择“BRCA1”，图形视图将聚焦于 BRCA1 权重最大的因子。稍微缩小视图后，可看到第二层节点 L2_8，该节点将此因子与其他相关因子连接起来。可比较生存信息和注释：点击 GO 术语 DNA repair 可高亮显示已注释的基因。添加一个蛋白质相互作用（PPI）窗口，以展示该因子中基因的网络互作关系。通过点击 添加窗口 按钮添加热图，可显示表达模式与生存之间的关联，提示 DNA 修复基因的表达升高可能与生存率降低相关。请点击此处查看此图的放大版本。

图 3：利用临床数据（生存数据）指导探索。 探索与卵巢癌生存最相关的首要因子（因子39），可揭示其邻近因子之间有趣的关联关系。在因子图中选择因子39并稍作缩小后，可见第二层中与因子39相连的因子还与其他五个因子相关联。附加的生存分析窗口可直接比较这些相关因子的生存差异。因子39和因子32均显示出与生存呈正相关，而因子52则无此相关性。所有蛋白质-蛋白质相互作用网络均结构清晰。通过链接至StringDB可进一步比较基因本体（GO）注释（图中未显示）：因子39与细胞因子信号通路网络相关，该通路涉及细胞毒性CD8+ T淋巴细胞的激活；而因子32则主要由MHC I类抗原呈递蛋白构成，这类蛋白可触发上述淋巴细胞的识别；然而，其邻近因子则主要由其他免疫系统组分（如CD4+辅助性T细胞）主导，且未显示出与生存的相关性。请点击此处查看该图的放大版本。

图4：探索最重要的生存因子提示潜在的治疗组合。 在主页菜单栏上的“数据集”链接可进入一个按p值排序的生存因子简明表格，并附有最主要的GO注释（未显示）。利用黑色素瘤的这一信息，因子171（免疫功能）与因子88（线粒体组织）的组合似乎具有互补性。该图并列展示了每个因子的注释窗口，以便进行对比。根据这两个因子单独或联合对患者进行分层的生存曲线表明，相较于任一单一因子，联合使用可增大生存差异。请点击此处查看此图的放大版本。

图5：搜索页面支持泛癌分析。可通过主页上的搜索链接，在所有数据集中查找基因或GO生物学过程术语。该图展示了对基因FLT1和GO术语“angiogenesis”（血管生成）的搜索结果。结果显示，在多种癌症中，被注释为“angiogenesis”（血管生成）的因子均含有FLT1。请点击此处查看该图的放大版本。

图6：热图可用于根据因子得分定性评估基因与样本之间的相关性。 当患者按热图中的因子得分排序时，高质量的基因表达关系表现为平滑的渐变。最左侧的因子18热图即为一个示例。图中模式也可能包含复杂的上下调表达特征，如中间因子11的大热图所示。低质量的模式有时会显示部分患者的基因表达出现突变，如右侧因子9的热图所示，或表现为高度噪声的相关性，如右下角因子161的热图所示。 请点击此处查看该图的高清版本。

讨论

我们介绍了CorExplorer网站，这是一个公开可访问的网络服务器，用于交互式探索由CorEx算法从肿瘤RNA-seq数据中学习到的最大相关基因表达因子。我们展示了如何利用该网站根据肿瘤基因表达对患者进行分层，并说明了这种分层与生物学功能及生存预后之间的对应关系。

已有其他用于RNA-seq分析的网络服务器。在cbioPortal中，可对肿瘤进行差异表达和共表达分析，并将其与其他类型的数据进行整合^19,20。GenePattern²¹、Mev²²和Morpheus²³等服务器集成了成熟的聚类技术，例如主成分分析（PCA）、k均值（kmeans）或自组织映射（SOMs）。更具创新性的尝试包括CamurWeb²⁴，其基于一种自动规则生成的分类器；以及TACCO²⁵，该工具实现了随机森林分类器和套索回归（lassos）。本文所采用的CorEx算法通过优化多变量信息，以发现能够解释数据模式的因子层级结构。相较于通过PCA获得的线性全局因子⁴，这种非线性且具有层级结构的因子学习方法在可解释性方面表现更优。此外，该技术对样本信号的细粒度解析能力，使得相较于常用的宽泛亚型分类，能够实现更为精确的肿瘤间比较。这种重叠性与层级性因子分析的结合，使CorExplorer区别于大多数其他方法，同时也需要开发新的可视化与总结工具。

CorExplorer 因子分析的一个关键部分是能够探索的因子不仅限于几个，而是超过 100 个具有信息性基因模式的因子，这些因子被置于一个重叠的层级结构中。CorExplorer 有助于挖掘这些大量因子与生物学和临床特征之间的关联，并实现对单个肿瘤的异常精细的表征。由于如此大量的因子是通过无监督学习获得的，因此并非所有因子都与疾病生物学相关。在这种情况下，必须利用注释信息或已知基因来提取感兴趣的因子，或搜索与生存等临床数据相关的因子。因此，CorExplorer 允许用户执行这一非常重要的筛选步骤。肿瘤中存在特定因子的基因模式，甚至可能提示个性化肿瘤治疗的策略。此外，每个肿瘤具有多个因子评分，这有助于发现潜在的有用治疗组合。

有时，与生存高度相关的因子可能未显示出显著的 GO 注释。这可能是由于数据噪声或采样不足所致，但也可能存在其他原因，例如聚类尺寸过小而无法获得显著的富集分数，或该组别是由来自不同通路且无明确生物学关联的单个基因组成的“混合类别”。此外，可能需要采用不同于 KEGG 和 GO 生物过程的注释类别，例如细胞组分。这些分析可通过本方案中所述链接至 StringDB 获取。目前 CorExplorer 网站上的基因本体富集分析尚未考虑因子中基因的权重，但这一问题有望在不久的将来得到解决。请注意，在“添加窗口”下提供了一个基因列表选项，可用于下载完整的因子基因列表，以便使用外部工具进行进一步分析。

为网站目的，CorEx 对每个数据集运行了五次，并保留产生最高总体总相关性的那次运行结果。对多次运行结果进行统计学表征可能更具信息量，这也是未来工作的目标。此外，服务器上可用的肿瘤类型集合相对较小，但预计将根据用户需求随时间逐步扩展。

如上所述，CorExplorer 可视化 CorEx RNA-seq 因子之间的关系，并整合临床和数据库信息，从而支持多种不同的分析模式。我们希望该工具能够推动进一步的研究，以利用 RNA-seq 分析的强大能力，在肿瘤学领域实现发现和临床应用。

材料

本文使用的材料清单

姓名	公司	目录编号	评论
CorExplorer 网站的公共服务器	USC	http://corex.isi.edu	Intel Xeon E5-2690 4 核 2.6 GHz，8GB 内存。后端架构为 LAMP：Linux、Apache、MySQL、PHP。
网页浏览器	Google/Apple	Chrome/Safari	已验证的网页浏览器。