胰腺癌非编码基因组中保守增强子突变的功能与临床意义研究

《Cancer Research Communications》：Mutational frequency in PDAC enhancers. A, Schematic of coding, promoter,... Open Access

本研究通过全基因组测序（WGS）系统分析胰腺导管腺癌（PDAC）编码区与非编码区（增强子、启动子等）的突变特征，发现增强子区域具有独特的突变谱（如SBS39富集），且保守增强子内的复发性体细胞突变可调控KLF5、TP63等胰腺发育相关转录因子表达，导致基因表达显著改变。该研究为PDAC非编码基因组研究提供了新视角。

　　

Abstract

胰腺癌的编码基因组已被详细表征，但非编码基因组的特征仍待探索。本研究利用全基因组测序技术，在两个独特的患者队列中研究了编码和非编码基因组（启动子、增强子和非编码非增强子区域）。研究发现，在所有四个基因组区域中，经过治疗的癌症比未经治疗的癌症具有显著更高的突变负荷。然而，尽管存在治疗诱导的遗传瓶颈，每个区域的突变相对比例得以保持。与其他非编码区域相比，增强子的突变数/Mb较低。增强子还具有独特的突变特征，表现为SBS39的富集。根据预测的人类和小鼠基因组同源区域重叠情况，将增强子序列分为保守和非保守区域，并对保守区域进行复发性体细胞突变筛查。研究发现，保守增强子区域中的复发性体细胞突变主要与已知在胰腺发育和癌症中起作用的转录因子相关，包括KLF5和TP63。基于RNA测序数据的转录表达分析显示，与没有增强子突变的癌症相比，具有增强子突变的癌症表现出显著不同的表达水平，最常见的是表达缺失，这表明了其功能性影响。这些发现扩展了我们对胰腺癌非编码基因组的认识，并指出了保守增强子突变在胰腺癌中未被探索的作用。

Introduction

尽管有关胰腺导管腺癌（PDAC）遗传学的数据丰富，但这种肿瘤类型的完整遗传特征谱仍有待发现。迄今为止，大规模全外显子组测序研究揭示了PDAC的复发性基因组特征，这些特征针对特定数量的核心通路。PDAC中的高频驱动基因突变包括KRAS、TP53、CDKN2A和SMAD4，它们在胰腺癌发生过程中的克隆扩增中起主要作用，并影响转移倾向。随后的亚克隆遗传事件对于肿瘤进展和转移也至关重要。新兴数据表明非编码基因组在癌症中也发挥作用。例如，长散在核元件1（LINE-1）活性增加与多种肿瘤类型相关，可能有助于已知驱动基因的体细胞结构变异。泛癌全基因组分析联盟对来自27种肿瘤类型的2,583个全基因组进行的全面重新分析也发现了非编码基因组的几种新的复发性体细胞改变。在PDAC中，Feigin及其同事发现启动子区域突变的影响，特别是与转录调控相关基因的顺式调控区域。尽管有越来越多的证据表明增强子区域的突变可能通过转录因子动力学在转录变化中发挥重要作用，但专门关注增强子突变的研究较少。在某些肿瘤类型中，如膀胱癌和乳腺癌，增强子突变被认为是真正的驱动突变，泛癌分析也支持这一观察结果。总体而言，虽然非编码基因组似乎在癌症中发挥作用，但这些改变比影响编码基因组的改变要少见得多。此外，组织特异性在识别复发性非编码体细胞改变中起作用；因此，泛癌研究中非编码改变的报道可能低估了特定肿瘤类型的改变。对于实体瘤克隆进展过程中非编码改变积累的时间点了解更少。为了解决PDAC这方面的知识空白，我们进行了多区域全基因组测序和克隆分析，以阐明非编码区域相对于编码基因组的特征性特征，包括它们与其他已知PDAC特征的进化关系。

Materials and Methods

Ethics statement

所有使用其组织的患者均获得了书面知情同意。该研究根据公认的伦理准则进行，并获得了约翰霍普金斯大学医学院和纪念斯隆凯特琳癌症中心机构审查委员会的批准。纪念斯隆凯特琳癌症中心进行的所有动物研究均经机构动物护理和使用委员会批准。

Tissue samples

