利用多种转录组分析方法开发用于预测三阴性乳腺癌预后的模型

《PLOS One》：Development of predictive models for the prognosis of triple-negative breast cancer using multiple transcriptomic analyses

　　摘要：三阴性乳腺癌（TNBC）是乳腺癌（BC）的一个亚型，约占所有BC病例的15-20%。这种亚型具有最具侵袭性的行为和最差的预后。在过去的几十年中，进行了大量研究来解决临床上可用的治疗选择有限的问题。特别是，针对有效治疗方案的潜在标记物已被积极研究。然而，这些努力受到TNBC复杂机制的阻碍，没有一项研究的预测性能超过0.85。本研究开发了预测性能超过0.94的TNBC预后模型。将这九个选定的标记物应用于五个独立数据集，证明了它们作为TNBC特异性预后标记物的潜力。这些基因中的大多数（包括GPR61、PZP、IGFL1和AHCTF1）与TNBC患者的总生存期（OS）相关。基于这些结果，这九个选定的基因可能作为TNBC患者OS的预后标记物。

引言：乳腺癌（BC）是全球女性中最常诊断出的癌症，并且具有最高的癌症相关死亡率[1,2]。它具有异质性，通常根据雌激素受体（ER）、孕激素受体（PR）和人表皮生长因子受体-2（HER2）的表达分为四种亚型[2,3]。三阴性乳腺癌（TNBC）是BC的四种亚型之一，占所有BC病例的约15-20%[2,4]。其定义为三种受体（ER、PR和HER2）均无表达[5]，因此对靶向激素受体和HER2的药物无反应。在解决这一问题方面取得了显著进展，特别是在针对TNBC患者的新药如多聚ADP-核糖聚合酶抑制剂和抗体-药物偶联物方面[6-8]。然而，这些研究主要关注TNBC患者的特定子集，如带有BRCA1/2种系突变或转移性TNBC的患者。尽管有这些进步，TNBC仍然是一种具有较差预后的疾病，与其他BC亚型相比。具体来说，原发性TNBC的5年总生存期（OS）率显著较低（77%），而其他BC亚型（Luminal A（LumA）、Luminal B（LumB）和HER2阳性）的OS率分别为93%[9-11]。鉴于这种不良预后以及大多数TNBC患者缺乏有效的靶向治疗，识别可靠的预后生物标志物具有重要的临床意义。这类生物标志物可以改善风险分层，指导治疗决策，并识别可能从新的治疗策略或强化监测中受益的患者[12]。因此，迫切需要额外的预后标记物和治疗靶点来进一步改善患者的预后，这是本研究的主要目的。在过去的几十年中，进行了大量研究，通过识别潜在的标记物来克服TNBC患者的治疗选择限制。新辅助化疗（NAC）旨在缩小肿瘤大小，已被确立为TNBC的标准治疗方法，对于实现病理完全缓解（pCR）的患者预后有所改善[2,13]。由于NAC的临床益处，许多研究集中在识别预测NAC后pCR的标记物，以便对患者进行更有效的NAC治疗[10,14,15]。研究还确定了TNBC对化疗的敏感性预测因子（如BCL2），包括术后辅助化疗[16,17]。此外，还进行了针对免疫检查点抑制剂（ICI）加化疗反应的预测因子研究，以识别新的治疗药物[18,19]。不幸的是，由于TNBC的复杂机制并非仅依赖于特定信号，这些努力未能成功，有效的治疗策略尚未开发出来[20]。除了化疗或ICI之外，先前的研究还试图识别预测TNBC患者预后（如复发、死亡、无病生存期（DFS）和OS）的标记物[5,21,22]。Campione等人开发了一个使用三种蛋白质特征（TrpRS、TSP1和DP）预测复发的模型，其曲线下面积（AUC）为0.82[5]。Xu等人使用临床病理数据开发了基于多种机器学习算法的死亡预测模型，具有最高的预测性能，AUC达到0.732[21]。Yang等人使用临床病理数据开发了一个预后模型来预测DFS和OS，该模型分别显示DFS的AUC为0.784，OS的AUC为0.783[22]。然而，目前没有模型的预测性能在验证数据上超过0.9，表明仍需开发更准确的模型。最近，利用LumA数据进行了研究，LumA表达两种激素受体，在BC亚型中预后最好。Choi等人构建了一个分子调控网络模型，将TNBC细胞重新编程为LumA细胞，并确定了BCL11A和HDAC1/2作为诱导向LumA细胞转化的最佳靶点[1]。Singhal等人建立了过表达SLFN12的TNBC细胞系驱动的人类BC异种移植物，导致更高水平的LumA标记物和HER2受体表达，最终获得更好的生存率[23]。在这项研究中，我们旨在通过使用多种机器学习算法并通过多队列转录组数据集验证，开发出在预测TNBC亚型预后方面表现出色的模型。

