基于机器学习筛选类风湿关节炎的诊断标志基因和免疫浸润分析

李玲琴, 周睿姣, 张燕妮, 贺泓霓, 袁心柱

类风湿关节炎(rheumatoid arthritis,RA)是一种慢性炎症性自身免疫疾病[1],主要影响滑膜关节,表现为滑膜微血管生成、滑膜衬里细胞增生,滑膜间质有大量免疫细胞浸润刺激局部炎症,若不及时治疗,可能会导致永久性关节损伤和残疾[2]。全世界大约每200名成年人中就有1人受RA的影响,且女性的发病率是男性的2～3倍,发病高峰期为50～59岁。据统计,RA患者的寿命较无RA者少6～11年[3]。虽然目前RA仍无法治愈,但早期诊断和及时治疗有助于减缓关节损害的进展,提高患者的生活质量[4]。RA的潜在发病机制尚未完全阐明,寻找与RA相关的发病关键基因并阐明其免疫机制有助于临床提高疾病的早期诊断能力[5]。目前,基于机器学习算法的免疫浸润和生物信息学分析可以挖掘新的诊断标志物,最近已应用于许多免疫相关疾病,包括溃疡性结肠炎、白癜风、骨关节炎和银屑病等[6]。鉴此,本研究通过生物信息学方法从RA数据集中筛选出差异表达基因(differentially expressed genes,DEGs),并进一步联合机器学习算法寻找最佳的RA诊断基因,并对RA可能的发病机制及免疫浸润机制进行探索,为RA的早期诊断以及靶向治疗研发提供参考。

1 资料与方法

