胰腺癌保守增强子突变：非编码基因组的新视角与功能影响

《Cancer Research Communications》：Conserved enhancer mutations and gene expression. A, Schematic of approac... Open Access

本研究通过全基因组测序揭示胰腺癌非编码基因组中保守增强子区域的复发性体细胞突变特征，发现增强子突变具有独特的突变特征（SBS39富集），并与KLF5、TP63等胰腺发育关键转录因子相关，其通过调控基因表达（多为表达缺失）影响肿瘤生物学行为，为胰腺癌发病机制提供了新见解。

　　

Abstract

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

Introduction

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

Materials and Methods

Ethics statement

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

Tissue samples

使用了来自MSKCC“Last Wish Program”的六个PDAC研究 autopsy 样本（MPAM01-06）。从冷冻切片上切取组织片并进行审查，以确定肿瘤细胞含量至少为20%且组织质量保存完好的样本。对正常样本进行审查以确认没有污染癌细胞。符合这些标准的样本在提取基因组DNA之前从连续未染色的切片中进行宏观解剖以富集肿瘤纯度。本研究还包括了先前通过约翰霍普金斯大学快速 autopsy 项目分析的四个未经治疗的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×和47.5×。

Filtering and annotation of variants

对所有患者，使用以下标准过滤体细胞变异：患者匹配的正常样本覆盖度≥10个读数，患者匹配的正常样本中的变异计数≤2，患者匹配的正常样本变异频率<0.02，肿瘤突变等位基因频率≥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

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

Pancreatic cancer enhancer regions and enhancer mutations

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

Mutational signature analysis

使用Palimpsest进行突变特征分析。每个样本在编码、增强子、启动子或非编码-非增强子区域中识别的所有突变用作输入。每个样本每个类别中的特征比例用作聚类分析的输入，以识别增强子区域突变特征的特征。执行Wilcoxon符号秩检验以识别组间的特征性特征。使用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

使用先前为autopsy队列和ICGC队列生成的RNA-seq数据，比较了有和没有保守增强子突变的样本之间的基因表达。从该分析中移除了具有非同义编码突变的样本，以避免对转录稳定性的潜在影响。使用Wilcoxon符号秩检验比较每个基因的基因表达变化。为了确认该分析的可靠性，对于每个基因，随机打乱样本间的突变类别，使用保守的胰腺癌增强子区域作为参考。使用Wilcox检验比较表达水平。该过程重复10,000次以生成参考P值分布。然后将原始P值与这些参考值进行比较，并计算小于原始P值的参考P值的比例，从而产生更新后的调整后P值。将autopsy和/或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期胰腺癌，2名为IIB期。后两名患者接受了以治愈为目的的胰十二指肠切除术，但随后出现了复发性转移性疾病。四名患者（PAM01-04）因诊断时体能状态差而未接受治疗，生存期较短（范围0.5-10个月），而六名患者（MPAM01-06）接受了化疗，生存期相对较长（范围9-49个月）。在 autopsy 时，所有10名患者至少有一个器官经病理学证实存在转移。10例中有9例为导管腺癌，1例（MPAM02）具有未分化癌伴破骨细胞样巨细胞的特征。两例导管腺癌PAM02和MPAM06显示局灶鳞状分化，第三例（PAM01）具有神经内分泌特征。

本研究共使用了来自这10名患者的86个PDAC基因组和10个正常基因组。我们在每个样本中识别出中位数为15,312个体细胞突变，对应于每个样本中位数为11,671个SNV、1,903个小插入和2,664个小缺失。检测到KRAS、TP53、CDKN2A、SMAD4、ATM、ARID1A、ARID2和PBRM1的编码驱动基因突变或拷贝数变化，与先前大型队列研究的结果一致。为了理解治疗与突变负荷之间的关系，我们根据接受标准护理化疗的临床病史独立评估病例。未经治疗的PDAC的突变数量显著低于经过治疗的PDAC。

为每位患者推导了样本系统发育，以确定体细胞变异发生的相对进化时间。常见的驱动基因在所有患者中均为主干起源。三名来自接受治疗患者的PDAC在MUC6、B2M或DNMT3A中存在亚克隆驱动基因突变。这些亚克隆驱动基因改变中的每一个都是患者分析的单个样本所特有的，最常见于原发肿瘤内，提示治疗诱导的遗传瓶颈。与这种解释一致，经过治疗的PDAC比未经治疗的PDAC具有更长的分支。在10例中的6例中，至少一个原发肿瘤样本与转移瘤关系最密切，表明原发肿瘤内存在亚克隆异质性和转移性亚克隆的形成。总的来说，这些发现表明该PDAC队列在驱动基因方面表现出预期的克隆动力学，并为解释晚期PDAC中的非编码遗传改变提供了基线。