材料与方法：数据来源：这项回顾性研究使用了六个数据集：四个RNA测序（RNA-seq）数据集和两个微阵列数据集。来自癌症基因组图谱乳腺浸润性癌症（TCGA-BRCA）队列的RNA-seq数据集作为模型开发的主要数据集；基因表达谱（Illumina HiSeq 2000）来自UCSC Xena Data Hub（https://xenabrowser.net/）。仅包括原发肿瘤以构建早期预后预测模型。亚型根据Microarray 50（PAM50）特征定义，包括143例TNBC（基底样）、386例LumA、186例LumB和69例HER2阳性样本。OS和死亡事件被用作主要结局。为了外部验证模型开发过程中识别的预后基因，我们使用了两个基因表达组（GEO）数据集（GSE65216和GSE215442）、两个细胞系数据集以及METABRIC队列。GSE65216是一个GPL570（Affymetrix U133 Plus 2.0）微阵列数据集，包含55例TNBC、29例LumA、30例LumB和39例HER2阳性样本。GSE215442是一个从MDA-MB-231 TNBC细胞生成的RNA-seq数据集，这些细胞过表达SLFN12，形成了预后较好的LumA样亚克隆，包括三个SLFN12过表达系和三个对照组[23]。细胞系数据集包括一个单细胞RNA-seq数据集（7,484例TNBC和4,599例LumA细胞）和一个批量RNA-seq数据集（31例TNBC和10例LumA细胞系）。METABRIC是一个包含320例TNBC肿瘤的微阵列队列，用于外部验证；在随访期间，168名患者（52%）死亡，中位OS为13.3年。

模型开发与评估：整个工作流程如图1所示。本研究的目标是利用Cox比例风险（CoxPH）回归基于事件时间数据开发TNBC的预后模型。包括了TCGA-BRCA队列中的所有143例TNBC样本，因为CoxPH模型可以在不需要预定义对照组的情况下估计队列内的相对风险。

数据集被随机分成70/30的训练集和测试集。在100次重采样中，在训练集上拟合了单变量Cox比例风险模型；在≥80/100次重采样中显著的基因被保留为候选基因。在训练集中基于AUC的逐步选择产生了最终的基因特征。对于每次重采样，都在训练集上训练该特征并在其对应的测试集上进行评估，并总结了所有100对的性能（AUC）。结局指标为5年和10年的总生存期（OS）。

模型开发与评估的流程如图1所示。数据集被随机分为训练集（70%）和测试集（30%），并按事件状态分层以保持死亡和删失病例的比例。对于每个基因，使用训练集拟合单变量CoxPH模型，p值<0.05的基因被视为显著基因。此过程重复100次，并在随机重采样中，至少在80次迭代中显著的基因被保留为候选预测因子。为了确定最佳的预测因子组合，我们应用了基于训练AUC的逐步选择程序。首先分别使用CoxPH拟合每个候选基因，并计算100次重采样的平均训练AUC。选择平均AUC最高的基因作为初始模型（Mk = Mk = 1）。在每个后续向前步骤中（k = k + 1），通过向当前模型（Mk-1）中添加一个之前未选中的预测因子来评估所有候选模型，选择平均AUC最高的模型作为Mk。只有当新模型的平均AUC超过Mk-1超过α = 0.005时，才接受新模型（Mk），这作为防止过度拟合的最低改进阈值；否则，保留Mk-1并停止该过程。在每个向后步骤中（k = k + 1），我们通过从Mk-1中移除一个预测因子来评估所有简化模型（一次删除一个预测因子；候选因素的数量等于Mk-1中的预测因子数量）。只有当简化模型的平均训练AUC超过Mk-1超过α = 0.005时，才接受新模型（Mk）；否则，保留Mk-1并返回到向前步骤。向前和向后步骤交替进行，直到在前向阶段不再观察到平均训练AUC的进一步增加。通过这种选择程序获得的最终预测因子用于开发三种算法的预后模型——CoxPH、随机生存森林（RSF）和生存支持向量机（Survival-SVM）。对于RSF和Survival-SVM，通过最大化训练AUC来优化超参数。使用相应的测试集评估了100次重采样的模型性能，以时间依赖的AUC、精确度-召回曲线下面积（AUPRC）和c指数作为5年和10年总生存结果的性能指标。

