骨关节炎(osteoarthritis，OA) 是一种以关节软骨损伤、骨赘形成为表现的退行性疾病，可导致关节疼痛、畸形和功能障碍。OA的致病因素繁多，增龄相关因素被认为是导致OA的重要原因。目前，临床上仍通过体格检查与影像学相结合的方式来诊断OA，但OA早期阶段症状隐匿，X线检查往往难以发现。临床上缺乏针对OA早期的有效筛查手段，使得OA无法得到有效的早期预防和治疗^[1]。因此，构建OA预测模型，对OA进行早期预测和干预尤为重要。

衰老是一种细胞损伤导致的不可逆的细胞周期永久停滞形式。衰老相关分泌表型(senescence-associated secretory phenotype，SASP) 由编码细胞因子、趋化因子和金属蛋白酶的基因共同组成，介导与衰老相关的生物学效应^[2]。软骨细胞维持着关节软骨的修复、再生、代谢平衡和结构完整性^[3]。许多研究已经证实了软骨细胞衰老参与了OA的进展^[4]，并且衰老的软骨细胞数量会随着年龄的增长而逐渐增多^[5]。当过量的软骨细胞发生衰老后，软骨细胞细胞外基质分解与合成的代谢平衡被打破，关节软骨遭到破坏，加速了OA的发生^[6]。

机器学习算法将大量的历史数据导入系统，通过特定的算法，识别历史数据中的模块特征和趋势，通过计算机在一定精度范围内对结局指标进行预测。近年来，随着计算机科学与医学领域的交叉融合，越来越多的算法被用于临床疾病的诊断、预测和预后^[7]。因此，本研究选择机器学习算法筛选OA的预测基因。然而，单独的算法可能在一定程度上增加数据的过拟合性，而多种机器学习算法联合使用能在一定程度上避免这种情况的发生。因此，本研究使用3种机器学习算法来避免数据的过拟合性，提高预测的精度。

1 材料与方法 1.1 基因表达综合(gene expression omnibus，GEO) 数据库获取OA数据集

从GEO数据库中获取OA微阵列数据集GSE48556，数据集样本来源皆为人外周血单核细胞。GSE48556中共有139个样本，其中OA患者外周血单核细胞样本(以下简称OA样本) 106个，正常人外周血单核细胞样本(以下简称正常样本) 33个。从文献[8]中收集125个SASP基因。

1.2 数据处理

在R语言4.1.3中对数据集进行背景校正、归一化和log₂转换处理。随后，筛选OA数据集中SASP相关基因及其表达值。

1.3 3种机器学习算法筛选候选预测标志物

使用最小绝对收缩和选择算子(least absolute shrinkage and selection operator，LASSO)、支持向量机递归特征消除(support vector machine-recursive feature elimination，SVM-RFE) 和随机森林(random forest，RF) 3种机器学习算法筛选OA候选预测标志物。LASSO是一种利用收缩方式的回归分析算法，能提高其生成的统计模型的预测准确性、可解释性和防止过拟合。SVM-RFE是一种广泛用于微阵列数据分析的技术，可以识别具有价值的特征，从而进行结局指标的分组。RF是一种基于递归划分构建二叉树的算法，由一组决策树组成，并考虑每个决策树的随机特征，可以有效地规避高维数据集小样本的过拟合。最后，为了提高特征的准确性，将3种机器学习算法分别筛选出的候选预测标志物取交集，得到共同基因(common genes，CGs)。

1.4 CIBERSORT免疫浸润

CIBERSORT可用于定量分析基因集中相关免疫细胞和功能。本研究采用CIBERSORT评估OA样本与正常样本的免疫浸润水平。使用“pheatmap”包将OA样本与正常样本的免疫浸润水平显示在热图中。免疫浸润分析中P < 0.05为差异有统计学意义。

1.5 OA预测模型构建

使用CGs构建OA预测模型，采用受试者操作特征(receiver operating characteristic curve，ROC) 曲线的曲线下面积(area under curve，AUC) 值评价模型的预测能力，并选取预测模型中最优基因(P < 0.001) 进行动物实验验证。

