骨肉瘤中铜死亡相关亚型鉴定、预后模型构建以及免疫细胞浸润分析

邵子晨 李华南 章晓云，3 孙伟康 袁启鹏 （.江西中医药大学临床医学院，南昌 330004； .江西中医药大学附属医院，南昌 330004； 3.广西中医药大学附属瑞康医院骨科，南宁 5300）

骨肉瘤（osteosarcoma，OS）是一种原发性骨恶性肿瘤，易发生局部侵袭和早期转移，常发生于肺部，具有双峰年龄分布特征，儿童和青少年发病率最高，其次为60岁以上人群［1-3］。目前临床常见的OS治疗策略包括手术切除和联合放化疗，极大改善了OS患者预后，但转移性或复发性OS预后未取得突破性进展［4-5］。调查显示，局限性OS患者五年生存率约为60%，但转移或复发患者五年生存率仅为20%［5-6］。尽管既往研究强调了铁死亡、细胞凋亡、细胞自噬和细胞焦亡在OS中的重要作用，但铜死亡在OS中的作用仍不清楚。因此，铜死亡相关基因探索可为OS治疗、预后及靶向药物研发提供有效参考，指导OS患者临床决策。

哈佛大学研究团队在人类细胞中证实了一种新型细胞死亡方式——铜死亡，同时揭示了铜死亡通过铜离子与线粒体呼吸的三羧酸循环中的脂酰化成分直接结合而发生，导致脂酰化蛋白聚集和随后的铁硫蛋白下调，从而导致蛋白毒性应激并最终导致细胞死亡［7-8］。这些发现有助于探索铜失调疾病（威尔逊病）及研发新型肿瘤治疗策略。铜可通过介导ROS积累、线粒体功能障碍、蛋白酶体抑制和抗血管生成等多种机制抑制肿瘤细胞生长及诱导肿瘤细胞死亡［9］。随着铜在癌症治疗中的地位提高，铜配合物因其较强的治疗诊断能力和较弱的副作用逐渐成为一种新的抗癌策略，但其长期稳定性和生物安全性有待进一步验证［10-11］。铜在肿瘤发生发展中的死亡机制已成为新的研究热点。

本研究联合TARGET与GEO数据库综合探讨OS中铜死亡基因相关亚型鉴定、预后和免疫浸润影响，强调了铜死亡相关基因（cuprotosis related genes，CRGs）在OS预后中的重要性，为铜配合物在OS中的治疗应用奠定基础。

1 材料与方法

1.1 数据采集 转录组数据包含88个OS样本，数据来自TARGET（https：//ocg.cancer.gov/programs/target）数据库［12］。GSE16091芯片包含34个OS样本，数据来自GEO（https：//www.ncbi.nlm.nih.gov/geo/）数据库［13］。本研究收集了19个CRGs，包括NFE2L2、NLRP3、ATP7B、ATP7A、SLC31A1、FDX1、LIAS、LIPT1、LIPT2、DLD、DLAT、PDHA1、PDHB、MTF1、GLS、CDKN2A、DBT、GCSH和DLST。

1.2 OS相关CRG的生存分析及预后 运用R软件将TARGET数据与GEO数据进行合并及批次矫正，提取CRG表达。结合TARGET与GEO临床数据对CRG表达进行COX分析与KM分析，选择最优的KM分析对P＜0.05的基因（即OS预后相关基因）进行生存曲线绘制，对OS相关CRG进行预后相关性分析，保存P＜0.000 1的结果，绘制预后网络图。

1.3 铜死亡分型 结合临床数据对OS样本进行铜死亡分型。对铜死亡分型进行生存分析、相关性热图分析和GSVA分析，同时运用SSGSEA对铜死亡分型进行免疫细胞差异分析。

1.4 铜死亡分型差异分析和GO、KEGG富集分析 设置log2FC=0.585和adj.Pvalue=0.05，对铜死亡分型进行差异分析，得出铜死亡分型的差异基因。设置Pvalue=0.05和qvalue=0.05，对铜死亡分型的差异基因进行GO及KEGG富集分析。

