文献解读:36个非癌组织体细胞突变

今天，小编为大家文献解读的是19年12月份发表在Genome Biology(IF≈14)杂志上的文章。体细胞突变在衰老和肿瘤发生中起着关键作用，为了系统地识别人体中的体细胞突变并了解其分布和功能影响，作者开发了一种方法，该方法可利用RNA携带的基因组信息来识别DNA体细胞突变，同时避免了大部分的假阳性。研究通过分析36个非癌组织中的体细胞突变，揭示了人体中组织特异性体细胞进化的基本模式。

 

The somatic mutation landscape of the human body

人体的体细胞突变谱

 

方法

原始RNA-seq检索和处理

原始的reads数据下载自dbGAP GTEx，过滤之后包含7584个样本，对应36个组织。同时，作者还处理了105个全血样本的DNA外显子组测序数据(有配对RNA-seq样本的)；然后，用STAR将reads比对到hg19参考基因组，为了避免种系突变的污染，dbSNP中的所有SNPs都被标注为Ns，并在下游处理中忽略；在比对之后，删除了重复的reads，以避免在文库准备中由PCR复制产生的偏差。

 

从RNA-seq数据检测DNA体细胞突变

为了过滤掉测序错误、RNA编辑事件、种系变异等造成的假阳性，作者将经典的DNA突变检测方法与其他过滤条件进行结合，自定义了一种新的从RNA-seq数据检测突变的方法。主要包括以下三个部分：

(1)识别有两个碱基响应的位置：在比对之后，分析bam文件以识别被拥有两个不同碱基响应的读段所覆盖的基因组位置。这一步，只考虑高覆盖度(c ≥ 40 reads)和高测序质量的位置(qs ≥ 30 in Phred scale)；最后，只保留最小等位基因得到至少6个reads支持的位置。除此之外，研究只考虑测序错误率p<0.0001的突变体。

(2)识别并消除种系突变：作者分别基于dbSNP和GATK’s HaplotypeCaller来消除常见及罕见的种系突变。

(3)过滤其他潜在的假突变：1、剔除了性染色体、未完成的染色体支架、线粒体基因组以及6号染色体上的HLA基因座等区域；2、RADAR和DARNED存储了编辑信息，作者去除了这两个数据库中描述的所有位置；3、排除了所有位于注释外显子末端7bp处以内的突变，以消除剪接造成的假阳性；4、过滤掉了有测序错误率大于0.01%的突变；5、为了消除读段位置、匹配质量、序列质量的偏向性以及链的偏差性，使用Mann-Whitney U检验过滤了p<0.05的突变(BCFtools软件)；6、作者删除了reads间替代等位基因位置的平均配对距离过高或过低的突变体(p<0.05)；7、剔除了等位基因突变频率(VAF) >0.7的突变体，以消除RNA特异性的等位频率偏差；8、为消除组织特异性突变的影响，作者删除了所有在某一组织至少40%的样本中被检测到的突变；9、此外，作者还剔除了在至少4%的总体样本中存在的突变；10、最后，研究删除了那些突变量比预期(基于测序深度和生物学因素)高出很多的样本。

 

方法验证

为了验证从RNA-seq识别DNA体细胞突变的准确性，作者将其与DNA外显子测序数据(下载自GTEx dbGaP)进行了比较，对105名随机选择的样本原始数据进行了相同的步骤，来识别突变。为了验证整个检测方法，研究将同一个体的RNA-seq和外显子DNA-seq样本进行了匹配，并分别计算了突变百分比的校正FDR值。

 

与生物和非生物因素有关的突变

为了识别和减弱非生物因素对每个样本突变数量的影响，作者根据组织将样本分类并利用非生物因素(包括测序深度、总缺血时间、RNA完整性数量等)作为解释变量对突变数做了线性回归，通过多元回归分析，根据每个系数的p值来评估它们与突变数的关联显著性。

为了发现与突变数相关的生物因子，研究取了先前线性回归的残差，并以生物因素(性别、年龄、BMI等)作为解释变量又进行了一次回归分析，同样通过多元回归计算了每个系数的显著性p值，并且分别计算了每个生物因素和通过线性回归得到的突变数与非生物因素残差之间的Spearman相关系数。

 

样本突变谱相似性分析(tSNE)

