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

范以东，秦刚，刘金富，苏国威，肖世富，刘俊良，李威材，吴广涛

1广西中医药大学研究生院，南宁 530222

2广西中医药大学第一附属医院骨病创伤骨科，南宁 530022

骨肉瘤（osteosarcoma，OS）是一种主要由间充质细胞引起的恶性肿瘤，主要发生在长管状骨中[1]。OS 是发病率仅次于多发性骨髓瘤的第二常见的原发性骨恶性肿瘤，占所有骨肿瘤的15%[2]。随着医疗水平的提升，OS 患者的5 年生存率可达到70%，但治疗效果不佳且转移的患者，其5 年生存率低于30%[3]。虽然新辅助化疗的出现提高了OS 患者的5 年生存率，但OS 患者的生存情况仍然较差。二代测序（next-generation sequencing，NGS）技术的出现彻底改变了临床研究、基础医学和应用医学，通过NGS 技术有可能发现新的预后生物标志物，为OS 患者的个体化治疗提供依据，有助于预测患者预后和提高治疗效果[4]。因此，确定新的靶点至关重要。铁死亡是一种铁依赖性，区别于坏死和凋亡的细胞死亡模式，其特征在于脂质活性氧（reactive oxygen species，ROS）的积累[5]，由Stockwell 等于2012 年首次提出[6]。铁死亡在细胞形态学和生物学方面各有特点，形态学方面表现为线粒体体积变小、双层膜密度增加、线粒体嵴减少或消失，生物学方面主要表现为谷胱甘肽过氧化物酶4（glutathione peroxidase 4，GPX4）失活，导致脂质过氧化物积累，Fe2+氧化脂质产生ROS，从而发生铁死亡[7]。研究表明，铁死亡可以影响多种肿瘤的发生发展，包括结直肠癌[8]、胶质母细胞瘤[9]、乳腺癌[10]、胃癌[11]等。本研究应用生物信息学方法通过公共数据库探讨铁死亡相关基因（ferroptosis related gene，FRG）在OS 中的表达及与患者预后的关系，并分析FRG 在OS 中的作用机制，构建FRG相关的OS 预后模型，旨在为OS 的诊断、治疗及预后评估提供参考依据，现报道如下。

1 材料与方法

1.1 数据来源与整理

从UCSC Xena 数据库（http://xena.ucsc.edu/）中获取88例OS患者的转录组测序数据和其中85例患者的临床资料（2 例患者生存状态不明，1 例患者总生存时间不明，故予以删除）。在基因型-组织表达（Genotype-Tissue Expression，GTEx）数据库（https://www.gtexportal.org/home/）中获取396例正常骨组织样本，对所有GTEx 数据进行log2（x+1）转换。使用“limma”包对数据进行合并和矫正。从FerrDb数据库（http://www.zhounan.org/ferrdb）中下载FRG，所有数据下载截止日期为2022年12月12日。

1.2 获取差异表达的FRG

以错误发现率（false discovery rate，FDR）＜0.05，|log2FC|≥2 为过滤条件筛选正常骨组织与OS组织中的差异表达基因。从筛选的差异基因列表中进一步提取差异表达的FRG。

1.3 基因富集分析

利用R 语言“org.Hs.eg.db”“enrichplot”包对差异表达的FRG 进行基因本体（Gene Ontology，GO）和京都基因与基因组百科全书（Kyoto Encyclopedia of Genes and Genomes，KEGG）富集分析，探索OS 中FRG 的生物学功能。

1.4 预后模型的构建

采用单因素Cox 回归模型分析差异表达的FRG 对OS 患者总生存期的影响，初步筛选与预后相关的基因。然后使用R 语言“glmnet”包进行最小绝对收缩和选择算子（least absolute shrinkage and selection operator，Lasso），最终依据最佳λ值构建一个基于FRG 的预后模型。每例患者的风险值公式：风险值=，其中xi为每个FRG 的表达量，Coefi为相应的基因系数，将其求和后即为该患者的风险评分。根据风险评分的中位值将患者分为高风险组和低风险组。

1.5 预后模型的分析

应用R 语言“survival”包绘制高风险组和低风险组患者的生存曲线，比较两组患者的生存情况，采用受试者工作特征（receiver operating characteristic，ROC）曲线分析总体风险评分对预后的预测价值。应用R 语言“pheatmap”包绘制两组患者的风险评分、生存状态分布图和模型基因表达热图。对两组患者进行主成分分析（principal component analysis，PCA），观察两组患者的分布是否具有差异。同时分析风险评分、转移情况和性别是否为OS 患者预后的独立影响因素。