使用了来自纪念斯隆凯特琳癌症中心“Last Wish”项目的六个PDAC研究尸检样本。从冷冻切片上切下组织切片并进行审查，以确定那些至少具有20%肿瘤细胞性且组织质量保存完好的样本。对正常样本进行审查以确认没有污染癌细胞。符合这些标准的样本在提取基因组DNA之前从连续未染色的切片中进行宏观解剖以富集肿瘤纯度。本研究还包括了先前通过约翰霍普金斯大学快速尸检项目分析的四个未经治疗的PDAC病例。

WGS

患者PAM01-04的WGS原始数据是先前生成的。从患者MPAM01-06的所有组织中提取基因组DNA。DNA定量、文库制备和WGS在MSK综合基因组学操作中心进行，体细胞变异的生物信息学分析由MSK生物信息学核心完成。MPAM01-06的测序、比对和分析按照PAM01-04的描述进行。简而言之，使用Illumina HiSeq 2000、HiSeq 2500、HiSeq 4000或NovaSeq 6000平台，目标覆盖度对于肿瘤样本为60X或80X，对于正常样本为60X或30X。对产生的测序读数进行计算机内分析，以评估质量、覆盖度以及使用Burrows-Wheeler Aligner与人类参考基因组的比对情况。在使用Picard套件和GATK版本3.1完成读数去重复、碱基质量重新校准和多重序列重新比对后，使用MuTect版本1.1.6和HaplotypeCaller版本2.4检测体细胞单核苷酸变异/插入-缺失。从系统发育分析中排除了低质量或比对不良的读数。在整个队列中获得了中位数覆盖度为76.5×（范围59-124）的肿瘤和中位数覆盖度为47.5×（范围28-76）的正常样本。

Filtering and annotation of variants

对所有患者，使用以下标准过滤体细胞变异：患者匹配的正常样本覆盖度≥10个读数，患者匹配的正常样本中的变异计数≤2，患者匹配的正常样本变异频率<0.02，肿瘤突变等位基因读数≥10，且至少在一個肿瘤样本中肿瘤变异等位基因频率≥0.05，条件是所有样本中肿瘤覆盖度≥20个读数。编码变异经过进一步的生物信息学注释，以确定致病性和来自全球健康人群的种系等位基因频率。通过FACETs推断拷贝数改变和全基因组复制。

RNA sequencing

RNA提取和测序按最近描述的方法进行。简而言之，使用TRIzol提取总RNA，然后使用RNeasy Plus Mini Kit。使用TruSeq Stranded Total RNA LT Kit进行文库制备后，对样本进行条形码编码，并在HiSeq 4000上以100 bp/100 bp或125/125 bp双末端运行。输出数据使用rnaStar比对器映射到目标基因组，并使用PICARD工具对输出的SAM文件进行后处理以添加读数组并将其转换为压缩的BAM格式。使用HTSeq确定来自比对读数的表达计数矩阵，并使用R/Bioconductor包DESeq2处理HTSeq生成的原始计数矩阵以标准化整个数据集。Log₂转换后的数据用作下游分析的标准化表达。

The Cancer Genome Atlas data set

从数据协调中心数据发布/国际癌症基因组联盟数据门户获取癌症基因组图谱胰腺癌WGS数据。从MuTect获得突变调用，并鉴定了总共2,271,144个独特突变。使用vcf2maf工具重新注释这些突变，并将新注释的文件作为下游分析的输入。通过FireBrowse下载TCGA胰腺癌RNA测序数据。从下载的RNA-seq数据计算每百万转录本数。TPM用于基因集富集分析，log₂转换后的TPM值用作相对mRNA表达。

Pancreatic cancer enhancer regions and enhancer mutations

增强子区域先前基于H3K27ac染色质免疫沉淀测序数据定义。从ChIP atlas中选择了23个胰腺癌H3K27ac ChIP数据集，基于峰值数量并结合集成基因组查看器的手动审查。使用MACS2合并和分析所有数据，并识别得分大于50的区域。使用GREAT将这些区域注释为最近和第二近的基因。包括位于转录起始点2至50 kb之间的区域，而排除距离转录起始点超过50 kb或小于2 kb的区域。总共有62,015个区域被纳入本研究作为“胰腺癌增强子区域”。使用Galaxy bedtools ClosestBed，将体细胞变异与PDAC增强子区域重叠的变异定义为增强子突变。根据先前的研究，使用微滴数字PCR验证了代表性增强子突变。

