基于铜死亡相关基因的胶质瘤预后模型构建

张博弘，张祺

胶质瘤是最常见的原发性中枢神经系统肿瘤，约占颅内恶性肿瘤的80%，根据组织和分子特征，可以分为WHO Ⅰ~Ⅳ级。胶质母细胞瘤是恶性程度最高的胶质瘤，尽管采取了包括最大限度的手术切除、放疗和化疗的标准治疗措施，其中位生存期仍只有大约16 个月。WHO Ⅱ~Ⅲ级的胶质瘤也具有侵袭性生长和进展到更高级别的恶性行为，大部分病人最终死于这种疾病［1，2］。探索新的预后模型，识别高危患者，从而精准治疗，是提高胶质瘤临床预后的有效措施，除了临床病理特征外，个体基因特征可望构建有效的预后模型。

铜死亡（cuproptosis），是2022 年3 月Tsvetkov等在《Science》上首次提出的一种铜依赖性的、全新的细胞死亡方式，不同于已知的细胞凋亡、坏死性凋亡、焦亡和铁死亡，它是通过铜离子与线粒体呼吸中的三羧酸循环中的脂酰化成分结合，使得脂酰化蛋白质异常聚集，铁硫簇蛋白下调，导致蛋白质毒性应激并最终引起细胞死亡［3］。此外，研究者还初步鉴定出12 个铜死亡相关基因（CRGs），包括：FDX1，LIAS，LIPT1，DLD，DLAT，PDHA1，PDHB，ATP7B，MTF1，GLS，CDKN2A 和SLC31A1［3］。近半年来，在多种肿瘤中报道了基于铜死亡相关基因构建的预后模型，准确预测患者预后、肿瘤免疫微环境以及对药物的治疗反应，其中包括：肾透明细胞癌［4］、肝细胞癌［5］、黑色素瘤［6］、头颈部鳞状细胞癌［7］、肺癌［8］和乳腺癌［9］。然而，铜死亡相关基因在预测胶质瘤患者预后中的作用仍不清楚。此外，最近，基于新型程序性细胞死亡（焦亡、铁死亡）相关基因构建的胶质瘤预后模型被相继报道，但，预测效能仍不能令人满意。铜死亡是最新报道的程序性细胞死亡方式，探索基于铜死亡相关基因的预后模型有望获得预测效能更高的预后模型。为探讨CRGs 对胶质瘤患者预后预测的价值，本研究拟通过分析公共数据库中胶质瘤患者的转录组表达谱和相应的临床数据，基于CRGs 构建并验证预后风险评分模型，进一步构建诺莫图提高模型临床适用性，为改善胶质瘤患者预后评估以及精准治疗提供新的思路。

1 资料与方法

1.1 数据收集

从GTEx（https：//www.gtexportal.org/）数据库中下载正常脑组织的转录组表达谱数据，共纳入1157 个正常脑组织样本。从TCGA（The Cancer Genome Atlas）（https：//protal.gdc.cancer.gov/）和GEO（Gene Expression Omnibus）（https：//www.ncbi.nl m.nih.gov/geo/）数据库中下载胶质瘤患者的转录组表达谱和相应的临床数据。获得了一个TCGA 胶质瘤队列和一个GEO 胶质瘤队列（GSE43378）用于后续研究。每个数据集采用“limma”包的normalize BetweenArrays 方法进行背景调整和分位数归一化处理，通过“Combat”算法消除批次效应后合并两个数据集。通过删除缺乏生存数据的样本，共纳入720 个胶质瘤样本进行研究。临床变量包括年龄、性别、分级、总生存时间、生存状态。

1.2 铜死亡相关基因分型