为了评估样本之间的突变谱相似性及其组织特异性，作者用tSNE对突变谱进行了聚类分析，然后基于轮廓系数打分。研究使用了包含1536种突变的图谱，由于样本特异的基因表达会造成背景序列的差异，作者使用了标准化的突变值(详细的标准化过程就不在此赘述了，感兴趣的小伙伴可以去原文查看)。

为了量化tSNE二维空间中不同分组之间的聚类，研究计算了每个样本的轮廓系数得分Si，然后计算了某一组中所有样本的平均得分s和95%的置信区间。最后，作者通过计算20个随机个体所有组织的轮廓系数平均得分，基于个体进行分组，并通过改变组织标签、重复计算得分得到了随机期望值。

 

血液和肺样本的细胞类型分解

研究使用CIBERSORT识别了GTEx数据中每个全血样本的细胞类型组成。使用的参数为默认参数和“absolute”模式。

 

基因表达关联性

为了避免人群对所有表达关联性分析的影响，研究只选择了白种人。对于每个组织，作者使用单个基因表达，通过线性回归来对突变值建模，为了评估SNP对表型的影响，计算了基因组中在至少20%的样本中TPM>1的每个基因的γ值；显着性p值是对γ进行非零t检验获得的，并进行Bonferroni校正。pi为突变负荷向量，将其定义为之前描述的非生物因素回归分析后的残差。在所有案例中，pi和γ都利用分位数进行了标准化。

 

基因本体分析

为了分析基因表达与突变负荷的关系，研究使用GOrilla进行了GO生物富集，并使用REVIGO去除冗余，输入的基因是根据它们显著表达的组织数进行排秩后的(p < 0.05)。

 

基因组不稳定基因的表达关联

研究选择了一组已知参与DNA修复或转录复制等通路的基因(图4)，并对每个组织的所有不同突变类型进行了表达与突变负荷之间的关联分析。作者使用了之前描述过的线性回归分析，并分别计算了单个基因–组织关联显著性、基因在所有组织的富集显著性以及单个组织中通路的关联显著性所对应的FDR值。

 

染色质分析

为了去除染色质对体细胞突变的影响，作者将GTEx和Roadmap Epigenomics Project中的组织进行了匹配，筛选出18个组织。对于每个组织，基于单个外显子，计算了H3K36me3, H3K4me1, H3K4me3, H3K27me3和H3K9me3(组蛋白修饰)的突变率和信号。然后以所有组蛋白修饰为特征，对突变率进行线性回归分析，计算了每个组蛋白修饰与突变率之间的关联显著性p值。此外，每个组蛋白修饰的个体效应大小和方差也使用线性回归进行了计算。

 

富集到COSMIC中的突变

作者从COSMIC下载了整套癌症突变数据，只保留了没有indel(插入或缺失)的单核苷酸变异。对每个样本，计算了它们的突变与COSMIC中突变的重叠率及其显著性(运用超几何分布)。

 

dN/dS分析

为了计算dN/dS(非同义替换率/同义替换率)比率，作者使用了泊松分布来对具有不同影响因素的突变数建模。对于非同义突变，用另外一个参数来表示dN/dS值，所用方法为R中的dndsloc。

 

癌症驱动基因的突变分析

基于TCGA中的数据，研究下载了已知包含至少一个癌症驱动突变的基因列表，并进一步移除了那些可能会使基因水平的分析产生偏差，但不影响全基因组分析的潜在假阳性突变响应。

为了评估这些基因的突变负荷及其显著性，对于每个组织，研究计算了这些基因的突变率，并分析了是否高于或低于这些基因在所有组织中的总体突变率(图5.b)。对于突变率高于预期的组织，使用二项分布的右侧积分来计算观察到n个突变的概率，相反，则使用左侧积分，并用Benjamini-Hochberg方法对p值进行FDR校正。此外，这些驱动基因对应突变的致癌状况用oncokb工具注释。

 

结果

547个个体的36个非疾病组织的体细胞突变

研究所使用的方法考虑了在RNA-seq读段中观察到两个等位的基因组位置，并评估了它们是否是真正的DNA体细胞突变(图1.a)。然后，根据多个过滤条件，消除了来源于生物和技术因素的假阳性(图1.b-c)。

