乳腺癌是全球排名第二的最常见恶性肿瘤,也是女性中发病率排名首位的恶性肿瘤。2022年全球有230万例乳腺癌新发病例,占女性癌症新发病例的25%,以及67万例乳腺癌死亡病例,占女性癌症死亡病例的15.5%
[1]。肿瘤的耐药性和肿瘤异质性仍是乳腺癌基础和临床研究面临的巨大挑战。乳腺癌已经取代肺癌成为全球第一大癌症。在我国,乳腺癌同样是女性恶性肿瘤中发病率最高的,占到了新发病例的16.72%,在15~44岁女性中乳腺癌导致了最多的癌症相关死亡。
三阴性乳腺癌(triple-negative breast cancer,TNBC)是指雌激素受体(estrogen receptor,ER)、孕激素受体(progesterone receptor,PR)及人表皮生长因子受体2(human epidermal growth factor receptor 2,HER2)表达均缺失的乳腺癌亚型,占所有乳腺癌的10%~15%,预后极差
[2]。对于TNBC,目前临床上治疗手段并不多,而且患者总生存期(overall survival,OS)短。现已有许多基因靶点的药物
[3],比如靶点基因为CDK4基因和CDK6基因的CDK4/6相关药物
[4]、HER2的靶向药物
[5]和抗体偶联药物(antibody drug conjugate,ADC)
[6]、PIK3CA基因编码PI3K的阿培利司
[7]、mTOR信号传导的mTORC1被依维莫司抑制
[8];针对BRCA基因突变TNBC的PARP抑制剂
[9],靶向TROP2的ADC药物等
[10],以及免疫治疗药物PD-1/PD-L1抑制剂
[11-12]。尽管乳腺癌的药物已取得了巨大成就,但仍面临诸多挑战。
随着基因测序,代谢通路等数据库的不断更新,为寻找预测乳腺癌治疗提供了新方向。细胞代谢的变化有助于转化和肿瘤进展,代谢在转移和治疗抵抗中发挥作用,改善代谢可能改善癌症治疗
[13]。代谢重编程在肿瘤进展过程中发挥重要作用,肿瘤的代谢重编程使代谢过程发生不同程度的改变,影响三羧酸循环、核苷酸和脂质代谢等
[14]。代谢途径的重新编程对于乳腺癌的适应和抗肿瘤治疗的耐药性至关重要,乳腺癌治疗后葡萄糖代谢,脂质重编程,氨基酸代谢改变,增强核苷酸代谢与乳腺癌耐药有关,同样免疫代谢影响肿瘤对免疫治疗的反应,可通过代谢成像监测乳腺癌治疗反应
[15]。仍需要关于乳腺癌的预后模型对乳腺癌进行分层治疗,为改善患者预后,提供更好的治疗方案。
基于此,代谢重编程,尤其是氧化磷酸化相关代谢异常,可能成为影响TNBC恶性进展及患者预后的重要生物学基础。然而,目前针对TNBC氧化磷酸化相关代谢基因的系统性筛选及其预后价值研究仍相对不足。为进一步探索TNBC代谢异常与患者预后的关系,本研究整合人类蛋白质图谱(The Human Protein Atlas,HPA)、国际癌症基因组联盟(International Cancer Genome Consortium,ICGC)、癌症基因组图谱-乳腺浸润性癌(The Cancer Genome Atlas-Breast Invasive Carcinoma,TCGA-BRCA)及单细胞测序等多组学数据库,筛选出与TNBC不良预后密切相关的氧化磷酸化代谢基因,构建并验证预后风险模型,同时结合肿瘤免疫微环境分析其潜在作用机制,以期为TNBC患者的风险分层、个体化治疗及潜在代谢治疗靶点的发现提供新的理论依据。
1 资料与方法
1.1 数据获取
TCGA-BRCA数据作为癌症信息数据库,保存包括基因组学、患者临床及生存数据等(
https://portal.gdc.cancer.gov)。下载已处理的基因表达数据,根据PAM50基因表达分型获取临床数据,将患者分为Basal-like(
n=184,多呈三阴性表型)、luminal A [
n=499,ER(+)/HER2(-),低增殖]、luminal B [
n=197,ER(+)/HER2(±),高增殖]及HER2-enriched [
n=78,HER2(+)/ER(-)]四个亚型。外部验证数据集标准化读取计数数据获取自ICGC BRCA-US数据集,包含192个肿瘤样本(
https://dcc.icgc.org)。在HPA数据库(
https://www.proteinatlas.org)乳腺癌数据库中下载“unfavorable genes analysis”数据,在数据“TCGA”栏筛选“unprognostic”指征的基因,以筛选出与乳腺癌患者不良预后显著相关的基因。GSE176078单细胞RNA测序数据集,从GEO数据库(
https://www.ncbi.nlm.nih.gov/geo)下载。
1.2 功能富集分析
使用STRING数据库(
https://string-db.org)对差异表达基因构建蛋白质-蛋白质相互作用(PPI)网络,以分析其相互作用关系。利用“cluster Profiler”包对差异基因进行基因本体论(Gene Ontology,GO)、京都基因与基因组百科全书(Kyoto Encyclopedia of Genes and Genomes,KEGG)富集分析,使用“ReactomePA”包进行Reactome通路富集分析。
1.3 预后模型构建
进行预后特征的筛选,从候选基因中识别出独立的预后因子,对前述筛选得到的与乳腺癌预后相关的基因进行单因素、多因素Cox回归分析构建风险评分公式,根据预后模型对TCGA-BRCA的Basal亚型队列进行评分,并通过热图展示高风险组和低风险组的相关基因表达,利用散点图展示风险评分与生存信息的关联;通过“survival包”的Kaplan-Meier方法绘制高风险组与低风险组的OS曲线;利用“survival ROC”包绘制时间依赖性受试者工作特征曲线(time-dependent ROC curve)。使用R语言中的“timeROC”包(v0.4),评估模型预测患者1、3、5、10年OS率的区分能力,并计算相应的曲线下面积(AUC)。此外,使用“survcomp”包(v1.40.0)计算模型的一致性指数(C-index)及其95%置信区间(confidence interval,CI),以评估模型在整个随访期间的整体预测准确率。
1.4 列线图模型
模型结合临床病理特征对乳腺癌的预测能力,纳入患者年龄、肿瘤分级、淋巴结转移、远处转移、风险评分进行单因素分析、LASSO回归、多因素分析,构建列线图并给出95% CI、风险比和双侧P值。计算模型的一致性指数并绘制了ROC曲线。
1.5 独立队列验证
将训练队列中确定的最终风险评分公式应用于ICGC乳腺癌三阴性队列的标准化表达数据,计算每例患者的风险评分。风险评分后分组进行基因表达热图,风险评分散点图,OS分析,绘制ROC曲线并计算一致性指数。
1.6 免疫浸润分析
从HPA数据库下载免疫相关基因或标志物的蛋白表达数据,结合转录组水平验证关键免疫相关蛋白在高、低风险组间的差异。采用ESTIMATE包
[16](v1.0.3)评估每个样本的Stromal score(基质评分)、Immune score(免疫评分)及ESTIMATE score(综合评分),并计算肿瘤纯度(tumor purity=1-ESTIMATE score)。
同时,利用CIBERSORT反卷积算法
[17],以默认的LM22基因特征矩阵(包含22种免疫细胞亚型的特征基因)计算每个样本中各免疫细胞亚型的相对丰度。运行1 000次置换检验,仅保留
P<0.05的样本用于后续分析。
1.7 统计学处理
本研究的所有统计分析及可视化均在R软件环境(版本4.2.0)和联川生物云平台
[18]中完成。比较采用非参数的Wilcoxon秩和检验。除非特别说明,所有统计学检验均为双侧检验,
P<0.05为差异有统计学意义。
2 结 果
2.1 乳腺癌风险基因与通路富集分析
为了寻找最具风险因素的泛乳腺癌基因特征,首先从HPA数据库中初筛获得562个与乳腺癌不良预后显著相关的基因(
图1A)。将筛选出562个相关基因进行KEGG功能富集分析,结果显示不良预后基因集主要富集于代谢通路(metabolic pathways)和氧化磷酸化通路(oxidative phosphorylation)等(
图1B)。最显著富集的代谢通路,包含63个不良预后基因。scCancer单细胞数据检索发现大部分基因表达主要分布于上皮来源的细胞,其中脱氧胞苷三磷酸焦磷酸酶1(dCTP pyrophosphatase 1,DCTPP1)、细胞色素C1(cytochrome C1,CYC1)等16个基因在肿瘤细胞中显著高表达(
图2A)。进一步利用COSMIC细胞系数据描述63个不良预后基因之间的相关性和乳腺癌亚型表达分布。DCTD、MAT2A等基因在TNBC中相对高表达;相反CMBL,MCEE和ALG3等基因在TNBC中低表达(补充
图1)。63个基因表达相关性集中于与ATP能量代谢相关(
图2B),后续利用PPI分析证明了这一点。数据显示从KEGG通路,Reactome通路和GO功能富集三个层面这63个基因在氧化磷酸化介导的ATP能量代谢途径上显著富集(
图2C-D)。以上,初步筛选找到63个乳腺癌风险基因,该基因集功能主要富集于氧化磷酸化介导的能量代谢途径。
2.2 TNBC代谢相关风险基因单细胞数据分析
分析GSE176078TNBC单细胞数据集,获得了1 537个在肿瘤细胞群中高表达的基因(
图3A);这些高表达基因的KEGG功能富集分析同样显著富集于代谢通路(metabolic pathways)和氧化磷酸化通路(oxidative phosphorylation)等(
图3B);与之前的63个乳腺癌风险基因取交集发现16个与TNBC相关的风险代谢基因(
图3C);PPI分析表明这16个基因主要参与电子呼吸传递链过程(
图3D)。
2.3 利用TCGA乳腺癌队列构建风险模型
进一步对上述16个相关基因进行单因素、多因素Cox回归分析。在单因素分析和LASSO回归分析中找到DCTPP1、CYC1、GPAA1为显著风险因素(补充表1)(补充
图2A),其中DCTPP1的风险比最高;在多因素分析中观察到DCTPP1、CYC1是独立的预后风险基因;进行多重共线性诊断,方差膨胀因子(variance inflation factor,VIF)共线性诊断CYC1与DCTPP1的VIF值均为1.133,不存在多重共线性并构建风险模型:风险评分=exp(DCTPP1)×0.729 2+exp(CYC1)×0.460 0。利用风险评分公式进行评分,使用中位风险评分作为阈值划分点,将TCGA-BRCA Basal队列中184个样本分为高危组(高风险评分,
n=92)和低危组(低风险评分,
n=92),用热图展示了DCTPP1、CYC1在高危组中均为高表达(
图4A);利用样本生存信息构建散点图,发现高风险评分的OS较短(
图4B)。Kaplan-Meier分析结果显示高危组的OS显著低于低危组(Log-rank
P=6.305e-03)(
图4C)。此外,ROC曲线显示了该风险模型的1、3、5、10年OS的预测能力,AUC值分别为0.664、0.655、0.662和0.775(
图4D),通过一致性指数评估为0.634(95%
CI=0.501 9~0.760 1,
P=0.037 2)进一步表明模型的高敏感度和特异度(
图4E)。不同亚型分层ROC分析表明,在TCGA-BRCA队列中,该风险评分模型对TNBC 10年OS的预测能力较稳定,优于其他亚型(
图4F)(补充
图2B)。
2.4 列线图模型构建
利用TCGA-BRCA Basal队列的年龄、TNM分期和风险评分进行单因素和多因素分析,最后把有意义的变量,年龄、远处转移(M分期)和风险评分构建列线图,可以直观地算出患者的OS率(
图5A)。评估结果显示,模型预测1、3、5、10年OS的AUC分别为0.664、0.655、0.630、0.702(
图5B)。总的区分度指标一致性指数为0.726(95%
CI=0.543~0.910,
P=0.016)(
图5C),提示预测效果好于随机猜测。Bootstrap重采样(1 000次)绘制的校准曲线,该列线图模型具有良好的校准度,其预测的OS率与实际观测结果无明显差异(
图5D)。
2.5 风险模型的验证及机制分析
为验证上述构建的风险模型,将其应用于ICGC乳腺癌队列(
n=192)进行验证。根据风险评分的中位值,将所有样本划分为高危组(
n=96)与低危组(
n=96),热图结果显示,风险模型中的关键基因在高风险组中同样显著上调(补充
图3A)。风险评分分组下的风险评分散点图,如补充
图3B;与上述结果一致,通过生存信息和风险评分的散点图,发现高风险组的预后明显低于低风险组(
图6A);与上述结果一致,高风险组的OS相对于低风险组较差(log-rank
P=3.952×10
-4)(
图6B)。进一步通过ROC分析验证评估模型的预测效能,预测的OS风险模型在1、3、5、10年的AUC分别为0.630、0.656、0.656、0.700(
图6C);模型评估,一致性指数为0.627 7(
P=3.522×10
-5)(
图6D)。综上所述,这些结果共同证实了该风险模型在ICGC乳腺癌队列中具有良好的预后预测能力。HPA免疫组化结果显示:相对于正常乳腺组织,DCTPP1在乳腺癌组织中表达显著上调;CYC1在正常组织中为低表达,在乳腺癌组织中呈异常高表达(
图6E)。利用ESTIMATE算法计算各样本的免疫、基质及综合评分再进行Spearman相关分析,结果显示,风险评分与ESTIMATE score(
R=-0.23,
P=0.001 4)、Immune score(
R=-0.16,
P=0.025)及Stromal score(
R=-0.26,
P=0.000 33)均呈明显负相关,与tumorPurity呈明显正相关(
R=0.23,
P=0.001 2)(
图6F)。其次,通过CIBERSORT对22种免疫细胞亚型的丰度进行高、低风险组间比较,虽然大多数细胞亚型未显著差异,但CD8
+ T细胞在高风险组中浸润比例明显降低(
图6G),这提示高风险TNBC预后差可能与抗肿瘤免疫浸润低有关。
3 讨 论
乳腺癌发病率持续增长。从1990—2021年,全球乳腺癌负担逐渐增加,年龄标准化发病率从每10万人16.42例上升到26.88例
[19];且中国乳腺癌2015—2030年间,死亡人数预计将分别增加13.1%和40.9%
[20]。
乳腺癌药物治疗面临诸多挑战。在早期TNBC的新辅助治疗中,PD-L1抑制剂阿特珠单抗联合白蛋白紫杉醇及蒽环类药物化疗,可将病理完全缓解率提高至69%
[21]。在NeoPanDa03 Ⅱ期研究中,PD-1抑制剂总病理完全缓解率(tpCR)达67.6%,总体缓解率达94.1%
[22],BRCA的PARP抑制剂新辅助化疗联合奥拉帕利无事件生存期提高16.3%
[23];KEYNOTE-522的Ⅲ期试验结果显示在化疗新辅助+辅助治疗的基础上,加用帕博利珠单抗使患者的5年OS率从81.7%提升至86.6%,绝对获益为4.9个百分点
[24];而对于晚期TNBC治疗,全球性Ⅲ期ASCENT研究,SG组的中位无进展生存期为5.6个月,中位OS为12.1个月
[10],ASCENT-04(Ⅲ期):戈沙妥珠单抗联合帕博利珠单抗一线治疗PD-L1阳性患者,中位无进展生存期为11.2个月
[25]。虽有较好的治疗的重要突破,但对于TNBC预后在临床治疗中,仍有较大的进步空间。部分患者初始治疗效果良好,但随着治疗时间的延长,肿瘤细胞可能通过突变或信号通路重编程产生耐药性,导致疾病复发或进展。患者分层治疗,实现对预后更精准的预测与干预,仍是待突破的瓶颈。
随着对肿瘤的不断研究,笔者越发关注肿瘤代谢基因重编程。为了维持持续增殖和转移,肿瘤细胞经历多种代谢适应以应对营养缺乏的微环境
[26];肿瘤细胞的代谢途径发生了显著变化,以支持肿瘤的发生发展,并且代谢过程不是独立发生的,是相互影响并共同促进肿瘤发生、进展和转移
[27]。代谢重编程在原发性和转移性乳腺癌治疗耐药性发展中发挥重要作用,代谢适应不仅推动肿瘤进展,还直接导致肿瘤化疗耐药;转移性病灶的代谢可塑性使癌细胞能够在不同的器官环境中茁壮成长,并逃避治疗
[28]。在乳腺癌预后风险基因的筛选中,笔者同样发现乳腺癌预后风险基因集中在代谢通路上聚集,并且笔者发现了DCTPP1、CYC1基因与乳腺癌预后相关。
DCTPP1影响乳腺癌的生长、迁移和预后。基于两个数据集的分析,从风险分组与基因表达可见高风险组DCTPP1基因具有较高的表达,并且在预后分析中,同样印证了此结果,DCTPP1高表达的乳腺癌患者预后较差,且通过免疫组化染色发现乳腺癌组织的表达高于正常乳腺癌组织。其他研究有类似的结果,DCTPP1可通过对dNTP的高水解活性来维持dNTP池稳态和异常核苷酸净化有助于基因组稳定性,但在癌细胞核中积累,并在癌组织中显著过表达,促进肿瘤进展,并与预后不良和免疫、化疗治疗耐药性相关
[29]。通过DCTPP1缺陷的MCF-7细胞和DCTPP1过表达的MDA-MB-231细胞,发现DCTPP1主要控制5-甲基-dCTP代谢和整体DNA低甲基化参与促进乳腺癌细胞增殖和干性维持
[30]。有研究
[31]发现,唐氏综合征细胞黏附分子反义RNA 1(DSCAM-AS1)可促进DCTPP1的过表达,促进MCF7和T47D细胞的生长、迁移和侵袭,同样敲低DSCAM-AS1降低DCTPP1的表达,可抑制luminal型乳腺癌的生长、迁移和侵袭。
CYC1在乳腺癌组织中表达上调,与肿瘤的侵袭性进展和不良预后密切相关。CYC1是线粒体复合体Ⅲ的重要亚基,与氧化磷酸化和ATP产生相关,其在TNBC的研究并不多。笔者通过两个三阴性乳腺癌数据集,发现CYC1与TNBC预后息息相关,并在高风险组有较高的表达;Han等
[32]使用在线数据库和工具,同样发现乳腺癌患者的CYC1表达较高与预后较差相关,利用两种高度转移的TNBC细胞系,发现CYC1抑制AMPK的活化促进肿瘤转移,并通过促进ATP的产生促进肿瘤生长。Chishiki等
[33]发现40%的乳腺导管癌病例有CYC1免疫反应性,且免疫组化CYC1状态与肿瘤大小、核级、粉刺坏死、范奈斯分类及Ki-67标记指数显著相关,这与数据库验证的免疫组化,乳腺癌组织的表达高于正常乳腺癌组织相关。有研究发现CYC1等多个与线粒体功能和生物合成密切相关的基因表达,可抑制ER阳性乳腺癌细胞的生长
[34]。
总之,本研究聚焦于代谢风险基因与预后差的TNBC,通过数据库筛选乳腺癌预后相关的基因,发现DCTPP1、CYC1代谢基因影响乳腺癌的预后,并建立风险模型,进一步与临床变量的关联构建了列线图模型,经过验证集评估了该模型预测工具有良好区分能力和稳定性。本研究的结果可能有助于TNBC预后风险评估,为乳腺癌分层治疗和个体化精准治疗提供依据,以及可能发现乳腺癌治疗的新型分子生物标志物。
湖南省科技创新计划资助(2022RC1019)
湖南省科技创新计划资助(2022SK2041)
吴阶平医学基金会科研专项资助(320.6750.2025-20-4)
湖南省自然科学基金项目资助(2023JJ40986)