从文献报道中提取出12 个铜死亡相关基因，根据预后相关CRGs 的表达谱，通过一致性聚类分析把胶质瘤样本分为不同的铜死亡亚型：CRGcluster A 和CRGcluster B。一致性聚类是一种无监督聚类方法，利用“ConsensusClusterPlus”R 包完成分析，本研究使用的聚类算法是K-means 聚类算法，距离是基于欧氏距离（euclidean）。聚类的标准如下：第一，累积分布函数（CDF）曲线逐渐平缓增加；第二，没有一个组的样本量很小；最后，聚类后，组内相关性增加，组间相关性降低。接着，利用prcomp 函数进行主成分分析（PCA）。然后，采用R 包“survival”和“survminer”，通过Kaplan-Meier法对不同亚型进行生存分析。利用R 包“limma”，筛选CRGcluster A 和CRGcluster B 的差异表达基因（DEGs），筛选条件为FC>1.5 和adjustedP<0.05，进一步通过单因素Cox 回归分析鉴定出1831 个预后相关的DEGs。

1.3 构建预后风险评分模型

首先，按照1∶1 比例，把胶质瘤样本随机分为训练集（n=360）和测试集（n=360）。在训练集中，基于预后相关的DEGs，使用R 包“glmnet”，通过Lasso Cox 回归分析建立预后模型，并根据最小偏似然偏差确定预后模型。采用多因素Cox 分析选择候选基因，得到风险评分计算公式如下：

风险评分=Σ（Expi*Coefi）

其中，Expi 和Coefi 分别表示各基因的表达量和风险系数。根据风险评分中位数将患者分为低风险组和高风险组，并对两组患者进行Kaplan-Meier 生存分析。利用R 包“timeROC”建立时间依赖ROC 曲线以评估预后模型的预测性能。

1.4 构建诺莫图评分系统

采用R 包“rms”，根据临床特征和风险评分构建预测预后的诺莫图。在诺莫图评分系统中，每个变量都对应一个分数，通过将所有变量的分数相加得到总分。校准图用于描述1、3 和5 年生存率的预测值与实际观察到的结果之间的差异。

1.5 统计学分析

所有统计分析均采用R（Version 4.1.3）进行。如果未特别注明，P<0.05 被认为具有统计学意义。

2 结 果

2.1 胶质瘤中铜死亡相关基因的差异表达和预后分析

从TCGA、GEO 和GTEx 数据库中下载胶质瘤和正常对照样本的转录组数据和相应的临床数据，共获得1877 个样本，其中胶质瘤样本720 个，正常对照样本1157个。本研究共纳入12个CRGs，包括：FDX1，LIAS，LIPT1，DLD，DLAT，PDHA1，PDHB，ATP7B，MTF1，GLS，CDKN2A 和SLC31A1。通过比较胶质瘤和正常组织的CRGs 的mRNA表达水平差异，发现10 个CRGs 在胶质瘤显著上调，2 个CRGs 显著下调，差异都具有统计学意义（图1）。进一步通过单因素Cox回归和Kaplan-Meier分析发现，包括ATP7B、SLC31A1、FDX1、LIAS、LIPT1、DLD、DLAT 和PDHA1 在内的8 个CRGs 具有预后意义（表1）。提示CRGs 在胶质瘤发生、发展中具有潜在的作用。

表1 胶质瘤中8 个具有预后意义的CRGs

图1 胶质瘤铜死亡相关基因差异表达分析 胶质瘤和正常组织的CRGs 的mRNA 表达水平，比较采用Wilcoxon rank sum test，****P<0.0001 vs.正常组

2.2 胶质瘤铜死亡亚型的鉴定

根据8 个预后相关CRGs 的表达谱，采用一致性聚类算法，对胶质瘤样本进行分型。结果表明，K=2 似乎是最佳选择，可以把整个队列分成CRGcluster A 和CRGcluster B（图2A-B），且主成分分析（principal component analysis，PCA）表明两个亚型之间具有较好的区分度（图2C）。进一步生存分析，Kaplan-Meier 曲线显示CRGcluster A 的胶质瘤患者具有较长的总体生存时间（OS）（图2D）。采用R 包“limma”比较不同铜死亡亚型的基因表达谱，得到差异表达基因（DEGs），进一步通过单因素Cox 回归分析获得1831 个预后相关的DEGs。

