前列腺肿瘤甲基化现象的概况

《Cancer Discovery》：The Landscape of Prostate Tumor Methylation

开放型图形查看器

摘要

前列腺癌的特点是临床和分子异质性极为显著。尽管其基因组异质性已被充分研究，但其表观基因组异质性仍了解较少。因此，我们创建了一个包含3,001个多祖先前列腺甲基组的汇编，这些甲基组涵盖了从正常组织到所有级别的局部疾病以及多灶性转移疾病。其中884个样本进行了多组学DNA和/或RNA分析。我们确定了几种表观基因组亚型，这些亚型可以分层患者的风险，并反映了不同的进化轨迹。我们发现了DNA拷贝数与甲基化之间的广泛调控相互作用，这种相互作用在不同基因和疾病阶段表现出不同的转录后果。我们为15个重要的临床分子特征定义了表观遗传失调的标志，为每个特征创建了预测模型。例如，我们识别出一些特定的表观遗传特征，这些特征可以预测患者的预后，并与临床预后特征具有协同作用。这些结果定义了肿瘤遗传学与表观遗传学之间的复杂相互作用，这种相互作用部分通过调节表观遗传老化来改变基因表达程序和临床表现。

重要性：我们创建了迄今为止最大的前列腺癌甲基组资源，揭示了四种表观基因组亚型，这些亚型根据癌症的遗传学、进化和临床表型进行分层。这些数据表明，基因组-表观基因组的相互作用取决于疾病阶段，从而将DNA甲基化定位为前列腺癌异质性的一个核心、具有临床信息价值的驱动因素。

引言

2023年，美国有超过280,000名男性被诊断出患有前列腺癌；这是男性癌症死亡的第二大原因（1）。前列腺癌具有高度遗传性，约60%的发病率变异可以由罕见和常见的种系遗传因素解释（2）。虽然大多数新诊断的病例在疾病处于局部阶段时被发现，因此有可能治愈，但治疗结果差异很大（3）。当前的临床预后因素只能部分解释对一线局部和系统性治疗的临床结果的异质性（4）。大多数局部疾病具有相对缓慢的临床进程，但一旦前列腺癌发生转移并对抗雄激素疗法产生耐药性，就变得无法治愈（5）。这种临床异质性与广泛的基因组异质性相匹配。特别是局部疾病表现出显著的肿瘤内和肿瘤间变异，这与不同的临床结果相关（6-10）。体细胞遗传异质性与广泛的表观基因组和基因表达变异同时存在（7, 11-16）。CpG岛（CGI）在基因启动子区域的异常甲基化与肿瘤抑制基因的转录失活、癌基因激活和基因组不稳定性有关（17-20）。在包括前列腺癌在内的多种癌症中鉴定出一种高甲基化亚型（21-23）。数百个位点的甲基化与不良临床结果相关（8, 24-27），并且与临床使用的生物标志物具有叠加或协同作用（28）。RNA甲基化驱动异常基因表达和前列腺癌的进展（29）。尽管这些研究提供了有价值的信息，但它们在样本量和临床数据方面存在局限性，并未涵盖疾病的全部自然史。因此，我们对DNA甲基化（DNAm）如何与肿瘤遗传学相互作用以在疾病不同阶段塑造临床分子表型知之甚少。为了填补这一空白，我们创建了一个包含3,001份前列腺标本的DNAm汇编：2,355份来自癌症样本，646份来自非癌性前列腺组织样本。这些肿瘤代表了局部和转移性疾病中的所有风险组，并来自具有不同祖先背景的个体。我们确定了与疾病阶段和进化相关的不同DNAm亚型。这使我们能够量化DNAm和拷贝数如何共同调控基因表达，以及这种调控在疾病进展过程中的演变。我们全面定义了体细胞驱动突变的表观基因组和转录组后果，并证明这些可以准确预测临床结果。由于肿瘤进展不仅反映了基因组改变，还反映了通路状态，我们系统地检查了不同甲基化亚型中的通路活性，包括在PTEN丢失等重复驱动事件发生的肿瘤中。特别是，我们关注了缺氧及相关过程（如血管生成和糖酵解），因为这些过程与肿瘤进展和治疗耐药性密切相关。总之，这些分析表明，遗传和表观遗传变异之间的相互作用合理化了前列腺癌的临床异质性。

结果

前列腺癌的甲基化亚型

我们创建了一个包含3,001个前列腺组织甲基组的队列，这些数据来自14个已发表和1个未发表的队列（7, 8, 13, 15, 25, 30-38）。该汇编涵盖了从正常组织（22%；n = 646）到局部癌症（70%；n = 2,102）再到去势抵抗性转移疾病（8%；n = 253）的全谱。局部癌症代表了国际泌尿病理学家学会（ISUP）的所有五个等级组（GG），涉及117至403名患者（补充图S1A-S1E；补充表S1-S3）。我们首先在最大的单一样本队列“The Cancer Genome Atlas Prostate Adenocarcinoma”（TCGA-PRAD）中确定了全面的甲基化亚型，该队列包括498个肿瘤和50个正常前列腺样本，然后在其余14个队列中进行了验证。在之前的小型亚型研究（7, 32, 39, 40）的基础上，我们定义了四种甲基化亚型，分别称为MS-1至MS-4（图1A；补充图S1F）。为了方便将这些亚型应用于新样本，我们确定了5,486个定义亚型的CpG，并创建了一个用于单样本亚型分配的简化的质心分类器（补充图S1G和S2；补充表S4）。我们在一个新的开源R包PrCaMethy中实现了这一点，以便于应用和扩展这里报告的所有方法和标志。图1。

