利用机器学习鉴定并验证骨关节炎中潜在的双硫死亡相关基因

Identification and Validation of Disulfidptosis-Related Genes in Osteoarthritis by Machine Learning

摘要: 骨关节炎(Osteoarthritis, OA)是一种普遍影响中老年人群的慢性炎性疾病，其典型表现为软骨退行性变、关节炎症以及功能受损。面对目前缺乏有效治疗手段的现状，迫切需要更为全面的理论指导。双硫死亡，作为一种新近发现的细胞死亡方式，在多种疾病中均有涉及，其与骨关节炎的关联尚未明确。本研究基于已报道的双硫死亡相关基因，运用GEO数据库中的三组数据集(GSE55235, GSE55457, GSE82107)，采用多元机器学习技术将骨关节炎分为两种亚型。研究识别出九个潜在的双硫死亡相关关键差异表达基因(hubDEGs)，其中APOD、EDNRB、FOXC2、JUN、LRCH1和MAFF在骨关节炎患者中呈下调趋势，而NUDT1、PNMAL1和ZNF668则呈上调趋势。深入探究发现，这些hubDEGs与骨关节炎中免疫细胞浸润的变化密切相关，主要表现为嗜酸细胞和Th2细胞浸润的减少，以及iDC、巨噬细胞、NK细胞、Treg细胞、Tfh细胞和Th1细胞浸润的增加。利用GSE89408作为独立验证数据集，结合qRT-PCR实验验证了这九个基因在双硫死亡与骨关节炎关联中的关键作用。此外，本研究还构建了一个TF-mRNA-miRNA网络，鉴定出可能的干预靶点，包括CTCF、SP1、hsa-miR-29和hsa-miR-424。这项研究首次揭示了骨关节炎与双硫死亡之间的潜在联系，并识别出参与骨关节炎中双硫死亡过程的关键基因和靶点，为骨关节炎的临床诊断和治疗提供了新颖视角。

Abstract: Osteoarthritis (OA) is a chronic inflammatory disease that afflicts a large population of middle-aged people, manifesting as cartilage degeneration, joint inflammation, and functional impairment. Currently, effective therapeutic interventions are lacking, which calls for more comprehensive theoretical guidance. Disulfidptosis, a newly discovered form of cell death implicated in various diseases, may be related to OA, but this has not been explored before. Based on previously reported disulfidptosis-related genes, this study used three datasets from the GEO database (GSE55235, GSE55457, GSE82107) and applied various machine learning methods to classify OA into two subtypes. Nine potential disulfidptosis-related hub of Differentially Expressed Genes (hubDEGs) were identified, among which APOD, EDNRB, FOXC2, JUN, LRCH1, and MAFF were down-regulated, and NUDT1, PNMAL1, and ZNF668 were up-regulated. Further exploration revealed that hubDEGs were associated with altered immune cell infiltration in OA, marked by decreased infiltration of eosinophils and Th2 cells, and increased infiltration of iDC, macrophages, NK cells, Treg cells, Tfh cells, and Th1 cells. GSE89408, as an independent dataset, and qRT-PCR confirmed the critical role of these nine genes in the association between disulfidptosis and OA. Moreover, a TF-mRNA-miRNA network was constructed, identifying potential intervention targets, such as CTCF, SP1, hsa-miR-29, and hsa-miR-424. This study is the first to reveal a possible link between OA and disulfidptosis, and to uncover key genes and targets involved in disulfidptosis in OA, providing novel insights for the clinical diagnosis and treatment of OA.

1. 引言

骨关节炎(Osteoarthritis, OA)是一种以软骨退化、骨赘形成、滑膜炎症等为病理特征并导致关节结构和功能的损害，引起关节疼痛、僵硬、活动受限和畸形的退行性疾病。其受多种因素影响，如年龄、性别、遗传、体重、运动、创伤等 [1] [2] 。2019年，全球有5.28亿骨关节炎患者，其中73%超过55岁 [3] 。随着人口老龄化，骨关节炎将对更多老年人的生活质量和社会医疗服务造成更大负担。然而现有的治疗方法主要是缓解症状、减轻炎症、保护关节和改善功能，不能逆转或阻止软骨的损失 [4] 。骨关节炎的发生与炎症环境下以软骨细胞为代表的多种细胞损伤和死亡密切相关，然而骨关节炎中细胞的死亡方式及其与免疫的关系尚不明确，迫切需要完善的骨关节炎病理机制来指引潜在的治疗方向。