1.5 基因分型 结合临床数据对OS样本进行基因分型。同时对样本数据进行单因素COX分析，得出单因素显著基因表达及预后情况，并对基因分型进行生存分析和相关性热图分析，分析基因分型在CRG中的差异表达。

1.6 预后模型构建及风险评估分析 将样品组数据分为训练组和测试组，进行Lasso回归分析及交叉验证。依据得到的特征基因构建预后模型，得到样品组、训练组和测试组风险评分及属性，以桑基图进行可视化展示。对样品组、训练组和测试组风险评分进行生存分析及风险曲线绘制。运用ROC曲线评估模型预测的准确性。

1.7 列线图构建及风险差异分析 结合临床性状绘制列线图预测患者生存期。对铜死亡分型与基因分型进行风险评分差异分析，同时分析CRG在高、低风险组的差异表达。

1.8 免疫细胞浸润与肿瘤微环境分析 对OS样本进行免疫细胞浸润分析；对参与模型构建的基因进行免疫细胞相关性分析。对OS样本进行肿瘤微环境分析，得出各OS样本基质细胞、免疫细胞、肿瘤纯度评分（其中基质细胞和免疫细胞综合评分即为ESTIMATE评分）。分析各肿瘤微环境评分类型在高、低风险组的差异表达。

1.9 药物敏感性分析 运用R软件“pRRophetic”包对OS样本进行药物敏感性分析。设置P=0.001作为过滤条件，对P＜0.001的结果进行可视化展示。

1.10 统计学分析 采用线性回归分析探索风险评分与免疫细胞间的相关性。K-W检验比较组间差异，所有参数分析均基于双尾检验，P＜0.05为差异有统计学意义，所有统计分析均使用R4.1.2版本。

2 结果

2.1 OS相关CRG生存分析及预后结果 设置P＜0.05作为过滤条件，KM分析筛选出10个OS预后相关基因，LIPT1与FDX1最为显著（图1）。DLST、LIPT1、NLRP3、SLC31A1、NFE2L2和LIAS基因高表达组预后优于低表达组预后；ATP7B、FDX1、DLAT和CDKN2A基因低表达组预后优于高表达组预后。对17个OS相关CRG进行预后相关性分析，设置P＜0.000 1为过滤条件，筛选出13个OS相关CRGs绘制预后网络图（图2），表明FDX1（为主）、GLS、DLAT和PDHB为预后的高风险基因，具有致癌特征。SLC31A1、ATP7A、NFE2L2、PDHA1、DBT、DLST、MTF1、DLD和LIAS为预后的低风险基因，具有抑癌特征。

图1 10个OS预后相关基因的生存曲线Fig.1 Survival curve of 10 OS prognosis related genes

图2 13个OS相关CRGs的预后网络图Fig.2 Prognostic network of 13 OS related CRGs

2.2 铜死亡分型 结合临床数据对OS样本进行铜死亡分型（图3），可将OS分为：CRGclusterA与CRG-clusterB。对铜死亡分型进行生存分析（图4A），CRGclusterA与CRGclusterB间OS患者生存差异显著（P＜0.001），且CRGclusterB生存率高于CRGclusterA，表明CRGclusterB预后优于CRGclusterA预后。对铜死亡分型进行相关性热图分析（图4B），大部分CRG在CRGclusterA中呈高表达。对铜死亡分型进行GSVA分析（图4C），造血细胞系、细胞因子受体相互作用和神经活性配体受体相互作用等在CRG-clusterB中呈高表达。p53信号通路、Mismatch_Repair信号通路、DNA复制、核苷酸切除修复、同源重组、卵母细胞减数分裂、细胞周期、剪接体、基础转录因子、碱基切除修复、嘧啶代谢、丙酮酸代谢、半胱氨酸及蛋氨酸代谢、蛋白质输出、硒氨基酸代谢和一个碳池等在CRGclusterA中呈高表达。运用ssGSEA对铜死亡分型进行免疫细胞差异分析（图4D），CRGclusterB中活化B细胞、NK细胞、骨髓来源抑制性细胞（MDSCs）、巨噬细胞、单核细胞和Th1细胞表达高于CRGclusterA，而Th2细胞在CRG-clusterA中的表达高于CRGclusterB。

