肿瘤信号通路生信分析如何发12分的NC

今天和大家分享的是2020年2月发表在Nature Communications(IF=12.121)杂志上的一篇文章Pathway and network analysis of more than 2500 whole cancer genomes。文章作者完成了2500多个完整癌症基因组的信号和网络分析。

Pathway and network analysis of more than 2500 whole cancer genomes

2500多个完整癌症基因组的信号途径和网络分析

（分享者：科研菌-桑葚）

一. 研究背景

在过去的十年中，蛋白质编码基因中癌症驱动程序突变的目录已大大扩展。然而，非编码的癌症驱动基因突变的特征较差，仅报道了少数复发性非编码基因突变，最著名的是TERT启动子突变。PCAWG汇总了来自38种肿瘤类型的2,658例癌症的全基因组测序数据，作者对PCAWG项目编辑的来自27种肿瘤类型的2,583个整个癌症基因组中的非编码区域突变进行了多方面的途径和网络分析，并成功进行了途径和网络分析而确定了蛋白质编码基因中的稀有突变的动机。

二. 分析流程

三. 结果解读

1. 编码和非编码的驱动基因突变

作者分析了来自27种肿瘤类型的2,538个ICGC PCAWG肿瘤样品中的单核苷酸变异(SNV)和短插入和缺失(indels)的基因，这些基因通过全基因组测序确定。作者的途径和网络分析是关注2,252个肿瘤的子集(这些子集排除了黑色素瘤和淋巴瘤，因为它们在调控区域的突变非典型分布)。

作者利用了PCAWG驱动程序和功能解释工作组分析中单个蛋白质编码和非编码元素的泛癌驱动基因p值，包括外显子，启动子，非翻译区(5'UTR和3'UTR)和增强子。该分析整合了来自16种驱动程序发现方法的预测，从而得出了编码和非编码元素的一致驱动基因p值。

在泛癌队列的蛋白质编码驱动基因p值中，有75个基因具有统计学意义，另外还有7个基因处于接近显著水平。这些数字与先前报道的驱动基因的长尾一致，很少有高度突变的基因，而很多基因在癌症类型之间很少发生突变。非编码突变表现出相似的长尾分布，甚至有更少的显著基因。没有一个单一的基因同时具有显著或接近显著的编码和非编码驱动程序p值，这表明非编码突变靶向互补的基因，成为编码突变。

2. 潜在驱动基因突变的途径和网络分析

作者使用由PCAWG驱动程序和功能解释工作组计算的单元素驱动基因p值，对癌症驱动程序进行了全面的途径和网络分析。作者应用了七种不同的途径和网络方法(ActivePathways，CanIsoNet，Hotarchical HotNet，超几何分析(Vazquez)，诱导子网分析，NBDI和SSA-ME)，它们各自利用了来自以下方面的信息分子途径或蛋白质相互作用网络(图1，Method部分)来放大单元素分析中的微弱信号。所有方法都根据随机数据进行了校准。

如图1，作者通过汇总PCAWG驱动程序预测中各个元素的驱动程序p值，包括注释的编码和非编码元素(启动子，5'UTR，3'UTR和增强子)。这些基因分数输入到五个网络分析算法(CanIsoNet，Hotarchical HotNet，诱导子网分析，NBDI和SSA-ME)中，这些算法利用了多个蛋白质-蛋白质相互作用网络，并通过两个途径分析算法(ActivePathways19和超几何分析(Vazquez))，利用多个途径/基因组数据库。

作者定义了一种非编码增值(NCVA)程序，以确定其非编码得分对组合编码和非编码分析的结果有重大贡献的基因，其中方法的NCVA结果会增加其在非编码数据上的结果。作者还定义了一种共识程序，将这七个算法确定的重要路径和网络结合起来。87个具有编码变异的途径相关驱动基因(pathway-implicated driver genes with coding variants ，PID-C)是大多数方法根据编码数据报告的一组基因。93个具有非编码变异的途径相关驱动基因(pathway-implicated driver genes with noncoding variants，PID-N)是大多数方法在非编码数据或其NCVA结果中报告的一组基因。只有五个基因(CTNNB1，DDX3X，SF3B1，TGFBR2和TP53)同时是PID-C和PID-N基因。