目前，已发现细胞焦亡 [5] [6] 、铁死亡 [7] [8] 、铜死亡 [9] 等细胞死亡与骨关节炎的相关性，他们主要依赖于ROS来影响骨关节炎。然而，最近有报道表明乳酸盐作为代谢物通过将葡萄糖代谢分流到戊糖磷酸途径来增加还原型烟酰胺腺嘌呤二核苷酸磷酸(NADPH)水平，并且作为信号传导分子通过人软骨细胞中的受体HCARl激活PI3K/Akt信号传导途径来上调NADPH氧化酶4 (NOX4)从而促进软骨细胞凋亡及软骨基质降解 [10] 。还发现胱氨酸水平与骨关节炎的严重程度呈正相关 [11] ，其通过半胱氨酸/GSH/GPX 4 (磷脂过氧化氢谷胱甘肽过氧化物酶)轴诱导软骨损伤 [8] 。S. Lambrecht, M. Pharm等人还在骨关节炎患者的关节软骨细胞中观察到了波形蛋白细胞骨架的实质性改变 [12] 。这些特征提示了骨关节炎与一种新的细胞死亡方式——双硫死亡的联系。

双硫死亡是一种在葡萄糖匮乏状态下迅速耗竭，导致无法被还原为半胱氨酸的胱氨酸等二硫化物异常积累 [13] ，从而诱发二硫化物应激导致细胞死亡的新方式，其机制与肌动蛋白细胞骨架、细胞伪足等高度相关 [14] 。然而，关于骨关节炎中的双硫死亡理论尚为空白，需要进一步的探索。

在本研究中，我们首次利用机器学习的方法筛选出骨关节炎中与双硫死亡相关的关键生物标志物并结合结果探讨了其影响骨关节炎的可能途径，也初步探索了这些生物标志物同免疫浸润的可能联系。同时利用验证集及qRT-PCR初步证明了结果的可靠性，还通过构建TF-mRNA-miRNA网络预测了潜在的调控靶点。这些结论可能揭示了骨关节炎与双硫死亡之间的联系，并将为骨关节炎的临床诊疗提供一种新的独特视角。

2. 材料与方法

2.1. 相关基因数据来源

从Gene Expression Omnibus (GEO)数据库中下载原始数据GSE55235 [15] 、GSE55457 [15] 、GSE82107 [16] 、GSE89408 [17] 。其中，GSE55235芯片数据包含10个骨关节炎样本和10个对照样本，GSE55457芯片数据包含10个骨关节炎样本和10个对照样本，GSE82107芯片数据包含10个骨关节炎样本和7个对照样本。GSE55235和GSE55457数据集在GPL96平台上测序，GSE82107则在GPL570平台上测序(表1)。对上述三个数据集使用R包“sva” [18] 和“limma” [19] 进行合并并去除批次效应，得到了包含30个骨关节炎样本和27个对照样本的整合数据集(NewGSE)。GSE89408中包含22个骨关节炎样本和28个对照样本，其将作为验证集以验证分析结果的可靠性。

Table 1. Original data sets

表1. 原始数据集

双硫死亡相关基因(Disulfidptosis-Related Genes, DSRG)由Liu等人 [14] 在近期提出，在我们的研究中，纳入了12个DSRG，包括SLC7A11、SLC3A2、RPN1、NCKAP1、WASF2、RAC1、NUBPL、NDUFA11、LRPPRC、OXSM、NDUFS1、GYS1。

2.2. 主成分分析(Principal Components Analysis, PCA) [20]

主成分分析是一种多变量数据降维和可视化方法，用于揭示数据中的模式和关系。我们分别使用了合并前后的数据集矩阵用于执行PCA。PCA是通过R包“FactoMineR” [21] 来进行的，该分析计算了数据中的主成分。生成PCA图时使用了R包“factoextra”的函数，帮助我们理解样本之间的相似性和差异性。

2.3. 基于DSRGs的骨关节炎亚型鉴定