图3 k=2～9的CRG共识矩阵Fig.3 Consensus matrix of CRG with k=2～9

图4 铜死亡分型的相关性分析Fig.4 Correlation analysis of copper death typing

2.3 铜死亡分型差异分析和GO及KEGG富集分析 设置log2FC=0.585和adj.Pvalue=0.05，对铜死亡分型进行差异分析，得到铜死亡分型差异基因246个。设置P=0.05和qvalue=0.05，对246个差异基因进行GO及KEGG富集分析。对BP、CC、MF排名前10的结果进行可视化（图5A）。GO富集显示主要涉及骨化、肾脏发育、肾系统发育、泌尿生殖系统发育、成骨细胞分化、细胞外基质组织、细胞外结构组织、外部封装结构组织和肾上皮发育等322条生物过程。KEGG富集显示主要参与人乳头瘤病毒感染、PI3K-Akt信号通路、Hippo信号通路、Wnt信号通路、TGF-β信号通路、p53信号通路、MAPK信号通路、库欣综合征、焦点粘连、乳癌、细胞黏附分子、癌症中的蛋白多糖、ECM-受体相互作用、调节干细胞多能性的信号通路、基底细胞癌、黑素生成、细胞周期、胃癌等（图5B）。

图5 DEG富集分析Fig.5 Enrichment analysis of DEG

2.4 基因分型 结合临床数据对OS样本进行基因分型（图6），可将OS分为两型：geneClusterA与gene-ClusterB。同时对样本数据进行单因素COX分析，得出单因素显著基因表达及预后情况，并对基因分型进行生存分析（图7A），geneClusterA与geneClusterB间OS患者生存具有显著差异（P＜0.001），且geneClusterB生存率高于geneClusterA，表明gene-ClusterB预后优于geneClusterA预后，与铜死亡分型情况一致。对基因分型进行相关性热图分析（图7B），大部分基因在geneClusterA中呈高表达，与铜死亡分型情况一致。分析基因分型在CRG中的差异表达（图7C），结果显示，FDX1、DLAT、PDHB在gene-ClusterA中的表达高于geneClusterB，而PDHA1在geneClusterB中的表达高于geneClusterA。

图6 k=2～9的246个差异基因共识矩阵Fig.6 Consensus matrix of 246 differential genes with k=2～9

图7 基因分型的相关性分析Fig.7 Correlation analysis of genotyping

2.5 预后模型构建及风险评估分析 将样品组数据分为训练组和测试组，进行Lasso回归分析（图8A），通过交叉验证（图8B）得到6个特征基因MYO6、IRX5、PRKX、CLDN11、MAGEA11、BAIAP2L2。依据6个特征基因构建预后模型，得到样品组、训练组和测试组风险评分及属性，以桑基图（图8C）进行可视化展示。对样品组（图9A）、训练组（图9B）和测试组（图9C）风险评分进行生存分析，结果表明高、低风险组间患者生存具有显著差异（P＜0.05），说明构建的模型可区分高低风险组患者，且低风险组患者预后优于高风险组患者。绘制样品组（图9A）、训练组（图9B）和测试组（图9C）风险曲线以及绘制1年、3年、5年ROC曲线评估模型预测准确性，结果显示样品组、训练组和测试组3年或5年AUC更大，构建的模型预测OS患者生存期准确性更高。因此构建的铜死亡相关风险特征与OS生存率显著相关。

图8 预后模型构建Fig.8 Construction of prognosis model

图9 风险评估分析Fig.9 Risk assessment analysis