Mutational characteristics of the PDAC genome

所有突变被细分为发生在编码区、启动子、增强子或非编码-NOS区域的突变。总体而言，增强子区域存在的突变最少，而非编码-NOS区域最多。当我们考虑每百万碱基的突变频率时，我们发现编码、启动子、增强子和非编码-NOS区域的中位数突变频率。为了确定经过治疗的PDAC与未经治疗的PDAC的差异程度，我们计算了每组中每个基因组区域的mut/Mb。在未经治疗的PDAC中，我们发现编码、启动子、增强子和非编码NOS区域的克隆mut/Mb中位数值，以及亚克隆mut/Mb中位数值，表明大多数突变发生在最近共同祖先形成之前。我们在经过治疗的PDAC中发现了类似的关系。总的来说，这些发现提示了非编码基因组的两个特征：第一，大多数突变发生在最近共同祖先形成之前；第二，尽管存在治疗诱导的遗传瓶颈，编码、启动子、增强子和非编码-NOS区域的突变比例得以保持。

Mutational signatures in PDAC

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

Orthogonal validation in the ICGC dataset

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

Mutations in conserved enhancer regions in PDAC

由于启动子区域的复发性突变先前已有报道，并且因为本研究中注意到的增强子的独特特征，我们专注于增强子突变在PDAC中复发的程度。从物种进化的角度来看，增强子区域分为两类：保守和非保守区域。因此，我们依赖先前生成的ChIP-seq数据，该数据基于H3K27ac激活的组蛋白标记定义了PDAC增强子区域，并根据小鼠和人类基因组中预测的同源增强子的重叠将这些区域分类为保守和非保守类别。该策略识别了8,966个保守增强子区域。在ICGC和 autopsy 队列中，保守增强子区域的突变频率均低于非保守增强子区域。

我们注释了这些保守增强子区域以识别它们最近和第二近的基因。在这5,618个基因中，我们专注于那些在其相关的保守增强子区域中存在复发性突变的基因。在 autopsy 队列中，有162个基因在三个或更多独特样本中发生突变，其中32个基因在两名或更多患者中受到影响。在ICGC队列中，我们发现了135个基因，其中三名或更多患者发生突变，其中18个基因在五名或更多患者中复发性突变。使用这些阈值（三个 autopsy 样本和三个ICGC患者）在两个队列中识别出11个基因，其中几个与癌症生物学或PDAC特异性相关。这11个基因编码的蛋白质产物功能为转录因子或辅助转录因子、缺氧反应、细胞外基质调节和肌动蛋白重塑。诸如DCBLD2和TPRG1等基因在癌症中研究较少，尽管后者进一步暗示了TP63轴。

Gene expression in association with mutated enhancers

保守增强子区域中的体细胞改变并不表明是否可能发生基因表达的变化。因此，对每个增强子最近的两个基因进行了转录表达评估。具体来说，对于每个带有突变的保守增强子，我们比较了同一组织中相关基因的基因表达。对于后续分析，仅考虑 autopsy 和ICGC队列中WGS-RNA-seq匹配的样本。

我们在合并的 autopsy 和ICGC队列中识别出92个基因，其基因表达在有和没有增强子突变的样本之间存在显著差异。基因表达最常在具有相关增强子突变的样本中降低。下调的基因再次包括 autopsy 队列中的转录因子KLF5和POU5F1B，或ICGC队列中的ELF3和RUNX1等。TP63也被再次识别，并且显示与增强子突变相关的基因表达降低了10倍以上。相比之下，AXIN2因附近增强子的体细胞突变而表达增加而引人注目。

虽然单个基因为了解增强子突变所针对的表型提供了一些线索，但我们结合了与 autopsy 和ICGC队列中增强子突变相关的显著改变表达的数据，并进行了基因本体分析以进一步了解。下调基因高度富集了参与DNA模板转录调控的基因。相比之下，没有发现与上调基因表达相关的通路。

为了确定所识别的复发性增强子突变影响已知或预测的增强子基序区域的程度，我们使用Hypergeometric Optimization of Motif EnRichment来识别在突变的保守增强子区域相对于野生型序列特别富集或缺乏的调控元件。增强子序列根据突变计数和基因表达变化进行分类。先前分析识别的基因均未被识别。然而，该分析揭示了两个已知的与ELF5和Isl1相关的共有区域显著富集。我们接下来进行了de novo分析以评估预测的基序区域，揭示了与Smad2、Nanog和Sox18相关的基序序列。

