基于焦亡相关基因构建乳腺癌预后风险模型

陈睿鹏，臧洪婧，孙 意

（中南大学湘雅二医院病理科，湖南 长沙 410011）

乳腺癌（breast cancer，BC）是女性最常见的恶性肿瘤，是一组高度异质性的疾病[1]。随着生活方式的改变，自20 世纪70 年代以来，乳腺癌在世界范围内的发病率显著上升[2]。焦亡（pyroptosis）也称细胞炎性死亡，定义为炎性介质介导的程序性死亡[3]。细胞焦亡的特点是gasdermin 被启动形成细胞膜孔，这与细胞凋亡机制不同[4]。在促炎因子作用下，完整的细胞膜被破坏，最终导致细胞死亡[5]。研究表明[6]，焦亡在人类恶性肿瘤的发生发展中起不可或缺的作用。焦亡与癌症的作用机制是复杂的[7，8]：一方面，介导肿瘤细胞焦亡以抑制其生长；另一方面，焦亡为肿瘤的生长提供了营养和免疫环境支持[9，10]。已有研究发现[11-13]，焦亡相关基因与卵巢癌、肺腺癌、胃癌有关。这些研究表明，焦亡基因可作为一种新的基因标记预测患者的生存预后。然而，焦亡基因在乳腺癌中的表达水平、与临床预后特征和免疫浸润状态的关系尚未明确。为此，本研究利用癌症基因组图谱（The Cancer Genome Atlas，TCGA）的基因表达资料构建预后风险模型，拟评估乳腺癌患者的生存状态，并使用基因表达综合（Gene Expression Omnibus，GE0）数据库进行验证。

1 资料与方法

1.1 资料来源 TCGA-乳腺癌队列的基因表达RNA测序数据及临床信息从genome data Commons data Portal（https://portal.gdc.cancer.gov/）获得。验证集GSE20685 和GSE42568 的数据来自GEO 数据门户网站（https://www.ncbi.nlm.nih.gov/geo/）。在去除没有临床信息的数据后，本研究的训练集——TCGA 队列纳入了1045 个乳腺癌样本和113 个相邻正常样本。验证集GEO 队列纳入了431 个乳腺癌样本。此外，179 份正常人类乳腺组织样本的基因表达RNA-seq 数据从GTEx 数据集（https://xenabrowser.net/datapages/）中获得。TCGA 数据库中的体细胞突变图谱来源于genome Data Commons Data Portal（https://portal.gdc.cancer.gov/）。

1.2 焦亡相关差异表达基因及通路的关键基因 使用TCGA 数据进行差异基因表达分析，发现乳腺癌样本（1045 个）与癌旁正常组织样本（113 个）的数量级存在差异，这可能导致统计误差。因此，GTEX 数据集中正常乳腺组织样本（179 个）也被纳入研究，以确保结果的正确性。该网站（https://xenabrowser.net/datapages/）已对数据预处理，去除了批次效应，可直接运用。此外，在线数据库STRING（v11.0,http://www.string-db.org/）用于可视化焦亡相关基因的蛋白质相互作用网络（PPI），PPI 分析所需的最小交互评分设置为0.4（默认推荐值）[14]。将结果文件导入Cytoscape（https://cytoscape.org/）通过MNC 算法鉴定焦亡通路的关键基因（hubgene）[15，16]。

1.3 构建与验证风险评分模型 使用Kaplan-Meier曲线和Log-rank 检验从上述步骤得到的焦亡相关差异表达基因集中筛选出具有预后价值的基因（P＜0.05）。然后，基于具有显著预后价值的焦亡相关基因进一步分析。应用最小绝对收缩和选择操作数（LASSO）惩罚方法构建预后模型。通过10 倍交叉验证确定最佳参数λ。根据“survminer”R 包计算的最佳截断值，将TCGA 训练集分为高、低风险组，通过“timeROC”软件包绘制依赖于时间的受试者工作特征（ROC）曲线，评估模型在TCGA 和GEO 队列中的特异性和敏感性。

1.4 风险评分与临床分子特征及免疫浸润的相关性使用R 语言，采用Wilcoxon 秩和检验（分类变量为2 个）及Kruskal-Wallis 检验（分类变量大于2 个）探究焦亡风险评分与临床分子特征的相关性。随后，使用TIMER 数据库（https://cistrome.shiny apps.io/timer/）阐明TCGA 队列中高、低风险组的免疫浸润水平差异[17]。

1.5 基于风险亚组的差异基因进行通路富集及基因组分析 在R 语言中通过Wilcoxon 秩和检验评估高、低风险亚组的差异表达基因。阈值设为错误发现率（FDR）＜0.05，样本质检表达量的差异倍数（FC）≥2。随后，利用WebGestalt（https://www.webgestalt.org/option.php）对焦亡高、低风险组的差异表达基因进行功能和通路富集分析，包括基因本体论（GO）和京都基因与基因组百科全书（KEGG）[18]。将目标基因输入WebGestalt 网页，应用“ggplot2”R 包绘制结果。此外，应用R 包通过“maftools”R 包处理以突变注释格式（MAF）形式排序的体细胞突变谱，探究高、低风险亚组的基因组改变情况。

