超过2500个全癌基因组的通路和网络分析

Pathway and Network Analysis of More Than 2500 Whole Cancer Genomes

Nature Communications（IF：12.121），2020.02.05

导语

GUIDE ╲

在过去的十年中，癌症基因组测序，如TCGA，已经确定了数以百万计的体细胞畸变（aberrations）；然而，这些畸变的注释和解释仍然是一个主要的挑战。具体来说，虽然某些体细胞畸变经常发生在特定类型的癌症中，但在中等规模的患者群体中，稀有畸变存在“长尾巴”分布情况，难以与随机passenger畸变区分。在许多癌症中，相当比例的患者在蛋白质编码区没有已知的driver突变，说明另外的driver突变仍未被发现。绝大多数已知的driver突变影响蛋白质编码区域，只有少数反复出现的非编码driver基因突变，最明显的是TERT启动子突变。在其他研究中，全基因组分析发现了几个调控元件中的频发突变（recurrent mutations），而表达定量性状遗传位点（eQTLs）分析发现了一些癌症类型中与基因表达变化相关的非编码体细胞突变。

背景介绍

癌症driver突变通过改变hallmark通路的活性来解锁细胞的致癌特性。因此，癌症基因已经被证明聚集在少量的细胞通路和相互作用的子网。因此，通路和网络分析已被证明通过基于不经常突变的基因的通路成员和与频发突变形成的物理或调控互作来揭示不经常突变的基因作为癌症基因是有用的。然而，在已知或新的通路中，编码和非编码driver突变之间的相互作用尚未被系统地探索。本工作对 Pan-Cancer Analysis of Whole Genomes（PCAWG）的一部分，即27种肿瘤类型的2583个肿瘤的编码和非编码体细胞突变进行了通路和网络分析。PCAWG联盟整理了38肿瘤类型的2658例样本的全基因组测序数据。这项工作提供了迄今为止最大的经过统一处理的癌症全基因组集合，这些全基因组包括生殖系和体细胞变异，这些变异来自与人类基因组 （reference build hs37d5）对齐的重新分析的测序数据，使用的是标准化和高精度的方法。最近来自ICGC的PCAWG项目的工作显示，在对单个基因和调控区域的分析中发现很少有频发的非编码driver。本工作采用了七种不同的通路和网络分析方法，并从这些方法的预测中获得了与通路相关driver（PID，pathway-implicated driver）基因的一致基因集。本工作检查了编码和非编码突变对改变生物过程的相对贡献。

数据介绍

1. 样本：PCAWG的一部分，即27种肿瘤类型的2583个肿瘤样本数据。

2. 通路和突变数据：使用PCAWG Drivers and Functional Interpretation Working Group分析的gene scores，结合本工作通路和网络分析中的几个通路和互作网络分析。

术语“pathway methods”是指利用相关基因集进行分析的方法。

术语“network methods”是指利用基因和/或基因产物之间的成对相互作用的方法。

3. 体细胞突变数据

使用PCAWG（无皮肤黑色素瘤和淋巴样本集）Drivers and Functional Interpretation Working Group分析的编码和非编码（core promoter， 5′ UTR， 3′ UTR， enhancers）的consensus driver P值。纳入了≤5个基因靶标的增强子，覆盖了89%的PCAWG Drivers and Functional Interpretation Working Group分析的增强子元件。在PCAWG Drivers and Functional Interpretation Working Group分析报告的同一基因组元件有多个P值的情况下，对该元件使用最小的报告P值。

4. 通路和网络数据库

（1）通路方法使用6个数据库的基因集：CORUM、GO、InterPro、KEGG、NCI Nature和Reactome

（2）网络方法使用3个互作网络数据：

①ReactomeFI 2015互作网络的最大连接子网，把它当作无定向的。

②iRefIndex14交互网络的最大连接子网，用KEGG通路的相互作用来增强它。

③BioGRID互作网络，也被用于评价和注释结果。

方法介绍

1. 推断gene scores

通路数据库和基因互作网络通常记录单个基因水平上的信息。因此，本工作通过结合编码和/或非编码基因组元件中的PCAWG driver P值，形成编码和非编码基因评分。

对于分析个体突变的通路和网络方法，使用了来自PCAWG MAF的与PCAWG Drivers and Functional Interpretation Working Group 分析相同的基因组元件（elements，即promoter、 5′ UTR、3′ UTR和enhancers）。

2. 个体通路和网络算法

