构建和验证基于铁死亡相关基因的骨肉瘤预后评估模型

罗辉,孙恒昌（中山大学附属第三医院,广东 广州 510630）

骨肉瘤（Osteosarcoma）是儿童和青少年最常见的原发性的恶性骨肿瘤[1-2]。尽管骨肉瘤的治疗方法如化疗、新辅助化疗、放疗及免疫治疗等取得长足的进步，但是骨肉瘤患者的5年生存率依然处在较低水平。主要是由于大部分骨肉瘤会发生转移[3]，有随访数据显示，未发生肺转移的骨肉瘤患者5年生存率约为70%，而发生肺部转移的骨肉瘤患者5年生存率仅有30%[4]。骨肉瘤发生转移的机制并不十分明确，且缺乏早期诊断和预后评估的生物标志。目前，其依然是青少年发生癌症相关死亡的第二大因素。因此，对骨肉瘤患者预后进行早期评估对骨肉瘤患者的预后以及治疗显得尤为重要。

铁死亡（Ferroptosis）是新发现的细胞程序性死亡形式之一，主要特点是铁依赖的脂质过氧化物损伤诱导细胞的死亡[5]。研究证实铁死亡参与了肿瘤的发生发展、侵袭、转移、耐药性与肿瘤免疫等各个方面。Harris[6]等学者研究表明通过铁死亡激活剂能够促进小细胞肺癌、结肠癌、乳腺癌、淋巴瘤等肿瘤细胞的死亡。此外，也有越来越多的研究利用公共数据库中的数据结合生物信息分析构建基于铁死亡相关基因预后评估模型[7-10]。例如Huang X[10]等学者基于12个铁死亡相关基因构建了急性髓细胞白血病的预后评估模型，其1年、3年、5年总体生存率ROC曲线下面积（AUC）均可达到0.8以上。然而，目前铁死亡相关基因在骨肉瘤预后评估中的作用并不明确。

本文从公共数据库中获取了骨肉瘤患者的组织基因表达数据，探讨铁死亡相关基因对于骨肉瘤患者预后评估的潜在价值。

1 材料与方法

1.1 数据来源和处理 从公共数据库（https://xenabrowser.net/datapages/）中下载88例骨肉瘤患者高通量测序基因表达数据（GDC-TARGET-OS，HTSeq-FPKM（n=88）），同时下载相关的临床信息数据。用R语言软件进行表达基因名转换以及整合临床信息，删除生存时间小于30天的样本，最终获得84例骨肉瘤患者数据。从数据库（http://www.datjar.com:40013/bt2104/）中获取铁死亡相关的基因（Ferroptosis-related genes，FRGs，共259个基因），并提取骨肉瘤患者铁死亡相关基因的表达信息。

1.2 获取预后相关的铁死亡基因 结合骨肉瘤患者铁死亡相关基因表达数据和生存数据，用“survival”包对FRGs进行单因素COX生存分析，以P＜0.05为条件筛选与骨肉瘤患者预后相关的FRGs。

1.3 预后模型构建 将骨肉瘤患者的基因表达数据随机抽取80%（n=67）为模型训练集，剩余的20%（n=17）为避免单因素COX分析筛选的预后相关FRG出现过度拟合，采用R软件“glmet”包进行进一步的LASSO回归分析，进一步筛选预后基因。LASSO回归中采用十则交叉验证确定λ值，选择偏似然偏差最小的λ（Lamda.min）为最佳λ。用最终筛选的基因构建预后模型，预后风险评分计算公式如下：Riskscore=C1*Expr1+C2* Expr2+……Cn*Exprn（n代表基因个数，C表示LASSO回归系数，Expr表示基因表达量）。根据上述公式计算出骨肉瘤患者的风险评分，风险评分大于中位数的骨肉瘤患者为高风险组，小于中位数为低风险组。绘制患者风险评分的分布图以及风险评分与生存时间的关系图。同时比较纳入模型的几个铁死亡相关基因在高低风险两组间的表达水平。

1.4 预后模型的预测效果验证 分别采用R软件的“survival”包和“survivalROC”包对高低风险两组进行K-M生存分析以及ROC曲线分析，计算1年、3年、5年总体生存率，对模型的预测效果进行评估。同时采用验证集数据做K-M生存分析和ROC曲线分析，评估预后模型效果。

1.5 统计学处理方法 采用R3.5.1进行数据处理和统计分析。组间比较用独立样本t检验；采用单因素COX分析筛选预后因子；分别采用“glmnet”“survival”“survivalROC”等R包进行LASSO回归分析、生存分析、及时间依赖的ROC曲线绘制。生存资料单因素及多因素分析采用COX比例风险模型。当P＜0.05时，表示差异具有统计学意义。