前列腺癌的表观基因组景观。A，来自15个队列的2,355个前列腺肿瘤和646个正常样本的亚型甲基化谱，包括临床治疗前的变量。热图包含了IMPACC选定的5,486个对亚型划分重要的CpG。CpG使用完全连锁和欧几里得距离进行分层聚类，然后分成五个大小相等的组进行可视化。缺失值以白色显示，而正常样本的临床变量（年龄除外）没有显示。队列标签对应于以下参考文献：TCGA-PRAD（7），Zhao及其同事（25），Li及其同事（30），ICGC PRAD-CA（8），Berglund及其同事（33），Dairo及其同事（35），Kirby及其同事（31），ICGC PRAD-FR（32），Ramakrishnan及其同事（34），Zhao及其同事（15），Ylitalo及其同事（36），Mundbjerg及其同事（38），Toth及其同事（37），以及Brockses及其同事（13）；墨尔本队列之前未发表。B-E，对于每个队列（垂直条形），按局部前列腺肿瘤样本（B）、转移样本（C）、正常样本（D）以及合并所有队列（E）显示各亚型的样本百分比。每个队列的总样本数显示在条形上方。L代表局部疾病；M代表转移疾病；N代表正常样本。F，使用Kaplan–Meier曲线分别估计每个亚型的无BCR概率，合并所有具有可用BCR数据的患者（n = 882）。右上角的P值来自对数秩检验。在左上图中，使用Cox比例风险模型估计HR值及其95%置信区间（CI）和Wald检验P值，以MS-1作为参考组，并调整了队列因素。图表下方的表格给出了0-5年内有风险的患者数量。右下角的森林图是对七个队列的BCR相关结果的荟萃分析（详情见“方法”部分）。`“总体”行包含了比较MS-4和MS-1的总体荟萃分析HR值。点表示由相邻线条给出的估计HR值及其95% CI。G，点图总结了亚型与四个临床特征（年龄、T分类、ISUP GG和PSA）的总体关联，合并所有患者时（总体）以及按队列分别统计。每个临床变量均按照补充表S1中的分类方式进行分类，并使用Fisher精确检验来检测其与四种甲基化亚型的关联。点的大小反映了校正偏差后的Cramér’s V效应大小，背景颜色表示q值，该值分别调整了每个临床变量的多重检验。单元格中的X表示该特征在该队列中不可用。

H-J，使用Fisher精确检验和OR值比较了每条染色体内定义亚型的“背景”CpGs与非定义亚型的“背景”CpGs的比率（H），在不同类型的CpG中（岛屿、开阔海域、大陆架和海岸；I）以及在基因体、基因间和基因启动子区域中（如“方法”中所定义；J）。图1。

四种甲基化亚型与疾病的临床特征紧密相关。例如，亚型MS-2是“正常样”亚型。它主要出现在低级别局部疾病中（64%的ISUP等级1或2；图1B），在转移疾病中极为罕见（1%；图1C），并在非恶性前列腺组织中占主导地位（95%；图1D）。亚型MS-4则表现出相反的模式，在侵袭性疾病中最常见：非恶性组织中占1%，局部疾病中占34%，转移疾病中占80%。亚型不受肿瘤纯度、自我鉴定种族（SIR）或遗传祖先的影响（补充图S3A；补充信息S1），并且在不同疾病阶段之间存在显著差异（P < 10^-45；图1E；补充图S3B）。对于局部疾病患者，亚型MS-1和MS-2的生化复发（BCR）风险最低，而MS-4的风险最高（HR = 2.4；P = 3.9 × 10^-5；图1F；补充图S3C-S3J）。即使调整了队列、年龄、病理T分类、ISUP GG、治疗前血清前列腺特异性抗原（PSA）和总拷贝数改变（CNA）负担后，MS-4仍然具有独立的预后预测价值（调整后的HR = 2.01；P = 0.0078）。这项汇总分析通过七个具有BCR相关结果的队列的荟萃分析得到了验证：在1,534名患者中，MS-4的复发风险是MS-1的2.7倍[图1F（插图）]，并且没有证据表明各队列之间的效应大小存在异质性（P = 0.52）。甲基化亚型与治疗前的预后特征相关，包括肿瘤大小和范围（T类别）、ISUP GG和血清PSA水平，这些效应大小在各队列以及整个研究中都是一致的（q值<0.05，图1G；补充图S3K–S3N）。这些MS亚型与之前在TCGA-PRAD队列中识别出的CGI甲基化表型（CIMP）显示出显著的一致性（Qannari检验P < 0.001；调整后的Rand指数（ARI）= 0.14；95%置信区间（CI），0.13–0.15；补充图S4A；参考文献23）。CIMP?样本主要被归类为MS-1或MS-2（84%；78/93），而CIMP+样本主要被归类为MS-3或MS-4（88%；52/59）。然后我们计算了所有3,001个样本的CIMP得分（补充表S4）。在整个队列中，即使调整了平台和队列因素后，CIMP得分仍然显著依赖于MS亚型（Kruskal–Wallis检验，P = 9.9 × 10?277；ε2 = 0.43；95%置信区间，0.40–0.46）。MS-3和MS-4的CIMP得分较高，而MS-2的得分最低（成对Mann–Whitney U检验，q < 0.05；补充图S4B）。几个临床分子特征，包括ISUP GG和T类别，与CIMP得分相关（Kruskal–Wallis检验，q < 0.05；补充图S4C–S4I）。几个关键的驱动突变与CIMP增加有关，但值得注意的是TMPRSS2::ERG融合事件则没有这种关联（补充图S4J）。绝大多数定义亚型的CpGs（5,407/5,486；98.6%）在肿瘤和正常组织之间显示出显著的差异（U检验q值<0.05；补充图S5A）。定义亚型的CpGs在肿瘤和正常样本之间的差异比非定义亚型的CpGs更明显（相对风险=1.16；P = 7.2 × 10?259）。同样，96.3%的定义亚型的CpGs在局部位点和转移性位点之间的差异也显著（相对风险=1.19；P = 9 × 10?257；补充图S5B）。在四种甲基化亚型中，这些CpGs也显示出不同的β值分布，当使用正常样本表示MS-2时，所有成对比较都显著（补充图S5C）。定义亚型的CpGs在五个染色体上的富集程度较高，在七个染色体上则较少（图1H），并且更倾向于出现在开放海域而不是岛屿或海岸区域（图1I），以及基因间区域而不是基因启动子或基因本体（图1J）。肿瘤的拷贝数异常（CNA）驱动着区域性的甲基化和基因表达。

这些定义亚型的CpGs在基因组中的位置偏差表明拷贝数状态和甲基化状态之间存在关联（8, 41–43）。因此，我们考虑了我们的研究汇编中具有匹配CNA分析的870个肿瘤。甲基化与拷贝数状态主要呈负相关，86%的基因显示出负相关（图2A和B）。鉴于基于阵列的技术无法量化等位基因特异性的甲基化状态，目前尚不清楚每个等位基因的平均甲基化程度是否随着拷贝数的不同而变化。启动子区域的CpG的甲基化与转录本丰度的降低显著相关，92%的基因显示负相关（图2C）。图2。前列腺癌的调控景观。A轴表示在基因组的该位置具有拷贝数增加（红色条形）或减少（蓝色条形）的患者百分比（x轴）。垂直的洋红色条形表示在F、B-D中后续显示的基因。B-C项表示每个基因内甲基化与CNA（B）、甲基化与mRNA（C）以及CNA与mRNA（D）之间的Spearman相关性。正相关性用蓝色表示，负相关性用橙色表示。每个患者的基因水平甲基化是通过基因启动子区域内CpG的中位β值来测量的。E项表示二元甲基化（过度甲基化与低甲基化；中位数二分法）和二元CNA状态（如果丢失是更常见的突变则为丢失，否则为增加与中性的）之间的线性回归交互系数。因变量是RNA丰度（z分数–转换后的log2 TPM）。显示具有显著交互效应的基因（q < 0.05），以及基因BCL2L10（q = 0.055）。F项展示了六个基因的甲基化与CNA对mRNA丰度的影响示例。检验交互效应的q值由qint表示。每个患者的基因水平甲基化是通过基因启动子区域内CpG的中位β值来测量的，然后中位数二分化为过度甲基化与低甲基化。G项表示，当CNA状态在不同样本之间变化时（更常见的为增加或减少），使用主要的CNA状态。然后，比较了具有CNA的患者与拷贝数正常的患者之间基因启动子区域内CpG的中位甲基化水平（log2倍数变化；y轴）。结果按肿瘤类型以及已知癌症驱动基因与其他基因进行了比较（x轴）。大于0的值表示在拷贝数改变的患者中甲基化程度较高的基因；小于0的值表示在CNA组中甲基化程度较低的基因。接近0的值表示CNA对这些基因的甲基化影响较小。小提琴上方的字母表示组间的成对差异：具有相同字母的两个组之间没有显著差异（Wilcoxon秩和检验q > 0.05），而具有不同字母的两个组之间有差异（q < 0.05）。对于增加（顶部）和减少（底部）的图表，应分别解释这些字母。H-J项表示局部位点和转移性样本中基因内甲基化、CNA和RNA的Spearman相关性分布。基因甲基化是通过基因启动子区域内CpG的中位β值来计算的。图2。鉴于CNA在前列腺癌中是常见且早期的驱动事件，我们假设拷贝数变化可能与甲基化相互作用共同调节RNA丰度。我们首先验证了CNA状态与大多数基因的mRNA丰度呈正相关（86%；图2D）。接下来，我们询问CNA对mRNA的影响是否受该位点甲基化状态的调节。我们通过线性回归中的交互项量化了这些CNA–甲基化对mRNA的影响，重点是具有所有三种数据类型的725个肿瘤。几乎一半的基因表现出显著的CNA–甲基化相互作用（10,919个基因中的4,950个 = 45.3%；q < 0.05；图2E）。因此，对于大约一半的基因组，要完全理解肿瘤CNA变化对基因表达的影响，需要量化该基因组区域的甲基化状态。这些CNA–甲基化相互作用既有协同作用也有拮抗作用，在某些情况下甲基化变化抵消了CNA的影响，在其他情况下则增强了CNA的影响。图2F给出了六个显著例子的相互作用，包括驱动基因PTEN。在每种情况下，CNA增加或减少的患者与拷贝数正常的患者相比，RNA丰度的差异要大得多。PTEN的相互作用与T类别和ISUP GG无关（补充图S5D和S5E）。为了确定这些CNA–甲基化变化是否优先影响癌症驱动基因，我们考虑了先前研究中反复识别出的91个前列腺癌驱动基因（补充表S5；参考文献8, 44–46）。我们按CNA状态和疾病阶段（局部位点 vs. 转移性）对分析进行了分层。我们首先检查了14个已知的反复出现的驱动突变事件（补充图S6）。在转移性样本中，无论驱动事件的状态如何，甲基化–CNA分布相似，而局部位点样本通常在有无驱动突变的患者之间显示出相关分布的中心和范围的移动。接下来，我们评估了91个驱动事件中CNA和甲基化之间的总体关联。在局部位点疾病中，驱动基因和乘客基因的变化平均来说无法区分：驱动基因和非驱动基因的差异甲基化量相似：无论是减少还是增加都是如此（图2G）。相比之下，在转移性疾病中，CNA变化与更大的甲基化变化相关，并且在驱动基因中这种变化显著增强（图2G）。这反映了一个普遍模式：增加的拷贝数与局部位点疾病中的甲基化减少显著相关（图2H；补充图S7A）。同样，mRNA的增加与局部位点疾病中的甲基化减少显著相关（图2I；补充图S7B）。相比之下，拷贝数与mRNA丰度之间的关联在局部位点疾病和转移性疾病之间几乎没有区别（图2J；补充图S7C）。这些数据表明，甲基化在局部位点前列腺癌中的调节作用比在转移性前列腺癌中更为重要（补充图S7D）。为了证实这一结果，我们评估了数据类型对之间的互信息（MI），这可以捕捉非线性关系。与我们的简单相关性分析一致，MI显示局部位点和转移性前列腺癌之间存在显著差异，尤其是在甲基化–RNA和甲基化–CNA关联方面。相比之下，CNA对RNA丰度的影响在局部位点和转移性疾病中相对一致（补充图S8A和S8B）。鉴于CNA–甲基化相互作用在肿瘤间和阶段间的异质性，我们假设在肿瘤间甲基化状态高度可变的基因可以识别出具有预后关联的位点。在甲基化变化最大的前500个基因中，有30个基因与BCR显著相关（q < 0.10；补充图S9A）。与非HVM基因相比，HVM基因在与CNA相关的基因组百分比（PGA）的变化大小和方向上表现出更大的变化，30%的HVM基因与PGA呈正相关，而非HVM基因中只有9%（P = 4.5 × 10?78；补充图S9B）。基因内甲基化与mRNA的相关性也表现出类似的模式（P = 0.01；补充图S9C）。HVM基因与CNA的变化率没有显著差异（P > 0.05；补充图S9D），也没有根据染色体位置表现出富集（补充图）。S9E），并且没有显示出对任何特定生物途径的富集。相反，与BCR显著相关的30个HVM基因显示出强烈的CNA-甲基化相关性，表明它们反映了个体异常事件（补充图S9F和S9G）。这些数据揭示了前列腺癌中拷贝数变化、甲基化和下游基因表达之间的复杂相互作用。拷贝数变化对甲基化具有阶段特异性影响，使得单一体细胞突变能够随着疾病的发展以不同的方式影响基因表达。拷贝数变化的表观基因组后果在肿瘤间的异质性直接影响了驱动基因的表达。

甲基化亚型反映了肿瘤的遗传学和微环境特征。

CNA、甲基化和mRNA丰度之间的复杂关联促使我们探索甲基化亚型的基因组特征。我们使用了序数逻辑回归比较了每个CpG的CNA频率，并调整了PGA和队列因素。对于约20%包含亚型定义CpGs的基因，不同甲基化亚型之间的CNA频率存在显著差异（图3A）。MS-3和MS-4在整个基因组中CNA增多的情况更为普遍，而MS-4在2号、5号、6号和13号染色体上的CNA丢失样本比例更高（图3B；补充图S10）。这些亚型在总体CNA负担（Kruskal–Wallis P = 2.02 × 10?55；图3C）、染色体级别的CNA负担（补充图S11A和S11B）以及驱动CNA负担（补充图S12A和S12B）方面存在显著差异。MS-2肿瘤中的CNA数量最少，这与大多数正常样本（95%）属于该亚型的观察结果一致。综合这些数据表明，定义亚型的CpGs可以区分肿瘤组织和正常组织，并且不同亚型的CNA谱也有所不同。MS-3和MS-4具有显著的PGA，这与它们较高的BCR风险相一致。

甲基化亚型具有不同的CNA和RNA特征：A、包含定义亚型CpGs的1,653个基因的CNA情况。每个基因（热图的一行）通过每个亚型中拷贝数增加或减少的患者百分比进行总结。行按照层次聚类和欧几里得距离进行排序。序数逻辑回归被用来模拟基因的拷贝数状态（丢失、中性、增加）作为因变量，亚型、PGA和队列作为自变量。计算了亚型的似然比检验P值，并在右侧条形图中显示了多重检验调整后的q值。B、y轴表示基因组（x轴）该位置上拷贝数增加（红色条形）或减少（蓝色条形）的患者百分比。结果按甲基化亚型进行分层。C、使用Kruskal–Wallis秩和检验P值比较了不同甲基化亚型之间的PGA。共享相同字母的组之间没有显著差异（Mann–Whitney U检验q > 0.05），而不同字母的组之间q < 0.05。D、对于3,516个mRNA转录本（已根据“方法”部分调整PGA和队列），它们在不同亚型之间存在显著差异[Kruskal–Wallis秩和检验q < 0.05且亚型中值的标准化差（SD）>0.4]，这些转录本的亚型中值从1排到4（从小到大）。按列显示了给定亚型每种排名的时间百分比。E、使用GSVA（根据“方法”部分调整PGA和队列）为50个标志途径计算了患者特定的途径评分。对于每个途径，热图中显示了亚型中值调整后的mRNA丰度z分数。Kruskal–Wallis检验的q值显示在右侧条形图中，途径按照亚型中值的标准化差降序排列。显示了最重要的20个途径。对于缺氧响应，使用了Buffa基因集。F、对于每个途径，热图显示了PTEN野生型（WT）和PTEN缺失肿瘤中GSVA途径z分数的中值。右侧条形图报告了Wilcoxon秩和检验比较PTEN WT与PTEN缺失的FDR调整后的q值，仅显示q < 0.05的途径。G、在每个甲基化亚型内，PTEN CNA状态与GSVA途径评分之间的关联。点表示效应大小，定义为具有与没有PTEN CNA丢失的肿瘤之间GSVA途径评分的中值差异。背景颜色表示FDR调整后的q值。H、在每个甲基化亚型内，反复出现的CNA驱动事件与GSVA缺氧评分之间的关联。点表示效应大小，定义为具有与没有该改变的肿瘤之间缺氧评分的中值差异。背景颜色表示FDR调整后的q值。

为了确定基因组-表观基因组相互作用是否与下游基因表达相关，我们接着考虑了738个具有匹配的甲基化和RNA分析的前列腺肿瘤。在评估亚型差异时，53,078个mRNA转录本中有40,496个（76%）在不同亚型之间存在显著差异（q < 0.05；调整了PGA和队列差异）。我们选择关注3,516个效应大小较大的显著转录本子集（补充表S6）。这些转录本在不同甲基化亚型之间显示出不同的表达谱（补充图S13A），其中MS-1和MS-2的表达量倾向于更高（图3D）。为了将这些丰度差异放在具体背景下，我们接下来使用基因集变异分析（GSVA）总结了标志途径的活性。由于GSVA产生的富集分数是基于排名的，我们将这些值解释为样本间的相对途径活性而非绝对水平。在途径层面，这些差异反映了广泛的癌症标志途径（图3E），其中MYC靶标在不同亚型之间的差异最大（48, 49）。鉴于在PTEN缺陷肿瘤中观察到的显著CNA-甲基化相互作用（图2F），我们接下来评估了PTEN缺失病例中的亚型特异性途径活性。MS-1中的血管生成途径活性升高，而MS-3和MS-4中没有显著变化。相比之下，MS-4中的缺氧和糖酵解评分显著增加（MS-3中呈现一致但不显著的趋势），而在MS-1中没有显著变化（图3F和G）。鉴于各亚型之间途径的显著差异，我们接下来关注了与前列腺癌进展和进化动态密切相关的缺氧，并且其与预后和治疗反应有明确的关联。我们在四种甲基化亚型中比较了一种经过验证的缺氧RNA特征（53）。MS-2（类似正常的亚型）的缺氧程度显著低于其他亚型（补充图S13B）。定义亚型的CpGs与缺氧的相关性是非定义亚型的CpGs的2.9倍（95%置信区间，2.92–2.95；P < 2.2 × 10?16；补充图S13C）。整个基因组中的CpGs与缺氧大多呈负相关（69%；q < 0.05；补充图S13D）。基因启动子区域（CGIs）的甲基化与1,046个基因中的缺氧评分显著相关（补充图S13E；补充表S6）。最后，我们检查了缺氧评分与每个亚型中反复出现的驱动改变之间的关联（存在于≥20%的肿瘤样本中）。在MS-3中，具有MYC增加的肿瘤显示出显著较低的缺氧评分，而具有NKX3-1丢失的肿瘤显示出显著更高的评分。在MS-4中，缺氧评分随着RB1和PTEN的丢失而显著升高（图3H）。总体而言，不同甲基化亚型之间的途径活性有显著差异，反复出现的驱动改变与缺氧之间的关联因亚型而异，这表明存在亚型特异的进化背景。

前列腺癌是一种C类肿瘤，具有较少的反复出现的单核苷酸变异（SNV），但有许多驱动CNA和基因组重排（GR）。我们评估了870个前列腺肿瘤中全局甲基化与28个反复出现的驱动异常的关联（8），以及SNVs、CNA和GR的突变负担（图4A）。与我们的基因表达分析一致，PGA对甲基化水平的影响最大（87.5%的CpG显著变化，q < 0.05），其次是全局SNV负担、MYC拷贝数增加和GR负担（图4B）。这些CpG变化分布在整个基因组中，影响CGI和非CGI位点（补充图S14A和S14B），并影响了下游mRNA的丰度（补充图S14C）。

体细胞驱动突变塑造了前列腺癌的表观基因组。前列腺癌是一种C类肿瘤，具有较少的反复出现的单核苷酸变异（SNV），但有许多驱动CNA和基因组重排（GR）。我们评估了870个前列腺肿瘤中全局甲基化与28个反复出现的驱动异常的关联，以及SNVs、CNA和GR的突变负担（图4A）。与我们的基因表达分析一致，PGA对甲基化水平的影响最大（87.5%的CpG显著变化，q < 0.05），其次是全局SNV负担、MYC拷贝数增加和GR负担（图4B）。这些CpG变化分布在整个基因组中，影响到CGI和非CGI位点（补充图S14A和S14B），并影响了下游mRNA的丰度（补充图S14C）。D部分：使用逻辑回归来检测驱动事件与亚型之间的关联，同时调整了队列和PGA（基因组复制数异常）的影响。当以MS-2作为参考组时，显示了每个亚型（列）发生特定突变的概率。PGA和SNVs/Mbp（单核苷酸变异/百万碱基对）被中位数二分；显示了12个q值小于0.05的驱动因素。E部分：按照甲基化亚型绘制了驱动事件总数的箱线图，并给出了Kruskal-Wallis P值。使用Mann-Whitney U检验对组间进行了成对比较，不同字母的组之间存在显著差异（q < 0.05）。F部分：最终选定的模型在预测特征方面的表现（y轴标签）。点表示在测试数据集上的估计性能，相邻的线表示95%的自助法置信区间（bootstrap CIs）。对于数值和分类特征分别使用了皮尔逊相关系数（r）和平衡准确率（BA）。G部分：测试数据集在预测log2 SNVs/Mbp方面的表现。皮尔逊相关系数（r）及其95%的自助法置信区间（95% bootstrap CI）以括号形式给出。H部分：测试数据集在预测MYC拷贝数增加方面的性能。平衡准确率（BA）及其95%的自助法置信区间（95% bootstrap CI）以括号形式给出。每个单元格根据该行总数的百分比进行着色。I部分：训练了随机森林（RF，蓝色）和正则化Cox比例风险模型（rCoxPH，红色），使用不同的数据组合作为预测因子来预测至BCR（无并发症进展）的时间[x轴标签；详见(J)部分定义]。条形图的高度表示在训练和测试数据集中的估计C指数。C指数越高，表示性能越好。J部分：图(I)中使用的预测因子组合。对于给定的行，y轴标签表示数据组合，框表示包含的数据类型。K部分：使用(I)中的最佳模型在161名患者的测试数据集中预测BCR风险评分，并将风险分为高和低两类。Kaplan-Meier图显示了高风险和低风险组的BCR-free概率，以及从Cox比例风险模型中估计的HR（风险比）、95%置信区间（95% CI）和Wald检验P值。图下方的表格给出了0至5年内处于风险中的患者数量。L部分：使用(I)中的最佳模型预测52个正常样本、测试数据集中的161个局部前列腺癌样本和99个转移性样本的BCR风险评分。使用Kruskal-Wallis秩和检验P值比较了三种组织组之间的风险评分。不同字母的组具有显著不同的风险评分特征（成对Wilcoxon检验q < 0.05）。图4。

要确定单个基因上的不同体细胞突变是否会引起类似的甲基化效应，我们关注了肿瘤抑制基因TP53。我们比较了三组中的TP53失活情况：仅SNV（6%的患者）、仅拷贝数丢失（24%）以及SNV和拷贝数丢失同时存在（8%）。仅SNV的TP53类别与显著的甲基组变化相关联，其程度远大于TP53拷贝数丢失（56%的CpGs对比20%；图4C）。特定的CpGs部分重叠（补充图S14D），表明不同的突变类型会引起不同的甲基化变化。与这些结果一致，28个驱动突变中有12个与甲基化亚型显著相关（logistic回归调整了队列和PGA的影响），MS-3和MS-4相对于正常亚型MS-2具有更高的突变概率（图4D；补充图S14E和S14F）。同样，MS-4的驱动事件总数更多，而MS-2的驱动事件较少（图4E；补充图S14G；补充表S7），即使在控制了局部癌和转移性癌之间巨大的突变负荷和驱动突变频率差异后也是如此（补充图S15）。同时拥有足够的资源和组织来分析肿瘤DNA和甲基化情况是相对罕见的。鉴于甲基化与癌症的突变组成之间的强烈关联，我们尝试创建预测模型，从而从癌症的表观基因组推断出其体细胞基因组的特征。我们将1,436名患者随机分为70%用于训练（n = 1,005）和30%用于测试（n = 431）。我们关注了15个具有足够数据用于模型训练的重要临床分子特征：患者年龄和血清PSA（连续和离散值）、肿瘤级别和T分期、突变负荷（SNVs和CNAs）、七个肿瘤抑制基因和一个癌基因。我们在训练数据集中拟合了随机森林和glmnet模型，进行了标准特征选择、超参数优化和过拟合评估（“方法”部分；补充图S16和S17；补充表S7）。最终模型在最终的保留测试数据集上表现良好，除了年龄之外的所有特征的性能都超过了0.60（图4F；补充图S18）。例如，我们的甲基化模型非常好地预测了肿瘤SNV的突变负荷，我们的预测与实际DNA测序之间的皮尔逊相关系数为0.72（图4G）。同样，尽管只有34%的癌症具有MYC增加，但我们的MYC增加预测模型准确预测了83%的MYC中性病例和76%的MYC增加病例（图4H）。只有一个模型显示出队列特异性偏差：PSA，反映了局部癌和转移性癌之间的强烈临床差异（补充图S19）。所有预测模型都可以在PrCaMethy R包中找到（补充图S20A），我们定义了一组约7,000个在预测模型中重叠的CpGs用于靶向检测（补充图S20B）。

为了将这些临床分子甲基化特征扩展到患者风险分层，我们关注了774名具有至BCR时间信息和多组学数据的局部前列腺癌患者。我们再次将这些患者随机分为70%用于训练和30%用于测试，并使用了标准的模型和特征选择方法（“方法”部分；图4I和J；补充图S21A和S21B）。训练数据集中的最佳模型被转移到验证数据集中，并包括了两个临床特征（PSA和ISUP GG）、甲基化亚型、10个启动子甲基化事件以及七个基因的RNA丰度（补充图S21C）。该模型在测试数据集中的风险评分C指数为0.86（95%置信区间，0.77–0.91），而仅基于临床特征的模型的风险评分为0.81（95%置信区间，差异为0.0011–0.11；补充图S21D）。中位数二分的风险评分在保留的测试数据集中识别出一组BCR风险增加8.2倍的患者（95%置信区间，1.9–36.1；P = 0.0055；图4K）。模型风险评分在非恶性、局部和转移性前列腺样本之间有显著差异（图4L），即使在被分类为中等风险的保留测试患者子集中也具有预后意义（HR = 16.23；P = 0.013；补充图S21E）。总体而言，这些结果表明前列腺癌的表观基因组及其与转录组的相互作用可以作为非常准确的临床结果预测器。

在肿瘤甲基化中，我们无法准确预测的唯一临床分子特征是患者年龄（图4F）。这提示我们前列腺肿瘤的表观遗传年龄可能与患者的实际年龄有显著差异。我们使用10个广泛可用的甲基化时钟（55个）估计了102个正常样本的DNAm年龄。Horvath时钟获得了最高的相关性（皮尔逊相关系数r = 0.48；95%置信区间，0.31–0.61；q值 = 5.18 × 10^-6；图5A和B），与原始研究的结果一致（56, 57）。然后我们使用Horvath时钟分析了配对的正常和肿瘤样本。肿瘤样本的DNAm年龄比其正常对应样本显示出减缓趋势，表明肿瘤中的衰老路径与相邻正常组织中的衰老路径不同（图5C）。图5。

体细胞驱动突变影响前列腺癌的衰老过程。A部分：10个DNAm时钟的性能。皮尔逊相关系数（r）及其95%置信区间被报告。右侧的条形图中显示了经FDR调整后的q值。B部分：10个DNAm时钟之间的成对相关性（皮尔逊相关系数r）。C部分：正常和肿瘤样本之间的DNAm年龄的成对比较。虚线表示拟合的线性模型（斜率 = 0.53；截距 = 23.1）。D部分和E部分：使用Kruskal-Wallis秩和检验比较了不同甲基化亚型（D）和ISUP组（E）之间的相对年龄加速，同时调整了实际年龄。F部分和G部分：基于MYC拷贝数增加（F）和NKX3-1拷贝数丢失（G），使用Mann-Whitney U检验比较了相对年龄加速，同时调整了实际年龄。H，ICGC PRAD-CA、TCGA及合并队列中CNA驱动因素与相对年龄加速（调整后考虑了实际年龄）之间关联的总结。点代表效应大小，背景颜色表示q值。查看大图下载幻灯片。

体细胞驱动基因突变影响前列腺癌的衰老过程。A，10种DNA甲基化时钟的表现。帕尔逊相关系数（r）连同95%置信区间一同报告。右边的条形图中显示了经FDR校正后的q值。B，10种DNA甲基化时钟之间的成对相关性（帕尔逊r）。C，正常组织和肿瘤样本间DNA甲基化年龄的成对比较。虚线代表拟合的线性模型（斜率=0.53；截距=23.1）。D和E，使用Kruskal–Wallis秩和检验，基于经实际年龄校正的残差，比较了甲基化亚型和ISUP组间的相对年龄加速。F和G，使用Mann–Whitney U检验，基于经实际年龄校正的残差，比较了MYC CNA扩增（F）和NKX3-1 CNA缺失（G）与相对年龄加速的关系。H，ICGC PRAD-CA、TCGA及合并队列中CNA驱动因素与相对年龄加速（调整后考虑了实际年龄）之间关联的总结。点代表效应大小，背景颜色表示q值。关闭图5。

接下来，我们研究了临床分子特征与表观遗传年龄加速之间的关联，同时调整了实际年龄（“方法”部分）。甲基化亚型（P = 1.32 × 10^-6）和ISUP分级（P = 8.43 × 10^-3）均显示出与相对年龄加速的显著关联（图5D和E）。总体突变负担与相对DNA甲基化加速呈负相关（SNVs/Mbp：Spearman’s ρ = ?0.30；P = 2.82 × 10^-16；PGA：Spearman’s ρ = ?0.26；P = 2.53 × 10^-13），这与其他癌症中观察到的SNV负担与DNA甲基化年龄之间的关联一致（56, 57）。两个特定的驱动因素与DNA甲基化年龄的减缓有关：MYC拷贝数扩增（P = 2.47 × 10^-5）和NKX3-1丢失（P = 8 × 10^-4；图5F–H）。总体而言，这些数据表明前列腺肿瘤表现出年龄缩短的特征，并且这与临床表型、总体突变负担以及特定的体细胞驱动基因突变相关。

癌症是一种持续进行基因组演化的疾病。表观基因组在体细胞突变及其基因表达后果之间起着关键的中介作用，最终决定了诊断、预后和治疗反应等临床表型。大量证据表明，各种形式的表观基因组（失调）调控，如染色质修饰，可以驱动前列腺癌的表型变化（58–60）。在这项研究中，我们关注DNA甲基化，以最大化样本量，并推断其在疾病演化过程中如何变化以及如何影响表型。我们发现前列腺癌可以分为四种主要的甲基化亚型。出乎意料的是，每种亚型都可以从正常前列腺到转移性疾病的各个阶段都被识别出来，尽管频率不同（图6）。这表明甲基化亚型会随时间演变，但至少部分反映了肿瘤的早期启动状态。与此一致的是，这些亚型与CIMP和特定的体细胞突变、转录和临床特征相关。为了将这些甲基化亚型与我们之前基于全基因组测序（WGS）定义的整合分子亚型（IMS）进行对比，我们比较了重叠样本中的MS和IMS分类结果。MS亚型分布于多个IMS组中，两种分类之间存在关联，但一致性有限（Qannari检验P < 0.001；ARI = 0.09；95%置信区间，0.07–0.11；补充图S22）。这些数据表明，MS和IMS捕捉到了前列腺癌异质性的相关但不等同的方面。图6。查看大图下载幻灯片。

不同疾病状态下，甲基化亚型的分布各不相同。图6提供了疾病进展过程中甲基化亚型、临床事件和复发驱动事件分布的概览。原发肿瘤根据ISUP分级被分为低级别（GG1）、中级别（GG2–GG3）和高级别（GG4–GG5），转移性肿瘤单独显示。对于每种疾病状态，展示了至少10%的肿瘤中出现的MS亚型比例、5年无复发生存率（BCR）和复发驱动事件的比例。转移性肿瘤被标记为统一事件阳性。正常、局限性和转移性病例的总体MS亚型分布显示在图1E中。作为对比，我们之前研究中IMS（IMS-1–IMS-7）的比例也显示在底部。顶部插图总结了转移性肿瘤与局限性肿瘤之间DNA甲基化–CNA关联的变化，表示为Δ中位数ρ（转移性 vs. 局限性；约0.22）。关闭图6。

尽管全基因组测序（WGS）能够推断肿瘤亚克隆结构（10, 16），但将基于甲基化的亚型划分扩展到亚克隆水平仍然具有挑战性。来自Illumina 450K/850K芯片的批量甲基化谱代表了表观遗传状态的混合体，虽然全基因组亚硫酸盐测序（WGBS）提供了读段级别的信息，但其典型覆盖范围限制了低频亚克隆的可靠检测。因此，理解前列腺癌甲基组的单细胞异质性（及其亚型）对于进一步阐明DNA甲基化如何影响癌症演化是一个重要机会。已知最显著的前列腺肿瘤演化驱动因素是拷贝数改变（CNA），这些改变与甲基化变化密切相关（8, 13, 61）。我们发现这些变化在个体之间以及在不同疾病阶段之间存在显著差异。PTEN肿瘤抑制基因是一个很好的例子，其高甲基化和半合子缺失都可能导致基因表达下降，伴随特征性的DNA甲基化谱型（61），这可能导致PTEN单倍体不足的功能变化（62–64）。我们确定了30个与无复发生存（BCR）显著相关的HVM基因。这些与BCR相关的HVM基因没有显示出通路富集，但确实显示出强烈的CNA–甲基化相关性，突显了甲基化对基因调控和临床结果的复杂性和异质性。先前的研究显示，前列腺癌中的DNA甲基化异质性与拷贝数演化平行，且CNA–甲基化关系高度依赖于具体背景（13）。远端顺式调控元素上的甲基化可以影响长距离基因表达，简单的启动子级别CpG分配可能会稀释连贯的通路级信号（65, 66）。在与BCR相关的HVM基因中，有几个与前列腺癌生物学有已知联系，包括HNF1B（67, 68）、LDHB（69）和MT1E（70），这表明这种方法可能对其他类型的癌症也有价值。使用GSVA总结标志性通路活性（基于排名的富集，反映相对水平而非绝对水平），我们发现PTEN缺陷肿瘤内的通路活性在不同甲基化亚型之间有显著差异。在MS-1中，血管生成通路活性最为明显，这与证据一致，即血管生成程序可能在前列腺肿瘤发生早期就被诱导，并可能在PTEN缺失的情况下通过PI3K/AKT信号通路得到增强（71–73）。相比之下，MS-4中的糖酵解活性显著升高，这与报告相符，即糖酵解代谢和乳酸产生在晚期前列腺癌中更为明显，并且可以由PI3K/AKT信号通路促进（74, 75）。我们还观察到MS亚型与复发驱动改变和Buffa缺氧特征活性之间的特定关联。我们小组之前的工作显示，MYC扩增与较高的缺氧水平相关，使用的是中位数二值化的缺氧评分（基因的上半部分=+1，下半部分=-1，跨基因求和；参考文献50）。然而，在基于排名的框架中，MYC扩增与MS-3中的较低相对缺氧信号相关，表明这种转录程序的平衡发生了特定于该甲基化背景的变化。这些结果共同强调了复发驱动改变与通路特征活性之间的亚型特异性相互作用，这与亚型特定的演化和微环境限制一致。我们扩展了之前的相关性研究（76），定义了关键前列腺癌驱动事件的精确表观基因组特征，并证明这些特征是可以高度预测的。这些发现支持了甲基化捕捉稳定肿瘤状态的观点，这些状态反映了谱系程序，并受到基因组破坏和微环境压力的影响。随着循环肿瘤DNA的甲基化谱型分析变得越来越可行（77, 78），这为开发非侵入性技术提供了途径，以纵向估计个体驱动事件和整体突变特征。此外，由于不同甲基化亚型之间的驱动–通路耦合存在差异，将甲基化谱型分析与临床分子特征的预测相结合，可能有助于支持具有亚型意识的治疗假设和分层策略。最后，我们预计下一阶段的验证将在 unfavorable intermediate-risk 和 high-risk 局限性前列腺癌中提供更多信息，而不是在通过主动监测管理的低风险疾病中，因为改进的分子分层可以支持最佳治疗策略并指导研究设计。我们的DNA甲基化年龄分析再现了DNA甲基化年龄的减缓（56, 57），并发现了与临床分子特征的关联，包括ISUP分级和复发驱动事件，表明“更年轻”的DNA甲基化年龄并非纯随机现象。Horvath时钟是在非恶性组织上训练的，其在恶性甲基组中的时钟CpG可能受到直接影响。因此，明显的时钟逆转可能反映了肿瘤对这些位点的重塑，而不是真正的恢复活力。从机制上讲，这与细胞分化受损（NKX3-1拷贝数丢失）和/或诱导类似胚胎干细胞的状态（MYC拷贝数扩增；参考文献48, 79）一致。尽管这项研究为前列腺肿瘤甲基组景观提供了新的见解，但这些发现应考虑到几个限制。首先，有两个队列[Li及其同事（30）和Zhao及其同事（15）使用的是WGBS进行分析，而大多数其他队列使用的是Illumina甲基化芯片（“方法”；补充表S3）。为了减轻平台相关的影响，我们通过将分析限制在Illumina 450K芯片共有的CpG上实现了数据的一致性。重要的是，即使调整了队列和甲基化分析平台后，使用独立定义的CIMP评分进行的亚型–CIMP一致性仍然稳健。我们研究的另一个限制是，并非所有队列都提供了原始IDAT文件，因此无法在数据集中统一处理强度级别。TCGA亚型发现队列使用了SeSAMe处理的数据，这可以减少假阳性甲基化测量（包括与拷贝数丢失区域相关的伪影）。其他队列使用了一系列不同的处理流程。为了 mitigate 这一问题，我们在下游模型中加入了全局CNA负担（即PGA）的测量，以减少拷贝数相关甲基化偏倚对我们发现的影响。然而，PGA是一个全局性指标，并不能替代针对特定位点的CNA校正。有多种方法可用于表观遗传亚型发现，包括离散聚类、整合多组学模型（80, 81）和考虑混杂因素的去卷积方法（82, 83）。在这项研究中，我们应用了基于轮廓的共识聚类方法，因为它产生了可重复的、适用于分类器的亚型标签，新队列加入时可以轻松扩展。这类方法有两个一般性限制。首先，它使用单一模态，而已经提出了结合多个组学层和分子数据的整合亚型框架。实际上，在多组学丰富性和保持足够样本量以严格评估普遍性之间存在权衡。其次，我们的分析依赖于批量肿瘤谱型，这些谱型平均了恶性组分和微环境组分，压缩了肿瘤内部的异质性。这些考虑因素促使未来使用单细胞来明确建模混杂情况，并将批量亚型分配与恶性和微环境细胞状态的混合物相关联。除了这些技术限制外，我们的转移性队列规模较小，治疗注释也有限，这限制了我们完全表征在治疗压力下可能出现的表观遗传亚型的能力。更大规模的转移性数据集，如果具有详细的治疗历史，并且理想情况下是纵向采样的，对于更全面地定义治疗诱导的表观基因组状态至关重要。由于这个汇编整合了具有不同临床组成和技术特征的队列，肿瘤纯度和遗传背景可能是潜在的混淆来源。为了评估我们发现的稳健性，我们进行了广泛的敏感性分析和队列调整分析（“方法”；补充信息S1）。我们的主要结论保持不变，尽管对纯度估计算法进行进一步的系统基准测试，并在各个队列中协调遗传祖先注释将非常有价值。总体而言，我们的数据阐明了DNA甲基化（DNAm）在前列腺癌整个生命周期中的核心作用，从早期老化减缓开始，到局部疾病中的广泛基因调控效应，以及作为复发和转移等临床重要事件的强大生物标志物。尽管我们的数据集规模较大，但它缺乏对未治疗病例和淋巴结转移的覆盖，而这些信息可能有助于填补我们对前列腺癌从局部病变发展成致命转移性病变机制理解中的重要空白。

**患者队列**

我们使用了TCGA-PRAD队列（n = 498个肿瘤样本；参考文献7）来鉴定甲基化亚型，并在14个外部前列腺癌队列中验证了这些亚型（13个已发表，1个未发表——墨尔本；图1B–E；补充图S2；补充表S1和S2；参考文献7, 8, 13, 15, 25, 30–38）。TCGA、国际癌症基因组联盟（ICGC）PRAD-CA（加拿大）以及赵及其同事（15）的前列腺癌队列被纳入了甲基化、mRNA、拷贝数异常（CNA）和体细胞驱动突变分析。有关每个队列的样本量和所用数据类型的详细信息，请参见补充图S1A。所有ICGC样本均来自患有局部病变的患者的前列腺组织，并遵循大学健康网络研究伦理委员会批准的研究方案（UHN 06-0822-CE、UHN 11-0024-CE和CHUQ 2012-913:H12-03-192）进行了处理。

**病理学和组织处理**

ICGC PRAD-CA的新鲜冷冻组织样本由两位独立的泌尿生殖病理学家检查肿瘤含量、Gleason/ISUP分级和肿瘤细胞占比，如先前报道的那样（8）。至少含有70%肿瘤细胞的样本被手动大解剖并用于核酸提取，同样依照先前的方法（8, 9）。

**微阵列数据生成**

ICGC PRAD-CA的甲基化、SNP和mRNA微阵列数据生成方法如先前所述（8）。简要来说，使用Illumina Infinium HumanMethylation450 BeadChip试剂盒在麦吉尔大学和Genome Quebec创新中心评估了全局甲基化水平，输入的基因组DNA量为500 ng。SNP微阵列使用Affymetrix OncoScan FFPE Express 2.0和3.0芯片进行检测，所需DNA量为200 ng。在DNA量有限的情况下（17个样本），我们采用了全基因组扩增（WGA；WGA2，Sigma-Aldrich），并验证了WGA不会显著改变CNA谱型（9）。总RNA使用mirVana miRNA分离试剂盒（Life Technologies）从交替的相邻切片中提取（以减少相对于DNA样本的潜在异质性）。提取的总RNA（100–150 ng）在应用基因组学中心（The Hospital for Sick Children）和伦敦区域基因组学中心（Robarts Research Institute）的Affymetrix Human Gene 2.0 ST（HuGene 2.0 ST）和Human Transcriptome Array 2.0（HTA 2.0）平台上进行了检测。

**甲基化分析**

甲基化处理是在R统计环境中进行的（版本v3.4.0）。ICGC PRAD-CA队列的甲基化数据存储在Gene Expression Omnibus: GSE107298中，并使用Dasen方法（84, 85）进行了处理，如先前所述（8）。排除那些检测灵敏度低于背景噪声的探针、非CpG甲基化、其标记的次要等位基因频率超过5%且与最近的CpG位点距离最大为两个核苷酸的SNP，以及已知会在多个基因组位置杂交的SNP。探针的平均强度值是通过对技术重复样本进行计算得出的。探针通过IlluminaHumanMethylation450kanno.ilmn12.hg19注释包（版本v0.6.0）映射到基因符号。除非另有说明，所有甲基化分析均使用β值进行，该值代表甲基化强度与未甲基化强度的比率，反映了特定CpG位点的甲基化水平。肿瘤纯度是根据白细胞的去甲基化程度（LUMP）来估算的，方法如先前所述（86）。总体而言，LUMP估算的纯度较高（Q1 0.83；中位数0.88；Q3 0.94），没有样本的纯度低于0.3，并且没有设定特定的肿瘤纯度阈值。使用了10种算法对纯度进行了额外估算，以进行敏感性分析（补充信息S1；参考文献86–93）。以下队列的CpG–处理β值是通过Illumina Infinium Human Methylation 450K或850K BeadChip获得的：Kirby及其同事（31）、Zhao及其同事（25）、TCGA-PRAD（7, 94）、ICGC PRAD-FR（32）、Ramakrishnan及其同事（34）、Ylitalo及其同事（36）和Mundbjerg及其同事（38）。对于用于亚型发现的TCGA-PRAD队列（498个前列腺肿瘤和50个正常前列腺样本），我们使用了Genomic Data Commons数据发布34中的处理过的DNA甲基化数据集，该数据集是通过SeSAMe流程生成的（95）。有关数据处理的更多细节，请参阅相应的参考文献。在所有分析中，仅使用了Illumina 450K芯片中的CpG信息。对于墨尔本的甲基化数据，使用了Illumina的Infinium Human Methylation 450K芯片来分析所有澳大利亚前列腺癌样本的甲基化状态。原始信号数据使用minfi包（96）进行了处理和标准化，以消除芯片内和芯片间的变异。minfi实现了基于Illumina参考的标准化和Subset-quantile Within Array Normalization方法，分别用于芯片间和芯片内的标准化（校正探针类型偏差）。任何检测P值大于0.01的甲基化测量结果都被从后续分析中剔除。对于Brockses及其同事（13）、Tothes及其同事（37）和Dairo及其同事（35）的队列，原始的Illumina甲基化IDAT文件使用minfi R包（版本v1.52.1；参考文献96）进行了处理。检测P值大于0.05的甲基化值被排除。在样本中检测失败率超过10%的探针也被移除。此外，任何有超过10%的探针未能检测到的样本也被从后续分析中剔除。标准化使用wateRmelon R包（版本v2.12；参考文献84, 85）中的dasen方法完成。Li及其同事和Zhao及其同事的研究中使用了WGBS甲基化数据（15, 30）。Rockermeth R包（版本v.0.1.1；参考文献97）中的“illumina450k_hg38”数据集用于识别与HumanMethylation450 BeadChip对应的CpG位置，从而使我们能够使用与其他队列相同的CpG信息。在Zhao及其同事（15）的研究中，读取次数少于10次的CpG和患者级别的CpG调用次数少于10次的CpG被排除在分析之外，而在Li及其同事（30）的研究中则使用了5次的截止标准，因为覆盖率较低。对于每个CpG，β值是根据该CpG位置上的甲基化读取次数计算得出的：#甲基化读取次数 / (#甲基化读取次数 + #未甲基化读取次数)。

**遗传祖先**

我们已经为ICGC PRAD-CA（16）和TCGA-PRAD（98）队列获得了公开的遗传祖先群体分类。对于Zhao及其同事（15）的队列，我们按以下方式估算遗传祖先：根据1000 Genomes（1KG）参考谱系的超大陆祖先群体（欧洲、非洲、东亚、南亚和混合美洲）推断出分类祖先：首先，将基因型与1KG参考队列（第三阶段主要发布版本，基因组构建GRCh38，仅包含创始人）合并，然后使用bcftools提取相交的SNP集。合并后的基因型转换为PLINK2（100）格式，并过滤出缺失值最大为0.05、等位基因频率最小为0.01的样本，然后修剪掉连锁不平衡较高的变异。使用PLINK2的–pca标志计算基因型的主成分（PC）。将K最近邻居（KNN）应用于前五个PC，将样本分配到1KG注释的祖先簇中。连续祖先定义为来自1KG群体的混合比例。使用ADMIXTURE（101）工具在无监督模式下估计每个样本合并后SNP集中的混合比例。超参数K表示预期的群体数量，设置为5，对应于五个1KG超群体。根据每个已知参考样本集中贡献最大混合比例的群体标签对推断出的簇进行标记。有五个患者的KNN PC–基于的祖先分类与ADMIXTURE基于的分类不一致，这些患者被标记为“混合”并从涉及遗传祖先的进一步分析中剔除。使用Cramér的V和ε2效应量分别量化了SIR与遗传祖先之间的关联以及这些特征与关键临床分子特征之间的关联。

**拷贝数异常（CNA）**

ICGC PRAD-CA的CNA数据获取方法如先前所述（8）。Affymetrix OncoScan FFPE Express 3.0 SNP探针检测的结果使用OncoScan Console v1.1与自定义参考序列对齐（8）。随后使用BioDiscovery的Nexus Express for OncoScan软件和SNP-FASST2算法调用拷贝数片段。参数与默认设置相同，除了每个片段的最小探针数量设置为20个。在需要时使用Nexus Express软件进行二倍体重新定位。使用BEDTools，根据refGene注释（2014年7月15日；参考文献102）识别基因水平的CNA。PGA定义为PGA =（所有具有增益或丢失的基因组区域的碱基对长度）/ 32亿个碱基 × 100。Zhao及其同事（15）的CNA数据来自（103），TCGA CNA数据从TCGAAbiolinks R包v2.25.3（7, 94）下载，其中拷贝数增益为“中性”的视为2，小于2的视为“丢失”，大于2的视为“增益”。

**mRNA**

ICGC PRAD-CA的预处理mRNA数据从GSE107299下载并在R（版本v3.4.3）中进行了处理。简要来说，使用了affy（v1.48.0）包中的Robust Multi-Array Average算法，并修改了Entrez Gene Chip Definition Files（CDF）来预处理原始强度数据。探针使用自定义CDF（v18）映射到HTA 2.0和HuGene 2.0 ST芯片的Entrez Gene ID。sva包（v3.24.4）用于校正不同芯片之间的批量效应（106）。基于Entrez Gene ID合并了HuGene 2.0 ST和HTA 2.0的mRNA丰度水平。然后在重复的Entrez Gene ID之间平均mRNA丰度。接着根据UCSC Table Browser中的人类参考基因组GRCh37将Entrez Gene ID转换为基因符号和染色体位置（下载日期：2018年1月30日）。TCGA处理的mRNA丰度水平从TCGAAbiolinks R包v2.25.3（7, 94）下载。Zhao及其同事（15）的FASTQ文件使用STAR（v2.7.6a；参考文献107）与GRCh38对齐，并使用RSEM（v1.3.3；参考文献108）计算mRNA丰度。除非另有说明，所有mRNA丰度水平均以log2转录本每百万（log2 TPM）的形式进行分析。

**WGS**

ICGC PRAD-CA队列的数据处理和变异调用方法在（16）中有详细描述。使用ShatterProof（v0.14，默认设置；参考文献109）生成了染色体断裂评分。Kataegis使用固定宽度的滑动窗口来测试观察到的SNV三核苷酸含量和突变距离与预期值的偏差，使用精确的二项式检验（https://www.biorxiv.org/content/10.1101/287839v2）。Zhao及其同事队列（15, 103）的相应WGS数据作为BAM文件下载（参见“数据可用性”部分）。BAM文件被转换为FASTQ格式并进行重新对齐，然后使用Metapipeline-DNA（v5.0.0-rc.4 - 5.3.1；参考文献110）调用变异。通过BCFtools（v1.17；参考文献111）保留了在两个或多个调用工具中都存在的SNV，创建了共识SNV调用集。然后使用allow和deny列表过滤SNV。使用的deny列表包括AccuSNP、Complete Genomics、dbSNP（v156；参考文献112）、Duplicate Gene（v68；参考文献113）、Encode Data Analysis Center和Duke（114）、gnomAD（v3.1.2，带有“PASS”标志的变异；参考文献115）以及1000 Genome项目（99, 116）。allow列表包括Catalogue of Somatic Mutations in Cancer（v92；参考文献117）中的编码和非编码变异。使用Manta（v1.6.0；参考文献118）调用体细胞结构变异，并使用StructuralVariantAnnotation R包（v1.18.0；参考文献119）进行注释。

**一般统计考虑**

使用Benjamini–Hochberg方法计算多重检验调整后的P值，以控制假发现率（FDR）在0.05的水平（在整个文中称为“q”或“q值”；参考文献120）。在对大于2 × 2的列联表进行Fisher精确检验时，使用了stats R包的fisher.test函数，并进行了100,000次模拟复制。

**甲基化亚型**

在TCGA-PRAD队列（n = 498）中，使用Interpretable MiniPatch Adaptive Consensus Clustering（IMPACC v0.1.0；参考文献121）和默认设置、10,000个聚类划分以及CpG β值识别甲基化亚型。过滤掉前50,000个方差最大的CpG用于聚类。使用silhouette方法（122）选择最大化中位轮廓分的亚型数量，最终得到四个亚型。具体来说，亚型的总数从2个变化到6个，每次都会计算每个簇的中位数轮廓得分，以簇得分的中位数作为聚类准确性的整体指标。训练了一个最近邻收缩质心分类器（123）来对TCGA亚型进行分类，该分类器用于对前七个队列的样本进行分类，缺失值通过KNN插补法进行填补（使用impute R包v1.74.1，默认设置；参考文献124）。随着更多队列纳入项目，其中一些队列使用了Illumina 850K甲基化阵列进行分析，该阵列缺失了大约9%用于确定亚型的Illumina 450K阵列上的CpG位点。为了解决由此产生的缺失数据问题，训练了一个随机森林分类器（randomForestSRC R包v3.3.3）来处理TCGA亚型。选择这种方法是因为它具有内置的插补功能，即使CpG覆盖不完全，也能实现对剩余七个队列样本的稳健分类（125）。与至少一个CpG显著相关的驱动事件（q < 0.05）被汇总，以计算每位患者的显著驱动事件总数，然后使用Kruskal–Wallis检验和成对Mann–Whitney U检验在不同亚型之间进行比较。最后，比较了局部前列腺癌和转移性前列腺癌中每个驱动突变的分布。对于连续型驱动因素，计算了Mann–Whitney U检验的P值；对于分类型驱动因素，计算了Fisher精确检验的q值和OR值，并对总突变负担进行了调整（CNA驱动因素根据PGA进行调整，SNVs根据log2 SNVs/Mbp进行调整，GRs根据log2 GR计数进行调整）。表观遗传学预测因子与临床特征和驱动突变的关系。

共有1,436名患者被随机分为70%的训练数据集（n = 1,005）和30%的测试数据集（n = 431）。比较了两种机器学习模型在预测临床特征和驱动突变方面的表现：随机森林和正则化通用线性模型（glmnet；参考文献135, 136）。选择这两种模型是因为随机森林和正则化线性模型在多种应用中都表现出色（参考文献137, 138）。随机森林可以自动捕捉潜在的交互作用或非线性效应（参考文献139），但相比glmnet更难以解释。相比之下，glmnet默认假设线性且加性的效应（无交互作用），其解释方式与传统逻辑或多项式分类模型相似（估计每个预测因子的log OR值），但通过将系数缩小到零来自动选择特征，试图移除不重要的特征。随机森林使用randomForest R包（v4.7.1.1）进行训练，mtry参数通过caret R包（v6.0.94；参考文献140）进行调整，并通过五折交叉验证选择最佳值。glmnet使用glmnet R包（v4.1.8；参考文献135）进行训练，α参数设置为0.5（即“弹性网络”惩罚，介于ridge和LASSO之间；参考文献141, 142），主要惩罚参数lambda通过五折交叉验证确定。缺失值使用KNN插补法（impute R包v1.74.1，默认设置；参考文献124）进行填充。连续型特征的预测准确性通过皮尔逊相关性进行评估，分类型特征的预测准确性通过平衡准确性[BA = (敏感性 + 特异性)/2]进行评估。对于有三个或更多类别的分类型特征，BA分别对每个类别使用“一对所有”方法进行计算，然后取平均值。预测因子包括13,603个基因的前列腺癌基因水平甲基化，定义为基因 promoter区域中CGI的中位数β值。在训练数据集中进行了特征选择，以尝试提高性能并识别出更小的重要预测因子子集。具体来说，使用Kruskal–Wallis秩和检验评估每个基因甲基化预测因子与每个临床特征和驱动突变的单变量关联，对于连续型特征使用Spearman相关性检验。为了考虑队列差异，使用线性回归对每个基因与队列进行回归，并将残差用于上述单变量检验。然后，使用前10个、50个、100个、1,000个或5,000个最显著的基因以及考虑所有1,3603个基因的模型进行建模。为了鼓励更简洁的模型，选择预测给定特征的最终模型是预测因子数量最少的模型，其皮尔逊相关性或平衡准确性分别与最佳值相差不超过0.02。最后，使用测试数据集，通过替换重采样患者，并使用boot R包（v1.3.28.1；参考文献143）计算95%的自举校正和加速（BCa）置信区间（CI）来评估性能指标。

共有774名患有局部前列腺癌的患者，这些患者具有BCR时间数据和多组学数据（甲基化、RNA和CNA），并被随机分为70%的训练数据集和30%的测试数据集。比较了两种机器学习模型在预测BCR时间方面的表现：随机生存森林（random survival forest；参考文献144）和正则化Cox比例风险模型（regularized Cox proportional hazards model；参考文献145），预测准确性通过C指数（C-index；参考文献146）进行评估。随机森林使用randomForestSRC R包（v3.2.30；参考文献144）和默认设置进行训练，而正则化Cox模型使用glmnet R包（v4.1.8；参考文献145）进行训练，α参数设置为0.5（即“弹性网络”惩罚，介于ridge和LASSO之间；参考文献141, 142），主要惩罚参数lambda通过十折交叉验证确定，以在训练数据中最大化C指数。缺失值使用KNN插补法进行填充（impute R包v1.74.1，默认设置；参考文献124）。考虑了几种数据类型的组合进行模型训练（补充图S21A）。由于潜在特征数量众多，在训练数据集中进行了特征选择，以识别出更小的潜在预测因子子集。为了进行特征选择，分别对每个特征作为BCR时间的预测因子进行Cox模型训练，同时调整预处理后的临床特征。然后只保留Wald检验P值< α的特征。考虑了几个α值（0.001、0.005、0.01、0.05和1），最终选择α = 0.005以在训练和测试数据集中最大化随机森林的C指数（补充图S21D）。使用测试数据集，通过对患者进行替换重采样，并使用boot R包（v1.3.28.1；参考文献143）计算两种模型C指数差异的95%自举BCa CI。使用多组学数据预测BCR时间

在患有局部前列腺癌的774名患者中，这些患者具有BCR时间和多组学数据（甲基化、RNA和CNA），被随机分为70%的训练数据集和30%的测试数据集。比较了两种机器学习模型在预测BCR时间方面的表现：随机生存森林（random survival forest；参考文献144）和正则化Cox比例风险模型（regularized Cox proportional hazards model；参考文献145），预测准确性通过C指数（C-index；参考文献146）进行评估。随机森林使用randomForestSRC R包（v3.2.30；参考文献144）和1,000棵树以及默认设置进行训练，而正则化Cox模型使用glmnet R包（v4.1.8；参考文献145）进行训练，α参数设置为0.5（即“弹性网络”惩罚，介于ridge和LASSO之间；参考文献141, 142），主要惩罚参数lambda通过十折交叉验证确定，以在训练数据中最大化C指数。缺失值使用KNN插补法进行填充（impute R包v1.74.1，默认设置；参考文献124）。基于大量潜在特征，进行了特征选择。为了进行特征选择，对每个特征分别作为BCR时间的预测因子进行Cox模型训练，同时调整预处理后的临床特征。然后只保留Wald检验P值< α的特征。考虑了几个α值（0.001、0.005、0.01、0.05和1），最终选择α = 0.005以在训练和测试数据集中最大化随机森林的C指数（补充图S21D）。使用测试数据集，通过对患者进行替换重采样，并使用boot R包（v1.3.28.1；参考文献143）计算两种模型C指数差异的95%自举BCa CI。类似地，为表现最好的模型计算了95%自举BCa CI。最后，将最佳模型的预测BCR风险评分与52个非恶性TCGA样本、测试数据集中的161个局部前列腺癌样本和99个转移性样本进行了比较。鉴于最佳模型包含了PSA和ISUP等级作为预测因子，而这些指标在非恶性样本中未测量，因此使用randomForestSRC impute.rfsrc R函数（默认设置；参考文献125）对这些指标进行了插补。

DNAm年龄是使用methylclock R包（v1.4.0；参考文献55）计算的，适用于具有连续时间年龄数据的队列[TCGA、ICGC PRAD-CA和Kirby及其同事（31）]。排除了WGBS数据集，因为甲基化时钟是基于微阵列数据训练的。相对DNAm年龄加速定义为在假设癌症老年患者线性衰老的情况下捕捉与年龄相关的比例变化：Relative DNAm age acceleration = DNAm age - Chronologic age。

相对DNAm年龄加速在不同年龄类别、甲基化亚型和ISUP等级之间使用Kruskal–Wallis秩和检验进行比较，并考虑了ε2效应大小。在≥20%的肿瘤样本中存在的驱动突变通过Mann–Whitney U检验与相对DNAm年龄加速进行了分析。计算了FDR校正后的q值。为了考虑时间年龄，应用线性回归对相对DNAm年龄加速与时间年龄进行回归，并将残差用于上述单变量检验。

PrCaMethy R包提供了使用Illumina 450K和850K阵列生成的前列腺癌DNAm数据预测甲基化亚型、临床特征和体细胞驱动突变的工具。该R包可在GitHub上公开获取：https://github.com/TheBoutrosLab/package-PrCaMethy，旧版本可在https://uclahs-cds.github.io/package-PrCaMethy/获取。

数据可视化是在R统计环境（v4.3.3 - 4.5.0）中使用BPG包（v7.1.2；参考文献147）完成的。Venn图是使用VennDiagram包（v1.7.3；参考文献148）制作的。

每个队列的数据下载来源在补充表S3中给出。

N. Zeltser在研究期间报告了来自NCI的资助。R. Agrawal在研究期间报告了来自国家普通医学科学研究所（NIGMS）的资助。M. Fraser报告了来自AstraZeneca的额外支持。A.U. Kishan报告了来自Varian Medical Systems、Novartis和ViewRay Systems的资助，以及来自Janssen Pharmaceuticals的资助和其他支持；此外，A.U. Kishan报告了在MiraDx中的少量股权。M.L.K. Chua在研究期间报告了来自新加坡国家医学研究委员会、Duke-NUS医学院和新加坡国家癌症中心的资助，以及来自ArteraAI、Bayer、Decipher Biosciences和Veracyte的资助，还有来自BeOne Medicines、Astellas Pharma、Bayer、IQVIA、Janssen Pharmaceuticals、Novartis、Seagen、Top Alliance和Varian的其他支持。P.C. Boutros在研究期间报告了来自Prostate Cancer Canada、安大略省癌症研究机构、Movember、Terry Fox研究机构、加拿大卫生研究院、加拿大癌症协会、NIH/NCI、Prostate Cancer Foundation和国防部前列腺癌研究计划的资助，以及来自Intersect Diagnostics、Sage Bionetworks和BioSymetrics的其他支持；此外，P.C. Boutros还拥有待审批的临时专利。其他作者没有报告披露信息。

作者的贡献：

J. Arbet：概念化、数据整理、软件、正式分析、可视化、方法学、写作-初稿、项目管理、写作-审阅和编辑。T.N. Yamaguchi：概念化、资源获取、数据整理、正式分析、可视化、方法学、写作-初稿、写作-审阅和编辑。Y.-J. Shiah：概念化、资源获取、数据整理、软件、正式分析、可视化、方法学、写作-初稿、写作-审阅和编辑。R. Hugh-White：资源获取、数据整理、软件、写作-审阅和编辑。A. Wiggins：正式分析、可视化、写作-审阅和编辑。J. Oh：资源获取、数据整理、写作-审阅和编辑。N. Zeltser：正式分析、写作-审阅和编辑。T. Gebo：资源获取、数据整理、软件、写作-审阅和编辑。A. Foucal：资源获取、数据整理、软件、写作-审阅和编辑。R. Lesurf：资源获取、软件、写作-审阅和编辑。C.-H. Jung：资源获取、写作-审阅和编辑。R.M.A. Dang：软件、写作-审阅和编辑。R. Agrawal：软件、写作-审阅和编辑。J. Livingstone：资源获取、数据整理、软件、可视化、写作-审阅和编辑。A. Salcedo：资源获取、数据整理、软件、可视化、写作-审阅和编辑。C.Q. Yao：资源获取、软件、可视化、写作-审阅和编辑。S.M.G. Espiritu：资源获取、数据整理、软件、写作-审阅和编辑。K. Houlahan：资源获取、数据整理、软件、可视化、写作-审阅和编辑。F. Yousif：软件、写作-审阅和编辑。L.E. Heisler：资源获取、数据整理、软件、写作-审阅和编辑。A.T. Papenfuss：监督、资金获取、写作-审阅和编辑。M. Fraser：概念化、监督、项目管理。B.J. Pope：监督、资金获取、写作-审阅和编辑。A.U. Kishan：监督、资金获取、写作-审阅和编辑。A. Berlin：监督、写作-审阅和编辑。M.L.K. Chua：监督、写作-审阅和编辑。N.M. Corcoran：监督、资金获取、写作-审阅和编辑。T. van der Kwast：概念化、资金获取、写作-审阅和编辑。C.M. Hovens：监督、资金获取、写作-审阅和编辑。R.G. Bristow：概念化、监督、资金获取、写作-初稿、项目管理、写作-审阅和编辑。P.C. Boutros：概念化、监督、资金获取、调查、可视化、方法学、写作-初稿、项目管理、写作-审阅和编辑。

作者感谢Boutros实验室的所有成员，特别是Veronica Y. Sabelnykova提出的有益建议。这里描述的结果部分基于TCGA研究网络生成的数据。作者非常感谢West Coast Dream Team联盟分享转移性前列腺肿瘤活检的加工后的甲基化数据，以及Kirby及其同事（31）分享的临床数据。这项研究得到了Movember基金通过Prostate Cancer Canada的支持，并得到了安大略省癌症研究机构的额外支持（由安大略省政府资助）。这项工作得到了Prostate Cancer Canada的支持，并由Movember基金会自豪地资助——资助编号#RS2014-01给P.C. Boutros。P.C. Boutros还得到了Terry Fox研究机构的新研究者奖和加拿大卫生研究院的新研究者奖的资助。这项研究得到了加拿大癌症协会（资助编号#705649）和加拿大卫生研究院的项目资助给P.C. Boutros。R.G. Bristow得到了加拿大癌症协会的研究科学家奖和Cancer Research UK的核心研究者资助（C5759）。这项研究还得到了加拿大表观遗传学、环境和健康研究联盟对P.C. Boutros的资助。这项工作还得到了NIH通过奖项P30CA016042、U01CA214194、U2CCA271894和R01CA270108以及NIH NIGMS资助T32GM008042和T32GM152342的资助。本研究得到了前列腺癌基金会（Prostate Cancer Foundation）特别挑战奖（奖项编号：20CHAS01）的支持，该奖项的资金来源于Larry Ruvo的慷慨捐赠。此外，本研究还得到了美国国防部前列腺癌研究计划（Department of Defense Prostate Cancer Research Program）的资助，相关奖项分别为W81XWH2210247和W81XWH2210751。J. Oh的研究工作得到了加州大学洛杉矶分校Jonsson综合癌症中心（UCLA Jonsson Comprehensive Cancer Center）的奖学金支持。补充数据可查阅《Cancer Discovery Online》网站（http://cancerdiscovery.aacrjournals.org/）。