无监督聚类是一种用于发现数据内部模式和结构的分析方法。使用R中的“ConsensusClusterPlus”包 [22] 来执行无监督聚类分析，这个包通过多次对数据进行随机重采样和聚类来确定最佳的聚类数目，参数设置为maxK=9，reps=50，pItem=0.8，pFeature=1，clusterAlg=pam，distance=euclidean，seed=123456。最佳聚类结果基于DSRG，这些基因在聚类分析中用于识别样本之间的相似性和差异性。DSRGs的选择是为了更好地捕捉骨关节炎样本之间的分型特征，从而划分骨关节炎亚型。

2.4. 双硫死亡亚型间的差异表达基因(SubDEGs)鉴定

使用整合数据集中的骨关节炎两亚型的样本，通过limma包进行差异基因分析，以识别不同亚型之间的差异表达基因。在分析过程中，我们设定了以下阈值：|logFC| > 0.5，表示基因在不同条件之间的表达变化具有生物学上的显著性。调整后的p < 0.05，用于控制多重检验的错误率，确保鉴定出的SubDEGs具有统计显著性。这些阈值的选择有助于筛选出最具有生物学意义的差异表达基因。利用“ggplot2”包 [23] 绘制火山图和“pheatmap”包作出热图用于显示差异基因表达的结果。

2.5. 双硫死亡相关差异表达基因(DSR-DEGs)鉴定

具体分析代码参数同“SubDEG”。使用“limma”包进行，以确定整合数据集中骨关节炎样本和对照样本之间的差异表达基因(DEG)，以研究DSRGs的表达水平对骨关节炎的影响。再将DEGs与SubDEGs取交集即可得到DSR-DEGs。

2.6. 基因集变异分析(Gene Set Variation Analysis, GSVA) [24]

使用R包“GSVA”执行GSVA分析，用于评估基因集在不同样本或条件之间的变异性的方法，以揭示与生物学过程、通路或功能相关的差异。此处将聚类所得骨关节炎两亚型作为依据，结合KEGG及Reactome数据库 [25] 进行分析。|logFC| > 0.1且p < 0.05的差异通路被认为是显著的，再利用“pheatmap”包作出热图。

2.7. GO (Gene Ontology) [26] 功能富集分析

使用R包“clusterProfiler” [27] 来执行GO功能富集分析，是一种用于识别基因集合中显著富集的生物学过程(BP)、细胞组分(CC)和分子功能(MF)的方法。从公共数据库或文献中获取基因注释数据，用于将基因与其对应的GO术语相关联。p < 0.05被认为是显著富集的。

2.8. KEGG (Kyoto Encyclopedia of Genes and Genomes) [28] 通路富集分析

使用R包“clusterProfiler”来执行KEGG通路富集分析。这个包提供了函数来确定基因集合中是否富集了与KEGG通路相关的生物通路。p < 0.05被认为是显著富集。

2.9. 关键基因(hubDEGs)的筛选

特征基因的筛选是本研究中的关键步骤，用于识别与目标变量相关且最重要的基因。使用了三种不同的机器学习算法(LASSO回归、随机森林、SVM支持向量机)来进行特征基因的筛选。在LASSO回归中，使用了“kernlab”包 [29] ，通过L1范数正则化来减小模型中不重要的特征的系数；使用“randomForest”包 [30] 进行随机森林计算，通过评估特征的重要性得分来确定哪些特征对预测任务最关键；利用“kernlab”包判断SVM的特征重要性或支持向量来确定哪些特征对分类任务最关键。最后再将三者的结果基因集取交集，得到hubDEGs。使用受试者工作特征曲线(ROC)和曲线下面积(AUC)估计hubDEGs的诊断效力，且p < 0.05被认为是显著的。

2.10. 细胞培养

使用THP-1细胞专用培养基(Procell, CM-0233)在T25培养瓶中培养THP-1细胞(Procell, CL-0233)至密度为200万细胞/ml、体积为5 ml，将细胞接种于12孔板中，50万细胞/孔，共6孔(对照组及骨关节炎组各3孔)。12 h后加入150 nM的PMA (Solarbio, P6741)诱导24 h，使THP-1细胞分化为巨噬细胞，此后换液，在骨关节炎组中加入1 μg/ml LPS (Sigma, L2880)刺激24 h。收取细胞样品。