2 结果

2.1 与骨肉瘤预后相关的关键铁死亡相关基因的筛选 经过单因素COX生存分析筛选，共有11个基因与骨肉瘤的总体生存率相关（P＜0.05）（见表1）。他们分别是MUC1、VEGFA、MAP3K5、HILPDA、G6PD、ARNTL、SCD、BNIP3、PML、SOCS1、CBS。

表1 单因素COX生存分析筛选出的与骨肉瘤预后相关的铁死亡相关基因

2.2 LASSO回归分析进一步筛选预后相关关键基因 根据上述单因素COX分析筛选P＜0.05的11个基因（见表1）及训练集样本数据进一步进行LASSO回归分析，采用交叉验证迭代分析，结果显示当变量个数为9时，模型的均方误差最小（λ=0.032）（见图1AB）。构建基于9个铁死亡相关基因的骨肉瘤预后评估模型：Risk score=（1.496）*CBS表达量+（0.456）*MUC1表达量+（0.096）*VEGFA表达量+（0.483）*HILPDA表达量+（-0.226）*SOCS1表达量+（-0.263）*PML表达量+（-0.504）*MAP3K5表达量+（-1.288）*G6PD表达量+（-2.40）*ARNTL表达量。

2.3 预后模型的预测效果 根据上一步得到的模型计算训练数据集中的每个样本的风险值（Risk score），按风险值的中位数（-5.51）为界，将训练集数据分为高风险组和低风险组，同时绘制两组患者的生存状态分布图、基因表达热图（见图2A）。训练数据集中，1年、3年、5年总体生存率的ROC曲线下面积AUC分别为0.928、0.964、0.947（见图2B）。生存分析显示高风险组（risk score＞-5.51）骨肉瘤患者总体生存率显著低于低风险组（risk score＜-5.51）（P＜0.01）（图2C）。而在验证数据集中，预后模型1年、3年、5年总体生存率的ROC曲线下面积（AUC）比训练数据集中要低，但是AUC也都分别达到0.937、0.855、0.896（见图3A-C）。且在验证数据集中生存分析也显示高风险组的生存率明显低于低风险组（P＜0.01）。说明铁死亡相关的预后模型对骨肉瘤患者预后的预测效果较好。

2.4 预后基因在高低风险两组间的差异表达分析 根据模型计算的风险值，将骨肉瘤患者分为高风险组和低风险组，分别计算两组患者之间9个预后因子的基因表达量差异。结果高风险组中CBS、HILPDA、VEGFA、MUC1基因的表达量显著高于低风险组（P＜0.05）。而SOCS1、PML、ARNTL、G6PD、MAP3K5在高风险组的表达量则显著降低（P＜0.05）（见图4）。

3 讨论

恶性骨肿瘤具有病情变化快、病死率高等特点[11]。骨肉瘤的发病机制是一个复杂的多步骤和多因素的过程，其中涉及广泛的分子异常和肿瘤异质性[12]。由于恶性程度高、转移早、耐药，致残率高、死亡率高，骨肉瘤的预后极差。尽管人们已经尝试了很多新的治疗药物和方法，但骨肉瘤的治疗效果并不理想，患者5年生存率仍然较低[13]。目前骨肉瘤的治疗前景并不乐观，而预后情况的预测，能有效地掌握疾病的发展趋势，对临床用药具有一定的指导作用。自从2012年铁死亡的概念被提出以来，其在肿瘤及非肿瘤疾病中的作用越来越受到关注，其在骨肉瘤的发生发展及治疗等方面也有报道。

本文从公共数据库中下载了骨肉瘤患者的基因表达数据，并且提取出259个FRGs的表达数据。采用单因素COX对现有患者的生存数据分析，筛选出11个与骨肉瘤预后相关的FRGs。进一步用LASSO回归筛选，最终获得9个FRGs用于骨肉瘤患者预后模型构建。该预测模型经过训练集和验证集的验证，发现在预测骨肉瘤患者1年、3年、5年生存率中具有较高的AUC值，与Liu[13]等学者构建的基于5个缺氧相关基因的骨肉瘤预后模型的AUC相当，而与Huang[14]等学者基于m6A相关基因构建的预后模型相比，本模型AUC明显更大，提示本模型具有比其他模型更好的预后预测能力。