为了验证该方法的有效性，作者将从105个血液样本的RNA-seq数据和外显子DNA测序数据中识别到的体细胞突变进行了比较。计算所得的FDR为29%，表明在RNA-seq中识别到的突变百分比得到相应外显子样本的支持(图1.d)。等位基因突变频率(VAF)较高的突变具有较低的FDR值，这使得整体突变相应的FDR值预估为34%，与之前的一项研究(40%)相当。

 图1. 从RNA-seq识别DNA体细胞突变的流程

 

最终，研究保留了来自36个不同组织、547个不同个体(未检测到癌症)的7584个样本。基于这些样本，共识别出280,843个不同的突变(样本中的平均频率为0.026%，个体中的平均频率为0.37%)，没有发现VAF接近0.5的突变富集，表明杂合子种系突变的影响不大。

在分析了影响每个样本和组织突变计数的因素后，研究发现影响最大的是测序深度，其中，突变数高于预期的多为接触环境诱变剂或细胞周转率高的组织(图2.a)。而且，作者观察到所有六种可能的点突变类型也都有类似的趋势，C>T和T>C频率最高，C>A次之。

 图2. 体细胞突变的跨组织分析

 

突变负荷受年龄、性别、种族和自然选择的影响

在控制了所有已知的技术因素后，研究人员发现年龄与大多数组织的突变负荷呈正相关，且血液中的突变与年龄相关性最显著，此外，皮肤中的C>T突变(与紫外辐射有关)是最显著的，除年龄以外，其他生物因素也会导致体细胞突变。例如，女性的乳腺突变负荷高于男性，白人种族中C>T突变与是否暴露于紫外线相关，而非裔人群则影响不大(图2.b-c)。

为了了解自然选择对体细胞突变的作用，研究计算了观察到某一突变的样本的VA F，发现同义突变的VAF显著高于错义和无义突变(图2.d)。

 

体细胞突变具有组织特异性且与染色质状态有关

为了可视化体细胞突变模式的组织特异性，作者将t-SNE应用于每个样本中的全部突变，并使用轮廓系数得分(SS)来量化样本在二维空间中的聚类，SS = 1表示所有样本聚在一个类中。结果表明，根据捐献者来源进行聚类时，平均SS = 0；根据组织进行聚类时，0.7 > SS > 0.1，表现出组织特异性；脾脏、血液、皮肤、肝脏和食道粘膜的得分最高，动脉、肺、胃和甲状腺则最低(图2.f)。作者假设同癌症中一样，染色质对这种组织特异性起着重要作用，发现在大多数组织（除了大脑）中，突变率与异染色质标记H3K9me3之间有很强的正相关关系，与H3K36me3显著负相关(图2.g)。这些结果表明，染色质与人体突变率的关联无处不在，并出现在癌症发展之前。

 

突变链的不对称性普遍存在且在个体间存在差异

癌症突变通常优先发生在一条DNA链上，由于该研究的突变来源于转录外显子，所以作者关注了转录突变链的不对称性。作者观察到C>A突变具有最明显的不对称性，它优先发生在除大脑以外的大多数组织中的转录链上，转录链上的C>A突变也显示出个体间的差异性(图3.a-b)，C>A突变在同一个体的不同组织间往往是不对称的(图3.c)。

链不对称可能是由RNA特异性编辑或转录错误引起的。因此，研究使用了配对的血液样本的外显子数据，观察DNA与RNA-seq中的不对称水平是否一致(图3.d)。此外，研究发现血液样本中C>T和T>C的不对称水平最高。然而，这些不对称性的方向性表明，样本的高度偏差可能是由A>I(即T>C)和C>U(即C>T)这两种RNA编辑造成的。值得关注的是，在每个血液样本中，静息NK细胞和CD8+ T细胞丰度最高，在应激条件下，这两种细胞中的两种RNA编辑水平也都有提高。这表明，虽然未在血液中观察到这两种突变类型FDR值的增加(图3.d)，但一些RNA编辑位点可能存在于体细胞突变中。

 图3. 体细胞突变的链不对称性

 

基因表达揭示了与体细胞突变负荷相关的通路

作者对那些表达值与全外显子的突变负荷相关的基因进行了组织水平的分析，寻找在非疾病组织中体细胞突变负荷增加的细胞因子。重点关注的是丰度最高的C>T突变，研究发现只有特定的1-3个组织关联性最强(图4.a)；在与多个组织的突变负荷正相关的基因中，作者发现它们富集到了核苷酸切除修复、细胞转运、细胞粘附、自噬和免疫反应等GO功能类(图4.b)，其中一些生物学过程已知与癌症或正常细胞的突变负荷有关。

 图4. 突变负荷与基因表达和通路有关

 