2.11. 实时荧光定量聚合酶链式反应(qRT-PCR)

从THP-1细胞中提取总RNA，使用RNAeasy™ Plus动物RNA分离试剂盒带旋转柱(Beyotime, R0032)，严格按照生产商的说明书操作。取0.5毫克RNA，利用PrimeScript RT Master Mix (TAKARA, RR036)进行逆转录，同样遵循生产商的指导。采用TB Green Premix Ex Taq II (TAKARA, RR820L)进行qRT-PCR，以检测mRNA表达。

2.12. 统计分析

研究中的所有统计分析均使用R软件(版本4.3.0)进行。对于所有图：*表示p < 0.05，**表示p < 0.01，***表示p < 0.001。

3. 结果

3.1. 利用DSRGs将骨关节炎分为A、B两种亚型

对整合数据集去批次前后的样本进行PCA主成分分析(图1A-B)。

利用DSRG，对30个骨关节炎样本进行无监督聚类，结果显示当骨关节炎被分为2种亚型时具有最佳的分类效果(图1C-D)，因此通过DSRG我们在骨关节炎中得到了与双硫死亡相关的A、B两种亚型。

对A、B两种亚型的骨关节炎，我们作出了A、B两亚型DSRG的表达差异可视化图，从结果中可以看到SLC3A2、RPN1和NCKAP1在B亚型中的表达显著高于A亚型；而WASF2和GYS1则是在A亚型中的表达显著高于B亚型(图1E)。

使用KEGG及Reactome数据库，对A、B两亚型的所有基因进行GSVA，这是一种非参数、无监督的算法，它可以反映出A、B两亚型之间基因表达模式更本质的不同，以帮助我们认识两亚型之间的生物学功能差异。在KEGG通路富集中，A亚型主要集中于丁酸、淀粉和蔗糖代谢，钙信号通路，下丘脑分泌的促性腺激素释放激素(GnRH)信号通路等；B亚型则主要看到在烟酸、烟酰胺以及细胞分裂相关功能等通路上的活跃(图1F)。同KEGG相映证的是，在Reactome中，我们发现了类似的现象，即在B亚型中，烟酸、烟酰胺等代谢活动以及细胞分裂等相关活动有显著的增强；而在A亚型中，离子通道活跃了起来，与之相关联的IL-2、TNF受体等的信号介导及转导也得到了富集(图1G)。

3.2. 各差异表达基因分析

为了筛选出与骨关节炎关系最为密切的，且与双硫死亡相关的关键基因。取A、B两亚型间及骨关节炎、对照间的差异表达基因交集。

(A) 批效应移除前的PCA分析。(B) 批效应移除后的PCA分析。(C-D) 结合DSRGs，使用无监督聚类方法将骨关节炎分为两种亚型。(E) 骨关节炎患者DSRGs的总表达箱线图：蓝色代表A亚型样本，黄色代表B亚型样本，水平轴表示基因，垂直轴表示基因表达水平。(F) 结合KEGG数据库生成的A和B亚型的GSVA差异热图，红色代表高表达，蓝色代表低表达。(G) 结合Reactome数据库生成的A和B类群的GSVA差异热图，红色代表高表达，蓝色代表低表达。

Figure 1. Using DSRGs to divide OA into two subtypes A and B, moreover analyzing the differences in DSRGs and GSVA between the two

图1. 利用DSRGs将骨关节炎分为A、B两种亚型，并分析两者DSRGs和GSVA的差异

首先，通过PCA主成分分析的结果(图2A)可以看出，此前，我们通过DSRG利用“ConsensusClusterPlus”包将骨关节炎分为2种亚型的结果是较为可靠的。查看火山图结果(图2B)可得，在SubDEG中取并集可以得到471个差异基因。

为了对骨关节炎及对照组进行差异基因分析，以便同上述结果进行交叉分析。我们作出火山图(图2C)，筛选出差异基因共3866个，再作出差异表达基因热图(图2D)展示出对照组与骨关节炎组前40个显著差异表达基因。