1.6 单基因富集分析

为探究预后相关基因在OS 进展过程中的潜在机制，采用单基因基因集富集分析（gene set enrichment analysis，GSEA）预测9 个预后相关基因的KEGG 下游途径。通过R 语言“clusterProiler”包对9 个预后相关基因进行单基因GSEA 分析，将基因表达矩阵中的正常样本删除，保留肿瘤样本作为GSEA 分析的输入文件。同时在MSigDB 数据库（https://www.gsea-msigdb.org/gsea/msigdb）中下载“c2.cp.kegg.v7.5.1.symbols.gmt”基因集，用以评估可能的相关途径。|NES|＞1、p-val＜0.05、q-val＜0.25 的通路被认为是显著富集的，展示富集最显著的前5 条通路。

1.7 预后相关基因的验证

1.7.1 预后相关基因的选择 通过查阅文献发现，酰基辅酶A 合成酶家族成员2（acyl-CoA synthetase family member 2，ACSF2）、脂肪酸去饱和酶2（fatty acid desaturase 2，FADS2）在OS 的相关研究中仍是空白，而参与预后模型构建的其他基因与OS 的关系已有了一定的研究。因此，本研究选择ACSF2、FADS2作为验证的目标基因。

1.7.2 仪器与设备 U2OS、MG63 细胞系购自武汉普诺赛生命科技有限公司，hFOB1.19（人SV40 转染成骨细胞）及hFOB1.19 细胞专用培养基均购自上海沪震生物科技有限公司，去基因组与逆转录一管化三代预混液（MR05101M）购自莫纳生物科技有限公司，荧光定量聚合酶链反应（quantitative polymerase chain reaction，qPCR）仪Light-Cycler®96 购自瑞士Roche 公司，梯度PCR 仪Biometra Tone 96 购自德国耶拿公司，超微量分光光度计NANODROP ONE 购自美国Thermo Fisher Scientific 公司。

1.7.3 细胞培养 以hFOB1.19 成骨细胞系作为对照组，OS 细胞系U2OS 和MG63 作为实验组。应用hFOB1.19 细胞专用培养基培养hFOB1.19 细胞，培养基包含DMEM/F12、0.3 g/L G418 及10%胎牛血清，将细胞置于34 ℃、5%CO2培养箱中孵育，隔2天更换1 次新的培养基。当瓶壁细胞密度达90%左右时，按1∶2 传代。U2OS、MG63 细胞在37 ℃、5%CO2培养箱中孵育。

1.7.4 qPCR 检测ACSF2、FADS2基因表达 当细胞培养至密度达90%左右时，应用Trizol 法提取样本总RNA，并采用MR05101M 合成互补DNA（complementary DNA，cDNA），然后进行qPCR。反应条件：95 ℃预变性10 min，95 ℃变性10 s，60 ℃退火30 s，共40 个循环。以β-肌动蛋白（β-actin）作为内参，采用2-△△Ct法计算ACSF2、FADS2mRNA 相对表达量。设置3 个复孔并计算均数和标准差，引物序列见表1。

表1 β-actin、FADS2、ACSF2 引物序列

1.8 统计学方法

采用R 语言4.1.1 版对生物信息学数据进行统计分析，采用Lasso 回归分析筛选预后相关基因并构建预后模型；采用Kaplan-Meier 法绘制生存曲线，组间比较采用Log-rank 检验；采用ROC 曲线评估预后模型对OS 患者1、3、5 年生存率的预测价值；采用单因素和多因素Cox 回归模型分析OS 患者总生存期的独立影响因素。采用GraphPad Prism 8.4.0 软件对qPCR 实验数据进行统计分析，符合正态分布的计量资料以均数±标准差（±s）表示，多组间比较采用单因素方差分析，多组间两两比较采用LSD-t检验。以P＜0.05 为差异有统计学意义。

2 结果

2.1 差异表达的FRG

通过差异分析共获得表达上调基因1404 个，下调基因1597 个。从FerrDb 数据库中共下载得到258 个FRG，将其与正常骨组织和OS 组织中的差异表达基因取交集，共得到57 个差异表达的FRG。（图1）

图1 差异表达基因的箱线图

2.2 GO 及KEGG 富集分析

GO 富集分析结果显示，在生物过程（biological process，BP）方面，FRG 主要与对营养水平的反应、缺氧反应、细胞对化学应激的反应等有关；在细胞成分（cellular component，CC）方面，FRG 主要与线粒体外膜、细胞器外膜、黑素体等有关；在分子功能（molecular function，MF）方面，FRG 主要与铁离子结合、泛素样蛋白连接酶结合、转氨酶活性等有关。KEGG 通路富集分析结果显示，FRG 主要参与多种疾病的神经退行性变、线粒体自噬、化学致癌-活性氧以及铁死亡等途径。（图2）