1.6 列线图的构建与评估 单因素及多因素Cox 回归分析评估风险评分及临床特征与生存预后的关系。整合风险评分及有预后价值的临床特征构建列线图（nomogram）。应用1、3 和5 年生存率的校准曲线来评估列线图与理想模型的偏差。构建ROC 曲线，计算曲线下面积（AUC），以确定列线图生存预后的预测能力。

1.7 统计学方法 所有数据通过R 语言（版本4.0.5）和R Bioconductor 包进行处理。本研究的主要终点是总生存期（OS），即从最初诊断到最终死亡的时间间隔。用Kaplan-Meier 曲线和Log-rank 检验比较风险评分组的总生存期，使用“rms”包构造列线图，C 指数用“survcomp”软件包计算。除非特别说明，所有统计检验均为双向检验，P＜0.05 为差异有统计学意义。

2 结果

2.1 焦亡基因的差异表达分析及关键基因 共收集35 个焦亡基因，检测出32 个为差异表达基因（P＜0.05），见图1A；其中18 个基因（CASP1、CASP4、CASP8、CASP9、ELANE、GPX4、GSDMB、GSDMD、GSDME、IL6、NLRP1、NLRP3、NOD1、PJVK、PLCG1、PRKACA、SCAF11 和TIRAP）在肿瘤组织中表达下调，14 个基因（AIM2、CASP3、CASP6、GSDMA、GSDMC、GZMA、IL18、IL1B、NLRC4、NLRP2、NLRP7、NOD2、PYCARD 和TNF）在肿瘤组织中表达上调。此外，35个基因集中PYCARD、AIM2、IL18、CASP1、CASP5、CASP8、NLRC4、IL6、NLRP3 和TNF 为焦亡通路的关键基因，见图1B。

图1 焦亡相关基因在乳腺癌中的表达水平及蛋白相互作用

2.2 焦亡风险评分模型构建 从32 个焦亡相关差异表达基因中鉴定出5 个具有预后价值的基因：CASP9、GSDMC、GZMA、IL18 和PYCARD，均与乳腺癌患者的预后有关，见图2。基于这些具有预后价值的基因，使用LASSO-Cox 回归分析构建评分模型（图3A），通过10 倍交叉验证（图3B），得到最小lambda 值=0.000 753 102 3，将相应系数代入，得到计算公式：风险评分=0.107 146 60×GSDMC 表达值-0.207 295 50×CASP9 表达值-0.205 471 10×IL18表达值-0.038 234 58×PYCARD 表达值-0.092 932 92×GZMA 表达值。由此可知，CASP9、IL18、PYCARD、GZMA 表示上调与预后良好相关，而GSDMC 表达上调与预后不良相关。使用最佳截断值（cut-off=0.056 399 17）将TCGA 队列分为高、低风险组，结果显示高、低风险组的KM 生存曲线比较，差异有统计学意义（P＜0.05），见图3C；高低风险组的生存死亡状况、人群分布趋势以及构成风险评分的基因表达热图见图3D；1、3、5 年的ROC 曲线的AUC 分别为0.65、0.61、0.62，见图3E，表明风险评分预测乳腺癌的预后有较高敏感性和准确性。

图2 有预后价值的焦亡基因KM 曲线

图3 乳腺癌焦亡基因预后模型的构建

2.3 外部数据集验证模型 以GSE20685（n=327）和GSE42568（n=104）作为外部验证数据集，使用TCGA 训练集的最佳临界值（cut-off=0.056 399 17）将GEO 验证队列分为高、低风险组，KM 生存曲线显示，高风险组预后差于低风险组，见图4A、图4B，其中两个GEO 队列的生存和死亡状态、人群分布和预后相关焦亡基因的表达水平见图4C、图4D。在验证集GSE20685 中，1、3、5 年的ROC 曲线的AUC 分别为0.77、0.68、0.69，见图4E；在验证集GSE42568中，1、3、5 年的ROC 曲线的AUC 分别为0.62、0.66、0.72，见图4F。

图4 GSE20685 和GSE42568 队列验证风险评分

2.4 风险评分与临床特征、分子亚型及免疫浸润的相关性 TCGA 队列中高、低风险组的患者分布趋势见图5A～图5D；男性、导管癌、肿瘤大于2 cm、三阴分子亚型、高MSI 和TMB 状态的患者风险评分较高（P＜0.05），见图5E～图5J。免疫浸润结果显示，低风险组活化CD4 T 细胞、活化CD8 T 细胞等效应细胞浸润相对较高，见图5K。