Mutational signature analysis

使用Palimpsest进行突变特征分析。每个样本在编码、增强子、启动子或非编码非增强子区域中识别的所有突变用作输入。每个样本每个类别中的特征比例用作聚类分析的输入，以识别增强子区域突变特征的特征。执行Wilcoxon符号秩检验以识别组间的特征性特征。使用Benjamini、Krieger和Yekutieli的FDR方法进行突变特征t检验，以识别克隆与亚克隆以及治疗和未治疗组之间的特征性特征。

Evolutionary analysis

通过使用Treeomics 1.7.9推导了每个样本集的系统发育。每个系统发育以匹配患者的正常样本为根，叶子代表肿瘤样本。Treeomics使用贝叶斯推理模型来解释容易出错的测序和变化的肿瘤细胞含量，以计算特定变异存在或不存在的概率。全局最优树基于混合整数线性规划。所有进化分析均基于WGS进行。存在于PDAC所有分析样本中的体细胞改变被认为是克隆性的，而存在于样本子集或单个样本中的则被认为是亚克隆性的。

Conserved pancreatic cancer enhancer regions and gene annotation

从先前的ENCODE研究中获得了人类和小鼠之间保守的增强子区域。与ENCODE保守增强子区域有任何重叠的胰腺癌增强子区域被定义为“保守的胰腺癌增强子区域”。在62,015个胰腺癌增强子区域中，有8,966个区域符合此标准。使用GREAT将每个保守增强子区域和包含在这些区域中的突变分配给500 kb内的第一和第二近的基因。

Conserved enhancer mutation and gene expression

使用先前为尸检队列和ICGC队列生成的RNA-seq数据，比较了具有和不具有保守增强子突变的样本之间的基因表达。从该分析中移除了具有非同义编码突变的样本，以避免对转录稳定性的潜在影响。使用Wilcoxon符号秩检验比较每个基因的基因表达变化。为了确认该分析的可靠性，对于每个基因，跨样本随机打乱突变类别，使用保守的胰腺癌增强子区域作为参考。使用Wilcox检验比较表达水平。该过程重复10,000次以生成参考P值分布。然后将原始P值与这些参考值进行比较，并计算小于原始的参考P值的比例，从而产生更新后的调整后P值。将尸检和/或TCGA队列中具有显著上调或下调变化的基因用于通过Enrichr进行基因本体分析。

Allele-specific expression analysis

使用cis-X进行等位基因特异性表达分析，遵循原始出版物中描述的方法。简而言之，我们组织来自PDAC样本的多模式基因组数据，使用cis-X计算框架进行整合分析。输入数据包括由Control-FREEC生成的拷贝数变异谱、由Strelka2调用的SNV和插入/缺失、使用STAR和featureCounts量化为FPKM值的RNA表达以及WGS比对文件。使用具有默认参数的cis-X Singularity容器进行ASE分析。

Motif analysis on enhancer sequences

使用原始突变汇总表提取在保守的胰腺癌增强子区域突变和野生型变异之间显示显著表达差异的基因。然后从原始突变数据中提取突变位置，并将其映射到上述原始保守增强子区域。使用Hypergeometric Optimization of Motif EnRichment中的findGenomeMotif.pl分析这些保守增强子区域，以识别每个区域中的已知基序以及突变发生在这些基序内的程度。

Statistics and reproducibility

所有统计和图表均使用GraphPad Prism和/或R进行和制作。通过双侧χ²检验比较参数分布，对于样本量<5的情况使用Fisher精确检验进行校正。使用Mann-Whitney U检验比较非参数分布，对于列联表分析，使用双侧Fisher精确检验。使用Kaplan-Meier方法进行总生存分析，并通过对数秩检验比较曲线。如果P值小于0.05则认为具有统计学显著性。基因集富集分析使用FDR q值。未使用统计方法预先确定样本量。只要文库和/或测序质量符合我们的标准，就没有数据被排除在分析之外。实验未进行随机化。除组织学切片审查外，研究人员在实验和结果评估期间对分配不知情。

Results

Clinicopathologic and genetic features of the patient cohort

