开发一种无偏的工具来评估体细胞变异体的等位基因表达 虽然有多种方法来描述某一特定突变是否被表达，但大部分工作都涉及基因特异性分析，而不是对癌症基因组的整体评估。为了全面评估肿瘤内突变的表达，我们量化了来自pdac患者来源的细胞系、患者来源的异种移植(Pdx)和从癌症基因组图谱(Tcga)中获得的数千个突变(误解、无意义和小插入/缺失)的等位基因表达。 24 ,25 ..由于突变可能以体细胞单核苷酸变异(SNV)或小插入/缺失(INDEL)的形式存在，因此发展了一种独特的方法来确定所有突变的等位基因表达。而不是依赖于肿瘤RNA-序列读取正确地排列到一个标准化的参考基因组(例如Hg 19)，Maxx(Maxx：突变等位基因E)。 x 压E x 建立了一种新的肿瘤选择性参考基因组，并专门提取突变等位基因的表达。Maxx方法利用来自dna测序的突变调用来生成肿瘤特异性参考基因组(如图所示)。 1A )。这个精确的参考基因组随后与tophat 2一起使用。 26 为了提高肿瘤RNA序列对突变等位基因的比对。Maxx方法可以计算DNA测序法定义的任何类型的RNA突变等位基因频率，从而能够对突变等位基因的表达进行全面、无偏的分析。
图1 开发一种无偏见的工具来评估体细胞变体的表达。( a )Maxx管道的流程图，最终确定突变等位基因和野生型等位基因的RNA读取计数。( b 突变表达群定位方法：V-ex(蓝色)、W-ex(绿色)和N-ex(红色)：V-ex，表达突变等位基因的突变；W-ex，仅表达野生型等位基因的突变；N-ex，不表达野生型等位基因或突变等位基因的突变。( c )与每个突变表达组相关的突变数目。这些突变来自19个PDAC患者来源的细胞系。( d )所有患者来源细胞株的每个突变表达组的比例。
肿瘤模型中不同等位基因表达模式的差异 为了初步了解Maxx管道的有效性，我们为19个PDAC患者衍生细胞系生成了突变表达谱，这是一个包含所有已识别突变的DNA和RNA等位基因频率的文件。这些细胞系相对于正常组织进行了RNA测序和外显子测序.细胞模型是评估Maxx的理想方法，因为它没有基质污染，大大降低了肿瘤纯度的混杂特征。 27 ,28 ..对19种突变表达谱的评估表明，肿瘤突变可以根据与野生型等位基因、突变等位基因或这两种等位基因相一致的读取数来分类(图1)。 1B )。为了防止突变基因与不表达突变等位基因的突变之间的错误分类，研究中放弃了总RNA读计数在3到9之间的突变。所有的非丢弃突变被分为以下三个突变表达组之一：1.表达变异等位基因的突变，标记为变异表达群(V-ex，Blue)2，其中只有野生型等位基因被表达，标记为野生型表达群(W-ex，Green)。3.不表达野生型或突变等位基因的突变，标记为非表达群(N-ex，Red)。这三个突变群在多发性骨髓瘤患者衍生细胞系中也得到了鉴定。 29 ..三个突变组的分布总结了所有的患者衍生细胞系(图一)。 1C )。N-ex组占所有变异的一半，而只有一小部分突变属于W-ex组。在所有单个模型中都观察到了这种总体分布，表明这种表达分布并不是特定于特定的测序运行或肿瘤(如图所示)。1D ).
Maxx管道精确映射indels，计算效率高。 先前的研究表明，rna序列比对器能够对snv进行校准，但是对于含有indel突变的读码，其校准精度很低。 30 ,31 ..然而，Maxx独特的工作流程允许对SNV和indel突变都有信心地计算等位基因表达(如图所示)。 2A )。增强rna序列数据的indel对齐，这是准确识别rna等位基因频率的一个关键方面。 20 ，Maxx根据来自dna测序的突变调用生成肿瘤特异性参考基因组。当比较bam文件中的rna等位基因频率时，bam文件使用Maxx或hg 19参考基因组与tophat 2对齐，并经过了重要的筛选步骤进行等位基因表达分析。 21 ,32 与Hg 19参考基因组相比，Maxx参考基因组显着地增强了indel突变的等位基因表达。 2B )。使用Maxx或Hg 19参考基因组，SNVRNA等位基因频率变化较小。
图2 Maxx管道精确映射indels，计算效率高。( a )整合的基因组观察可视化的原始RNA-seq读取与野生型等位基因或突变等位基因的SNV，插入和删除突变。与Hg 19参考基因组相比，非灰色的颜色代表了另一种核苷酸( b )比较利用Maxx产生的参考基因组或Hg 19参考基因组计算的RNA突变等位基因频率。( c )利用Maxx产生的参考基因组或Hg 19参考基因组识别的RNA突变等位基因频率与Maxx所创建的附加突变基因组之间的对比。
为每个肿瘤样本生成一个参考基因组将需要大量的计算处理和存储来进行大规模的研究。然而，Maxx产生的参考基因组明显小于Hg 19参考基因组。这种缩小的大小是由于肿瘤参考基因组只包含野生型和突变序列的每一个基因携带一个突变。与Hg 19参考基因组相比，使用只包含突变基因的参考基因组显着地减少了用于索引参考基因组和运行Tophat 2的存储空间和计算资源(补充图)。 1 )。为了确定是否由于高度截短的基因组而过度选择对读序列，我们比较了从Hg 19参考基因组中建立的RNA突变等位基因频率与附加的Maxx突变体基因组与Maxx生成的参考基因组之间的差异(图1)。 2C )。这两个参考基因组提供了几乎相同的RNA突变等位基因频率和突变的位置到表达群(补充图)。 2 )。总之，这表明一个简化的参考基因组为文本读取的差异比对提供了一个健壮的底物，并且可以自信地用于生成突变表达谱。
突变等位基因表达与dna突变等位基因频率有关，而与rna表达水平无关。 为了解决PDAC中RNA突变等位基因频率是否与DNA突变等位基因频率有关的问题，我们评估了三个突变表达组之间的双重关系(图一)。 3A )一般说来，V-ex突变的变异等位基因表达与其dna等位基因频率(R)相关。 2 =0.59)。在另一项使用ct 26、b16f10和4t1小鼠细胞系的研究中，也发现了这种相关性。 33 ..然而，在我们的数据集中观察到，Vex突变的多个基因有100%的RNA突变等位基因频率，尽管DNA突变等位基因频率约为50%。dna和rna突变等位基因频率的差异表明，这些基因的非突变等位基因被选择性地沉默，这是抑癌基因的常见现象。 34 ..N-ex和W-ex突变表达组未表达突变等位基因，其RNA突变等位基因频率为0(如图1所示)。 3B )。为了确定突变表达组是否偏向于特定的突变类型(误解、删除、插入、无义)，如无意义介导的衰变过程中所期望的那样，在突变表达组之间测量突变类型的比例(补充图)。 3 )。统计检验表明，突变类型与突变表达群之间没有相关性。有趣的是，突变表达组的变异等位基因频率确定了突变表达提供了突变选择性的洞察力。 3C )。Vex组的平均变异等位基因频率明显高于其他两组，表明Vex突变对肿瘤增殖的影响最大。对于N-ex组，大多数突变的变异等位基因频率在30%到40%之间，这表明这些突变中有许多虽然没有表达，但实际上在肿瘤中基本上是克隆的。相反，W-ex组的平均变异等位基因频率明显低于其他两个突变表达组。为了询问是否存在导致转录缺失的突变的选择性特征，我们将没有突变等位基因的样本的平均基因表达与含有突变等位基因的样本的平均基因表达进行了比较(图1)。 三维空间 )。如预期的那样，含有N-ex突变的基因的总体表达量大大低于Vex或W-ex突变的基因。这证实了Nex突变的基因很少被表达，因此这类突变对PDAC增殖的影响很小。相反，具有V-ex或W-ex突变的基因在肿瘤中的基因表达相对较高，表明Vex或W-ex突变基因对肿瘤的进展具有重要意义。各表达组间突变基因与未突变基因差异显著，但幅度变化不大。有趣的是，我们从基因型组织表达(GTEx)数据库中获得的193个正常胰腺组织中，在每个突变表达组中发现了相似的基因表达模式。 3E 进一步提示，与Vex和W-ex突变基因相比，N-ex突变基因与胰腺组织的相关性不大。然而，在组织特异性方面，没有观察到突变表达组含有很高比例的胰腺组织特异性基因(附图)。 4 )。为了确定基因表达是否增加，但没有突变等位基因表达的W-ex突变是由于RNA剪接事件，对突变和非突变外显子的外显子表达进行了量化(补充图)。 5 )。与转录本分析一样，突变型和非突变型外显子的表达差异不大。这些发现表明，在任何一个突变表达组中，没有针对转录或外显子表达的强选择。
图3 突变等位基因表达与DNA突变等位基因频率有关，而与RNA表达水平无关。( a )患者衍生细胞系DNA等位基因频率与其相应的RNA等位基因频率相对应。( b )V-ex突变、W-ex突变和N-ex突变的外显子测序和RNA测序的整合基因组观察。非灰色代表与Hg 19参考基因组相一致的相互冲突的核苷酸的存在。( c )基于双尾p值的双样本t检验，分析了3个突变表达组的dna等位基因频率分布及其统计学意义。( d )两个样本配对t检验，每个突变表达组的变异基因平均表达水平与无突变基因样本的平均基因表达水平之间存在一个双尾p值。( e )在193个正常胰腺组织中，病人衍生细胞系中已鉴定的突变基因的平均表达的盒形图。采用双尾p值的双样本t检验对突变表达组进行比较.
PDAC相关基因主要属于V-ex和W-ex两类。 为了更全面地理解突变表达谱的意义，我们评估了三个表达组中的突变基因。在PDAC肿瘤方面，众所周知有四个基因经常发生突变并被认为是疾病驱动因素：kras、tp 53、smad 4和ckn2A。 35 ..这四个基因被反复观察到含有V-ex突变，扭曲了V-ex组的突变频率分布，如小提琴图所示(图1)。 4A )。相比之下，突变基因在W-ex和N-ex群中的频率分布主要集中在一个群体，这意味着这些基因在PDAC中的频率分布很少。像kras这样的癌基因的dna突变等位基因频率一般为50%，因为只有一个等位基因被突变才能作为癌基因。 36 ..然而，在选择的情况下，只有KRAS突变等位基因的表达有特定的选择。有趣的是，在这些样本中，KRAS突变等位基因的100%表达主要是由于DNA突变等位基因频率的增加，而不是等位基因特异性基因的表达(图一)。 4B )。相比之下，抑癌基因TP 53、Smad 4和CDKN2A的突变等位基因频率接近100%，主要表达突变等位基因。就Smad 4和CDKN2A而言，这是意料之中的。然而，由于TP 53突变可以获得功能突变，这一数据表明，在胰腺内肿瘤发展过程中失去野生型等位基因的压力非常大。
图4 PDAC相关基因主要属于V-ex和W-ex两类.( a )所有V-ex、W-ex和N-ex突变的基因突变频率。( b 比较最著名的PDAC突变(KRAS、CDKN2A、Smad 4、TP 53)的DNA等位基因频率和RNA等位基因频率。( c )每个突变表达组与宇宙和坦博莱罗癌症基因数据集相关的部分突变。每个突变表达组之间采用两比例z检验进行统计学意义。( d )V-ex、W-ex和N-ex突变基因所产生的网络，这些变异基因存在于宇宙数据集或Tamborero数据集中。( e 前50位突变基因的分布情况由驱动突变检测方法2020+、Muffin和OncodriveFM输出。( f )对三种驱动突变检测方法进行了Fisher‘s精确测试，根据宇宙和坦博莱罗的50个顶级突变基因，对其识别癌症相关基因的能力进行了量化。计算了所有前50位突变基因或不含N-ex突变的前50位突变基因的癌症基因统计意义。
为了阐明对癌症的推测意义，我们将重点放在COSMIC.v80中的每一个突变表达组中突变的比例上。 10 或者坦博罗 等人 .37 癌症相关基因数据集(图1。 4C )。正如预期的那样，由于KRAS、TP 53、Smad 4和CDKN2A的频率，V-ex组中与癌症相关的突变基因所占比例最大。将W-EX组与N-ex组进行比较，发现W-EX组与癌症相关的突变基因比例明显高于Nex组。为了确定突变的癌症相关基因与突变表达组之间是否有关系，进行了网络分析(图1)。 4D )。对产生的网络的评价表明，具有W-ex突变的癌症相关基因与含有V-ex突变的癌症相关基因整合得很好。至于含有N-ex突变的基因，通常在网络外部发现，对网络结构几乎没有贡献。这些结果支持了这样的假设，即尽管W-ex基因缺乏克隆选择和突变等位基因的表达，但它们发生在对肿瘤很重要的基因中。相反，N-ex突变似乎存在于对肿瘤进展几乎没有影响的基因中。
大多数，如果不是全部，目前的司机突变检测方法没有整合等位基因表达信息时，预测有影响的癌症突变。因此，为了确定突变等位基因表达的掺入是否能提高肿瘤基因的特异性，我们首先采用了三种常用的驱动突变检测方法。选择的三种方法为2020+。 19 穆芬 17 ，以及OncodriveFM 15 ..之所以选择这些方法，是因为它们具有处理indel突变的能力和相对较小的样本大小。当比较2020年+、Muffin和OncodriveFM前50位突变基因的突变表达群分布时，V-ex和N-ex组之间存在一致的权衡(图一)。 4E )。驱动突变检测方法2020+的VEX突变基因比例最高，OncodriveFM突变基因的比例最低。为了测试每种方法从CosmicandTamborero数据集中识别癌症相关基因的能力，进行了Fisher‘s精确测试，假设每种方法预测的前50位突变基因都与癌症相关(图一)。 4F )。总的来说，Muffin的表现超过了2020+和OncodriveFM。当我们从前50位被认为是驱动基因的突变中去除N-ex突变时，这三种方法在预测癌症相关突变基因方面都有统计学意义。这表明，结合突变的等位基因表达有助于降低驾驶员突变检测方法的假阳性发现率，改善下游分析。
PDX肿瘤细胞突变表达特征的保存 为了确定Maxx管道在肿瘤突变表达谱中的有效性，采用PDX模型对Maxx管道进行了评价。19株患者衍生细胞系中有8株与PDX模型相匹配。 24 因此，这8种PDX模型与其相应的细胞系具有相同的肿瘤特异性参考基因组。然而，这些模型有自己的转录体序列。与细胞系数据相似的是，根据野生型和突变型等位基因读取计数来确定三个突变表达组，并根据DNA突变等位基因频率绘制RNA突变等位基因频率图(图1)。 5A )。PDX模型表明，突变表达群在两种基因中都存在。 体外 和 体内 ..RNA突变等位基因频率与细胞株相同，与DNA突变等位基因频率密切相关。在V-ex突变中，PDX与细胞系之间的表达特征有着非常高的共性(图一)。 5B )。这种保守性对于N-ex突变来说更为极端。 5C )。然而，W-ex突变在PDX与其相应的细胞系之间表现出更大的变异性，其中大多数独特的W-ex突变都在PDX样本中(图1)。 5D )。在评估突变等位基因频率与PDX和细胞系之间的共同和独特的联系时，V-ex变体显示了明显的关系(图1)。 5E )。虽然没有统计学意义，但W-ex突变显示了相反的关系.这一结果在一定程度上表明，突变表达缺乏保守性是特定变异体亚克隆特征的反映。
图5 从细胞系到PDX肿瘤突变表达特征的保存。( a )比较PDX各突变体的DNA等位基因频率和RNA等位基因频率。( b –d V-ex、W-ex和N-ex突变在相应细胞株和PDX样品之间的重叠.细胞株或PDX数据中被归类为“丢弃”的突变被排除在外。( e )对细胞株和PDX样品中存在的突变的DNA等位基因频率和对细胞株或PDX样本唯一突变的DNA等位基因频率进行了双尾Wilcoxon Mann Whitney检验。
原发肿瘤突变表达谱 原发性肿瘤常含有间质污染，这降低了突变呼叫的敏感性，减少了与肿瘤相关的rna序列读取。 38 ..因此，据报道基质污染在进行dna甲基化和亚型分析时会导致不准确的结果。 28 ,39 ..因此，为了限制肿瘤纯度的问题，我们重点研究了75例tcga pdac病例，这些病例的肿瘤纯度评分最高，由软件估计。 40 Vex突变率至少为10%.与细胞株和PDX模型相似，TCGA样本中存在三个突变表达组(图1)。 6A )。然而，肿瘤dna突变等位基因频率的分布大多在0-50%之间，而不是细胞株的0-100%(图5)。 6A v.V. 3A )。为了确定基质污染是否也改变了RNA等位基因频率，我们将重点放在KRAS、CDKN2A、TP 53和Smad 4突变上。 6B )。KRAS和CDKN2A的RNA等位基因频率范围可与细胞系数据相比较，TP 53和Smad 4的RNA等位基因频率范围有显著变化。提示Smad 4和TP 53在肿瘤标本和间质室均有表达。尽管如此复杂，这些典型的PDAC突变中的大多数突变被确定为在肿瘤中表达，如细胞株和PDX数据所示。
图6 原发肿瘤突变表达谱。( a )比较TCGA各突变体的DNA等位基因频率和RNA等位基因频率。( b TCGA中CDKN2A、KRAS、Smad 4和TP 53突变的DNA等位基因频率与RNA等位基因频率的相关性。( c )不包含突变基因的样本的平均第一外显子甲基化与含有突变基因的样本的平均第一外显子甲基化相比较。这是对具有第一外显子甲基化数据的非丢弃突变基因进行的。( d )对每个突变表达组之间的第一外显子甲基化进行双尾p值的双样本t检验。
利用来自TCGA样本的突变表达和甲基化特征，我们能够询问不同突变表达组的DNA甲基化水平。为测量含有突变基因的甲基化程度，计算了每个记录到的第一外显子内CpG的平均甲基化百分比。 41 ..与突变基因表达分析类似，没有突变的样本的第一外显子甲基化的平均值与含有突变基因的样本的平均第一外显子甲基化程度进行了比较(图一)。 6C )。决定系数的评价(R) 2 )对于每个表达组，野生型与突变基因第一外显子甲基化之间的差异不大。然而，当比较突变表达组间的平均第一外显子甲基化时，N-ex突变的平均第一外显子甲基化程度明显高于V-ex和W-ex突变。 6d )。这一发现表明，基因的过甲基化可能是N-ex突变基因缺乏表达的原因之一。当比较V-ex组和W-ex组第一外显子甲基化时，发现含有W-ex突变的基因第一外显子甲基化程度较高。为什么W-ex组的第一外显子甲基化程度低于N-ex组，而第一外显子甲基化程度却高于V-ex组，一个可能的解释是，这些基因易受等位基因特异性甲基化的影响。 42 ..带有W-ex突变的基因的等位基因特异性甲基化也解释了为什么这些基因的突变等位基因转录本是无法检测到的。
通过表达分离突变显著增强PDAC亚型 为了确定突变表达组是否提供预后相关性的洞察力，基于网络的分层(NBS) 43 对19株患者来源的细胞系和75例TCGA标本均采用管道。NBS使用蛋白质-蛋白质相互作用网络将样本划分为基于突变样本矩阵的预定义的子类型数。为了量化突变表达组的预后能力，只对TCGA样本进行生存日志秩统计-log 10(p-值)的计算。用3-8个预定义的亚型数计算所有可能的突变表达组的对数秩统计p值(图3~8)。 7A )。突变表达组对预后的统计预测能力最强的组合是V-ex和W-ex突变。无论预定义子类型的数量如何，这都是正确的。这些发现进一步证实了V-ex和W-ex突变与肿瘤生物学有关，并且发生在蛋白质-蛋白质相互作用网络的相似邻域内。相比之下，N-ex突变对肿瘤的进展影响不大，似乎是随机分布在蛋白-蛋白质相互作用网络中的。为了比较V-ex和W-ex突变的预后预测和流行的司机突变检测方法的基因排名，我们对2020+、Muffin和OncodriverFM的结果进行了类似的NBS分析。 7b )。为了在所有方法之间进行不偏不倚的比较，每种方法都使用相同的输入数据来产生一个突变基因的排序列表。在NBS分析中，使用相同数量的排列突变基因生成突变样本矩阵(n=1，349)。与每一种驱动突变检测方法相比，V-ex和W-ex突变的组合为每种预定义的亚型提供了最佳的预后预测。因此，这些结果意味着类似于图。 4F 目前的驾驶员突变检测方法报告了大量的假阳性，污染了检测结果，影响了下游分析。然而，通过根据它们的等位基因表达来分离突变，有可能删除对肿瘤进展几乎没有影响的大量突变子集，并以突变数据为重点显著提高预后预测。
图7 通过表达分离突变增强了PDAC亚型。( a )基于NBS结果对75个TCGA样本的对数秩统计量-log 10(p-值)。用对数秩统计量log 10(p-值)计算了3~8个预定义亚型突变表达群的所有配对组合的NBS结果。( b 对75个TCGA样本采用基于2020+、Muffin、OncodriveFM或V-ex和W-ex突变的分层样本进行对数秩-log 10(p-值)统计。( c )基于NBS结果的TCGA样本Kaplan Meier图，使用V-ex和W-ex突变，亚型数等于3。( d )网络，表示对亚型1和亚型2所特有的变异通路。( e )细胞株药物反应数据，通过Wilcoxon Mann Whitney检验，按已识别的亚型和统计意义分组。
以V-ex和W-ex组为基础的分层亚型中最具统计学意义的有3种。由这一结果生成的Kaplan-Meier图如图所示。 7C ..为了确定哪些基因是新建立的亚型的唯一基因，我们对这三种亚型的每一种亚型之间的nbs基因平滑评分进行了一个大于检验的Wilcoxon Mann Whitney( e .g ..亚型1与亚型2和亚型3)。在调整p-值后，使用benamini/yekutieli方法 44 ，每个亚型的网络都是使用含有调整后的p值小于0.01的基因生成的(如图所示)。 7D )。亚型1和亚型2分别以TP 53和Smad 4为中心，具有唯一定义的网络。GO定期分析发现，1亚型突变影响细胞凋亡途径，2亚型突变参与细胞分化途径。相比之下，被归类为第3亚型的样本只有少数几个具有统计意义的变异相互作用，因此无法生成结构良好的网络。有趣的是，TP 53或Smad 4的唯一突变并不足以构成本研究所用样本的亚型(补充图)。 6 )。这些数据支持使用路径和Maxx分析来根据突变数据破译新的预后亚型。
为了评估这种亚型是否与PDAC的治疗相关，使用了我们的19种患者衍生细胞系中的13种药物反应数据进行了统计分析(图1)。 7E )。在有药物敏感数据的细胞株中，5株为1亚型，4株为2亚型，4株为3亚型。为确定药物治疗在各亚型之间的意义，对3种亚型的曲线(AUC)值下的药物反应区(AUC)进行了Wilcoxon Mann Whitney试验。在亚型1和亚型3之间最一致的药物分别是极光激酶和EGFR抑制剂。与其他亚型相比，3亚型对mTOR抑制剂反应更好。出乎意料的是，亚型2对任何药物治疗都没有表现出选择性的敏感性。这些结果表明，Maxx衍生数据除了可用于预测预后外，还可用于药物敏感性的预测。