图1. 路径和网络分析方法概述

由于非编码体细胞突变对于癌症的优先顺序尚未解决，因此很难预先知道哪种分析方法最合适，作者通过汇总跨路径或网络的微弱信号来区分驱动基因和乘客基因。作者使用不同途径或网络信息来源和不同优先级标准的方法，每种方法指定的基因，以及具有可能的编码和非编码驱动程序突变的共有基因集，均定义为通过7种方法中的至少4种发现的基因，并使用术语“途径相关驱动程序”(pathway- implicated driver，PID)基因来描述这些候选驱动程序基因。通过补充实验得知，作者发现这7种方法的输出之间只有适度的重叠，这表明不必对共识进行加权以减轻冗余方法的影响。仅使用编码突变就确定了一组PID-C基因的87条途径相关的驱动基因。87个PID-C基因包括COSMIC癌症基因普查(CGC)数据库所分类的68个先前鉴定的癌症基因(见图2a, c)。

图2是途径和驱动程序分析可在驱动程序p值的长尾中识别驱动基因，以进行编码和非编码突变。图2a展示了250个基因的蛋白质编码元件上的驱动程序p值具有最重要的编码驱动程序p值；虚线分别表示FDR = 0.1和0.25。深绿色条是PID-C基因，浅绿色条是非PID-C基因。x轴下方的蓝色方块表示COSMIC癌症基因(CGC)人口普查基因。总共87个PID-C基因中的31个具有FDR>0.1。如图标记了几个PID-C基因，包括所有编码FDR> 0.1的CGC基因。PID-C和PID-N基因之间的重叠用星号表示。

图2a. 途径和网络方法可以识别出重要的编码驱动突变

图2c为根据用CGC基因编码驱动程序p值和PID-C基因，得到基因之间的重叠，排名靠前的具有统计学意义。绿色方块表示PID-C基因和CGC基因之间的重叠。

图2c. 编码驱动基因p值排名

根据图2a,c和补充实验还得知，PID-C基因的编码基因得分明显高于非PID-C基因，并且网络邻域分析表明，PID-C基因并非仅由其网络邻居牵连，而是它们自身通过途径和网络方法的重要贡献。87个PID-C基因还包括31个在PCAWG驱动程序和功能解释工作组分析中没有统计学意义的基因，这说明网络邻域可以提名具有罕见突变的基因，即长尾中的那些作为可能的驱动基因。

之后作者仅使用非编码突变，就可以通过作者的共有途径和网络分析来识别出62个基因，从而产生的基因数量要少于编码突变的基因。但是，当作者对编码和非编码突变进行联合分析时，作者发现了编码突变中更强的信号主导编码和非编码突变中的组合信号。为了提高检测非编码突变贡献的敏感性，作者设计了“非编码增值”(NCVA)程序(见图1)，NCVA程序询问编码突变是否增强了潜在的非编码驱动基因的发现超出了仅非编码突变所发现的范围。

该程序确定了另外31个基因，当与仅具有非编码突变的62个基因融合后，产生了93个具有非编码变异(PID-N)的途径相关的驱动基因。PID-N基因表现为健壮且具有生物学相关性的集合，不受反映特定致癌物或DNA损伤过程的任何特定突变过程的偏见。

如图2b，非编码元素(基因的启动子，5'UTR和3'UTR)上的驱动基因p值具有250个具有最重要的非编码驱动程序p值的基因；虚线表示FDR分别为0.1和0.25。深黄色条是PID-N基因，浅黄色条是非PID-N基因。蓝色方块为CGC基因。总共93个PID-N基因中的3个(TERT，HES1和TOB1)具有非编码驱动基因p值，其FDR≤0.1，而90个存在于FDR>0.1。如图标记了几个PID-N基因，包括在顺式基因表达中具有重要意义的PID-N基因(见图3)和所有非编码FDR>0.25的PID-N基因。PID-C和PID-N基因之间的重叠用星号表示。