(A) 将训练集中骨关节炎组样本分为A和B亚型后的PCA分析，A亚型为蓝色，B亚型为黄色。(B) SubDEGs火山图，横坐标为log2FoldChange，纵坐标为-log10 (调整后P值)。红色节点表示上调的DEGs，蓝色节点表示下调的DEGs，灰色节点表示没有显著差异表达的基因。(C) DEGs火山图，横坐标为log₂FoldChange，纵坐标为-log10 (调整后P值)。红色节点表示上调的DEGs，蓝色节点表示下调的DEGs，灰色节点表示没有显著差异表达的基因。(D) DEGs差异热图，红色代表高表达，蓝色代表低表达。

Figure 2. The analysis of differentially expressed genes between cluster A and B and between OA and Control

图2. 不同差异表达基因在A和B亚型间以及骨关节炎与对照组间的分析

3.3. DSR-DEGs分析

我们作出Veen图分别将DEG中上下调的基因同SubDEG作交集，获得48个上调基因及69个下调

(A) SubDEGs与上调DEGs的文氏图：蓝色代表SubDEGs，黄色代表上调DEGs。(B) SubDEGs与下调DEGs的文氏图：蓝色代表SubDEGs，黄色代表下调DEGs。(C-E) GSEA-GO分析：水平坐标是基因比例；垂直坐标显示GO术语；颜色表示-log10 (p值)；节点大小表示在GO术语中富集的基因数量。(F) GSEA-KEGG分析：水平坐标是基因比例；垂直坐标显示KEGG途径；节点大小表示途径中富集的基因数量；节点颜色表示-log10 (p值)。

Figure 3. Enrichment analysis of SubDEGs

图3. SubDEGs的富集分析

基因(图3A-B)。接下来，利用这117个基因组成的DSR-DEG完成了GO及KEGG的富集分析，以期找到骨关节炎中与双硫死亡相关的功能及通路差异。

从GO富集结果看，在生物学过程中的差异主要体现在细胞迁移趋化及cGMP介导的信号转导及其生物合成中(图3C)；在细胞组分中的差异则主要是神经肌肉接头、肌浆网、微纤维、单元后缘等与肌动

(A) 通过LASSO回归筛选差异表达的关键基因。(B) 骨关节炎患者SubDEGs特征的RF图。(C) 通过SVM筛选差异表达的关键基因。(D) hubDEGs的文氏图：蓝色代表LASSO，绿色代表RF，红色代表SVM。(E) 9个基因标志物在骨关节炎诊断中的ROC曲线。(F) 9个hubDEGs每对之间的关系，红色代表正相关，绿色代表负相关。

Figure 4. Screening and validation of hubDEGs from SubDEGs using machine learning

图4. 利用机器学习的方法，从SubDEGs中筛选并验证hubDEGs

蛋白细胞骨架相关的组分(图3D)；而在分子功能中，MAPK信号通路与酪氨酸的代谢及CXCR趋化因子受体结合的富集(图3E)。从KEGG富集分析结果上看，MAPK信号通路、Toll样受体信号通路、NOD样受体信号通路、IL-17信号通路、TNF信号通路、NF-κB信号通路、类风湿关节炎、醛固酮合成分泌等的富集(图3F)。

3.4. 利用机器学习筛选hubDEGs

为了筛选出骨关节炎中与双硫死亡相关最为关键的基因，以便为骨关节炎理论及临床干预提供更精确的指导，我们将上述117个基因通过特征选择的三种机器学习方法进行筛选，再将三种方法筛选出来的结果取交集，得到hubDEG。其结果如下，通过LASSO回归(图4A)、随机森林(图4B)、SVM支持向量机(图4C)，分别得到18个、20个(随机森林的前20结果，图中仅展示出前10)、117个，取交集，得到9个关键基因(图4D)，其中APOD、EDNRB、FOXC2、JUN、LRCH1和MAFF在骨关节炎患者中呈下调趋势，而NUDT1、PNMAL1和ZNF668则呈上调趋势。

我们还绘制了整合数据集中9个hubDEGs预测骨关节炎的ROC曲线(图4E)。结果，这9个hubDEGs在预测骨关节炎方面表现出色。

同时，我们对9个基因之间的相关性作出了可视化，以便探索hubDEGs之间的联系(图4F)。

3.5. HubDEGs的验证