2.6 列线图构建及风险差异分析 结合临床性状绘制列线图预测患者生存期（图10A）。其中1年、3年、5年校准曲线（图10B）离灰色线条很近，说明通过列线图预测1年、3年、5年生存期准确性较高。对铜死亡分型与基因分型分别进行风险评分差异分析（图10C、10D），结果显示CRGclusterA与CRG-clusterB风险评分具有显著差异（P＜0.001），且CRGclusterA评分高于CRGclusterB评分。geneClusterA与geneClusterB风险评分具有显著差异（P＜0.001），且geneClusterA评分高于geneClusterB评分。同时分析CRG在高低风险组差异表达（图10E），结果显示FDX1、LIPT1和DLAT在高风险组表达上调。

图10 列线图的构建及风险差异分析Fig.10 Construction of nomogram and risk difference analysis

2.7 免疫细胞浸润与肿瘤微环境分析 对OS样本进行免疫细胞浸润分析；对参与模型构建的基因（PRKX、MAGEA11和BAIAP2L2）进行免疫细胞相关性分析（图11A），结果显示γδ T细胞、静息肥大细胞和静息树突状细胞与风险评分呈正相关，而CD8 T细胞与风险评分呈负相关。说明在OS高风险中，γδ T细胞、静息肥大细胞和静息树突状细胞可能存在高表达；在OS低风险中，CD8 T细胞可能存在高表达。对OS样本进行肿瘤微环境分析，得出各OS样本基质细胞、免疫细胞、肿瘤纯度评分（其中基质细胞和免疫细胞的综合评分即为ESTIMATE评分）。分析各肿瘤微环境评分类型在高低风险组差异表达（图11B），表明基质细胞评分和ESTIMATE评分在高低风险组具有显著差异（P＜0.05）；ESTIMATE综合评分在高风险组中下调，说明高风险组肿瘤纯度更高。

图11 免疫细胞浸润与肿瘤微环境分析Fig.11 Analysis of immune cell infiltration and tumor microenvironment

2.8 药物敏感性分析 运用R软件“pRRophetic”包对OS样本进行药物敏感性分析，设置P=0.001作为过滤条件，对P＜0.001的结果进行可视化展示（图12），结果显示Elesclomol与GW.441756药物在高风险组IC50更低，说明高风险组患者对这两种药物敏感性更高。

图12 药物敏感性分析（Elesclomol与GW.441756）Fig.12 Sensitivity analysis of drugs （Elesclomol and GW.441756）

3 讨论

本研究在铜死亡基因基础上建立了一个OS预后模型，且基于铜死亡基因识别出CRGclusterA与CRGclusterB两种亚型，其中CRGclusterA与Th2细胞相关，CRGclusterB与Th1细胞相关。铜死亡是一种不同于任何已知细胞死亡机制的新型细胞死亡方式，依赖于线粒体呼吸［14］。这种不寻常的机制可能为OS治疗带来新的解决方案。

铁氧还蛋白是铜离子载体诱导细胞死亡的关键调节因子，也是蛋白脂酰化的上游调节器，可编码一种将Cu2+还原成Cu+的还原酶，是Elesclomol的直接靶点［15］。二氢硫辛酸转乙酰基酶是丙酮酸脱氢酶复合体的重要组成，其脂酰化依赖性寡聚由铜与脂酰化TCA循环蛋白结合导致［15-16］。谷氨酰胺酶（glutaminase，GLS）是一种调控谷氨酰胺代谢关键酶，在上调肿瘤生长的细胞代谢方面起重要作用。GLS与肿瘤生长率和恶性度相关，并受癌蛋白c-Myc调节，具有致癌特征［17-18］。丙酮酸脱氢酶B（pyruvate delydrogense B，PDHB）是一种位于线粒体内在氧化磷酸化中起重要作用的酶。PDHB作为糖酵解循环到柠檬酸循环的切入点用于能量生产，在多种人类恶性肿瘤进展和转移中发挥关键作用［19-21］。这些基因可能在不同癌症中分别发挥抑癌和促癌作用。本研究中，虽然FDX1、GLS、DLAT和PDHB作为OS预后高风险基因可能表现出致癌特征，但仍需要相关体外及动物模型实验明确其在OS中的作用特征。

