类风湿关节炎（rheumatoid arthritis，RA）是一种病因尚不明确的自身免疫性疾病，多表现为关节滑膜慢性炎症、关节畸形和功能丧失等^[1]。骨关节炎（osteoarthritis，OA）是由多种因素引起的关节软骨损伤及缺失所导致的关节疾病，多表现为关节疼痛、僵硬、肿大、无力及活动障碍等^[2]。近年来，微阵列技术与集成生物信息学分析被用于识别与各种疾病相关的新基因，这些基因可能作为疾病诊断和预后的生物标志物。本研究拟应用CIBERSORT算法探讨RA与OA二者存在差异的免疫细胞和诊断生物标志物。

1 材料与方法 1.1 数据收集

以“rheumatoid arthritis”和“osteoarthritis”为关键词在基因表达综合（gene expression omnibus，GEO）数据库（https://www.ncbi.nlm.nih.gov/geo/）中检索包含RA与OA的芯片，选取编号为GSE55235和GSE55457的2组芯片，共包含23例RA和20例OA数据。再根据GEO数据库中获得的探针注释文件，将每个数据集中的探针更改为基因符号。将这2个数据集合并至1个元数据队列中，并进行集成分析。当同一基因符号对应多个探针时，将探针的平均值作为该基因的表达值。选取包含10例RA和6例OA样本的GSE55584数据集用作验证队列。

1.2 差异表达基因（differentially expressed gene，DEG）分析

使用R软件中的“limma”包对元数据队列进行差异分析，采用“fdr”处理方法，并设置过滤条件为logFC filter > 1，校正P < 0.05，进行筛选得到DEG。使用R软件中“clusterProfiler”包，以P（P value filter） < 0.05为过滤条件，对DEG进行基因本体（gene ontology，GO）、京都基因和基因组数据库（kyoto encyclopedia of genes and genomes，KEGG）和疾病本体（disease ontology，DO）富集分析。

1.3 特征生物标志物的筛选

采用最小绝对收缩和选择算子（least absolute shrinkage and selection operator，LASSO）及支持向量机（support Vector machine，SVM）2种机器学习算法来预测特征生物标志物。首先，利用R软件中的“glmnet”包，对DEG用LASSO回归算法进行交叉验证并筛选特征基因。递归特征消除（recursive feature elimination，RFE）算法可避免基因的过度拟合。因此，通过R软件中的“e1071”“kernlab”和“caret”包，采用SVM递归特征消除（SVM-RFE）方法识别DEG中具有最高分辨能力的基因集。

1.4 特征生物标志物对RA的诊断与治疗价值预测

为测试1.3中所获特征基因价值，将验证组中包含的10例RA和6例OA样本数据导入R软件中，利用“limma”和“ggpubr”2个包绘制箱线图，以判断LASSO与SVM-RFE算法所得交集基因的价值，再通过R软件中的“pROC”包使用相同的数据绘制受试者操作特征（receiver operating characteristic，ROC）曲线，并用ROC曲线下面积（area under the ROC curve，AUC）值确定区分RA和OA的诊断有效性。

1.5 免疫细胞含量分析

采用CIBERSORT算法对浸润免疫细胞进行相关分析和可视化。用R软件中“corrplot”和“vioplot”包，进一步分析DEG，并行可视化处理，得到RA与OA免疫细胞浸润的差异数据及代表差异大小的小提琴图。

1.6 已鉴定基因与浸润性免疫细胞的相关性分析

利用R软件中的Spearman等级相关分析，对已鉴定的基因生物标志物与免疫细胞渗入水平的相关性进行研究，再用图表函数和R软件中的“ggplot2”包对二者的相关性进行可视化处理。

2 结果 2.1 获取DEG

通过对2个数据集（GSE55235、GSE55457）中基因表达数据的分析，获取DEG共410个，其中上调基因180个，下调基因230个，包含显著下调基因25个，显著上调基因35个。

2.2 富集分析

对DEG进行GO、KEGG和DO富集分析，并绘制柱状图。GO结果表明，DEG主要集中在白细胞迁移、趋化因子介导的信号通路等生物学进程和质膜外侧、含有胶原蛋白的细胞外基质等细胞组分及趋化因子活性、趋化因子受体结合等分子功能。KEGG富集分析显示，DEG多富集在细胞因子-细胞因子受体相互作用、病毒蛋白与细胞因子和细胞因子受体的相互作用和趋化因子信号通路等。DO富集结果显示，DEG多富集在原发性细菌性传染病、肺病等疾病。

2.3 诊断特征生物标志物的识别和验证