该队列的临床特征总结在补充表中。10名患者中有8名在临床IV期胰腺癌诊断，两名为IIB期。后两名患者接受了具有治愈意图的胰十二指肠切除术，但随后发生了复发性转移性疾病。四名患者因诊断时体能状态差而未接受治疗，生存期短，而六名患者接受了化疗，生存期相对较长。在尸检时，所有10名患者至少在一个器官中经病理学证实有转移。10例中有9例是导管腺癌，1例具有未分化癌伴破骨细胞样巨细胞的特征。两例导管腺癌显示局灶鳞状分化，第三例具有神经内分泌特征。本研究共使用了来自这10名患者的86个PDAC基因组和10个正常基因组。我们在每个样本中鉴定出中位数为15,312个体细胞突变，对应于每个样本中位数为11,671个SNV、1,903个小插入和2,664个小缺失。检测到编码驱动基因突变或拷贝数变化，与先前大型队列研究的结果一致。为了理解治疗与突变负荷之间的关系，我们根据接受标准护理化疗的临床病史独立评估病例。未经治疗的PDAC比经过治疗的PDAC具有显著更少的突变数。为每位患者推导了样本系统发育以确定体细胞变异发生的相对进化时间。常见的驱动基因在所有患者中都是主干起源。三名来自接受治疗患者的PDAC在驱动基因中存在亚克隆突变。每个亚克隆驱动基因改变对所分析患者的单个样本是私有的，最常见于原发肿瘤内，提示治疗诱导的遗传瓶颈。与这种解释一致，接受治疗的PDAC比未经治疗的PDAC具有更长的分支。在10个病例中的6个中，至少一个原发肿瘤样本与转移瘤最密切相关，表明原发肿瘤内存在亚克隆异质性和转移亚克隆的形成。这些发现共同证明该PDAC队列表现出驱动基因预期的克隆动力学，并为解释晚期PDAC中的非编码遗传改变提供了基线。

Mutational characteristics of the PDAC genome

所有突变被细分为发生在编码区、启动子、增强子或非编码非增强子区域的突变。总体而言，增强子区域存在的突变最少，非编码非增强子区域最多。当我们考虑每百万碱基的突变频率时，我们发现编码、启动子、增强子和非编码非增强子区域的中位数突变频率。为了确定接受治疗的PDAC与未经治疗的PDAC在多大程度上不同，我们计算了每组中每个基因组区域的mut/Mb。在未经治疗的PDAC中，我们发现编码、启动子、增强子和非编码非增强子区域的中位数mut/Mb值。在接受治疗的PDAC中，每个区域的中位数mut/Mb至少高出三倍。接受治疗的PDAC与未经治疗的PDAC突变频率的成对比较表明，与未经治疗的PDAC相比，基因组所有区域的突变负荷显著更大。我们接下来确定了在每個PDAC的进化生命史过程中突变积累的相对时间点。在未经治疗的PDAC中，我们发现编码、启动子、增强子和非编码非增强子区域的克隆mut/Mb和亚克隆mut/Mb中位数值，表明大多数突变发生在最近共同祖先形成之前。我们在接受治疗的PDAC中发现了类似的关系。这些发现共同提示了非编码基因组的两个特征：第一，大多数突变发生在最近共同祖先形成之前；第二，尽管存在治疗诱导的遗传瓶颈，编码、启动子、增强子和非编码非增强子区域的突变比例得以保持。

Mutational signatures in PDAC

为了识别每个基因组区域突变特征的一般特征，我们使用了Palimpsest并基于每个队列的突变比例创建了热图。基于组合突变特征的无监督聚类分析识别出两个组，其中四个基因组区域非随机分布。具体而言，与基因组的其他区域相比，增强子区域表现出不同的特征特征。这伴随着SBS39特征的显著富集以及SBS1、SBS8、SBS40和SBS41特征的减少。为了解释样本量效应，我们使用相同数量的突变对所有类别进行子抽样分析，并确认聚类结果相似。基于突变特征的主成分分析也在第一主成分中将增强子改变识别为独特的。我们接下来确定了这些突变特征在多大程度上与PDAC的其他特征相关，例如突变的克隆性或治疗。SBS1和SBS18是唯一两个在克隆性突变部分中比在亚克隆进化过程中出现的突变中更普遍的特征。当将增强子改变细分为克隆性或亚克隆性时，它们在PC1中仍然保持独特。SBS1和SBS5在未经治疗的病例中也显著更普遍，与治疗期间或由治疗 specifically 引起的替代突变过程的富集一致。总体而言，这些发现表明突变过程在这四个基因组区域中以环境依赖的方式不同，增强子突变具有与其他区域相比最独特的特征。

