糖尿病骨质疏松症在分子水平上表现出异质性。铁死亡是一种由脂质过氧化积累引起的可控细胞死亡形式，是多种疾病发病和发展的原因之一。本研究旨在从分子水平探讨糖尿病骨质疏松症中与铁死亡相关的分子亚型，并进一步阐明其潜在的分子机制。
图1 研究设计和分析流程图
1. 鉴定在糖尿病骨质疏松症中特异表达的 FRGs
使用 "limma "软件包对 GSE35958 中 4 个对照样本和 5 个骨质疏松症样本的表达谱数据进行归一化处理（图 2A、B）。根据P-Value<0.05和log2倍变化（FC）≥1，通过差异分析共鉴定出1102个DEGs，包括677个上调基因和425个下调基因（图2C）。作者从 CTD 数据库和 GeneCards 数据库中分别获得了 38253 和 14818 个糖尿病相关基因。259 个铁死亡相关基因（FRGs）来自 FerrDb 数据库。作者将 DEGs、糖尿病相关基因和铁质疏松相关基因重叠，得到了 15 个重叠基因，即与糖尿病骨质疏松症相关的 FRGs，如维恩图所示（图 2D）。通过绘制热图（图 2E），作者可以观察到与糖尿病骨质疏松症相关的 15 个 FRGs 在骨质疏松症样本和对照样本之间有显著的差异表达。为了明确这 15 个 FRGs 之间的关系，作者采用了 Spearman 相关性分析（图 2F）。此外，15 个 FRGs 在染色体上的定位也显示在循环图中（图 2G）。
图2 糖尿病骨质疏松症中的差异 FRGs 筛查
2. 确定骨质疏松症中的铁死亡亚群
作者利用无监督聚类分析了与糖尿病骨质疏松症相关的 15 种 FRGs 在 42 个骨质疏松症样本中的表达情况，从而探索了骨质疏松症中的铁质疏松亚群。当共识矩阵的 k = 2 时，亚型的数量最为稳定，代表两个定义明确的亚型群（图 3A）。如图 3B 所示，k = 2 时的 CDF 曲线在一致性指数 0-1.0 范围内波动最小。CDF 图显示了 k 值变化时面积的相对变化（图 3C）。主成分分析（PCA）进一步证实了两个群组差异很大的结论（图 3D）。
图3 骨质疏松症中基于 FRGs 的分子聚类识别
3. 铁死亡亚群之间的差异
为了更好地理解两个铁质疏松亚群之间的区别，作者分析了两个亚群中15个FRGs的表达差异以及通路和生物活性的变化。两个亚型中 15 个 FRGs 的表达明显区别于对照组和骨质疏松症样本（图 4A）。亚型1中FBXW7、G6PD、MAPK3、PML PGD、SLC1A5、SQSTM1、TP53和YWHAE的表达水平较高，而亚型1中ALOX5、BAP1、BRD4、CDKN1A、EGFR和NNMT的表达水平较高（图4B）。GSVA 分析的生物功能结果表明，Cluster1 中含有核蛋白复合体、对创伤的反应和伤口愈合的表达下调，而 Cluster2 中神经发生、细胞对生物刺激的反应和端粒酶活性调控的表达上调（图 4C）。此外，Cluster1 的富集途径主要是上调，如结直肠癌、甲状腺癌和小细胞肺癌；而 Cluster2 则主要与下调途径相关，如基底细胞癌、亨廷顿氏病和肌萎缩侧索硬化症（图 4D）。这些结果表明，骨质疏松症患者的铁死亡亚型之间在 15 个 FRGs 的表达、富集通路和生物学作用方面存在显著差异。对于不同的铁死亡亚群，需要采取特定的治疗方法。
图4 两个铁死亡之间的差异分析
4. 铁死亡群组之间的差异基因分析
根据基因表达谱，通过 WGCNA 算法构建了与铁死亡亚群联系最紧密的可能模块。如图 5A 所示，序列号为 GSM1369716 的样本被排除在外。保持无标度拓扑网络的理想软阈值被确定为 6（R2 = 0.85）（图 5B）。根据相关性聚类，对 15 个特征模块进行了分类，并赋予不同的颜色标签（图 5C）。在这些模块中，蓝色模块（4 591 个基因）与 Cluster1（R = -0.89）和 Cluster2（R = -0.89）的相关性最强（图 5D）。作者观察到蓝色模块与模块相关基因之间存在明显的相关性（cor = 0.91）（图 5E）。随后，作者使用 adj. P-Value < 0.05 和 |log2fold change (FC)| ≥ 1 作为临界值，鉴定了铁死亡亚群的 DEGs。共发现 1,376 个 DEGs，其中 265 个出现上调，1,111 个出现下调（图 5F-H ）。
图5 铁死亡亚群间 DEGs 的鉴定
5. 与糖尿病骨质疏松症相关的铁死亡亚群的 FRGs 综合分析
结合数据库和数据集中的基因，作者共获得了 17 个与糖尿病骨质疏松症铁死亡亚群相关的 FRGs（图 6A）。PCA 结果显示了 17 个与糖尿病骨质疏松症相关的 FRGs，这些 FRGs 有效地区分了两个骨质疏松症亚群（图 6B）。此外，与糖尿病骨质疏松症相关的 17 个 FRGs 的关系网络图显示，它们之间存在显著的正相关，有助于全面分析基因之间的相互关系（图 6C）。同时，如图 6D 所示，除 BNIP3 外，所有基因都在集群 1 中高表达。为了进一步研究与糖尿病骨质疏松症相关的 17 个 FRGs 的可能生物学功能和通路活性，作者进行了 GO 和 KEGG 富集分析。GO 富集分析的重要结果显示，17 个 FRGs 主要与细胞对外界刺激的反应、对氧化应激的反应和神经元凋亡过程有关（图 6E）。此外，根据 KEGG 富集分析，17 个 FRGs 主要参与各种经典信号通路，包括脂质和动脉粥样硬化、有丝分裂动物和内分泌抵抗（图 6F）。
图6 铁死亡亚群之间 FRGs 的综合分析
6. 构建预测模型和确定关键基因
基于整个数据集，作者使用四种成熟的机器学习方法（LASSO、SVM RFE、Boruta 和 XGBoost）从 17 个 FRGs 中找到了与糖尿病骨质疏松症相关的重要基因。这些算法分别得到了 4、7、16 和 6 个基因（图 7A-E）。然后，作者利用 GSE56815 作为外部数据集，通过 ROC 曲线验证了四种机器学习算法的效率。四种算法的曲线下面积（AUC）值均大于 0.8，作者认为预测模型的结果是可靠的（图 7F）。IDH1 是所有四种算法的共同基因（图 7G）。考虑到所鉴定基因的准确性，作者用外部数据集 GSE56815 绘制了 ROC 曲线，结果显示预测效率很高（图 7H）。最终，通过四个高效预测模型，作者从 17 个 FRGs 中鉴定出 IDH1 作为糖尿病骨质疏松症亚型的前瞻性指标。
图7 构建预测模型并确定关键基因
总结
总之，作者在糖尿病骨质疏松症中发现了两个铁死亡亚群，并确认了各自的显著特征。基于 17 个铁死亡反应基因的四个不同的机器学习预测模型（LASSO、XGBoost、Boruta 和 SVM）发现了能够区分糖尿病骨质疏松症亚型的铁死亡反应调节因子。最终，作为一种经外部数据集验证的铁死亡调节因子，IDH1有能力精确区分糖尿病骨质疏松症的分子亚型，这可能会为糖尿病骨质疏松症临床症状和预后异质性的病理生理学提供新的见解。