1.6 微RNA (microRNA，miRNA)-转录因子(transcription factor，TF)-mRNA调控网络预测

利用miRTarBase、Starbase和TargetScan数据库预测CGs靶向miRNA。为了提高预测的准确性，只保留3个数据库中共同预测的miRNA。使用Enrichr数据库预测CGs的TF，以P < 0.05为阈值。在获得miRNA-TF-mRNA的调控关系后，使用Cytoscape可视化miRNA-TF-mRNA的调控网络。

1.7 动物实验

将12只SD大鼠适应性喂养7 d后，随机分为正常组和OA组，每组6只。OA组采用前交叉韧带切断法构建OA模型，2%戊巴比妥钠0.02 mL/kg腹腔注射麻醉，麻醉后对大鼠右膝关节进行备皮、碘伏消毒。在膝关节内侧做长约1 cm的切口，分离皮下筋膜层并暴露关节囊和髌韧带。剪断髌韧带与肌肉连接后，将髌骨和髌韧带向外侧剥离，暴露关节间隙。剪断前交叉韧带，行前抽屉实验，阳性表示已剪断前交叉韧带。术后每只大鼠腹腔注射青霉素(80 000 U/d)，连续3 d，以预防感染。30 d后处死大鼠，腹主动脉采血，切断腹主动脉，造成急性失血死亡。切开右膝关节，肉眼观察关节软骨表面退变情况，刮取关节面软骨组织至冻凝管，液氮快速冷冻。本研究获得辽宁中医药大学伦理委员会批准(21000042023054)。

1.8 实时定量PCR (real-time quantitative PCR，RT-qPCR)

取正常组(n = 6) 和OA组(n = 6) 大鼠各50 mg膝关节软骨组织，通过TRIzol法提取膝关节软骨组织总RNA，使用SYBR法进行RT-qPCR，引物采用Primer5软件设计。引物序列：TNFRSF1A，正向5’-CCTCCTCAGTGGGTTTCT-3’，反向5’-CGCCTTTCTATGCTTGTCC-3’；GAPDH，正向5’-TGCGACTTCAACAGCAACTC-3’，反向5’-ATGTAGGCCATGAGGTCCAC-3’。以TNFRSF1A与GAPDH的相对比值作为其表达量，采用2^-ΔΔCt法计算相对比值。为确保实验准确性，每组每只大鼠作3个复孔。

2 结果 2.1 SASP相关OA基因

首先，对GSE48556数据集进行预处理，共得到19 613个OA相关基因。其次，将OA基因与SASP基因相结合，分离出125个与SASP相关的OA基因。

2.2 候选基因筛选结果

LASSO和SVM-RFE均使用10折交叉验证。LASSO筛选出31个基因作为OA候选预测标志物(图 1)。SVM-RFE在特征基因数量为30个时得到最佳SVM-RFE模型(图 2)。RF选取重要性排名前10位的基因作为候选预测标志物(图 3)。将3种机器学习算法分别筛选出的候选预测标志物取交集，共得到7个CGs。

2.3 CIBERSORT免疫浸润分析

使用CIBERSORT分析正常样本与OA样本免疫浸润水平差异(图 4)，结果显示，正常样本与OA样本间浆细胞浸润水平存在显著差异(P = 0.001 3)。

2.4 OA预测模型构建

使用CGs构建OA预测模型，生成列线图并得到AUC值。列线图中，NAP1L4 (P = 0.005 479)、TNFRSF1A(P = 0.000 875)、白细胞介素(interleukin，IL)-1β (P = 0.012 166)、IL-32 (P = 0.150 475)、CD55 (P = 0.150 475) 和肿瘤坏死因子(tumor necrosis factor，TNF) (P = 0.616 024) 高表达与OA的发生率呈正相关，CXCL8低表达(P = 0.00 169) 与OA的发生率呈正相关(图 5)。其中，TNFRSF1A为预测模型中最优基因。预测模型的AUC值为0.891，说明模型具有良好的预测能力。