这些关键差异表达基因(hubDEGs)在骨关节炎中与双硫死亡至关重要。为了验证它们与骨关节炎的相关性，我们使用了一个外部数据集(GSE89408)作为验证集。9个基因的ROC曲线分析表明其对骨关节炎预测能力较为出色，揭示了它们作为临床生物标志物诊断骨关节炎的潜力(图5A)。我们还诱导了一个巨噬细胞(来自THP-1的分化)炎症细胞模型，qRT-PCR结果显示，与对照组相比，hubDEGs的mRNA表达水平在统计上显示出显著且与先前的分析结果一致的差异(图5B)。

(A) 在骨关节炎诊断中9个基因标志物的ROC曲线。(B) 在150nM PMA诱导后，THP-1细胞在LPS刚性作用24小时后，通过qRT-PCR测量的hubDEGs的mRNA相对表达水平(每组n = 3)。

Figure 5. Validation of hubDEGs

图5. 验证hubDEGs

4. 讨论

骨关节炎是一种多发于中老年人群的以软骨退化、关节炎症和关节功能障碍为特征的炎性疾病，其严重影响患者的生活质量，且为社会带来较大负担 [3] 。骨关节炎表现多样且缺乏有效的临床治疗手段，因而迫切需要更为完善的理论机制 [4] 。近期，发现了一种与二硫化物应激导致的新的细胞死亡方式称为双硫死亡，其与诸多疾病进展密切相关且与炎性疾病存在联系 [31] ，但目前，双硫死亡与炎性疾病的关系尚不明确，在骨关节炎中更是一片空白。

为了阐明双硫死亡与骨关节炎的可能联系，及其与免疫浸润间的关系，首先我们利用近期发现的双硫死亡基因(SLC7A11, SLC3A2, RPN1, NCKAP1, WASF2, RAC1, NUBPL, NDUFA11, LRPPRC, OXSM, NDUFS1, GYS1)将骨关节炎分为A、B两种亚型，发现SLC3A2、RPN1和NCKAP1在B亚型中的表达显著高于A亚型；WASF2和GYS1在A亚型中的表达显著高于B亚型(图2E)。SLC3A2编码的跨膜蛋白在转运氨基酸的同时，通过调节胞内钙水平影响软骨细胞的生存和凋亡 [32] ；尽管目前还没有直接证据表明RPN1基因与骨关节炎有关，但由于RPN1参与了内质网应激，它可能通过未折叠蛋白反应(Unfolded Protein Response, UPR)影响到软骨细胞的功能和代谢 [33] [34] ；此外，SLC3A2和RPN1都参与了蛋白质N-糖基化的过程，这是一种重要的翻译后修饰方式，影响蛋白质的稳定、折叠和功能。NCKAP1与WASF2作为WAVE复合物的组成部分，通过与Arp2/3复合物相互作用调节伪足形成，进而影响细胞的形态、运动和功能，值得注意的是，这种伪足的形成可以促进双硫死亡 [35] ；GYS1基因编码糖原合成酶，其活性受其氨基(N)和羧基(C)末端的磷酸化水平的调控，当6-磷酸葡萄糖发生变构时，这种磷酸化抑制被解除 [36] [37] 。

以上述分型结果为前提，取A、B两亚型间及骨关节炎、对照间的差异表达基因交集，得到了117个基因(其中48个上调，69个下调)。通过富集分析，我们发现在骨关节炎中许多与双硫死亡相关的功能富集。涉及到的有，葡糖基转移酶活性、烟酸和烟酰胺代谢、肌动蛋白及细胞骨架的变化、细胞迁移黏附趋化及伪足形成等。说明双硫死亡的确参与了骨关节炎的进展过程。

此外，富集结果显示MAPK信号通路、Toll样受体信号通路、NOD样受体信号通路、IL-17信号通路、TNF信号通路、NF-κB信号通路富集。