对gene scores和突变数据应用了七种通路和网络方法。使用了两种通路方法：ActivePathways（这个方法我们公众号之前的文章有介绍过，①ActivePathways整合多维组学通路分析，②多维组学通路分析R包ActivePathways的使用方法及Cytoscape绘制网络图的实用教程）和超几何分析（Vazquez）。使用了五种网络方法：CanIsoNet、Hierarchical HotNe、an induced subnetwork analysis （Reyna and Raphael， in preparation）、NBDI22和SSA-ME。表1显示了每种方法使用的通路和网络数据库。

利用这些通路和网络数据库，在GS-C、GS-N和GS-CN基因评分上运行每种方法，以确定三个相应的基因列表。

3. 非编码增值（NCVA，non-coding value-added）过程

GS-CN结果同时利用了编码和非编码突变数据，改进了弱通路和网络信号的检测。设计了一个非编码增值程序来分离编码信号和非编码信号，从而得到一组NCVA基因，其中非编码突变数据对GS-CN结果的发现做出了显著贡献。使用permutation test来评估GS-CN结果中基因的统计学意义。

4. 通路和网络方法的一致结果

为每一组结果定义了一组一致的基因：GS-C结果、GS-N结果、GS-CN结果和GS-N结合NCVA的结果（通过7种通路和网络方法）。具体来说，如果该基因被绝大多数通路和网络方法发现（≥4/7），那么将该基因定义为一致性基因（consensus gene）。在本工作的分析中，侧重于一致的GS-C结果，称之为编码变异的通路牵连driver基因（PID-C，pathway-implicated driver genes with coding variants），以及一致的GS-N结果与NCVA结果结合，称之为非编码变异的通路牵连driver基因（PID-N，pathway-implicated driver genes with noncoding variants）。将GS-C结果中87个基因定义为PID-C，将GS-N结果中93个基因定义为PID-N。进行了几项分析来评估PID-C和PID-N基因的生物学相关性。

5. 识别PID基因突变特征

对PCAWG的mutation signatures analysis中的突变特征进行了基于permutation的富集分析。确定了PID-N基因中每个非编码突变最有可能的突变特征，并将其与随机选择的non-PID-N中的非编码突变进行比较。

6. 改进的PID基因的网络邻居分数

为了评估PID基因的分数在多大程度上有助于通过通路和网络方法检测PID基因，考虑了每个PID基因的分数对其在BioGRID互作网络中的网络邻居得分的贡献。对于每个PID基因g，使用Fisher方法将g的一阶网络邻居的基因分数与g本身的分数合并。其中，对于基因g，设p（g）为g的基因得分，N（g）为g的网络邻域。

①当包含基因g时，g的网络邻域得分：

②当不包含基因g时，g的网络邻域得分：

如果g的网络邻域有g时的P值小于没有g时的P值，那么基因g提高网络邻居的得分，这表明基因g得分通过通路和网络方法对其进行检测；如果g的网络邻域有g时的P值大于没有g时的P值，那么基因g弱化网络邻居的得分，这表明基因g的网络邻居的分数主要是负责通过通路和网络方法对g检测。

7. PID基因表达分析

评估了是否每个PID基因突变状态与RNA表达相关。使用PCAWG-3基因表达数据，其平均值来自TopHat2和STAR-based的比对，并进行FPKM-UQ归一化。知道肿瘤类型和拷贝数畸变是基因表达的协变量，因此考虑了肿瘤类型并注释拷贝数变异。

（1）使用以下程序来评估个体肿瘤的表达相关性：

①只考虑具有至少三个突变样本和三个非突变样本的情况，以将分析限制在具有足够统计效力的情况下。

②对于每个PID-C基因或PID-N基因中的每个非编码元件，我将带有表达数据的样本划分为元件发生突变的A组样本和元件没有突变的B组样本。

③使用Wilcoxon rank-sum检验两组表达差异，并用Benjamini–Hochberg校正。

（2）使用以下程序来评估肿瘤类型之间的表达相关性：

①只考虑具有至少一个突变样本和一个非突变样本的案例，以限制分析具有足够的统计效力的案例。

②对于每个PID-C基因和一个PID-N基因中的每个非编码元件，将具有表达数据的样本划分为素发生突变的c样本群中突变样本的Ac集和未发生突变的c样本群中的Bc集。

③使用c样本群中未突变样本的表达将，表达值转换为 z-scores，Wilcoxon rank-sum检验基因集

的表达差异，C是所有样本中包含元件突变样本本集。Benjamini–Hochberg校正显著性。

结果解析

01

编码和非编码driver突变的“长尾”