大量研究表明肿瘤微环境在OS发生发展中发挥重要作用，其中肿瘤相关巨噬细胞（tumor associated macrophages，TAMs）广泛参与OS生长、血管生成、转移和干细胞样表型［22-23］。巨噬细胞被诱导分化为M1和M2两种类型，M1型参与Th1型免疫应答并发挥促炎作用，M2型参与Th2型免疫应答并发挥抑炎作用［24］。肿瘤发展过程中，TAM倾向于极化为M2型TAM［25］。激活Th2型免疫应答，可导致免疫逃逸发生［24］。目前，在包括肝癌、肺癌、胃癌在内的多种恶性肿瘤中发现高度浸润的M2型TAM［26-27］。REN等［24］研究发现单核细胞-巨噬细胞通过与OS S180细胞共培养被证明极化为M2型TAM。本研究中，大部分OS铜死亡基因在CRGclusterA中呈高表达。CRGclusterA与Th2细胞相关且M2型巨噬细胞参与Th2型免疫应答，表明在OS中，巨噬细胞可能更倾向于极化为M2型，进而可能增加Th2型免疫应答，导致肿瘤细胞生长和免疫逃逸，与REN等［24］研究结果相符。今后可能提供一种干扰巨噬细胞极化抑制肿瘤的方法，其中抑制M2巨噬细胞表型或阻断M1巨噬细胞向M2极化的治疗在多项研究中具有积极结果［23］。

目前γδ T细胞的抗肿瘤作用已在体内外许多研究中得到证明，主要基于其高细胞毒性和γ干扰素形成，在主要组织化合物中具有无限溶解度活性，与各种肿瘤具有高度相容性［28-29］。最近研究显示，小鼠和人类产生IL-17的γδ T细胞被发现具有意外的肿瘤促进作用［30］。可见γδ T细胞对肿瘤生长具有双重作用。已知CD8 T细胞可直接杀死癌细胞［31］。肥大细胞的高浸润与许多癌症预后不良、生存率低和转移增加有关。肥大细胞可通过生成血管生成因子和产生蛋白酶促进肿瘤生长［32-33］。研究表明，暴露于肿瘤抗原的树突状细胞通过诱导肿瘤特异性T细胞反应在抗肿瘤免疫应答中发挥关键作用［34］。以树突状细胞为基础研制治疗黑色素瘤和肾细胞癌的疫苗已在免疫和临床反应中得到了证实［35］。本研究发现高水平的CD8 T细胞与OS良好预后相关，而高水平的γδ T细胞、静息肥大细胞和静息树突状细胞与OS不良预后相关，与上述结论相矛盾［34］，即高水平的树突状细胞与OS不良预后相关。需要更多实验证明上述免疫细胞在骨肿瘤中具有特定的抗肿瘤或促肿瘤作用。基于免疫细胞的免疫治疗可能促进OS患者临床预后改善，尤其是针对OS复发者和转移者，其安全性和有效性有待进一步临床试验。

Elesclomol是一种强效铜离子载体，可促进特定肿瘤细胞凋亡，且作为潜在抗癌药物进入了临床试验，但具体分子作用机制尚不清楚［8，36］。编码Elesclomol靶蛋白的FDX1基因是促进铜死亡的基因［8］。本研究药物敏感性分析显示，OS细胞对Elesclomol与GW.441756更敏感。铁氧还蛋白与OS预后显著相关。因此，Elesclomol可治疗OS患者，可能与OS细胞更易受铜死亡影响有关。

综上，本研究系统分析了OS中CRG相关亚型鉴定、预后和免疫浸润的影响，识别出CRGclusterA与CRGclusterB两种亚型，鉴定出4个预后高风险基因FDX1、GLS、DLAT和PDHB，为OS预后评估和免疫治疗候选靶点提供了新见解。