图2b. 途径和网络方法可识别罕见的非编码驱动程序突变

如图2d为根据用CGC基因编码驱动程序p值和PID-N基因，得到基因之间的重叠，排名靠前的具有统计学意义。驱动程序FDR阈值0.1和0.25突出显示。黄色正方形表示PID-N基因和CGC基因之间的重叠。

图2d. 非编码驱动基因p值排名

在PCAWG驱动程序和功能解释工作组分析中，经过一系列分析比对，作者对途径和网络分析的共识程序显著增加了PCAWG驱动程序和功能解释工作组中的重要突变元素：PID-C和PID-N基因比PCAWG驱动程序和功能解释工作组分析中发现的基因多包含121个基因，包括90种可能的新非编码驱动程序。总体而言，与PID-C基因相比，PID-N基因中的非编码突变涵盖了151个样本。

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

非编码突变可以通过改变转录因子结合位点或其他类型的调节位点来起作用。因此，作者评估了PID-N基因中的非编码突变是否与同一基因中顺式表达的变化相关，并发现可以使用RNA-Seq数据测试的90个中的5个PID-N基因(FDR <0.3)显示出显著的顺式表达相关性(见图3) 。相反，在87个PID-C基因中，有34个在顺式表达相关性中表现出显著或接近显著性(FDR <0.3)。

由补充实验得知，PID-N基因在顺式表达中最重要的是TERT启动子突变与表达增加之间的相关性。作者发现其他四个PID-N基因：TP53，TLE4，TCF4和DUSP22在顺式表达中也很重要。

如图3，用PhyloP估计的基因组元素的进化保守性显示为灰色特征。在GM12878 HapMap B淋巴细胞细胞系中测得的H3组蛋白赖氨酸4个三甲基化位点在绿色轨迹中突出显示，表明转录起始位点附近的活性启动子区域。箱线图中心线显示中位数，箱线图边界显示第一个四分位数Q1和第三个四分位数Q3。

• 图a为TP53启动子，说明了具有TP53启动子区域放大图的TP53编码和非编码基因组位点。TP53启动子突变与显著降低的TP53基因表达有关， TP53启动子区域中拷贝数增加和减少的样品分别用红色和蓝色标注。六个TP53启动子突变中的两个与转录因子结合位点重叠(一个突变匹配三个基序)。
• 图b为TLE4启动子。说明了TLE4启动子区域放大图的TLE4编码和非编码基因组位点。肝脏-HCC样品中的TLE4启动子突变(三个突变)与较低的TLE4基因表达相关。拷贝数损益的样本分别用红色和蓝色标注。三种TLE4启动子突变之一具有ZNF263的转录因子结合位点。
• 图c为TCF4启动子。说明了具有TCF4启动子区域放大图的TCF4编码和非编码基因组位点。肺SCC样品中的TCF4启动子突变(三个突变)与较低的TCF4基因表达相关。拷贝数损益的样本分别用红色和蓝色标注。三个TCF4启动子突变之一具有ZEB1的转录因子结合位点。

图3. 基因表达的变化与PID-N基因的突变的相关性

最后，在五个患有DUSP22 3'UTR突变且显著过表达的肺腺癌患者中，DUSP22的表达显著降低；在3名患有DUSP22 5'UTR突变的肺腺癌患者中。这些UTR突变是互斥的。DUSP22编码磷酸酶信号蛋白，最近被提出是淋巴瘤中的肿瘤抑制因子。

此外，具有相关基因表达变化的PID-N基因数量少是由于PID-N基因突变的样品数量少，整个肿瘤类型中表达数据的不均匀可用性以及启动子区域的序列覆盖率降低所致 。这些问题进一步减少了具有非编码突变和RNA表达的样品数量，从而限制了顺式基因表达关联分析的功效。

4. 编码和非编码突变的模块化组织