首先分析了来自27种肿瘤类型的2538例ICGC PCAWG肿瘤样本中，全基因组测序发现的单核苷酸变异（SNVs）和短插入和缺失（indels）靶向的基因。通路和网络分析集中在2252个肿瘤子集上，这些肿瘤排除了黑色素瘤和淋巴瘤，因为它们在调控区域的非典型突变分布。利用PCAWG Drivers and Functional Interpretation Working Group分析的单个蛋白编码和非编码元件的泛癌driver P值，包括外显子、启动子、未翻译区域（5 ' UTR和3 ' UTR）和增强子。该分析整合了16种driver发现的预测方法，使得编码和非编码elements的driver P值一致。根据阳性选择、突变的功能影响、区域突变率、突变过程和特征等不同的方法，个体基因和非编码元件的P值表明它们作为driver的统计显著性。在蛋白编码driver P值的泛癌样本中，有75个基因具有统计学意义。这些数字与之前的报告一致，即在癌症类型之间，driver基因存在“长尾”，其中高度突变的基因很少，而突变罕见的基因很多。非编码突变表现出类似的长尾分布，显著基因更少。没有一个基因同时具有显著或接近显著的编码和非编码driverP值，这表明非编码突变以一组互补的基因为靶点，即编码突变。早期的研究表明，含有编码driver突变的蛋白质在分子通路和网络中相互作用的频率明显高。本工作观察到显著突变的编码和/或非编码元件之间存在大量的相互作用，这表明该通路和网络方法可能是有用的。

02

潜在driver突变的通路和网络分析

本工作使用PCAWG Drivers and Functional Interpretation Working Group计算的单元件（single-element）driver P值作为输入，对癌症drivers进行了全面的通路和网络分析。采用了七种不同的通路和网络方法，每一种都利用了来自分子通路或蛋白质相互作用网络的信息，来放大单元件分析中的微弱信号（图1）。

由于癌症中非编码体细胞突变的优先排序尚未解决，因此很难预先知道哪种分析方法最适合通过汇集不同通路或网络的微弱信号来区分drivers和passengers。因此，本工作对多种方法形成了consensus，遵循机器学习的“群体智慧”集成方法来提高结果的特异性。纳入了使用不同通路或网络信息来源和不同优先级标准的方法。使用术语通路牵连driver（PID）基因来描述候选driver基因。本工作的通路或网络分析方法使用不同的先验知识来源（例如pathway databases or interaction networks）以及输入数据（如driver p-values， point mutations， and/or gene expression），并依赖不同的技术来集成这些数据源。发现这7个方法的输出只有少量的重叠，说明采用非均匀加权来减轻冗余方法是不必要的。

仅使用编码突变，识别了87个PID-C基因，其中包括68个由COSMIC Cancer Gene Census（CGC）数据库编目的癌症基因（图2a， c）。PID-C基因编码基因得分显著高于non-PIDC基因，87个PID-C基因均高于其网络邻域得分。这一网络邻域分析表明，PID-C基因不仅与其网络邻域有关，而且其自身对其通路和网络方法的发现做出了重要贡献。87个PID-C基因还包括31个在PCAWG Drivers and Functional Interpretation Working Group分析上没有统计学意义的基因（图2a， c），说明网络邻域可以指定突变不频繁的基因，即那些在“长尾”的基因，可能是driver基因。有趣的是，根据CGC数据库，这31个FDR > 0.1基因中有13个也是已知的driver基因。因此，共识通路和网络分析重新获得了许多已知的蛋白质编码driver突变，并确定了其他可能的driver，这些driver不经常突变，因此低于基因特异性driver分析的统计显著性阈值。

仅使用非编码突变，使用一致通路和网络分析确定了一组62个基因。当对编码突变和非编码突变进行联合分析时，发现编码突变有更强的信号，并支配着编码突变和非编码突变的联合信号。为了提高检测非编码突变贡献的敏感性，设计了一种“非编码增值”（NCVA）程序（图1）。NCVA程序询问的是编码突变是否会增强对潜在非编码driver基因的发现，而非编码突变是否会增强对潜在非编码driver基因的发现。该程序识别了另外31个基因，与62个单独发现的非编码突变基因合并时，产生了93个非编码变异的通路相关driver基因（PID-N）。PID-N基因表现为一个健壮的以及和生物学相关的集合，不受任何特定突变过程的影响，反映了特定的致癌物质或DNA损伤过程。93个PID-N基因包括19个CGC的癌基因（图2 b， d），比预期的3.2个基因显著富集。排除的八个PCAWG Drivers and Functional Interpretation Working Group分析基因中有单个显著性的非编码元件的基因，即是PID-N又是CGC基因，超过3.1预期基因显著丰富，说明非编码突变可能会改变部分样本中重复编码或结构变异的基因。PID-N基因的非编码基因得分显著高于non-PIDN基因，92/93个PID-N基因（除HIST1H2BO外）提高其网络邻域得分。这表明PID-N基因不仅仅与它们的网络邻居有关。绝大多数PID-N基因（90/93，包括19个CGC基因）与PCAWG Drivers and Functional Interpretation Working Group分析（图2b）的不同，只有三个共同基因：TERT、HES1和TOB1。在这三个基因中，只有TERT在CGC数据库中记录为已知的癌症基因。