为了进一步探讨不同DNA修复通路对突变的影响，作者重点研究了修复通路成员的表达与突变负荷之间的关联，具体分析了参与双链断裂修复(D SB)、错配修复(MM R)、核苷酸切除修复(NER)、碱基切除修复(BER)、DNA脱氨酶APOBEC3A/B和反式聚合酶POLH的基因。最终，识别出了14个FDR< 20%、8个FDR< 10%和2个FDR< 5%的关联关系(图4.c)。接下来，分析了不同组织间，这些通路基因与突变负荷关联关系的一致性。NER、MMR和BER的结果是显著一致的，并且这些关联与预期的一样，都是基于每个基因已知的信息。

此外，研究分析了每个组织中整个通路或功能相关基因集的表达情况，以查找可能影响突变负荷的因素。在多个组织中，MMR, BER和NER都与突变负荷强相关，但它们的重叠度很低(图4.e)，因此，每个组织中的关联主要由一个或两个通路主导。

 

癌症驱动基因的突变富集和在健康组织中的正向选择

为了研究癌症相关突变在癌前的状态，作者计算了突变谱中所有COSMIC癌点突变的富集。许多样本是高度富集的，其中一些与COSMIC突变有超过20%的重叠，暴露于阳光的皮肤组织其重叠度最高，免照于阳光的皮肤和骨骼肌次之(图5.a)，七个脑区、主动脉和脾脏的重叠最低。

在53个已知的癌症驱动基因中，有31个识别到了突变。一些组织在这些基因上表现出显著更高或更低的突变率(图5.b)；MAP2K1和RRAS2在肌肉中高度突变，而IDH2和PPP2R1A在心脏和骨骼肌中都高度突变(图5.c)；

 图5. 癌症驱动基因的突变富集和正向选择

 

然后，观察到的癌症驱动基因突变根据其对蛋白质序列的影响进行了分类。结果发现，一些基因如IDH2主要是错义突变，而另一些基因如MAP2K1则是大量的的同义突变。

为了评估选择压力对这些基因的影响，研究使用dN/dS计算了非同义突变与同义突变的比率，同时考虑了整个基因组的序列组成和变突率。dN/dS接近1表示存在中性选择，dN/dS > 1表示正选择，dN/dS < 1则表示纯化选择。该研究的分析表明，大多数癌症驱动基因都是在正选择效应下进化的(图5.e)。

为进一步探讨选择对癌症相关突变的影响，作者进行了VAF分析。根据COSMIC中的癌症点突变，创建了两个体细胞突变列表：已经在癌症样本中观察到的，以及那些与已知癌症突变位置完全相同、但碱基变化不同的突变。然后，将这两个子集分别同义、错义和无义突变三种。在这三种情况中，癌症相关突变的VAFs明显高于相应的非癌症相关突变(图5.f)，表明癌症相关突变会增加细胞的增殖率。有趣的是，癌症/非癌症突变VAFs的比值在同义突变类型中最低，错义突变居中，在无义突变中最高，表明了对许多极端突变的显著正向选择。此外，与其他突变相比，癌症相关突变在更多的个体中存在，特别是在非脑组织中，这些结果与dN/dS分析一致(图5.e)，同样表明了对癌症驱动基因中无意义突变的强正向选择。

最后，研究筛选出了在OncoKB中被注释为致癌或可能致癌的驱动突变。结果发现，所有观察到的NOTCH1中的无义突变都被标注为致癌，RHOA和RAC1在它们的Ras结构域也有致癌突变(图5.g)。以上所有结果表明，癌症突变在非疾病组织中富集，它们不仅聚集在驱动基因中，而且大多是在正向选择下进化的，表明这些突变在癌症发生之前就会增加细胞的增殖。

总的来说，这篇文章描绘了一副复杂的人体体细胞突变谱，突出了它们的组织特异性分布和功能关联关系。癌症突变的普遍存在和癌驱动基因在非疾病组织中的正向选择表明了一种稳定的癌前状态的可能性。最后，该研究从RNA-seq数据中识别体细胞突变的方法对于研究体细胞进化及其在衰老和疾病中的作用具有重要意义。

 Best wishes for your holiday!