基于表达谱的拓扑数据分析识别癌相关的遗传变异

Identification of relevant genetic alterations in cancer using topological data analysis

Nature Communications（IF：12.121），2020.07.30

肿瘤靶向治疗的关键基础是确定肿瘤发展和维持所必需的分子机制。大规模的剖面的癌症分子研究，如The Cancer Genome Atlas （TCGA）和International Cancer Genome Consortium（ICGC）通过系统地汇编许多肿瘤的遗传变异（genetic alterations）来实现这种识别。出现复发性基因或通路改变的肿瘤可能是由常见的分子机制驱动的。通过利用计算方法寻找正向选择的特征，这些研究已经产生了大量频繁突变的癌症相关基因，这些研究也揭示了大多数癌症突变发生的频率很低（<10%的样本），包括潜在的可操作的治疗靶点。

背景介绍

使用recurrence-based方法来识别低频率的癌症相关突变需要大量的样本来实现统计能力，以及建模背景突变率固有复杂性。一种癌症类型的中性突变率在患者、基因组区域或突变类型之间存在显著差异，这限制了recurrence-based方法的效力。据估计，就目前正在进行的剖面研究的规模而言（通常包括少于1000名患者），只有发生在中高频率（>15%）的癌症相关突变可以用recurrence-based方法完全获得。recurrence-based方法是必要的，它可以模拟罕见事件。有一种识别癌症相关基因的方法不是基于对突变率的建模，而是整合来自肿瘤的其他类型的数据。如果突变伴随着拷贝数、基因表达和/或甲基化的一致变化，就有可能利用这些变化将突变事件与癌症进展联系起来。一些研究已经利用突变基因的拷贝数、表达或甲基化（cis-effects，顺式效应）的变化来识别新的癌症相关突变。然而，基于突变基因以外的基因变化（trans-effects，反式效应）来识别与癌症相关的突变更具挑战性，通常需要提供已知基因-基因关系的信息，以减少误报的数量。DriverNet、OncoIMPACT、CaMoDi和Xseq利用基因表达中的反式效应来识别癌症相关突变。为了限制表达空间的维度，这些方法使用表达模块（一组共同表达的基因、功能相关的基因或基因网络）。然而，目前缺乏一种充分考虑表达空间复杂性的方法。

目前使用表达数据识别癌症相关突变的方法对几种混淆效应很敏感。染色质开放的基因组区域更容易被DNA修复酶进入，导致基因表达水平与突变率之间的关系不明显。此外，具有不同表达特征的肿瘤，如基因组不稳定肿瘤，可能有不同的突变率。这些影响导致突变和表达特征之间的虚假相关性，是当前算法的一个误报源。

为了解决这些问题，本工作设计了一种方法，从多个肿瘤的表达数据来识别癌症相关的突变基因。利用拓扑数据分析（topological data analysis，TDA）来重建表达空间的结构，并在评估突变基因的显著性时考虑到上述伪效应。将其应用于12种肿瘤类型的4476例患者的突变和表达数据，识别出95个突变的癌症基因，其中38个为既往未报道的低频（low-prevalence）基因（在同一组肿瘤样本群中突变平均频率为5%）。因此，提出了一种对recurrence-based方法的补充方法，能够识别难以捉摸但可能与临床相关的突变癌症基因。

数据介绍

TCGA的12种肿瘤类型的基因表达和体细胞突变数据（表1）。只考虑了两种数据都可用的患者。RNA-seq表达数据恢复到RSEM 格式，并根据公式r = log2（1 + 106×x）转换各基因的估计相对丰度（x）。从Broad Institute TCGA GDACFirehose Portal检索到组织体细胞突变数据。

方法介绍

1. 拓扑表现（Topological representations）

使用Ayasdi软件 （https://www.ayasdi.com/platform/）中的Mapper算法，来构建每种癌症样本RNA-seq数据的拓扑表现。Mapper建立在任何降维算法（也称为“过滤器函数”）上，以产生一种新的低维网络表现，其中保留了局部关系。为此，Mapper用重叠bins覆盖低维表现，并对高维空间中的点进行单链接聚类。bins的数量及其重叠分别由“resolution”和“gain”参数指定，每个bin中的clusters数量由Singh et al.描述的方法确定。然后通过给每个cluster分配一个节点来构建低维网络，如果一个样本出现在两个节点中，则通过一条边将它们连接起来。

2. 统计分析

使用了Rizvi et al.介绍的拓扑关联的概念，确定与表型空间局部区域相关的特征。在本工作中，测试的特征是肿瘤样本中的体细胞突变，表型空间是肿瘤样本的表达空间。对于样本中每个突变基因g，考虑以下分数:

Γ表示拓扑表现中的节点集，A_ij 是它的邻接矩阵，N是拓扑表现中的节点数量，e_i（g）节点i中样本g的非同义突变的平均频率。C（g）评分是网络各边的总和，其中每边的贡献与每边连接的两个节点中含有突变基因的肿瘤比例的乘积成正比。定性地说，这种非参数评分评估了携带特定突变基因的样本在基因表达上的相似性程度。特别是，由于它是基于图的距离，它可以安全地用于非线性空间，如肿瘤群体跨越的基因表达空间。它在减少数据的稀疏性的简化网络表现（如 Mapper representations）中的应用是明确的。此外，作为一种局部方法，它对异常值相对不敏感。

为了比较不同患病率的突变基因的得分，为每个基因进行了permutation test。将患者id随机分布在外显子组数据中，建立C（g）的零分布，并根据其零分布为每个基因g的得分赋一个p值。此外，将每次分析的基因数量限制在非同义突变和突变总数比例最高的350个基因。

3. 对与表达的虚假关联进行控制

高突变肿瘤通常有独特的表达特征。在这种情况下，表达空间的某些局部区域就会出现突变率较高的肿瘤。这些局部区域将隐藏着passenger突变的积累，可能会干扰本工作方法。为了控制整体表达模式和肿瘤突变率之间的关联，使用2. 统计分析一段中所述的相同方法评估了拓扑表现上突变肿瘤负荷的定位。如果突变负荷的定位显著（p< 0.05），则手动设置突变负荷的阈值，将样本分为高突变和非高突变肿瘤。从每一个高突变肿瘤中随机二次采样突变，以便在二次采样后样本中高突变肿瘤的中位突变负荷与非高突变肿瘤的中位突变负荷相等。然后重新评估显著性，使用下降采样（down sampled）数据定位突变负荷。如果突变肿瘤负荷的定位程度不显著，继续使用向下采样的突变数据进行分析。如果再次采样后定位程度仍然显著，则不考虑该样本。

为了控制同一基因内的表达和突变率之间的关联，评估了拓扑表征上体细胞突变和mRNA表达谱之间的相似性，计算拓扑表现中每个基因的表达和突变谱之间的Jensen–Shannon度。为了控制由于不同突变 calling centers的差异产生的批次效应，使用与2. 统计分析相同的方法评估了在拓扑表现中批次的本地化程度。

结果解析

01

低级别胶质瘤的拓扑表现

肿瘤的表达谱可以精确地描述为高维表达空间中的一个点，其中每个维度代表一个基因的mRNA水平，空间的维数由表达基因的数量给出。在这个空间中彼此靠近的点对应于具有相似表达谱的肿瘤。一种癌症类型的所有可能肿瘤的集合跨越了表达空间的子空间。在剖面研究中测量个体肿瘤的表达概况相当于从这个子空间中采样有限的点集。

本工作分析了TCGA的513例原发性低级别胶质瘤（LGG），其中RNA-seq和全外显子组DNAseq数据均可用。为了从这些RNA-seq数据中推断LGG的表达空间结构，使用了一种拓扑方法。拓扑学是研究空间的不同部分如何相互连接的数学领域。TDA将拓扑的一些概念概括为点集和成对距离。因此，TDA的目的是在给定有限样本点的情况下，推断和总结空间的拓扑结构。TDA最近被用于研究病毒重组、人类重组、细胞分化、乳腺癌和其他复杂的遗传疾病。本工作使用TDA算法Mapper，通过利用TCGA的表达数据构建LGG表达空间的低维表现（图1a）。Mapper生成表达空间的网络表现，其中每个节点对应一组具有相似表达谱的肿瘤。一个给定的肿瘤可以出现在多个节点上，如果两个节点有一个或多个共同的肿瘤，它们就通过一条边连接起来。与其他降维方法（如主成分分析和多维标度）相反，由Mapper产生的拓扑表现保持高维表达空间的局部关系。拓扑表示中任何相邻的两个肿瘤（以连接两个肿瘤的最短路径所包含的边数来衡量）都保证在原始的高维表达空间中相邻。使用Pearson’s相关性作为单个肿瘤表达谱之间相似性的衡量。LGG表达空间的拓扑表示由三个区域组成（图1a），与聚类分析中发现的表达亚型一致。然而，这些区域在拓扑表现上被稀疏结构（thin structures）连接，表明一些肿瘤具有多表达亚型的表达谱特征。

02

识别LGG癌相关基因

假设如果一个突变基因出现在表达空间的局部，它与肿瘤子集一致的整体表达模式相关，因此，它是肿瘤进展的一个候选驱动因素（图1b）。另一方面，在相同的基因组中，如果一个基因的突变作为阳性选择突变，它是clonally expanded的，但与癌症无关，它们将随机地分散在表达子空间中（图1b）。