综合来看，与PCAWG Drivers and Functional Interpretation Working Group的分析相比，PID-C和PID-N基因多包含121个基因，90个新的可能的非编码driver。总的来说，PID-N基因的非编码突变比PID-C基因多151个样本，大多数PID-N基因的编码突变和大多数PID-N基因的非编码突变是克隆的。此外，绝大多数的PID-N基因与 PID-C基因，只有5个共同基因：CTNNB1、DDX3X、SF3B1、TGFBR2和TP53。虽然这表明编码和非编码driver突变发生在很大程度上不同的癌症基因组中，但下文显示，这两种类型的突变影响了许多具有相同标志的癌症过程的基因。

03

非编码突变对基因表达的影响

编码突变可以通过改变转录因子结合位点或其他类型的调控位点来起作用。因此，接下来评估了PID-N基因的非编码突变是否与同一基因的顺式表达变化相关。发现，在可使用RNA-Seq数据检测的90个PID-N基因中，有5个PID-N基因在顺式表达相关性方面表现显著（图3a）。毫不奇怪，在PID-N基因的顺式表达相关性中，最重要的是TERT启动子突变和表达增加之间的相关性，在11个ThyAdenoCA肿瘤、11 CNS-Oligo肿瘤和22 CNS-GBM肿瘤中发现了这种相关性。值得注意的是，这些关联是有限的，因为一些有TERT突变的样本无法获得RNA表达数据，而且启动子区域的低测序覆盖率限制了TERT启动子突变的检测。

发现其他四个PID-N基因的顺式表达相关性显著：TP53、TLE4、TCF4和DUSP22（图3）。在6种不同类型的肿瘤中，TP53启动子突变的6种肿瘤中，TP53表达显著降低（图3a）。突变样本表达的降低与众所周知的TP53作为抑癌基因的作用是一致的，并且TP53启动子的甲基化和表达两者之间存在联系在之前已经被研究过。本研究还通过PCAWG Drivers and Functional Interpretation Working Group解释描述了这种表达变化。在3例伴有TLE4启动子突变的肝细胞癌中，TLE4的表达显著降低（图3b）。TLE4是一个转录共抑制因子，它结合到几个转录因子，并且在急性髓系淋巴瘤中，TLE4通过与Wnt信号的相互作用作为一个肿瘤抑制基因。在急性髓系淋巴瘤细胞系中，TLE4沉默增加了细胞分裂率，而强制表达TLE4诱导细胞凋亡。然而，TLE4在实体肿瘤中的作用尚不清楚。在三个具有TCF4启动子突变的肺鳞状细胞癌中，TCF4表达显著降低（图3c）。TCF4是TCF4/ β-catenin复合体的一部分，并编码Wnt信号通路下游的转录因子。在肺鳞状细胞癌中观察到TCF4低表达。最后，DUSP22在5例伴有DUSP22 3’UTR突变的肺腺癌患者中显著降低表达，在3例伴有DUSP22 5’UTR突变的肺腺癌患者中显著过表达。这些UTR突变是相互排斥的。DUSP22编码一种磷酸酶信号蛋白，最近被认为是淋巴瘤的肿瘤抑制因子。

虽然这些基因表达相关性为PID-N基因的子集提供了额外的支持，但突变的变异等位基因频率和基因的拷贝数是基因表达的额外协变量。本工作发现，这些协变量在确定的相关中并没有发挥作用：每个PID基因的大多数突变都是克隆的，拷贝数的变化并不影响上述5个PID-N基因的表达相关性（图3）。此外，出现相关基因表达变化的PID-N基因数量较少，这是因为PIDN基因突变的样本数量较少，不同肿瘤类型的表达数据可用性不均匀，以及启动子区域序列覆盖的减少。这些问题进一步减少了非编码突变和RNA表达的样本数量，限制了顺式基因表达相关性分析的能力。