本研究中建立的预后模型主要是由9个（PML、MAP3K5、ARNTL、CBS、HILPDA、VEGFA、MUC1、SOCS1、G6PD）铁死亡相关基因构成。其中PML（promyelocytic leukemia protein）是一种锌指转录因子，最早发现于急性早幼粒细胞白血病（APL），在维持基因组稳定性中发挥重要作用[15]。研究发现，PML的缺失会导致ROS水平增加[16]。PML还是一种促凋亡的基因，在非小细胞肺癌中，可被基于顺铂的全身化疗间接抑制[17]。此外，PML在骨肉瘤细胞株中是抑制基因，已被证明与癌蛋白MDM2相互作用，从而调控肿瘤细胞的生物学行为[18]。MAP3K5（又称为ASK1）是丝裂原活化蛋白激酶，属于激酶（MAP3K）家族的成员之一，参与癌症、神经退行性变、炎症和糖尿病的发病机制[19]。据报道[20-21]，在MAP3K5通路激活后能够促进骨肉瘤细胞的死亡。ARNTL是细胞昼夜节律的调控分子，作为生物钟基因调节着癌症的发生发展和化疗的耐受性。昼夜节律的失调让心血管疾病、免疫系统疾病和肿瘤细胞发生的风险增加[22]。据报道[23]，在乳腺癌和恶性血液病中检测到AENTL启动子的高甲基化，进一步加强了恶性细胞昼夜节律的紊乱。而本文研究结果中，低风险组患者的PML、MAP3K5、ARNTL几个基因的表达量均显著高于高风险组患者，这与前人的研究结果相符合。上述基因的表达上调，提示骨肉瘤患者预后较好。

胱硫醚β合酶基因（cystathionine-beta-synthase，CBS）是编码胱硫醚-β-合成酶的基因，它参与同型半胱氨酸的转硫途径的第一步反应，将同型半胱氨酸不可逆的转化为半胱氨酸和α酮丁酸。有研究显示抑制CBS能够诱导肿瘤细胞发生铁死亡，从而影响肿瘤的预后及治疗效果[24]。MUC1是一种Ⅰ型跨膜糖蛋白，一般情况下在胰腺、乳腺、肺和胃肠道等上皮细胞的顶端表面低水平表达[25]。MUC1表达量上调，除了影响肿瘤侵袭以及不良预后外，还与肿瘤血管生成和化疗药物的耐药性相关[26]。研究发现[25]，MUC1/xCT信号通路可抑制铁死亡，从而促进肿瘤细胞存活。VEGFA是一种介导血管生成的因子，黏附于血管内皮上，能够促进细胞增殖以及增加血管通透性，在肿瘤血管生成中发挥重要作用[27]。VEGFA基因上调，影响着骨肉瘤的发展与转移，测量VEGFA的循环水平对于评估预后有一定的价值意义[26]。HILPDA（hypoxia-inducible lipid droplet-associated）在细胞限制氧气供应时，通过代谢重编程调节脂解，使每个消耗氧气的分子产生三磷酸腺苷（ATP）的比例更高，以减少活性氧的生成[28-29]。在肿瘤细胞中，HILPDA基因的过表达，通常预示患者预后不良[30]。本研究中，上述基因的过量表达均提示着患者的预后质量极差，治疗后的总体生存率也显著降低。

细胞因子信号抑制因子1（Suppressor Of Cytokine Signaling-1，SOCS1）是SOCS家族的一员，在DNA损伤部位中连接P53，通过介导P53信号通路促进细胞凋亡。SOCS1基因随着DNA甲基化介导的沉默，在人类癌症中的表达量减少[31]。葡萄糖-6-磷酸脱氢酶（Glucose-6-phosphate dehydrogenase，G6PD）是磷酸戊糖途径的限速酶，生成的NADPH对维持细胞氧化还原稳态和还原生物合成至关重要，癌细胞NADPH的生成和糖酵解都依赖于G6PD，在肾细胞癌[32]、膀胱癌[33]、胃癌[34]等多种人类癌症中都观察到G6PD活性的升高。除此之外，异常激活G6PD还可引起ROS升高，导致细胞过度生长和分化失控[35]。G6PD基因的过表达，密切联系着恶性肿瘤患者的病情进展，已经作为患者预后不良的指标[34]。本研究发现高风险组中SOCS1和G6PD的表达量显著降低，根据前人研究提示前者与预后差相关，后者反而提示良好的预后，可能不同肿瘤中G6PD发挥的作用并不相同，值得进一步研究。

本文也存在一些局限：构建的预后评估模型训练数据病例数较少。其次，该模型缺乏其他患者人群数据的外部验证。因此，模型的外推性还有待进一步验证。

综上所述，本研究建立了基于9个铁死亡相关基因的骨肉瘤预后评估模型，为骨肉瘤患者的个性化治疗提供预测因子，研究的结果可能有助于骨肉瘤患者的预后评估，同时有助于发现骨肉瘤潜在的治疗靶点。