为了验证这一假设，本工作实施了一种计算方法来评估肿瘤表达空间中非同义体细胞突变的定位。为了控制突变率和肿瘤表达谱之间存在的虚假相关性，评估了重建表达空间中的肿瘤突变负荷（每个肿瘤中的体细胞突变总数）的定位。基于此分析，对LGG中两个高突变肿瘤的突变进行了取样（n_mut > 102.5）。此外还评估了在重建表达空间中每个个体基因的表达和突变谱之间的相似性。在对这些伪相关进行校正后，在重构的LGG表达空间中显著定位16个突变基因（图1c）。这些基因包括众所周知的高患病率（>的15%）驱动基因，如IDH1、TP53、ATRX和CIC，以及一些低频突变基因，如NIPBL（在4%的肿瘤中发生突变）和ZNF292（在3%的肿瘤中发生突变），这些基因最近在一个大样本的胶质瘤分析中有报道。总的来说，16个显著突变基因中有15个是被报道过的，SYNE1基因是（2%的肿瘤发生突变）唯一的新候选者。没有观察到显著基因的显著性和患病率之间的显著相关。特别是，一些最显著的基因，如FUBP1、NOTCH1、PTEN、EGFR和NF1，在该肿瘤类型的患者中只有不到10%发生突变（图1c），这表明这些基因的突变与整体表达变化密切相关。这些结果在Mapper算法的参数空间中是稳定的（图1c）。

重构的LGG表达空间中显著基因的位置与Ceccarelli， M. et al分析的已知的成人弥漫性胶质瘤分子亚型一致（图1d）。特别值得注意的是，IDH2突变的肿瘤定位在少突神经胶质瘤的表达空间内，与IDH1突变的少突胶质细胞瘤表达谱不同（图1d）。这一观察结果与最近一项基于基因组变异的研究一致。神经标志（neuronal markers）表达在恶性（III/IV级）胶质瘤中有报道，除了典型的间变性神经节胶质瘤。在本工作的样本中，表达典型神经细胞标志物如神经丝（NEFL、NEFM和NEFH）和突触素（SYP）的肿瘤被显著定位在少突胶质细胞瘤的表达空间内。除了星形细胞胶质瘤特有的分子改变，如TP53和ATRX突变外，这些肿瘤还伴有arm 19q染色体的频繁缺失。虽然该组的平均肿瘤纯度估计值明显低于其他少突胶质细胞瘤表达组，但肿瘤含量估计值在很多情况下大于98%，表明神经标志物的表达并非由于肿瘤纯度差。

03

计算基准

为了评估通过本工作的方法确定的肿瘤相关基因的数量与样本大小的函数关系，在更小的样本集中重复了同样的分析，这些样本集由随机抽取原始LGG队列中的样本产生（图1e）。还通过生成随机数据集来评估假阳性的数量，在表达数据中对患者的标签进行了排列。观察到，本工作的方法需要至少100个肿瘤样本大小。对于较大的队列，预期的假阳性数在1到2之间（图1e）。

接下来，将本工作的结果与目前使用表达数据识别癌症相关基因的算法进行比较。为此使用最近发布的算法Xseq分析了相同的LGG样本集，Xseq是一个层次贝叶斯统计模型，通过预先计算的‘influence graph’来量化体细胞突变对表达谱的影响，该‘influence graph’编码两个已知的基因是否在功能上相关的。通过使用Xseq对LGG进行分析，只得到两个显著基因，其中只有一个（PTEN）在LGG中有报道。说明本工作的拓扑方法比最先进的算法的高灵敏度。除了Xseq之外，还将本工作的结果与MutSig2CV的结果进行了比较（图1f）。MutSig2CV建立了中性背景突变率模型，考虑了由于表达水平和复制时间的差异而引起的基因组变异。观察到MutSig2CV和本工作方法的结果有明显的重叠，23个突变基因中有15个根据MutSig2CV是显著的，根据本工作的方法也是显著的。通过MutSig2CV基于recurrence识别的一些最显著的癌症基因，如PIK3R1（4%的肿瘤发生突变），没有在本工作基于表达的方法识别到，这突出了基于recurrence和基于表达的方法的独立性。结合MutSig2CV和本工作综合拓扑方法的结果（图1f），筛选出新的low-prevalence突变基因，如NOTCH2，作为LGG中肿瘤进展的潜在驱动因素。