作者确定了特定的PPI子网络和生物学途径，这些子网和途径被编码突变，非编码突变或两种突变的组合所改变。作者发现，使用节点级保留置换测试偶然发现的PID-C基因之间的相互作用显著更多，接近大量的PID-N基因之间的相互作用，并且PID-C和PID-N基因之间的相互作用显著更多，证明了物理PPI网络上编码和非编码突变之间的相互作用。

作者将涉及PID-C和PID-N基因的相互作用子网组织为五个生物过程：核心驱动程序，染色质组织，细胞增殖，发育和RNA剪接(见图4a)。

如图4a，图4是包含PID-C和PID-N基因的途径和网络模块。其中，图4a的节点代表PID-C和PID-N基因，边缘显示来自ReactomeFI网络的功能相互作用(灰色)，来自BioGRID网络的物理PPI(蓝色)或两个网络中记录的相互作用(紫色)。以及表示PID-C基因(绿色)，PID-N基因(黄色)或PID-C和PID-N基因(橙色)；节点大小与基因得分成正比；每个节点中的饼图表示与相应基因相关的编码和非编码突变的相对比例。虚线轮廓表示在染色质组织和细胞增殖中起作用的基因簇，主要包含PID-C基因。发育包括相当数量的PID-C和PID-N基因；RNA剪接包含PID-N基因，还显示了具有许多已知驱动因子的核心基因簇。

图4a. PID-C和PID-N基因之间的功能相互作用网络

作者使用g：Profiller Web服务器进一步表征了在PID-C和PID-N中富集的分子途径(图4b)。总体而言，富集了PID-C基因的63条途径，富集了PID-N基因的13条途径。作者的基因优先方法使用途径数据库和相互作用网络作为先验知识，因此PID-C和PID-N基因富含多种分子途径。

如图4b，矩阵中的每一行对应于PID-C或PID-N基因，矩阵中的每一列对应于富含PID-C和/或PID-N基因的途径模块。一个已填写的条目表示一个基因(行)属于一个或多个通道(列)，这些通道根据PID-C基因(绿色)，PID-N基因(黄色)或PID-C和PID-N的基因成员进行着色基因(橙色)。深色条目表示PID-C或PID-N基因分别属于显著富集PID-C或PID-N基因的途径，浅色条目表示PID-C或PID-N基因属于显著富集的路径，用于PID-C和PID-N基因的并集，但不适用于PID-C或PID-N基因。丰富度由每个途径模块名称和PID基因名称附近的圆圈概述。空心圆表示途径模块包含的途径比单独的PID-C和PID-N结果显著丰富了PID-C和PID-N基因的并集。丰富的模块和PID基因被分为四个生物学过程：染色质，发育，增殖和RNA剪接。

图4b. 包含PID-C和PID-N基因的途径模块

作者发现染色质和细胞增殖过程中的途径(包括染色质重塑，组蛋白修饰，凋亡信号转导，信号转导，Ras信号转导和细胞生长)主要是由PID-C基因中的编码突变所改变的，因为这些途径包含许多癌症基因，同时，作者还发现多个信号传导途径包括大量的PID-C和PID-N基因，这表明非编码突变为破坏这些途径提供了编码突变的替代方法。

作者还发现，PID-C和PID-N基因的显著数目改变了几个发育过程。细胞命运确定主要受PID-N基因DUSP6，MEF2C，JAG1，SOX2，HES1，ACL1，ID2，SUFU和KLF4中非编码突变的影响，还可以通过编码PID-C基因BRAF，GATA3和NOTCH1 / 2中的突变来实现。神经系统发育有关的途径丰富了在癌症中起已知作用的PID-N基因ASCL1，CTNNB1，ID2，SUFU和TERT，补充了PID-C基因NOTCH1，PTEN和RHOA也具有癌症作用。

模式指定过程也受到编码和非编码突变的影响，包括PID-N基因ASCL1，SUFU和RELN以及PID-C基因ATM和SMAD4。在这些情况下，非编码突变可补充破坏这些途径的编码突变，从而覆盖其他患者。