统计分析：使用双侧Wilcoxon秩和检验测试基因表达的组间差异。对于生存分析，使用最大选择秩统计量（MaxStat）[24]确定每个基因的最佳表达截止值，然后将样本分为高表达组和低表达组。生存差异使用双侧对数秩检验进行比较。表达组（高 vs. 低）与总生存状态之间的关联使用Fisher的确切检验进行评估。统计显著性设定为p值<0.05。

结果：在TCGA-BRCA队列的16,336个蛋白质编码基因中，我们最初筛选了与OS相关的候选预测因子（基因）。在队列中的784名BC患者中，有143名（18%）被确定为TNBC（S1表）。其中，142名患者的NAC历史可用，他们都没有这样的历史；一名患者的数据缺失。在随访期间，18名TNBC患者（13%）死亡，中位OS为20.4年。为了防止在模型开发过程中过拟合，总数据集被随机分为7:3比例的训练集和测试集，并且仅使用训练集进行从筛选到模型开发的过程。使用CoxPH模型筛选候选预测因子。选择显著性水平为5%的基因。为了避免选择特定数据集依赖的预测因子，数据被随机分割100次以选择显著基因。预期这可以减少由于随机分割导致的选择偏差。经过100次迭代后，有53个候选变量被选为最终候选变量超过80次。我们使用这些53个候选变量进行了基于训练AUC的逐步选择。5年OS和10年OS（作为之前开发的TNBC预后预测模型的主要结局[22,25]）被用作响应变量。对于5年OS，超过90%的测试集的测试AUC超过0.9，平均AUC和AUPRC值分别为0.9459和0.8027（图2）。对于10年OS，超过98%的测试集的测试AUC超过0.9，平均AUC和AUPRC值分别为0.9570和0.9070（图2）。ROC曲线和精确度-召回率曲线的平均AUC和平均AUPRC已被分别报告。https://doi.org/10.1371/journal.pone.0348414.g002 除了作为本研究主要模型的CoxPH模型外，还应用了其他机器学习（ML）算法（RSF和Survival-SVM）到100个训练集和测试集上。CoxPH和其他ML算法在AUC、AUPRC和C指数方面的预测性能均超过了0.8（表1）。下载：PNG（大图像）TIFF（原始图像）表1. 每个响应变量的最终选定预测模型。https://doi.org/10.1371/journal.pone.0348414.t001 对OS有贡献的系数的方向和显著性在100个训练集和整个数据集中呈现了一致的趋势。根据整个数据集的拟合CoxPH模型（表2），CELF6、IGFL1、GPR61和TTLL2的系数为正，表明随着表达水平的增加，生存时间缩短。相比之下，其他预测因子（包括TMEM14B和CREB5）的系数为负，表明随着表达水平的降低，生存时间缩短。这些发现与Kaplan-Meier分析结果一致（图3）。此外，Fisher精确检验显示幸存者和非幸存者之间的基因表达组（高表达组与低表达组）存在显著差异（表3），进一步支持了基因表达与生存结果之间的关联。标准化系数表明CELF6和IGFL1是对OS最具有影响力的预测因子（图4）。下载：PNG（大图像）TIFF（原始图像）表2. 每个响应变量的CoxPH模型拟合结果（多元分析）。https://doi.org/10.1371/journal.pone.0348414.t002 下载：PNG（大图像）TIFF（原始图像）表3. 由基因特异性阈值定义的基因表达组（高表达组 vs. 低表达组）与生存状态（随访时存活 vs. 死亡）之间的关联。进行了Fisher精确检验以评估幸存者和非幸存者之间表达组的差异。https://doi.org/10.1371/journal.pone.0348414.t003 下载：PNG（大图像）TIFF（原始图像）图3. TNBC中九个基因的Kaplan-Meier生存曲线（n = 143）。对于每个基因，样本根据MaxStat得出的最佳临界值被分为高表达组（红色）和低表达组（蓝色）。使用对数秩检验评估生存差异。https://doi.org/10.1371/journal.pone.0348414.g003 下载：PNG（大图像）TIFF（原始图像）图4. OS预测模型中的标准化基因系数。对于每个结果（5年OS和10年OS），显示了针对整个数据集拟合的最终模型的系数。所有变量在拟合前都经过了z标准化，以便直接比较基因效应。https://doi.org/10.1371/journal.pone.0348414.g004 为了确认这九个选定的TNBC预后标志物的潜力，我们在多个队列数据集中对其进行了验证。我们假设如果表达水平的增加会导致预后恶化，那么其在TNBC中的表达水平将高于LumA中的表达水平，后者在BC亚型中预后较好。为了验证我们的假设，我们比较了四个数据集（包括TCGA-BRCA队列）中每个基因在TNBC和其他BC亚型之间的表达模式。在这九个基因中，TTLL2和GPR61在TCGA-BRCA和GSE65216数据集中都表现出与这一假设一致的趋势（图5A和5B）。具体来说，TTLL2和GPR61在TNBC中的表达水平显著高于LumA（TCGA-BRCA中的TTLL2和GPR61的p值分别为<2.2E-16和<4.1E-11）。此外，这两个基因在TNBC中的表达水平也高于其他BC亚型。这一趋势在GSE65216数据集中也是一致的。下载：PNG（大图像）TIFF（原始图像）图5. 四个独立队列中的外部验证。(a) TCGA-BRCA；(b) GSE65216；(c) GSE215442；(d) 细胞系批量RNA-seq。对于每个队列，箱线图显示了表达符合假设方向的九个基因。倍数变化（FC）是TNBC中的平均表达水平除以LumA中的表达水平。使用双侧Wilcoxon秩和检验比较了TNBC与LumA、LumB和HER2。https://doi.org/10.1371/journal.pone.0348414.g005 与TCGA-BRCA和GSE65216数据集不同，GSE215442数据集旨在从TNBC细胞系（MDA-MB-231）生成SLFN12过表达的异种移植瘤，从而产生预后更好的LumA样TNBC细胞系。GSE215442数据集包含两组RNA-seq数据：原始TNBC和LumA样TNBC细胞系。对GSE215442数据集的分析确定了三个基因，在LumA样TNBC组中的表达水平趋势与选为预测因子的九个基因的趋势一致（图5C）。GPR61基因在TCGA队列中显示出随着表达水平的增加，生存时间缩短的趋势，并且在原始TNBC组中的表达水平高于LumA样TNBC组。相比之下，TMEM14B和PZP在TCGA队列中显示出随着表达水平的降低，生存时间缩短的趋势，并且在原始TNBC组中的表达水平低于LumA样TNBC组。GSE65216数据集和基于细胞系的批量RNA-seq数据集显示TNBC中的IGFL1表达水平高于LumA中的表达水平，但差异不具有统计学意义（图5B和5D）。通过多队列数据的验证分析，我们确认了五个基因（TTLL2、GPR61、TMEM14B、PZP和IGFL1）在至少一个独立数据集中得到了验证。有趣的是，尽管GPR61（G蛋白偶联受体61）与TNBC预后之间的关联尚未明确，但GPR61在三个数据集中得到了验证。GPR161（属于同一受体家族）在TNBC中过表达，并且在敲低实验中会抑制TNBC细胞系的增殖[26]。考虑到GPR161是一个潜在的药物靶点，预计GPR61也是如此。讨论我们开发的TNBC预后模型在测试集中的AUC超过了0.94，优于之前报道的OS预测模型，后者的AUC通常低于0.85[22,25,27]。也有类似的高性能模型被报道，包括一个10基因的早期TNBC特征集[28]、一个基于干性的预后模型[29]和一个EMT相关基因特征集[30]。这些研究共同表明，紧凑的基因集可以实现临床上有意义的风险分层。我们的模型通过结合基于AUC的特征选择的系统重采样来扩展这种方法，从而提高了模型的泛化能力。除了技术性能外，长期生存预测具有重要的临床相关性。虽然治疗反应预测主要用于初始治疗决策，但长期生存预测（5年和10年OS）提供了互补但不同的临床价值。TNBC是一种高度侵袭性和异质性的疾病，即使在具有相似临床特征的患者中，生存结果也有很大的变异性。这种异质性需要个性化的预后评估和准确的长期风险分层[31]。尽管进行了大量的研究，但TNBC的稳健预后工具仍然有限，现有的临床标志物如pCR只能提供部分预后信息[32]。生存预测通过实现风险分层和支持长期管理决策（包括治疗强度和随访计划）提供了超出治疗反应的临床可操作性见解[31,33]。此外，相当一部分TNBC复发和死亡发生在诊断后的五年之后，即使最初取得良好反应（例如pCR）的患者也可能出现晚期复发[34,35]。这些发现表明，治疗反应和长期生存捕捉到了疾病进展的相关但不同的方面，突显了长期生存预测在TNBC中的重要性。在评估的算法中，CoxPH模型表现出最佳的预测性能（表1）。尽管先进的ML方法可以有类似的表现，但先前的研究表明，当样本量适中且关系大致线性时，基于CoxPH的模型仍然具有竞争力，并且通常更优秀[36,37]。与这些报告一致，CoxPH模型在TCGA-BRCA数据集中实现了最高的区分度。鉴于我们数据集中的类别不平衡，AUPRC提供了一个补充的性能指标。在142名患者中只有18人（13%）发生了事件，相应地，随机分类器的基线AUPRC为0.127。尽管存在这种不平衡，我们的模型分别实现了5年和10年OS的AUPRC值为0.8027和0.9070，显著超过了基线，证明了其对少数类的强大预测能力。相应的AUC值分别为0.9459和0.9570。由于AUC和AUPRC具有不同的基线和尺度，直接进行数值比较并不合适；然而，两个指标的一致高值支持了我们的模型在类别不平衡下的稳健性。使用METABRIC队列进行的外部验证显示出较低的预测性能，这可能是由于METABRIC（微阵列）和TCGA-BRCA（RNA-seq）之间的临床组成和检测平台差异所致。类似的地跨平台差异已被报道，像EMBER这样的框架已经证明统计协调可以改善数据集之间的整合[38]。应用这样的方法可能会进一步提高跨队列的可重复性。Kaplan-Meier分析显示TCGA-BRCA和METABRIC之间的生存趋势一致（S1图），支持了所识别标志物的生物学合理性。在选定的九个基因中，有五个（PZP、GPR61、TTLL2、TMEM14B和IGFL1）在各个队列中显示出一致的表达模式和效应方向，这与将它们与TNBC增殖和生存联系起来的先前的研究一致[1,23,26,39,40]。其余四个基因显示出不一致的系数，但表达模式相似（S2图），表明可能存在亚型特异性效应[41]。最近的研究确定了免疫相关和B细胞相关的特征是早期TNBC的强大预后决定因素[42]，这表明将免疫相关特征纳入我们的模型可能会进一步提高预测性能。除了OS，我们还使用相同的框架开发了无进展生存（PFS）和DFS的模型（S2表）。由于这些终点反映了不同的生物学过程，因此预计模型性能会有所不同[43,44]。这些发现突显了我们框架在多个预后结果上的潜在适用性。总之，我们提出了一个可重复且性能优异的TNBC预后模型，其性能超过了之前的基准，并与新兴文献一致。尽管在外部验证和实验确认方面存在局限性，但模型在多个数据集中的强可重复性和与最近预后研究的一致性[28–30,38,45]强调了其稳健性。未来的研究将集中在实验验证、跨平台协调以及多组学和免疫特征的整合上，以进一步推进TNBC的精确预后。支持信息

参考文献