图2 OS中FRG基因的GO和KEGG富集分析

2.3 预后模型的构建

采用单因素Cox 回归模型分析差异表达的FRG对患者总生存期的影响，共得到12个预后相关基因，绘制森林图（图3A）。随后将上述12 个基因纳入Lasso 回归模型进一步筛选，回归模型的类型为Cox，采用10 倍交叉验证，最大迭代次数设为1000，筛选最佳lambda 值来确定最优模型，最终得出预后最相关的9 个FRG（图3B、3C）。其中，B 细胞淋巴瘤/白血病-2 相互作用蛋白3（B-cell lymphoma/leukemia-2 interacting protein 3，BNIP3）、FADS2、血管内皮生长因子A（vascular endothelial growth factor A，VEGFA）为高风险基因，ACSF2、芳香烃受体核转运因子样蛋白（aryl hydrocarbon receptor nuclear translocator like，ARNTL）、葡萄糖-6-磷酸脱氢酶（glucose-6-phosphate dehydrogenase，G6PD）、细胞因子信号转导抑制因子1（suppressor of cytokine signaling 1，SOCS1）、转化生长因子β受体1（transforming growth factor beta receptor 1，TGFBR1）、磷酸葡萄糖酸脱氢酶（phosphogluconate dehydrogenase，PGD）为低风险基因，分别计算每例患者的风险评分，按照风险评分的中位值将患者分为高风险组（42 例）和低风险组（43 例）。

图3 预后模型基因的筛选

2.4 预后模型分析

将高风险组和低风险组患者进行PCA，从图中点的分布可以看出两组患者的分布具有显著差异，以PC1 为分割线，低风险组患者主要分布在图中的左侧，高风险组患者主要分布在右侧（图4A）。生存分析显示，低风险组患者的生存情况明显优于高风险组，差异有统计学意义（P＜0.01）（图4B）。使用ROC 曲线验证风险评分对OS 患者1、3、5 年生存率的预测价值，结果显示，风险评分预测OS 患者1、3、5 年生存率的曲线下面积分别为0.815、0.845、0.838（图4C）。在风险曲线中，随着风险评分的增大，死亡患者逐渐增加，生存时间逐渐缩短（图5A、5B）。风险热图显示，随着风险评分的增大，BNIP3、FADS2、VEGFA的表达量增加，ACSF2、ARNTL、G6PD、SOCS1、TGFBR1、PGD的表达量降低（图5C）。

图4 高风险组（n=42）和低风险组（n=42）OS患者的生存分析与预后预测

图5 预后模型中OS患者生存风险状态图

2.5 预后影响因素分析

单因素和多因素分析结果显示，转移和风险评分均是OS 患者预后的独立影响因素（P＜0.01），性别不是OS患者预后的独立影响因素（P＞0.05）。（图6）

图6 OS患者预后的影响因素分析

2.6 单基因GSEA 分析

GSEA 结果显示，这些预后相关基因主要与原发性免疫缺陷、趋化因子信号通路、嗅觉转导、自噬的调控、维甲酸诱导基因-Ⅰ（retinoic acid inducible gene-Ⅰ，RIG-Ⅰ）样受体信号通路、剪接体、细胞质DNA 感应通路等密切相关。其中，PGD与VEGFA未发现符合条件的显著富集通路。（图7）

图7 OS患者预后相关基因的单基因GSEA分析

2.7 ACSF2、FADS2 mRNA 相对表达量的比较

U2OS 和MG63 细胞中FADS2mRNA 相对表达量均高于hFOB1.19细胞，差异均有统计学意义（P＜0.05）；MG63 细胞中ACSF2mRNA 相对表达量高于hFOB1.19 细胞，差异有统计学意义（P＜0.05）。hFOB1.19 和U2OS 细胞中ACSF2mRNA 相对表达量比较，差异无统计学意义（P＞0.05）。（表2）

表2 不同细胞中ACSF2、FADS2 mRNA相对表达量的比较（±s）

表2 不同细胞中ACSF2、FADS2 mRNA相对表达量的比较（±s）

注：*与hFOB1.19细胞比较，P＜0.05

细胞hFOB1.19 U2OS MG63 F值P值ACSF2 mRNA相对表达量1.000±0.000 2.349±1.285 10.416±0.972*89.79＜0.01 FADS2 mRNA相对表达量1.000±0.000 5.242±1.229*2.836±0.491*23.25＜0.01