用LASSO算法缩小DEG的范围，最终确定了14个变量作为鉴别RA与OA发病的诊断生物标志物。用SVM -RFE算法确定了DEG中40个特征子集。最终选择了这2种算法之间的9个重叠特征基因（COPZ2、FAH、IL15RA、LTC4S、SCRG1、SFRP1、SLAMF8、SLC18A2、UNC5C）。利用GSE55584数据集验证这9个特征基因的表达水平，样品中COPZ2、FAH、IL15RA、LTC4S、SCRG1、SFRP1和SLAMF8的表达水平在OA与RA中差异显著（均P < 0.01），见图 1。因此，这7个显著差异基因被用来通过logistic回归算法建立诊断模型。

2.4 特征基因的诊断有效性

如图 2所示，7种生物标志物对区分RA与OA有显著的统计学差异。COPZ2、FAH、IL15RA、LTC4S、SCRG1、SFRP1、SLAMF8的AUC值分别为0.985（95%CI：0.948~1.000）、0.987（95%CI：0.954~1.000）、0.983（95%CI：0.935~1.000）、0.978（95%CI：0.935~1.000）、0.972（95%CI：0.922~1.000）、0.965（95%CI：0.902~1.000）、0.989（95%CI：0.963~1.000）。

2.5 免疫细胞浸润

利用R软件中“preprocessCore”“ e1071”包提取RA和OA样品中的免疫细胞，再利用R软件中的“vioplot”包进行免疫细胞差异分析，结果显示，幼稚B细胞、浆细胞、CD8⁺T细胞、记忆激活的CD4⁺T细胞、卵泡辅助细胞、T细胞调节（T调节）、T细胞γδ在RA中的比例显著高于OA，巨噬细胞M1、巨噬细胞M2、静息树突状细胞、静息肥大细胞在RA中的比例显著低于OA，见图 3。再利用R软件中的“corrplot”包对免疫细胞进行相关性分析，如图 4所示，红色越深代表二者正相关性越高，蓝色越深代表二者负相关性越高，负相关性较高的组合有浆细胞与静息记忆CD4⁺T细胞、静息NK细胞与巨噬细胞M1等，正相关性较高的组合有幼稚B细胞和T滤泡辅助细胞、静息树突状细胞和静息肥大细胞等。

2.6 7种生物标志物与浸润性免疫细胞的相关性分析

结果显示，COPZ2与静息肥大细胞、T细胞调节、静息树突状细胞呈正相关，与幼稚CD4⁺T细胞、T细胞γδ、记忆激活的CD4⁺T细胞、幼稚B细胞、CD8⁺T细胞、浆细胞、滤泡辅助性T细胞、M1巨噬细胞呈负相关。FAH与静息肥大细胞、T细胞调节、静息树突状细胞、M2巨噬细胞呈正相关，与T细胞γδ、幼稚B细胞、CD8⁺T细胞、滤泡辅助性T细胞、浆细胞、巨噬细胞M1呈负相关。IL15RA与浆细胞、巨噬细胞M1、CD8⁺T细胞、幼稚B细胞、滤泡辅助性T细胞、记忆激活的CD4⁺T细胞、T细胞γδ呈正相关，与巨噬细胞M2、静息树突状细胞、T细胞调节、静息肥大细胞呈负相关。LTC4S与静息肥大细胞、T细胞调节、静息树突状细胞呈正相关，与中性粒细胞、CD8⁺T细胞、幼稚B细胞、滤泡辅助性T细胞、浆细胞、巨噬细胞M1呈负相关。SCRG1与静息肥大细胞、T细胞调节、巨噬细胞M2、静息树突状细胞呈正相关，与幼稚B细胞、CD8⁺T细胞、滤泡辅助性T细胞、浆细胞、巨噬细胞M1呈负相关。SFRP1与静息肥大细胞、T细胞调节、静息树突状细胞、静息NK细胞呈正相关，与幼稚B细胞、T细胞γδ、CD8⁺T细胞、记忆激活的CD4⁺T细胞、浆细胞、巨噬细胞M1、滤泡辅助性T细胞呈负相关。SLAMF8与巨噬细胞M1、浆细胞、记忆激活的CD8⁺T细胞、滤泡辅助性T细胞、CD4⁺T细胞呈正相关，与静息树突状细胞、T细胞调节、静息肥大细胞呈负相关。

3 讨论

RA是一种自身免疫性疾病，主要病理改变为骨侵蚀和软骨破坏，严重影响关节功能^[3]。OA是一种常见的关节疾病，其主要病理改变为软骨退行性变、骨赘形成以及半月板损伤等^[4]。二者病因均尚未明确，也尚无有效的治疗方法。本研究欲通过探究RA与OA相关分子机制及分析免疫细胞浸润，加深对二者的认识，为相关疾病的研究提供新的思路。

本研究首次通过挖掘多个GEO数据集回溯性研究确定RA与OA患者免疫细胞渗透相关诊断生物标志物。从GEO数据集中选取并综合分析了2个数据集，共鉴定出DEG 410个，通过GO、KEGG和DO富集分析发现多种细胞组分、通路以及疾病均有DEG的参与。研究^[5]发现，RA患者的关节滑液中可检测到白细胞迁移抑制因子。炎症由趋化因子网络调控，研究^[6]表明，可通过提高抗趋化因子活性使趋化因子网络失效，从而减轻炎症反应。