Orthogonal validation in the ICGC dataset

为了确定我们中等规模尸检队列中的发现在多大程度上具有普遍性，我们评估了来自ICGC的正交集中的259个PDAC的突变负荷和特征。我们在每个样本中鉴定出中位数为5,899个体细胞突变，对应于每个区域中位数的编码区、启动子、增强子和非编码非增强子区域突变。mut/Mb揭示了编码区、启动子、增强子和非编码非增强子区域的中位数突变。基于突变特征的主成分分析也将增强子改变识别为独特的，尽管程度不如尸检队列中明显。最后，我们确定了最能区分每个基因组区域的突变特征。与尸检队列相似，我们发现SBS39在增强子区域富集，而SBS1、SBS5、SBS8和SBS40在增强子区域中代表性不足。总体而言，我们得出结论，这个大型早期PDAC队列中的非编码基因组特征与尸检队列中观察到的相似。

Mutations in conserved enhancer regions in PDAC

启动子区域的复发性突变先前已有报道。因此，并且由于本研究中注意到的增强子的独特特征，我们专注于增强子突变在PDAC中复发的程度。从物种进化的角度来看，增强子区域分为两类：保守和非保守区域。因此，我们依赖先前生成的ChIP-seq数据，该数据基于H3K27激活的组蛋白标记定义了PDAC增强子区域，并根据小鼠和人类基因组中预测的同源增强子的重叠将这些区域分类为保守和非保守类别。该策略识别了8,966个保守的增强子区域。在ICGC和尸检队列中，保守增强子区域的突变频率低于非保守增强子区域。我们注释了这些保守的增强子区域以识别它们最近和第二近的基因。在这些基因中，我们专注于那些在其相关的保守增强子区域中具有复发性突变的基因。在尸检队列中，有162个基因在三个或更多独特样本中发生突变，其中32个基因在两个或更多患者中受影响。通过ddPCR确认和验证了代表性的保守增强子突变。在ICGC队列中，我们发现了135个基因，其中三个或更多患者发生突变，其中18个在五个或更多患者中复发性突变。使用这些阈值在两个队列中识别出11个基因，其中几个与癌症生物学或PDAC specifically 相关。这11个基因编码的蛋白质产物功能为转录因子或辅助转录因子、对缺氧的反应、细胞外基质的调节和肌动蛋白重塑。诸如DCBLD2和TPRG1等基因在癌症中研究较少，尽管后者进一步暗示了TP63轴。

Gene expression in association with mutated enhancers

保守增强子区域中的体细胞改变并不表明是否可能发生基因表达的变化。因此，对每个增强子最近的两个基因的转录表达进行了评估。具体而言，对于每个具有突变的保守增强子，我们比较了同一组织中相关基因的基因表达。对于后续分析，仅在尸检和ICGC队列中考虑WGS-RNA-seq匹配的样本。我们在合并的尸检和ICGC队列中识别出92个基因，这些基因在具有和不具有增强子突变的样本之间存在显著的基因表达差异。与增强子突变相关的基因表达最常见的是降低。下调的基因再次包括尸检队列中的转录因子KLF5和POU5F1B，或ICGC队列中的ELF3和RUNX1等。TP63也被再次识别，并且与增强子突变相关的基因表达降低了10倍以上。相比之下，AXIN2因附近增强子的体细胞突变而表达增加而引人注目。虽然单个基因提供了一些关于增强子突变所针对表型的线索，但我们结合了ICGC和尸检队列中与增强子突变相关的显著改变表达的数据，并进行了基因本体分析以进一步了解。下调基因高度富集于参与DNA模板转录调控的基因。相比之下，没有发现与上调基因表达相关的通路。为了确定所识别的复发性增强子突变影响已知或预测的增强子基序区域的程度，我们使用Hypergeometric Optimization of Motif EnRichment来识别在突变保守增强子区域相对于野生型序列特别富集或缺乏的调控元件。根据突变计数和基因表达变化对增强子序列进行分类。先前分析识别的基因均未被识别。然而，该分析揭示了两个与ELF5和Isl1相关的已知共有区域显著富集。我们接下来进行了de novo分析以评估预测的基序区域，揭示了与Smad2、Nanog和Sox18相关的基序序列。