3 讨论

虽然新辅助化疗的出现提高了OS 患者的无复发生存率和总生存率，但随后的几十年间，OS 的治疗效果并没有得到太大提升，似乎到了瓶颈期[12]。铁死亡自2012 年被提出后，10 余年间相关研究不断增多，逐渐成了医学领域的研究热点，在多种类型的肿瘤中均有相应研究。虽然有研究表明铁死亡与OS 的发生发展关系密切，但其作用机制尚不清楚。通过生物信息学探索与预后有关的FRG 有望成为一种突破OS 治疗瓶颈的新方法。

本研究筛选出57 个差异表达的FRG，通过GO与KEGG 富集分析发现，这些FRG 与缺氧反应、线粒体外膜、细胞器外膜、铁离子结合等有关，FRG主要参与多种疾病的神经退行性变、线粒体自噬、化学致癌-活性氧以及铁死亡等途径，这些均与肿瘤的发生发展密切相关。缺氧反应在体内肿瘤微环境中具有重要地位。有研究表明，缺氧条件促进了OS 细胞的侵袭和迁移[13]。

本研究筛选了9 个基因进行预后模型的构建。ACSF2 主要位于线粒体内，可以通过与辅酶A（coenzyme A，CoA）形成硫酯来催化脂肪酸代谢中的初始反应[14]。研究发现，ACSF2 在溃疡性结肠炎中表达显著下调，且与免疫相关通路显著相关，铁死亡抑制剂可逆转ACSF2的表达，证明ACSF2在介导铁死亡过程中具有重要作用[15]。FADS2是铁死亡的限速因子，具有催化脂肪酸去饱和的作用，可诱导宿主细胞中的铁死亡，并阻碍丙型肝炎病毒的复制[16]。ARNTL是哺乳动物生物钟的核心组成部分，可通过抑制egl-9家族缺氧诱导因子2（egl-9 family hypoxia inducible factor 2，EGLN2）的转录来抑制铁死亡，而阻断ARNTL降解或抑制EGLN2活化可以诱导铁死亡发生[17]。G6PD 被认为是肝细胞癌患者预后的危险因素，可通过抑制细胞色素p450 氧化还原酶（cytochrome p450 oxidoreductase，POR）的表达促进肿瘤细胞生长、侵袭和转移，同时减少铁死亡[18]。PGD可以减少细胞中的糖酵解，增加戊糖磷酸途径（pentose phosphate pathway，PPP），从而保护细胞免受ROS 的侵害[19]。铁死亡与ROS 关系密切，未来需要更多研究探索PGD 与铁死亡之间的关系。SOCS1是细胞因子信号转导抑制因子家族成员，在细胞生长和分化过程中发挥重要作用，可通过与p53 相互作用抑制肿瘤发生发展[20]。溶质载体家族7成员11（solute carrier family 7 member 11，SLC7A11）和亚精胺/精胺N1-乙酰转移酶1（spermidine/spermine N1-acetyltransferase 1，SAT1）在铁死亡过程中具有关键作用，且SOCS1 能够抑制谷胱甘肽（glutathione，GSH）的表达，证明SOCS1 能够提高细胞对铁死亡的敏感性[21]。TGFBR1又称激活素受体样激酶（activin receptor-like kinase，ALK）4/5，在细胞中主要参与氧化应激反应[22]。在肾近端肾小管上皮细胞中，ALK4/5 信号通路已被证明与铁死亡相关，阻断ALK4/5 信号通路可抑制铁死亡[23]。VEGFA 被认为是一种重要的血管生长因子，研究发现，当子宫内膜基质细胞发生铁死亡时，p38-促分裂原活化的蛋白激酶（mitogen-activated protein kinase，MAPK）/信号转导及转录激活因子6（signal transduction and activator of transcription 6，STAT6）信号转导会抑制VEGFA 和白细胞介素-8（interleukin-8，IL-8）的表达[24]。因此，铁死亡可能通过旁分泌VEGFA和IL-8对邻近病变产生血管生成作用，从而在子宫内膜异位症的进展过程中发挥关键作用。

本研究成功构建了基于FRG 的骨肉瘤风险预后模型，发现ACSF2、ARNTL、G6PD、PGD、FADS2、SOCS1、BNIP3、TGFBR1、VEGFA等9 个FRG 能够作为OS 患者的预后生物标志物，并采用qPCR 技术验证了ACSF2和FADS2在OS 细胞与正常成骨细胞中的表达情况，为OS 患者的临床治疗和预后评估提供了参考。然而本研究也存在一定的局限性，需要更多的样本及基础研究进一步验证。