2.5 miRNA-TF-mRNA调控网络构建

利用miRTarBase、Starbase和TargetScan数据库进行CGs-miRNA预测，结果显示共有29个交集miRNA，其次在Enrichr数据库中进行CGs-TF预测，结果显示共有35个TF。通过mRNA-miRNA和mRNA-TF预测后，得到71个miRNA-TF-mRNA调控关系。通过Cyto-scape建立调控网络(图 6)，包括29个miRNA，35个TF，7个mRNA。

2.6 正常组与OA组大鼠膝关节软骨组织中TNFRSF1A的表达

通过RT-qPCR检测大鼠膝关节软骨组织中TNFRSF1A的表达，正常组和OA组TNFRSF1A表达水平分别为1.00±0.14和4.08±1.21，OA组TNFRSF1A表达水平明显高于正常组(P < 0.0001)，与OA预测模型分析结果一致，表明TNFRSF1A对OA具有潜在的预测价值。

3 讨论

OA是一种高致残率的疾病，早期预测对限制其致残率至关重要。许多研究^[9]发现，在OA早期，炎症通路的激活会加速SASP的分泌，诱导软骨细胞衰老。因此，SASP相关基因能否作为潜在的OA预测标志物，是本研究主要关注的问题。本研究使用3种机器学习算法和预测模型筛选出1个最优基因TNFRSF1A，并利用CIBERSORT探究正常样本与OA样本间免疫浸润水平的差异。

随着对OA研究的深入，许多研究开始关注免疫调控在预防和治疗OA中的作用^[10]。然而，对于OA相关的免疫调控仍存在许多未解之谜。因此，本研究通过CIBERSORT进行免疫浸润分析，结果显示，在正常样本和OA样本中浆细胞浸润水平存在显著差异。许多研究表明，关节软骨细胞在受损或应激后，会释放蛋白酶和胶原酶，这些降解酶在OA软骨细胞中表达失调，其表达和活性的增加是OA发生、发展过程中软骨降解的主要因素。关节软骨在损伤后巨噬细胞被激活成M1型巨噬细胞并浸润到损伤局部，分泌大量促炎性因子，如TNF-α、IL-1β、IL-6。其中IL-6作为炎症的重要调节因子，不仅可以与其他炎性因子相互作用，加速软骨降解与细胞外基质的破坏，更可以诱导B细胞分化为浆细胞^[11]。浆细胞是一种合成与储存抗体的细胞，可以产生抗瓜氨酸化的蛋白抗体，这种抗体会直接刺激破骨细胞，促进破骨细胞的生成，进一步加重OA^[12]。

机器学习算法和列线图结果显示，TNFRSF1A表现出最优秀的预测精度。现阶段的研究^[13]表明，炎症通路激活是OA的发病基础，随着炎症浸润程度的加深，软骨细胞发生不可逆的炎症性衰老，最终发展为OA。TNFRSF1A作为诱导核因子κB (nucler factor-κB，NF-κB)通路活性的主要介质^[14]，也是细胞因子TNF-α的关键细胞表面受体。TNFRSF1A通过结合TNF-α，诱导转录因子刺激NF-κB通路的激活。该通路的持续激活会导致体内炎症水平升高，并参与软骨退化、软骨下硬化等OA的早期病理过程^[15]。本研究通过动物实验验证了TNFRSF1A在OA大鼠膝关节软骨组织中高表达。因此，TNFRSF1A极有可能成为潜在的OA预测标志物。

本研究存在一定的局限性。由于数据库资源和样本量的限制，预测模型的准确性可能存在偏差。此外，本研究预测相关的miRNA-TF-mRNA调控网络未来需要深入的实验验证。

综上所述，本研究基于机器学习算法构建了一个包含7个预测候选标志物的OA预测模型。机器学习算法、免疫浸润和miRNA-TF-mRNA调控网络，可能为研究OA提供新的方向。