基于富集分析结果结合相关文献，我们得知NADPH是还原胱氨酸所必需的辅因子，也是双硫死亡的关键因素 [13] [14] 。NADPH主要来源于葡萄糖–戊糖磷酸途径，而这个途径受到多种信号通路的调节。例如，MAPK信号通路可以通过激活转录因子AP-1，促进葡萄糖-6-磷酸脱氢酶(G6PD)的表达，G6PD是葡萄糖–戊糖磷酸途径的限速酶 [38] [39] 。Toll样受体信号通路和TNF信号通路也可以通过激活NF-κB信号通路，增加G6PD的表达 [40] [41] [42] [43] [44] 。IL-17信号通路可以通过激活STAT3信号通路 [45] [46] ，增加6-磷酸葡萄糖脱氢酶(6PGD)的表达 [47] ，6PGD是葡萄糖–戊糖磷酸途径的另一个关键酶。因此，这些信号通路可能通过增加NADPH的产生，抵抗双硫死亡。

此外，ROS作为生物体内常见且重要的氧化物，会导致NADPH等还原剂的减少，从而导致二硫化物积累，可能是双硫死亡的重要诱因 [13] [14] 。ROS主要来源于线粒体呼吸链和NADPH氧化酶(NOX)等氧化应激源。而这些氧化应激源受到多种信号通路的调节。例如，MAPK信号通路可以通过激活c-Jun和c-Fos等转录因子，促进NOX的表达和活性 [39] [48] 。Toll样受体信号通路和TNF信号通路也可以通过激活NF-κB信号通路，增加NOX的表达和活性 [48] [49] 。IL-17信号通路可以通过激活STAT3信号通路，增加NOX2和NOX4的表达 [50] [51] 。NOD样受体信号通路可以通过形成NLRP3炎性小体，激活caspase-1和IL-1β等细胞因子，增加线粒体ROS的产生 [52] 。因此，这些信号通路可能通过增加ROS的产生，促进双硫死亡。

这些信号通路之间可能存在一种平衡或拮抗的关系，一方面通过增加NADPH的产生来抵抗双硫死亡，另一方面通过增加ROS的产生来促进双硫死亡。当这种平衡被打破时，例如在葡萄糖匮乏或ROS过量的情况下，双硫死亡就会发生。不过，这些都仅仅是我们通过分析结果所作的推测，尚且需要更为可靠的实验验证。

为了确定骨关节炎中更为关键的与双硫死亡相关的潜在基因，我们通过三种特征选择的机器学习方法LASSO回归、SVM、RF，筛选出了9个hubDEG，分别是APOD、EDNRB、FOXC2、JUN、LRCH1、MAFF、NUDT1、PNMAL1、ZNF668。他们是潜在的通过调控双硫死亡而影响骨关节炎进展的基因。

APOD是一种载脂蛋白，主要参与脂质运输和细胞应激反应 [53] [54] ，影响细胞死亡的方式。而FOXC2是一种转录因子，能够促进细胞的迁移变化，与伪足形成相关 [55] [56] 。EDNRB是一种内皮素受体，主要参与血管生成和血管平滑肌收缩。EDNRB的表达受到JUN的调控 [57] ，而JUN能够影响GLUT1的表达，从而影响细胞对葡萄糖的摄取和利用 [58] [59] 。LRCH1编码一种钙调蛋白同源蛋白，主要参与细胞骨架的重塑和细胞极性的建立 [60] [61] 。LRCH1和NCKAP1之间可能存在一定的联系，因为二者都是与肌动蛋白相关的蛋白，可以调节肌动蛋白丝的组装和稳定。而NUDT1是一种核苷酸水解酶，能够水解NADPH和其他氧化型核苷酸 [62] [63] [64] 。至于PNMAL1和ZNF668，可能需要更多的研究来明确它们与双硫死亡的确切联系。

不光是验证集的ROC分析，我们的qRT-PCR结果也通过THP-1细胞建立的模型证实了hubDEGs表达水平与对照组存在着显著的统计学差异，这初步证明了9个基因在双硫死亡与骨关节炎联系中的关键性。

这些结论提供了一个理论框架，为未来更深入地研究这些基因及免疫特征乃至潜在的调控靶点在双硫死亡和参与骨关节炎进展中的作用提供了方向。未来需要更多的实验研究，以验证这些结论并揭示这些基因及免疫浸润特征、靶点通过双硫死亡影响骨关节炎的确切作用和机制。

基金项目

这项研究得到了国家自然科学基金(项目编号U23A20447，82171010)的支持。