图2 胶质瘤铜死亡亚型鉴定 A、B：一致性聚类分析K=2 时的一致性矩阵热图，同时一致性累积分布函数图提示K=2 时CDF 下降坡度小，最佳的K 值为2；C：主成分分析显示两个铜死亡亚型之间转录组差异显著；D：KM 曲线，CRGcluster A 类患者的OS 时间更长，比较采用chi-square test：χ2=38.04，P<0.001

2.3 构建预后风险评分模型

采用R 包“caret package”把样本按1∶1 随机分成训练集（n=360）和测试集（n=360）。在训练组集中，对1831 个预后相关的DEGs 进行LASSO Cox回归分析，根据最小偏似然偏差确定了8 个基因的预后模型（图3A~B）：风险评分=MGME1×0.4304+GPR156×0.2167+ARNTL×0.3229+RAB38×0.2424+IGFBP5 × 0.1675 + SERPINA5 × 0.1579 + FABP5 ×0.1809+COL22A1×0.1415。计算训练集中所有患者的风险评分，根据中位风险评分，把患者分为高风险组（180 例）和低风险组（180 例），风险评分的分布图显示，随着风险评分的增加，生存时间减少，死亡率增加（图3C），Kaplan-Meier 生存曲线显示，低风险组患者总体生存期较长（图3D）。此外，接受者操作特征曲线（ROC）显示，曲线下面积（AUC）分别为：0.872（1 年）、0.933（3 年）、0.889（5 年），表明该模型具有良好的预测性能（图3E）。

图3 在训练集中构建预后风险评分模型 A~B：LASSO 回归系数分布图和10 倍交叉验证选取最优λ 值；C：训练集中患者风险评分、总体生存时间、生存状态的分布；D：Kaplan-Meier 生存曲线；E：模型的ROC 曲线及AUC 值

2.4 验证风险评分模型

计算测试集中所有样本的风险评分，根据中位风险评分，把测试集样本分为高风险组（183 例）和低风险组（177 例），与训练集中结果一致，风险评分的分布图显示，随着风险评分的增加，生存时间减少，死亡率增加（图4A），Kaplan-Meier 生存曲线显示，低风险组患者总体生存期较长（图4B），此外，ROC 曲线显示AUC 分别为：0.881（1 年）、0.894（3 年）、0.874（5 年），表明该模型对测试集胶质瘤患者的同样具有良好的预测能力（图4C）。

图4 在测试集中验证风险评分模型 A：测试集中患者风险评分、总体生存时间、生存状态的分布；B：Kaplan-Meier 生存曲线；C：模型的ROC 曲线及AUC 值

2.5 构建诺莫图预测胶质瘤患者生存率

考虑到风险评分在临床应用上缺乏便利性，我们建立了一个包含风险评分和临床参数的诺莫图，用于预测胶质瘤患者1 年、3 年和5 年的生存率（图5A），其中预测因素包括患者的性别、年龄和风险评分，同时校准图表明，本诺莫图具有理想模型相似的预测性能（图5B）。

3 讨 论

越来越多的证据显示，同样WHO 级别的胶质瘤，其生物学行为和临床预后却相差甚远，有的患者能长期存活或具有较长的生存期，而有的患者却表现为高度恶性的转归，因此准确的预后评估并选择性地给予不同的治疗策略是十分必要的。目前临床应用的胶质瘤预后风险分层标准主要有美国的RTOG9802 标准［10］和欧洲的EORTC22844标准［11］。然而，这两个标准仍存在不足的地方，有时候甚至出现较大的偏差。随着胶质瘤基因组测序的普及，海量的生物数据提供了更多的预后信息，补充甚至模糊了传统的WHO 分级标准。近几年来，许多研究者借助公共数据库，基于不同的基因表达谱，通过多种分析方法，在不同肿瘤中构建了预后模型。Chen 等［12］基于11 个铁死亡相关基因表达谱构建了一个胶质瘤预后风险评分模型，其1 年的AUC 值是0.879。Zhang 等［13］基于10 个焦亡相关基因表达谱构建了一个胶质瘤预后风险评分模型，其1年的AUC值是0.872。Feng等［14］基于6个自噬相关基因表达谱构建了一个低级别胶质瘤预后风险评分模型，其1 年的AUC 值是0.840。