Discussion

近期的大规模努力极大地促进了泛癌或特定肿瘤类型（包括PDAC）中非编码癌症驱动因素的发现和注释。我们现在将PDAC中发生的遗传事件范围扩展到保守增强子区域中的复发性突变。本研究中识别的任何基因均未被泛癌全基因组分析所识别，这与作者的结论一致，即肿瘤类型特异性研究将揭示与该组织细胞起源相关的其他改变。与这一预测一致，我们识别了与蛋白质产物有助于胰腺生物学和稳态的基因相关的复发性增强子突变，例如KLF5等。通过基序分析识别了其他潜在靶点，例如TGFΒ信号传导的经典介质Smad2。虽然需要额外的工作来完全理解复发性增强子突变对PDAC生物学的贡献程度，但基因表达与许多这些增强子突变显著相关，而其他突变与预测的调控基序区域相关，支持了它们的功能作用。

一个特别重要的发现是在两个研究的PDAC队列中，影响TP63及其靶基因TPRG1的复发性增强子突变。在ICGC队列中，与这些突变相关的TP63基因表达也显著降低，而TPRG1已被确定为B细胞淋巴瘤顺式调控区域中体细胞突变的复发性靶点。TP63是增强子重编程的已知主调控因子，其在PDAC中的核表达与基底样或鳞状特征高度相关。我们注意到 autopsy 队列中有两名患者具有基底样特征，但这两种PDAC或具有这种组织学的特定区域均未发现TP63附近增强子的突变。因此，可以想象这些调控元件中的突变可能负向调控TP63表达并强化经典型PDAC表型，特别是因为GATA6过表达或扩增是经典型表型和对吉西他滨/白蛋白结合型紫杉醇反应的不完美标志物。至少，在PDAC中 interrogation 这个保守增强子区域可能会扩大该疾病一线化疗的预测性生物标志物列表。

本研究中发现的总体突变率与文献一致，并且也突出了治疗的预期效果，即经过治疗的PDAC的总体mut/Mb比未经治疗的PDAC高3.8倍。然而，我们现在发现突变率因基因组区域而异。具体而言，我们发现增强子的突变率显著低于其他非编码区域，并且具有独特的突变特征。这些差异可能对应于许多因素，包括但不限于增强子（常染色质）和内含子DNA（异染色质）之间染色质结构的差异，以及基因组不同区域DNA修复效率的差异。可以想象，增强子区域内低突变率也可能反映了与其他非编码基因组区域相比的相对负选择。传统上，正选择或负选择的推断考虑了蛋白质编码基因中非同义与同义变异的比率。然而，这种方法不仅对非编码DNA不充分，而且没有考虑局部突变率的变异性。鉴于数据显示增强子DNase I超敏感位点内的突变比在内含子区域中发现的相同序列更不常见，我们的发现最有可能反映了突变积累的总体速率较低，而不是负选择，尽管需要专门的研究来正式证明这种可能性。

我们的研究仍然有局限性。例如，虽然 autopsy 队列允许评估终末期PDAC中的增强子突变，但它包含的患者数量较少。我们用来自ICGC数据集的大型队列补充了这一点，尽管该后一队列中并非所有患者都有RNA-seq数据来最好地确定保守增强子突变对基因表达的影响。我们还承认增强子可能影响比我们在本研究中关注的最近两个基因远得多的基因；因此，我们的发现可能相对保守地代表了受增强子突变影响的基因。最后，虽然我们在两个队列之间观察到一些重叠，但大多数与复发性增强子突变相关的基因是ICGC或 autopsy 队列所独有的。为了进一步评估功能性，我们使用cis-X在 autopsy 队列中进行了等位基因特异性表达分析，识别出了一组不同的候选基因。这是意料之中的，因为我们的研究专注于保守的增强子突变，而cis-X旨在检测所有增强子中的等位基因特异性效应。这些互补的方法自然会产生不同的候选基因，共同提供了增强子介导调控的更广泛视角。使用长读长测序将增强子突变与转录本等位基因进行分相将是有价值的，但这超出了本研究的范围，并作为未来工作的一个限制。总之，这些发现揭示了PDAC中基因组改变的额外层面，值得进一步研究，并拓宽了我们对癌症中增强子介导调控的理解。