04

编码和编码突变的模块化组织

接下来确定了编码突变、非编码突变或两种突变组合改变的特定蛋白质-蛋白质相互作用子网络和生物通路。利用置换检验发现在PID-C基因之间有显著多的相互作用，PID-N基因间的相互作用数量接近显著，在PID-C和PID-N基因之间的相互作用显著多。将涉及PID-C和PID-N基因的相互作用的子网络组织成五个生物学过程：核心drivers、染色质组织、细胞增殖、发育和RNA剪接（图4a）。而PID-C和PID-N 之间的分子相互作用频率预料之中的很高，这是由于这种被用作通路和网络方法的信号，这些相互作用的组织说明了个体子网络中编码和非编码突变的相对贡献。

进一步的，使用g：Profiler web server描述了PID-C和PID-N中富集的分子通路（图4b）。总的来说， PID-C基因显著富集了63条通路，PID-N基因显著富集了 13条通路。进一步利用通路富集图中标注的通路之间的重叠，将这些分子通路分组为29个模块。对于每个富集模块，检查了PID-C、PID-N或两种基因是否与观察到的富集有关。这将模块和PID基因分成四个生物过程：染色质组织、细胞增殖、发育和RNA剪接（图4b）。发现染色质和细胞增殖过程中的通路——包括染色质重塑和组织、组蛋白修饰、凋亡信号、信号转导、Ras信号和细胞生长——主要是由于PID-C 基因的编码突变而改变的。同时，发现多种信号通路同时包含大量的PID-C和PID-N基因，这表明非编码突变在干扰这些通路时提供了一种替代编码突变。特别是Wnt信号通路，它主要由编码基因靶向，也有一些PID-N基因的非编码突变。

由于没有发现剪接相关的PID-N基因非编码突变和表达改变之间存在显著的顺式关联，因此探索了这些基因非编码突变和其他基因表达之间的潜在反式效应。发现，剪接相关的PID-N基因的非编码突变在很大程度上概括了最近发表的一项TCGA PanCanAtlas分析结果，是几种剪接因子的编码突变与47条通路的差异表达之间的联系（图5）。特别是，对通路之间的差异表达模式进行层次聚类，在RNA剪接基因中发现了三组突变（C1、C2和C3；图5a，b）。利用t分布随机近邻嵌入（t-SNE， t-distributed stochastic neighbor embedding），发现了高度重叠的聚类集（图5a，顶部注释栏），表明该聚类对聚类方法的选择具有较强的鲁棒性。通过 silhouette scores和bootstrapping法，进一步证实了聚类的稳健性。这些簇中每一个在剪接基因SF3B1、FUBP1和RBM10中都包含至少一个编码突变，正如之前报道的，同时在剪接相关的 PID-N基因中也包含非编码突变，表明这两种突变导致了相似的基因表达特征。剪接因子编码和非编码突变的联合分析也发现了富集的通路中的两组（图5a）。一组（P1）的特征是免疫细胞signatures，另一组（P2）主要反映细胞周期、DDR和必要的细胞机制的细胞自主基因signatures。几种PID-N剪接因子非编码突变的基因表达特征与之前报道的剪接因子基因编码突变的基因表达特征之间的相似性支持了剪接相关的PID-N基因在改变类似基因表达程序中的功能作用。除了上述模块外，还发现转录因子在PID-C和 PID-N基因中都有很好的体现。共有9个PID-C基因是转录因子，19个PID-N基因为转录因子。这一观察结果表明，非编码突变可能影响转录调控网络。

本工作采用了七种不同的通路和网络分析方法，并从这些方法的预测中获得了与通路相关driver基因的一致基因集。使用非编码变异和编码变异识别通路相关driver基因。这两组PID基因，都包含很少的与已知的癌症基因相互作用的突变基因，121个新的PID-N和PID-C基因被揭示为有希望的候选基因，扩大了癌症中driver突变的范围。还检查了编码和非编码突变对改变生物过程的相对贡献。分析表明，在非翻译区和顺式调控区的体细胞非编码突变在编码突变方面构成一组互补的遗传扰动，影响多种生物通路和分子相互作用网络，应该进一步研究它们在癌症发生和进展中的作用。

引用：

Reyna MA， Haan D， Paczkowska M， et al. Pathway and network analysis of more than 2500 whole cancer genomes. Nat Commun. 2020;11（1）：729. Published 2020 Feb 5. doi：10.1038/s41467-020-14367-0