为了与现有方法进行更系统的比较，对多种肿瘤类型进行了一项与Bertrand et al.35相似的研究。基于前15个预测与癌症相关基因与金标准列表的重叠，估计了本工作的综合拓扑方法、Xseq、MutSig2CV、 OncodriveFML和 20/20+方法的精度、召回率和F1评分。除了本工作LGG样本，还分析了来自TCGA的208个结直肠腺癌（COAD）和930个乳腺浸润癌（BRCA）肿瘤样本。在三种癌型的样本中，本工作的综合拓扑方法的精度、召回率和F1评分在5种算法中最高或第二高，突出显示了它用于识别突变的癌症相关基因的效用。

04

识别12种癌症的癌相关基因

基于以上结果，决定将本工作的分析扩展到其他肿瘤类型。考虑了TCGA中12种肿瘤类型，因其有足够的样本（n > 140）和可用的RNA-seq和全外显子组数据（表1）。本工作分析的完整结果可以通过在线数据库访问。总的来说，本工作的方法识别了95个突变的癌症基因，其中16个基因在两种或两种以上的肿瘤类型中具有显著性（图2a）。一些最常见的基因有TP53、KRAS、HRAS、PIK3CA、ATRX、EGFR和NF1。各肿瘤类型中显著基因的数量与样本大小相关（图2b）。61%的显著基因出现在 Cancer Gene Census或OncoKB数据库的癌基因中（图2c）。75%的患者有一个显著基因发生突变，其中24%携带有被批准药物可操作基因突变。在整合拓扑分析的95个显著基因中，有38个根据MutSig2CV不显著（图2a），然而，这些推测的难以捉摸的癌症基因在MutSig2CV分析中确实经常显示出显著性倾向，这可能反映了样本大小的局限性。它们也与Cancer Gene Census和OncoKB有很大的重叠（图2d），以及参与发育过程的基因。结合本工作的分析结果和MutSig2CV的结果，可以对结肠腺癌中突变基因的研究进行优先排序，其中根据MutSig2CV的显著基因数量太多（n = 1698个基因）。特别地，本工作的分析强调了ARHGAP5和ARFGEF1是之前未报道的这类癌症中肿瘤进展的假定驱动基因。

05

肺腺癌中ADAMTS12的截断突变

利用TCGA生存数据，发现在以前未报道的癌症相关基因中，ADAMTS12的失活突变与较差的生存率相关（图3a）。ADAMTS12是一种血小板反应蛋白基的金属蛋白酶（metalloproteinase with thrombospondin motif），可阻断Ras-MAPK信号通路激活。当给免疫缺陷小鼠注射过表达ADAMTS12的A549肺腺癌（LUAD）细胞，与亲代A549细胞形成的肿瘤相比，其肿瘤生长不足。ADAMTS12基因位于染色体arm 5p中，该基因在60%以上的肺腺癌肿瘤中完全扩增。有研究表明，编码端粒酶催化亚基的TERT基因可能是扩增的靶点。与建议的ADAMTS12的抗肿瘤特性一致，本工作观察到与无染色体5p扩增者相比，LUAD患者染色体5p扩增且无ADAMTS12改变会有具有更好的总生存率（图3a）。相反，与ADAMTS12中有染色体5p扩增和截断突变的患者相比，有染色体5p扩增而没有突变的患者的生存率降低（图3a）。ADAMTS12的截断突变往往与染色体5p扩增同时发生（图3a）。

为了验证ADAMTS12失活是肺癌进展的驱动因素，利用shRNA在肺癌细胞系 LL/2-luc-M38中表达ADAMTS12，研究了ADAMTS12沉默的影响。体外增殖和侵袭实验显示，与对照细胞相比，转染shRNA质粒的细胞具有显著的增殖和侵袭潜能（图3b，c）。接下来生成ADAMTS12−/−小鼠，与对照组小鼠相比，ADAMTS12基因敲除小鼠的肺肿瘤数量增加了5倍（图3d）。与本工作的计算分析结果一致，ADAMTS12在肺癌中具有抑癌作用。

小编总结

本工作开发了一种方法，可以从多个肿瘤的表达数据来识别癌症相关的突变基因，通过利用拓扑数据分析来重建表达空间的结构。本工作对12种癌型的体细胞突变和表达数据进行了综合分析，发现肿瘤基因突变通常伴有表达的实质性变化，并发现38个新的的癌症相关候选基因，包括肺腺癌中金属蛋白酶ADAMTS12失活突变。小鼠实验显示，ADAMTS12−/−的肺肿瘤的易感性增加了5倍，说明ADAMTS12有肿瘤抑制作用。

引用：

Rabadán R, Mohamedi Y, Rubin U, et al. Identification of relevant genetic alterations in cancer using topological data analysis. Nat Commun. 2020;11(1):3808. Published 2020 Jul 30. doi:10.1038/s41467-020-17659-7