基于综合生物信息学和机器学习算法构建衰老相关分泌表型的骨关节炎预测模型

刘孝生，魏东升，2，何信用，方策

（1.辽宁中医药大学研究生学院，沈阳 110847；2.辽宁中医药大学中医脏象理论及应用教育部重点实验室，沈阳 110847；3.抚顺市中医院骨伤一科，辽宁 抚顺 113008）

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

衰老是一种细胞损伤导致的不可逆的细胞周期永久停滞形式。衰老相关分泌表型（senescenceassociated 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中对数据集进行背景校正、归一化和log2转换处理。随后，筛选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，RTqPCR）

取正常组（n=6）和OA组（n=6）大鼠各50 mg膝关节软骨组织，通过TRIzol法提取膝关节软骨组织总RNA，使用SYBR法进行RT-qPCR，引物采用Primer5软件设计。引物序列：TNFRSF1A，正向5’-CCTCCTCAGTGGGTTTCT-3’，反向5’-CGCCTTTC TATGCTTGTCC-3’；GAPDH，正向5’-TGCGACTTCA ACAGCAACTC-3’，反向5’-ATGTAGGCCATGAGGTC CAC-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个时得到最佳SVMRFE模型（图2）。RF选取重要性排名前10位的基因作为候选预测标志物（图3）。将3种机器学习算法分别筛选出的候选预测标志物取交集，共得到7个CGs。

图1 LASSO模型中的特征选择Fig.1 Feature selection in least absolute shrinkage and selection operator（LASSO）model

图2 特征基因数量为30时得到最佳SVM-RFE模型Fig.2 The optimal support vector machines recursive feature elimination（SVM-RFE）model was obtained at 30 featured genes

图3 RF模型图和RF预测基因重要性排序图Fig.3 Random forest（RF）model diagram and importance ranking diagram

2.3 CIBERSORT免疫浸润分析

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

图4 正常样本与OA样本中免疫浸润水平差异性厢图Fig.4 Differences in immune infiltration levels between normal and osteoarthritis（OA）samples

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，说明模型具有良好的预测能力。

图5 OA预测模型列线图Fig.5 Nomogram for osteoarthritis（OA）prediction model

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调控关系。通过Cytoscape建立调控网络（图6），包括29个miRNA，35个TF，7个mRNA。

图6 miRNA-TF-mRNA调控网络图Fig.6 Diagram of miRNA-TF-mRNA regulatory network

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提供新的方向。