1.1数据资料检索 在基因表达综合(Gene Expression Omnibus,GEO) 数据库(https://www.ncbi.nlm.nih.gov/geo/)以“rheumatoid arthritis”为检索词进行检索。纳入标准:(1)包含RA患者和健康对照数据;(2)样本类别为关节滑膜组织;(3)芯片的平台信息明确。排除标准:(1)芯片平台上的基因探针ID无法转化为基因名称;(2)下载文件中的数据不完整或无效。根据上述纳入、排除标准从GEO数据库中选择了3个基因数据集:GSE55235、GSE77298和GSE55457。将GSE55457作为独立验证数据集,将GSE55235和GSE77298作为联合芯片训练集。各芯片数据文件以及相应的平台信息见表1。

表1 数据集信息

1.2数据预处理 使用“Limma”R包(版本)中的“normalizeBetweenArrays”函数来标准化GSE55235、GSE77298和GSE55457数据集,箱线图用于规范化处理后的可视化。

1.3DEGs的筛选 使用R统计软件中的Limma包研究联合芯片训练集的DEGs。筛选DEGs的截止值:adjustedP<0.05且|logFC|≥1。

1.4DEGs富集分析 应用R软件的Cluster Profiler数据包对DEGs进行基因本体论(Gene Ontology,GO)富集分析和京都基因与基因组百科全书(Kyoto Encyclopedia of Genes and Genomes,KEGG)富集分析,以P<0.05为标准筛选。

1.5机器学习筛选诊断基因 使用支持向量机-递归特征消除(support vector machines-recursive feature elimination,SVM-RFE)、最小绝对值收敛和选择算子(least absolute shrinkage and selection operator,LASSO)和随机森林(random forest,RF)三种机器学习算法筛选诊断基因,最后取三种方法获得诊断基因的交集。

1.6诊断基因的验证以及诊断效能分析 将GSE55457作为独立的验证数据集,采用成组t检验,以P<0.05判定差异有统计学意义,建立受试者工作特征(receiver operating characteristic,ROC)曲线,计算曲线下面积(area under the curve,AUC)值,评估研究指标诊断RA的效能。

1.7诊断模型的构建 整合诊断基因在验证集数据集中的表达矩阵,使用logistic回归分析构建诊断模型,应用R软件构建可视化的列线图来诊断RA患者,采用ROC曲线评估模型性能。

1.8免疫细胞浸润分析 应用xCell算法[7]计算淋巴细胞、髓系细胞、基质细胞、干细胞以及其他免疫细胞的相对表达量,并采用秩和检验比较RA组与健康对照组(HC组)的表达差异,以P<0.05为差异有统计学意义。

2 结果

2.1数据预处理结果 应用R软件Limma包对数据集GSE55235、GSE77298和GSE55457进行标化处理,批次校正后各数据集间的数据分布趋于一致。见图1。

2.2DEGs分析结果 通过分析联合芯片训练集中的数据,总共筛选出RA组和HC组的DEGs共704个,其中上调476个,下调228个,并将结果可视化为火山图(图2ⓐ)和热图(图2ⓑ)。

图2 DEGs筛选结果图

2.3DEGs的GO富集分析和KEGG富集分析结果 GO富集分析结果显示,DEGs主要富集于白细胞介导的免疫、免疫应答的激活、白细胞迁移、淋巴细胞介导的免疫。见图3ⓐ。KEGG富集分析结果显示,DEGs主要富集于趋化因子信号通路、利什曼病、RA、金黄色葡萄球菌感染等。见图3ⓑ。

图3 DEGs的GO富集分析和KEGG富集分析结果图

2.4机器学习筛选诊断基因结果 SVM-RFE算法挑选出27个诊断基因(见图4ⓐⓑ),构建LASSO回归模型并进行交叉验证,误差最小值对应17个特征基因(见图4ⓒⓓ)。RF算法鉴定了39个特征基因(见图4ⓔⓕ)。取交集得到4个诊断基因:趋化因子CXC配体13(C-X-C motif chemokine ligand 13,CXCL13)、富含亮氨酸重复序列结构域15(leucine rich repeat containing 15,LRRC15)、多配体蛋白聚糖-1(syndecan 1,SDC-1)和核酸结合蛋白3(Y-box binding protein 3,YBX3)。见图4ⓖ。

2.5诊断基因的验证以及诊断效能分析 使用GSE55457芯片数据集对筛选出的4个诊断基因进行外部验证,结果显示,相对于HC组,YBX3基因表达在RA患者中显著下调(P<0.05),而CXCL13、LRRC15和SDC-1基因表达在RA患者中显著上调(P<0.05)。见图5。这与训练集中的基因表达趋势相同。ROC曲线分析结果显示,机器算法筛选得到的4个诊断基因在验证数据集中对RA具有较高的诊断效能(AUC>0.8)。见图6。

图5 4个诊断基因在验证数据集中表达水平比较图

图6 验证数据集中4个诊断基因ROC曲线图

2.6诊断RA的列线图模型构建结果 基于验证数据集的诊断基因表达矩阵,使用logistic回归方法构建诊断预测模型,并通过列线图对模型进行可视化。见图7ⓐ。ROC曲线分析结果显示,相对于任意一个单独的诊断基因,组合模型对RA的诊断性能更高(AUC=0.985)。见图7ⓑ。

图7 诊断RA的列线图模型构建结果图

2.7免疫浸润分析结果 与HC组相比,RA组B细胞、CD4+T细胞、树突状细胞和单核细胞水平显著上调(P<0.05)。见图8。

ⓐ淋巴细胞;ⓑ髓系细胞;ⓒ基质细胞;ⓓ干细胞和其他细胞。蓝色为HC组,红色为RA组;*P<0.05

3 讨论

在本研究中,笔者对RA的基因表达芯片数据进行了生物信息学分析,共鉴定出了704个DEGs。GO和KEGG富集分析结果表明,DEGs在与免疫反应有关的生物过程中明显富集。随后联合机器学习筛选出RA的诊断基因,最后得到4个诊断基因,即CXCL13、LRRC15、SDC-1和YBX3。在验证数据集中,ROC曲线分析结果显示这些基因对RA的诊断效能较高,具有较好的临床应用前景。

3.1RA是一种以滑膜炎症和进行性关节破坏为特征的慢性全身免疫介导疾病[8]。B细胞在RA发生中发挥了重要作用。自身反应性B细胞可产生直接参与关节损伤的IgG自身抗体[9]。CXCL13是B细胞一种关键的趋化因子,可通过调节局部免疫反应和抗体产生而参与疾病发生。研究表明,在RA患者中CXCL13水平显著升高[10]。不仅如此,CXCL13还可反映RA疾病的严重程度。Meeuwisse等[11]研究表明,RA患者CXCL13基线水平越高,其关节破坏率越高。使用阿达木单抗控制RA病情后,CXCL13水平显著下降,而在治疗中断时,CXCL13的水平上升至治疗前水平[12]。动物实验显示,使用抗CXCL13的多克隆抗体治疗可以降低小鼠关节疾病的严重程度并抑制生发中心的形成[13]。

3.2LRRC15目前被认为是一个有应用前景的抗癌靶点。而在RA方面,LRRC15在RA患者的滑膜组织中存在过表达,并与滑膜细胞的增殖、迁移、侵袭和血管生成能力增加以及促炎细胞因子的加速释放有关。动物实验显示,关节炎的发展导致平衡型滑膜成纤维细胞萎缩,出现以LRRC15表达为标志的滑膜成纤维细胞轮廓,其功能是增强炎症反应和基质分解过程[14]。敲低LRRC15表达不仅能够在体外抑制滑膜细胞的侵袭性表型,而且在体内也能显著抑制胶原诱导的关节炎小鼠的滑膜增殖,减少骨侵袭和破坏[2]。He等[5]也证实了LRRC15在RA模型大鼠滑膜组织中表达显著增加。本研究结果也显示LRRC15是参与RA微环境动态平衡失调的关键因子。

3.3SDC-1也被称为CD138,它可与许多配体相互作用,引发与细胞黏附、血管生成、炎症和组织修复相关的生物事件[15]。血清SDC-1被认为是炎症活动的潜在标志物[16]。阻断白介素(interleukin,IL)-34/SDC-1通路可减轻胶原诱导的关节炎症和骨破坏,并加剧血管生成[17]。Deyab等[18]研究表明,RA患者血清C反应蛋白水平与SDC-1水平呈显著正相关,当RA患者病情得到控制后,血清SDC-1水平显著降低。在其他风湿性疾病中也观察到相似的现象,如SDC-1水平与系统性红斑狼疮活动指数及抗dsDNA抗体水平呈正相关,活动性狼疮肾炎患者的血清SDC-1水平也高于非活动性狼疮肾炎患者及非肾炎患者[19]。

3.4YBX3是一个转录因子,参与调节上皮形态发生和稳态调节[20]。尤其在肿瘤性疾病方面的研究较为集中。YBX3通过调节胆汁酸生物合成途径导致肝细胞癌转移[21]。此外,YBX3高表达与直肠癌的深度浸润相关,抑制YBX3表达可以减少体内肿瘤生长[22]。另外,抑制YBX3表达也可以提高结直肠癌细胞对化疗药物的敏感性[23]。然而,尽管YBX3的上调增加了癌细胞侵袭和肿瘤化疗耐药,但在某些肿瘤中也显示出抗癌作用。有学者发现YBX3在低级别肾透明细胞癌中的表达水平较高,而在高级别肾透明细胞癌中表达降低[24]。本研究发现YBX3在RA中扮演重要角色,值得进一步通过基础实验开展研究。

3.5RA是一种全身炎症性自身免疫性疾病,理想的治疗策略应该是在组织损伤之前重新诱导自我耐受[25]。因此,明确疾病的免疫机制可为治疗提供重要帮助。本研究通过xCell算法探讨了RA的免疫特征,发现B细胞、CD4+T细胞、树突状细胞和单核细胞呈高表达,这可能是RA发生发展的关键免疫细胞。在RA中,B细胞主要向CD4+T辅助细胞呈递自身抗原,外周血中的B细胞可分泌多种不同的细胞因子参与骨破坏,包括肿瘤坏死因子-α(tumor necrosis factor-α,TNF-α)、IL-6和IL-1β等[26]。目前,B细胞抑制剂利妥昔单抗在治疗RA中取得良好效果。滑膜组织中的特殊成分和体内产生的内源性物质可由树突状细胞作为自身抗原呈现,激活CD4+T细胞并导致炎症。Inamo等[27]认为可使用CD4+T细胞中的转录组数据对RA缓解和非缓解情况进行分类。T滤泡辅助细胞是CD4+T细胞的一种亚型,可以帮助B细胞调节抗体产生,从而进一步参与RA的发生[28]。在RA中,单核细胞离开血液并浸润发炎的滑膜组织时,可以分化成外周血来源树突状细胞,通过促进IL-17的产生和向强效破骨细胞的转化参与骨关节破坏;反之,用于控制RA炎症的生物疗法可调节单核细胞向树突状细胞的转化过程。识别和控制外周血来源树突状细胞分化的环境介质以及潜在的分子信号通路,这可能是RA新疗法开发的突破关键[29]。

综上所述,本研究基于机器学习方法筛选出CXCL13、LRRC15、SDC-1和YBX3等对RA诊断具有应用前景的因子,阐释了相关标志基因在RA中的生物学意义及可能的免疫机制,为RA的诊断和治疗靶点开发提供了参考。但本研究由于数据资料的限制,无法评估生物标志物或免疫细胞与RA患者临床特征的关联性,还需要更多的临床研究数据来进一步探索标志基因和相关免疫细胞在RA中的作用。