2022 年3 月，Tsvetkov 等首次在《Science》上提出了一种具有铜依赖性的、全新的细胞死亡方式-“铜死亡（cuproptosis）”，研究指出Cu2+直接与三羧酸循环中的脂酰化成分结合，导致脂酰化蛋白异常聚集和铁硫簇蛋白丢失，进而引发蛋白质毒性应激反应，最终致使细胞死亡［3］，并初步列出了12 个铜死亡相关基因（CRGs）。受此启发，基于CRGs 构建的预后模型在多个肿瘤类型相继被报道［4-9］。为构建基于CRGs 的胶质瘤预后模型。本研究首先探讨了12 个CRGs 在肿瘤和正常样本的差异表达情况和预后意义，基于8 个预后相关CRGs 的表达谱，对胶质瘤样本进行铜死亡分型（CRGcluster A 和CRGcluster B），然后进行铜死亡分型的差异基因分析和预后分析，得到1831 个预后相关差异基因，接着通过LASSO 回归分析，在训练集中构建了一个包含8 个基因的预后风险评分模型，其1 年的AUC 值是0.872，并在测试集中进行了验证，其1 年的AUC 值是0.881，优于已报道模型的预测性能。

本研究所构建的模型包含以下8 个基因：MGME1、GPR156、ARNTL、RAB38、IGFBP5、SERPINA5、FABP5、COL22A1。MGME1（线粒体基因组维持外切酶1）对线粒体DNA（mtDNA）的有效合成至关重要，其功能丧失性突变会影响mtDNA 的正常维持，并导致线粒体疾病［15］。ARNTL 是控制昼夜节律的时钟基因之一，其高度甲基化可以促进鼻咽癌的发生并抑制其对顺铂的敏感性［16］。RAB38 属于小GTP 结合蛋白家族，通过调控线粒体代谢以及线粒体介导的细胞死亡促进胶质母细胞瘤恶性进展［17］。IGFBP5 是胰岛素样生长因子结合蛋白家族的一员，参与线粒体功能的维持，沉默IGFBP5 的表达会导致神经母细胞瘤细胞线粒体凋亡的发生［18，19］。SERPINA5 属于丝氨酸蛋白酶抑制剂超家族，可以作为低级别胶质瘤的预后因子，补充当前基于组织学的肿瘤分类系统［20］。FABP5 是脂肪酸结合蛋白家族的重要成员，在脂肪酸代谢中起着重要作用，一方面，FABP5 的缺失将损害心脏细胞的线粒体功能［21］；另一方面，FABP5 参与氧化应激时α 突触核蛋白诱导的线粒体损害［22］。相比之下，关于GPR156 和COL22A1的研究相对较少，与肿瘤的关系还不是很清楚。模型的大部分基因都涉及线粒体功能的调控，一致的是，铜死亡的发生与线粒体密切相关，这也反映了本研究结果的可靠性。

值得注意的是，本研究存在一些局限性。首先，本研究在训练集构建模型，通过测试集进行验证，但这一模型并没有在临床队列中得到进一步的验证；其次，铜死亡与胶质瘤的关系需要进一步通过基础实验验证；最后，本诺莫图仅包含了风险评分和3 个临床参数（性别、年龄和胶质瘤WHO分级），如果能纳入更多的参数，将进一步提高预测性能，例如IDH1 突变情况，1p19q 缺失情况和MGMT 甲基化情况，这些都是与胶质瘤预后和治疗反应密切相关的参数，但是由于公共数据库并没有详细提供这些数据，导致没办法进行相关分析，下一步我们将收集自己医院的胶质瘤样本，详细收集这方面的数据，再进行相关的分析。