本研究基于LASSO与SVM 2种机器学习算法，识别出7个特征基因。其中，COPZ2虽然并无肿瘤抑制性，但与其所包含的微RNA152一起作为肿瘤抑制因子存在于RA患者体内^[7]。CHENG等^[8]的研究表明，FAH突变的儿童患有1型肝肾酪氨酸血症，通常导致肾功能障碍、肝功能衰竭、神经损伤和癌症，因此，FAH的蛋白质替代疗法在治愈1型肝肾酪氨酸血症方面具有极大的潜力。有学者^[9]通过敲除小鼠IL15RA发现，缺少IL15RA可导致骨矿化及成骨细胞活性降低。LTC4S基因中的-444A/C多态性是白种人和阿司匹林耐受人群哮喘的危险因素^[10]。SCRG1可通过丝裂原活化蛋白激酶（mitogen-activated protein kinase，MAPK）信号通路阻断对免疫细胞具有趋化性的CCL22的产生^[11]。SFRP1作为分泌型糖蛋白SFRP家族的成员，广泛存在于包括人成纤维滑膜细胞在内的细胞中，能通过参与细胞增殖、迁移和细胞焦磷酸凋亡等生物学行为影响疾病进展和预后^[12]。有学者^[13]对TNF-α刺激的人RA滑膜细胞系和人成纤维细胞样滑膜细胞的迁移和侵袭进行检测后发现，SLAMF8的缺失可导致促炎反应显著减弱。

本研究采用CIBERSOTR算法对RA和OA对照组的免疫细胞浸润类型进行评估，发现多种免疫细胞亚型与急性心肌梗死（acute myocardial infarction，AMI）的重要生物学过程密切相关。幼稚B细胞、浆细胞、CD8⁺T细胞、记忆激活的CD4⁺T细胞、卵泡辅助细胞、T细胞调节、T细胞γδ在RA中的比例显著高于OA，巨噬细胞M1、巨噬细胞M2、静息树突状细胞、静息肥大细胞在RA中的比例显著低于OA。此外，通过对COPZ2、FAH、IL15RA、LTC4S、SCRG1、SFRP1、SLAMF8和免疫细胞之间的相关性分析发现，COPZ2、FAH、IL15RA、LTC4S、SCRG1、SFRP1、SLAMF8均与巨噬细胞M1、CD8⁺T细胞、静息肥大细胞、浆细胞、静息树突状细胞等相关。巨噬细胞存在2种不同的亚群，即经典活化型和巨噬细胞M1，二者均有促炎性，可通过被脂多糖、Th1细胞因子等物质极化后，产生促炎性细胞因子[如白细胞介素-1β（interlukin 1β，IL-1β）、IL-6、IL-12、IL-23和TNF-α等]，从而发挥抗炎和免疫调节的作用^[14]。CD8⁺T细胞的细胞溶解功能可促进全身免疫，减少感染风险，近期有研究^[15]发现，自身反应性CD8⁺T细胞也可以通过释放增加靶细胞对细胞毒性易感性的细胞因子，或通过分泌吸引其他免疫细胞至自身免疫部位的趋化因子，促进自身免疫。肥大细胞可通过淋巴器官内的同源相互作用，微调T细胞和B细胞的活性，从而发挥调节局部和全身炎症反应的作用^[16]。抗原可激活B细胞，并使其分化成浆母细胞和浆细胞，最近的研究^[17]表明，除产生抗体外，浆母细胞和浆细胞还可产生细胞因子（如IL-17、IL-10和IL-35），且二者已被证明在自身免疫性脑脊髓炎中起调节作用。树突状细胞是免疫系统中最有效的抗原呈递细胞，在调节免疫力方面起关键作用，不仅可指导T淋巴细胞和B淋巴细胞，还能激活自然杀伤细胞并产生干扰素，从而连接先天性和适应性免疫系统^[18]。

本研究存在一定的局限性。首先，本研究是生物信息学研究，血液中的生物标志物和免疫细胞图谱是从2个数据集中获得的，其重现性有待进一步验证。未来将通过生物信息学分析7种生物标志物的功能与免疫细胞浸润，并进行前瞻性分析，扩大研究的样本量，以进一步验证本研究结论。

综上所述，本研究结果显示，COPZ2、FAH、IL15RA、LTC4S、SCRG1、SFRP1和SLAMF8这7个基因可作为鉴别RA与OA的生物标志物。巨噬细胞M1、CD8⁺T细胞、静息肥大细胞、浆细胞、静息树突状细胞等免疫细胞可能在RA与OA病程中发挥不同的作用。