此外，作者发现RNA剪接途径主要受非编码突变的影响。但由于作者未发现非编码突变与剪接相关PID-N基因表达变化之间的顺式关联有任何显著意义，因此作者探索了这些基因中非编码突变与表达之间的反式作用潜力其它基因，作者发现，与剪接相关的PID-N基因中的非编码突变在很大程度上概括了最近从TCGA PanCanAtlas分析中发表的关联，该关联在多个剪接因子中的编码突变与47条途径的差异表达之间(见图5)。作者使用跨这些途径的差异表达模式的分层聚类，鉴定了RNA剪接因子中的三个突变簇(C1，C2和C3，见图5a,b)。使用t分布随机邻居嵌入(见图5a中的顶部注释栏)发现了一组高度重叠的集群，这表明该集群对于选择集群方法具有鲁棒性。如先前所报道，这些簇中的每个簇在剪接因子SF3B1，FUBP1和RBM10中均包含至少一个编码突变，以及与剪接相关的PID-N基因中的非编码突变，表明这两种类型的突变均导致相似的基因表达标志。一组(P1)的特征是免疫细胞信号，另一组(P2)主要反映细胞周期。几个PID-N剪接因子中非编码突变的基因表达特征与先前报道的剪接因子中编码突变的特征之间的相似性支持剪接相关PID-N基因在改变相似基因表达程序中的功能性作用。

如图5，图5说明RNA剪接因子主要以非编码突变为靶标，并改变与剪接因子编码突变类似的途径的表达。如图5a，矩阵的列指示参考文献中报道的与剪接相关的PID-N基因的非编码突变和剪接基因的编码突变，矩阵的行指示47种精选的基因集。红色热图条目表示相对于非突变样本，突变样本中的途径上调，蓝色热图条目表示下调。第一栏注解表明突变簇的成员符合共同途径的规定，第二栏注释指示突变是PID-N基因中的非编码突变还是编码突变，第三栏注释指定畸变类型(启动子，5'UTR，3'UTR，错义或截短)，第四列注释指示编码突变的癌症类型。突变分为三类：C1，C2和C3。途径分为两类：P1和P2，其中P1包含免疫特征基因集，P2包含细胞自主基因集。

图5a. 基因集富集分析(GSEA)的热图归一化富集得分(NES)

如图5b，具有非编码突变的样品的基因表达特征聚集在与剪接相关的PID-N基因中，具有编码突变的基因表达特征在先前发表的剪接因子中。每个点的形状表示突变簇的分配(C1，C2或C3)，颜色表示相应的基因是具有非编码突变的PID-N基因还是具有编码突变的剪接因子基因。

图5b. 突变元素的t-SNE图

除上述模块外，作者还发现转录因子在PID-C和PID-N基因之间都得到了很好的表达。总共有9个PID-C基因是转录因子，而19个PID-N基因是转录因子。这一观察表明非编码突变可能会影响转录调控网络。

本文作者提出了一个整合的途径和网络分析方法，将具有可能的非编码驱动程序突变的基因列表扩展到很少突变的长尾中，发现同时具有编码或非编码突变的基因在途径和网络中重叠。作者整合的途径和网络分析确定了87个具有编码变体(PID-C)的与途径相关的驱动基因和93个具有非编码变体(PID-N)的与途径相关的驱动基因。重要的是，通过非编码突变数据的单元素测试，90个PID-N基因在统计学上不显著，并且这些基因是未来实验表征的关键候选者。其中，作者还发现TP53，TLE4和TCF4中的启动子突变与这些基因的表达降低有关；编码和非编码驱动基因突变主要针对不同的基因，并且对受癌症困扰的途径和网络做出了不同的贡献。特别是，作者在RNA剪接途径中发现了多个基因中的非编码突变，具有这些突变的样品表现出的基因表达特征与在具有编码突变剪接因子SF3B1，FUBP1和RBM1037的样品中观察到的基因表达变化一致。总而言之，这些结果表明，罕见的非编码突变可能导致与普通和互补生物学过程相似的扰动。

我们一起期待下一篇精彩的文献吧！

编辑：桑葚

校审：螺蛳粉，糯米饭