2.5 高低风险组富集途径和基因组改变分析 在TCGA 队列的高风险组中，562 个基因上调，210 个基因下调，见图6。GO 分析显示，高风险组上调基因富集于上皮细胞的生长及角化过程（图7A），高风险组的下调基因富集于个别信号通路及肌球蛋白的合成（图7B）。KEGG 分析显示，上调基因主要富集在炎症介质介导的信号通路（图7C），下调基因主要富集于生物代谢反应（图7D）。高低风险组最频繁突变的基因并不完全一致（图7E、图7F），这揭示了焦亡风险评分组基因组改变的差异。采用Fisher 检验对高低风险组不同突变基因进行探究，结果发现TP53处于首位（图7G），表明TP53 与乳腺癌发生细胞焦亡有较高的相关性且提示预后不良。在高低风险组共发生和互排斥的突变中，TP53-PIK3CA 存在互斥突变现象（图7H、图7I），提示两者的突变在乳腺癌中可能产生拮抗作用。

图6 差异基因的火山图

图7 高低风险亚组的富集途径和基因组改变分析

图7 高低风险亚组的富集途径和基因组改变分析（续）

2.6 构建列线图 单因素和多因素Cox 分析表明，焦亡风险评分可作为预测乳腺癌生存状态的独立预后指标，见表1。整合风险评分和预后相关的临床特征（年龄、远程转移分期、TNM 分期）以构建列线图（图8A）。时间依赖性ROC 曲线显示，列线图预测效果优于单独的风险评分（图8B）；列线图的KM 生存曲线结果显著（图8C）；校准曲线表明，与理想模型相比，列线图性能良好（C-index=0.7618，图8D）。

表1 TCGA 训练集的单变量和多变量分析

表1（续）

图8 基于TCGA 队列构建列线图

图8 基于TCGA 队列构建列线图（续）

3 讨论

细胞焦亡是一种独特的程序性细胞死亡，不同人体组织和遗传景观下焦亡对癌症有不同影响[7]。本研究从已发表的文献中收集了35 个焦亡基因，其中大多数为差异表达基因。与正常组织相比，肿瘤组织中18 个焦亡基因表达减少，14 个基因表达增加。PYCARD、AIM2、IL18、CASP1、CASP5、CASP8、NLRC4、IL6、NLRP3 和TNF 是焦亡通路的关键基因。本研究采用LASSO Cox 回归分析建立了一个包括5 个焦亡基因的风险评分，该评分在预测乳腺癌患者1、3、5 年总生存率方面具有较高的准确性和实用性。结果表明，CASP9、IL18、PYCARD、GZMA 是预后的保护因素，GSDMC 是预后的危险因素。

CASP9 是一种半胱氨酸蛋白酶，可触发细胞凋亡[19，20]。有研究表明[21]，在人类乳腺癌中，通过上调CASP9 表达抑制miR-182-5p 来诱导肿瘤细胞凋亡并抑制其增殖。IL18 参与了多种恶性肿瘤的发生发展，文献报道[22]其在胃肠道癌、乳腺癌和恶性黑色素瘤中表达上调，这与本研究的结果一致。IL18 具有促肿瘤和抗肿瘤的双重作用，在结肠炎相关结直肠癌中可通过修复上皮屏障对机体起保护作用[23]；而在肝细胞癌中可驱动癌细胞发生转移致预后不良[24]。促凋亡基因PYCARD 在焦亡通路中起重要作用，目前尚缺少其对乳腺癌患者预后影响的报道，但有研究表明相较于正常组织，乳腺癌患者中PYCARD 高表达。这很可能与PYCARD 基因的甲基化水平有关[25，26]。毒素颗粒酶A（GZMA）是一种通过Caspase 途径介导细胞凋亡的胰蛋白酶[27]。GZMA 在大多数胰腺癌（70%）中表达，在＜35%的乳腺癌、宫颈癌、肝癌、卵巢癌、前列腺癌、肾癌、胃癌、睾丸癌和尿路上皮癌中表达，在＜10%的淋巴瘤和黑色素瘤中表达。CTL/NK 细胞可通过表达上调GZMA 杀灭癌细胞[28]。GSDMC 属于gasdermin 家族中一员，有报道称GSDMC 参与乳腺癌细胞焦亡并在该病中高表达，同时与患者预后不良相关[29]。

肿瘤微环境在抗肿瘤分子治疗中起至关重要的作用。同在卵巢癌及胃癌观察到的现象类似[12，13]，本研究中低风险组的免疫细胞浸润程度高于高风险组，这进一步证实了焦亡在肿瘤免疫微环境中发挥着重要作用，而高风险组患者较差的预后可能是由于抗肿瘤免疫水平降低所致。比较高低风险组的基因组变化，提示高风险组与更具有侵略性的分子改变（如TP53 突变）显著相关，同时本研究观察到一种特殊的TP53-PIK3CA 互斥突变现象，这些基因组的改变与乳腺癌患者预后不良高度相关。本研究通过单变量和多变量Cox 分析发现，结合焦亡风险评分、年龄、远程转移分期、TNM 分期的列线图具有优秀的预测能力，可作为一项有意义的个性化医疗指标。

综上所述，本研究构建了焦亡相关基因的乳腺癌预后模型，该模型有利于临床医生制订个性化的治疗方案，提高乳腺癌患者的远期生存率。