Discussion

最近的大规模努力极大地促进了癌症非编码驱动因子的发现和注释。我们现在将PDAC中发生的遗传事件范围扩展到保守增强子区域中的复发性突变。本研究中识别的基因均未被泛癌全基因组非编码DNA分析所识别，与作者的结论一致，即肿瘤类型特异性研究将揭示与该组织细胞起源相关的其他改变。与这一预测一致，我们识别了具有与胰腺生物学和稳态相关的基因产物潜在联系的复发性增强子突变。通过基序分析识别了其他潜在靶点，例如TGFΒ信号传导的经典介质Smad2。虽然需要额外的工作来完全理解复发性增强子突变在多大程度上促进PDAC生物学，但基因表达与许多这些增强子突变相关而显著变化，而其他突变与预测的调控基序区域相关，支持了它们的功能作用。一个特别重要的发现是在两个研究的PDAC队列中影响TP63及其靶基因TPRG1的复发性增强子突变。在ICGC队列中，与这些突变相关的TP63基因表达也显著降低，而TPRG1已被识别为B细胞淋巴瘤顺式调控区域中体细胞突变的复发性靶点。TP63是增强子重编程的已知主调控因子，其在PDAC中的核表达与基底样或鳞状特征高度相关。我们注意到我们尸检队列中的两名患者具有基底样特征，但这两个PDAC或具有这种组织学的特定区域均未发现TP63附近增强子的突变。因此，可以设想这些调控元件中的突变可能负调控TP63表达并强化经典型PDAC表型，特别是因为GATA6过表达或扩增是经典表型和对吉西他滨/白蛋白结合型紫杉醇反应的不完美标志物。至少，在PDAC中探查这个保守的增强子区域可能会扩大该疾病一线化疗的预测性生物标志物列表。本研究中发现的总体突变率与文献一致，并且也突出了治疗的预期效果，因为接受治疗的PDAC的总体mut/Mb比未经治疗的PDAC高3.8倍。然而，我们现在发现突变率因基因组区域而异。具体而言，我们发现增强子的突变率显著低于其他非编码区域，并且具有独特的突变特征。这些差异可能对应于许多因素，包括但不限于增强子和内含子DNA之间染色质结构的差异，以及基因组不同区域DNA修复效率的差异。可以设想，增强子区域内发现的低突变率也可能反映了与其他非编码基因组区域相比的相对负选择。传统上，正选择或负选择的推断考虑了蛋白质编码基因中非同义与同义变异的比率。然而，这种方法不仅对于非编码DNA不足，而且没有考虑局部突变率的可变性。鉴于数据显示增强子DNase I超敏感位点内的突变比在内含子区域中发现的相同序列更少见，我们的发现最可能反映了较低的总体突变积累率，而不是负选择，尽管需要专门的研究来正式证明这种可能性。我们的研究仍然有局限性。例如，虽然尸检队列允许评估终末期PDAC中的增强子突变，但它包含的患者数量较少。我们用来自ICGC数据集的大量患者队列补充了这一点，尽管并非该后一队列中的所有患者都有RNA-seq数据来最佳地确定保守增强子突变对基因表达的影响。我们还承认增强子可能影响比我们在本研究中关注的最近两个基因更远的基因；因此，我们的发现可能相对保守地代表了受增强子突变影响的基因。最后，虽然我们在两个队列之间观察到一些重叠，但与复发性增强子突变相关的大多数基因对于ICGC或尸检队列是独特的。为了进一步评估功能性，我们使用cis-X在尸检队列中进行了等位基因特异性表达分析，该分析识别了一组不同的候选基因。这是可以预见的，因为我们的研究专注于保守的增强子突变，而cis-X旨在检测所有增强子中的等位基因特异性效应。这些互补的方法自然产生不同的候选基因，共同提供了增强子介导调控的更广泛视角。对增强子突变与转录等位基因进行定相的长读长测序将是有价值的，但这超出了本研究的范围，并被视为未来工作的一个限制。总之，这些发现揭示了PDAC中值得额外研究的一层额外的基因组改变，并拓宽了我们对癌症中增强子介导调控的理解。