基于生物信息与机器学习分析类风湿关节炎疾病特征基因与免疫细胞的关系①

宋世雷 陈跃平 陈 锋 （广西中医药大学附属瑞康医院，南宁 530011）

类风湿关节炎（rheumatoid arthritis，RA）是以关节滑膜慢性炎症为主要病理特征的一种自身免疫性疾病，其发展与自身抗体病理过程相关，发病初期可在血清中发现类风湿因子（rheumatoid factor，RF）和抗瓜氨酸蛋白抗体（anti-citrulline protein antibody，ACPA）水平升高［1］。随着各种免疫细胞与成纤维细胞和细胞因子相互作用，滑膜组织逐渐产生慢性炎症，并发生骨侵蚀和关节破坏，出现临床症状和损伤，进而发展为严重残疾，此外全身炎症可损害多个器官系统，如心脏组织、血管系统、肾脏、肺组织和神经系统，严重影响患者生活质量［2-3］。

RA 发病机制较为复杂，目前普遍认为其发病机制与遗传、表观遗传学、环境、代谢、免疫和微生物因素相互作用有关，尚未得到系统性阐述。临床研究和实验室研究均发现RA 疾病进展受多种免疫细胞影响，如T 细胞、B 细胞、单核细胞、巨噬细胞、中性粒细胞、肥大细胞、树突状细胞、Treg 细胞和纤维细胞样滑膜细胞等［4］。这些免疫细胞在疾病微环境中表现出不同作用，而作用不同取决于疾病环境与疾病特征基因表达不同，因此免疫细胞介导的滑膜炎症和软骨破坏机制可能在不同患者中有不同表现，导致RA 存在临床异质性，极大影响治疗的准确性［5］。因此需要寻找一种可通过检测不同患者疾病特征基因表达情况进行提前诊断与评估患者RA发生及发生后免疫细胞浸润情况的方法，实现早期诊断与精准治疗。

生物信息学作为一种新型研究手段，可高效系统地了解疾病分子作用方式，利用公共基因芯片数据库中的RA 患者与健康者滑膜组织进行合并分析，寻找RA 滑膜组织中的差异基因与分子通路，运用机械学习法筛选RA 疾病特征基因并探讨疾病特征基因与免疫浸润情况的相关性，为系统阐述发病机制及预防和提早诊断RA 提供更多依据，最终实现精准治疗。

1 资料与方法

1.1 数据来源 从GEO 数据库下载GSE12021、GSE55235 和 GSE77298 基因表达谱芯片。GSE12021 数据集包含12 例RA 患者和9 例正常关节患者滑膜组织，RA 滑膜组织由德国汉诺威临床化学研究所胡贝儿教授对该类患者行关节置换术或滑膜切除术后提取，正常滑膜组织则由创伤性关节骨折患者提取［6］。GSE55235数据集包括10例RA患者和10 例正常关节患者滑膜组织，RA 滑膜组织由WOETZEL 等［7］在进行关节置换术或滑膜切除术时收集，正常组滑膜组织来源于死亡患者正常关节或创伤性关节骨折患者；GSE77298数据集包括16例RA患者和7例正常关节患者滑膜组织，由Mathijs活检 所 得，基 于GLP570 平 台 完 成［8］；GSE12021、GSE55235数据集均基于GPL96平台完成。

1.2 数据处理及差异基因筛选 运用Perl 脚本对GSE12021、GSE55235 数据集进行整理，R 软件（版本4.1.2）“SVA”包对整理后的数据进行处理合并，获取基因批次信息，运用R软件“limma”包分别对两组合并后的批次信息进行常规矫正，输出矫正后的结果，筛选矫正后同时符合adj.P<0.05 和|log2FC|>2 的基因作为差异表达基因，分析得到差异基因，输出矫正后的差异基因表达和差异结果，“pheatmap”“ggplot2”软件包构建整合数据集中差异表达基因的分层聚类图与火山图。

1.3 差异表达基因GO和KEGG富集分析 为探求差异表达基因主要功能和途径，通过R 软件“Bioconductor”包对差异表达基因进行GO 和KEGG 富集分析，筛选有显著差异的基因（P<0.05）；GO 富集分析包含生物学过程（biological process，BP）、细胞成分（cellular constituent，CC）和分子功能（molecular function，MF），通过“ggplot2”包对富集结果进行可视化，分别绘制GO和KEGG圈图、气泡图和柱状图。

1.4 机器学习法获取疾病特征基因 运用机器学习法在差异基因中寻找疾病特征基因，首先运用LASSO回归法：R 软件“glmnet”包对差异基因表达进行过滤，通过交叉验证寻找误差最小的基因则为疾病特征基因。运用SVM-RFE 法筛选疾病特征基因：R 软件“e1071”“kernlab”“caret”包对差异基因表达进行过滤，通过交叉验证寻找误差最小的基因则为疾病特征基因，将两者筛选的疾病特征基因取交集，绘制VENN图，筛选最终的疾病特征基因，对筛选获得的基因进行ROC 曲线检测，GSE77298 与交集疾病特征基因进行验证，绘制ROC曲线。

1.5 芯片的免疫浸润分析 使用R 软件（4.1.2 版本）通过CIBERSORT 算法对以上获得的联合芯片的矫正基因表达矩阵进行分析，筛选P<0.05 的样品，并计算各样品中22种免疫细胞占比。利用“vioplot”软件包比较筛选的RA 样品与正常样品滑膜组织，分析两组间各免疫细胞浸润水平，P<0.05 代表两组间差异有统计学意义。

1.6 疾病特征基因与免疫细胞的关系 使用R 软件（4.1.2 版 本）“limma”“reshape2”“ggpubr”“ggExtra”包对疾病特征基因免疫细胞浸润结果与校正后的差异基因表达数据进行依次相关性分析，观察各疾病特征基因表达与免疫细胞的相关性，并对其统计结果进行可视化分析，绘制相关性棒棒糖图。

2 结果

2.1 差异基因筛选 在GEO 网站［Home-GEONCBI（nih.gov）］下载原始数据，运用Perl 脚本对各基因进行标准化处理，处理后的数据用R 软件进行再次处理，将GSE55235 和GSE12021 整理后的数据进行处理、合并、分析，共得到差异基因12 548 个，包含5 935 个上调基因和6 613 个下调基因，设定adj.P<0.05 和|log2FC|>2 为过滤条件，最终筛选出90 个差异基因，包含64 个上调基因和26 个下调基因，绘制差异基因热图和火山图（图1）。图1A 中正常组为Control 组，疾病组为Treat 组，红色表示高表达，蓝色为低表达，颜色越深代表程度越高，图1B红色为高表达，绿色为低表达。

图1 差异基因热图与火山图Fig.1 Heatmap and volcanic map of difference genes

2.2 差异基因GO与KEGG 富集分析 通过R软件（4.1.2）对筛选出的90 个差异基因进行GO 和KEGG 富集分析，得到GO 和KEGG 圈图。结果显示GO 分析共得到主要条目209 个，KEGG 富集分析共得到11 个项目，其中GO BP 190 个，CC 5 个，MF 14 个，主要涉及白细胞激活、淋巴细胞活化、B 细胞受体信号通路、吞噬细胞介导的免疫、淋巴细胞分化、体液免疫反应等（图2A）；KEGG 富集分析发现RA 相关通路有趋化因子信号通路、IL-17信号通路、Toll 样受体信号通路、PPAR 信号通路、原发性免疫缺陷等（图2B）。

图2 GO与KEGG圈图Fig.2 GO and KEGG circle graphs

2.3 机械学习法获取疾病特征基因 运用LASSO回归，通过R 软件“glmnet”包对差异基因表达进行过滤，取误差最小值的前10 个基因为疾病特征基因，绘制相关可视化图（图3A），纵坐标代表误差大小，横坐标代表基因个数，基因误差按从高到低排名 依 次 得 到IGHM、SLAMF8、CXCL10、FNDC4、AIM2、EGR1、AKR1B10、FKBP5、TLR8、LIF，发 现第10个基因LIF误差最小；运用SVM-RFE 法筛选疾病特征基因：通过R 软件“e1071”“kernlab”“caret”包对差异基因表达进行过滤，交叉验证寻找误差最小的疾病特征基因并进行可视化分析（图3B），纵坐标代表误差大小，横坐标代表基因个数，发现误差表达最小的点共90个基因，将两者各自得到的疾病特征基因取交集，绘制韦恩图（图4），得到IGHM、SLAMF8、CXCL10、FNDC4、AIM2、EGR1、AKR1B10、FKBP5、TLR8、LIF 10 个疾病特征基因，绘制10 个疾病特征基因的ROC 曲线，发现10 个疾病特征基因AUC 均>80%，运用测试基因芯片对上述10 个特征基 因 进 行 检 测，发 现IGHM、SLAMF8、CXCL10、AIM2、AKR1B10 AUC均>0.8，说明IGHM、SLAMF8、CXCL10、AIM2、AKR1B10为疾病特征基因。

图3 LASSO回归图与SVM-RFE图Fig.3 LASSO regression and SVM-RFE

图4 VENN图Fig.4 VENN diagram

2.4 芯片免疫浸润分析 CIBERSORT 算法分析矫正后的两个芯片基因表达，其中RA 组有22 个滑膜组织样本，正常组有19 个滑膜组织样本，同时筛选得到符合条件的22 个RA 与19 个正常组的22 个免疫细胞浸润存在差异，计算41个样品免疫细胞占比（图5A）。与正常滑膜相比，RA 滑膜中浆细胞、CD8 T 细胞、记忆静息T 细胞CD4、记忆激活T 细胞CD4、T细胞滤泡辅助器、激活的NK 细胞、单核细胞、巨噬细胞M1、活化的肥大细胞9 种免疫细胞与正常组织存在明显差异（图5B）。

图5 免疫细胞占比柱状图及小提琴图Fig.5 Proportion of immune cells columnar and vioplot

2.5 疾病特征基因与免疫细胞的相关性分析 使用 R 软件（4.1.2 版本）“limma”“reshape2”“ggpubr”“ggExtra”包对5个疾病特征基因的免疫细胞浸润结果与校正后的差异基因表达数据进行依次相关性分析，观察各疾病特征基因表达与免疫细胞的相关性，绘制可视化棒棒糖图（图6），从左到右依次对应SLAMF8、IGHM、CXCL10、AIM2、AKR1B10，最左侧纵坐标代表免疫细胞，最右侧为P值大小，圆圈越大代表P绝对值越小，发现CD8 T 细胞、浆细胞、滤泡辅助细胞、激活记忆T 细胞CD4 在SLAMF8、IGHM、AIM2、CXCL10基因中的表达较高，且该基因与这些免疫细胞存在显著相关性，活化的肥大细胞、记忆静息CD4 T 细胞、单核细胞、激活的NK 细胞在AKR1B10基因中表达较高，且该基因与这些免疫细胞存在显著相关性。SLAMF8还有M1巨噬细胞和B记忆细胞，IGHM 基因还有M1 巨噬细胞、M0 巨噬细胞和γδT 细胞，AIM2 还有M1 巨噬细胞与γδT 细胞，CXCL10还有M1巨噬细胞。

图6 疾病特征基因与免疫细胞相关性的棒棒糖图Fig.6 Lollipop chart of correlation between disease characteristic genes and immune cells

3 讨论

RA 是一种慢性、多系统的自身免疫病，其特征为持续性滑膜炎症，通常对称分布于关节周围，某些病例中出现全身症状。遗传、环境因素和自身免疫3个重要因素在RA发病机制中发挥作用。

通过公共基因芯片数据库网站下载RA 芯片数据，对其进行数据合并分析，设定筛选条件后共筛选出差异基因90 个，包含64 个上调基因和26 个下调基因，对上述基因进行GO 与KEGG 富集分析发现，GO 结果显示其主要富集于淋巴细胞活化、B 细胞受体信号通路、吞噬细胞介导的免疫、淋巴细胞分化、体液免疫反应，上述机体功能表达均与免疫息息相关［9］。淋巴细胞是白细胞中最小的一种免疫细胞，由淋巴器官产生，根据淋巴细胞来源、表面标志物与免疫功能主要分为T 淋巴细胞、B 淋巴细胞和NK 细胞，研究发现以上3 种免疫细胞在RA 发生机制上有各自的重要作用［10-12］。当T 淋巴细胞被激活，经历一系列转变得到致敏T淋巴细胞，发挥抗感染、抗肿瘤、抗异体细胞作用，T 淋巴细胞激活和分化是促进炎症发生发展的重要生理病理过程。研究发现RA 初始阶段涉及T 细胞和B 细胞激活，由T细胞与B细胞主导和组织促炎细胞因子参与炎症反应，这些炎症因子包括TNF-α、IL-1、IL-17等，最终刺激骨和软骨炎症及降解［13］。B 淋巴细胞主要通过产生并分泌免疫球蛋白和细胞内因子，包括抗原呈递和T细胞激活，与趋化因子相互作用，促进TNF-α形成，激活巨噬细胞，参与免疫调节［14］。另有研究认为B细胞在RA发病机制中发挥更重要作用，主要原因为B 细胞产生RF（抗IgG 自身抗体），且B 细胞介导的抗CCP 产生后，RA 诊断的特异度已提高至90%～98%［15］。KEGG 富集分析发现RA 相关通路有趋化因子信号通路、IL-17 信号通路、Toll 样受体信号通路、PPAR 信号通路、原发性免疫缺陷等。IL-17属于T 细胞子集，主要由CD4+Th17 细胞产生，也可由CD8+T 细胞、NKT 细胞、γδT 细胞、中性粒细胞和淋巴组织诱导剂样细胞产生，IL-17 家族由IL-17A～IL-17F 6 个成员组成，参与自身免疫病发生发展［16］。多个动物实验发现IL-17 信号通路与RA 显著相关，主要与炎症发生、软骨细胞破坏、骨质破坏有关［17-19］。有证据表明IL-17 信号通路是关节炎症发生的重要诱因，在RA 大鼠关节中注射IL-17 后发现中性粒细胞增殖与功能增强，可能与IL-17 信号通过诱导各种趋化因子和多种炎症介质分泌，包括IL-8、CXCL1、CXCL6、IL-1β、IL-6、TNF-α、GM-CSF、MIP-2、MCP-1 和G-CSF 等，促进中性粒细胞和单核细胞募集有关，IL-17 还发挥IL-17RA 和IL-17RC 介导的抗凋亡作用，促进滑膜组织增生，表明IL-17 通过刺激滑膜炎症和滑膜组织增生引发RA［20］。IL-17还可上调RANK 和M-CSF 因子，在RANKL 和M-CSF刺激下，破骨细胞分泌各种中间体和促炎介质，如诱导型一氧化氮合酶，促进骨丢失，可见IL-17 信号通路在RA 发生中占据重要位置［21］。Toll 样受体属于蛋白质识别受体，人类中Toll 样受体家族由10 个成员组成，即TLR1～TLR10，Toll 样受体通过感知保守的病原体相关微生物模式和细胞损伤后产生的危险相关分子模式变化引发固有免疫应答［22］。研究发现RA 中TLR2、TLR3、TLR4、TLR7 免疫活性在滑膜组织的内膜和下层表达上调，而在骨关节炎患者和健康供体滑膜组织中低表达，早期RA 滑膜中主要表现为TLR3 和TLR4 表达上调，而后期RA 滑膜则表现为TLR2、TLR3 和TLR4 上调，Toll 样受体特异性高表达说明其在RA 发病中的重要作用［22］。Toll 样受体信号通路对RA 的作用机制主要通过与外源性或内源性配体结合被激活，其中髓系分化原代反应蛋白88（MyD88）和诱导含TLR 结构域的适配器蛋白IFN-β（TRIF）是TLR 信号转导的两个主要配体，TLRs与MyD88结合后刺激IL-1受体相关激酶（IRAK）磷酸化，导致TRAF（TNF 受体相关因子）募集，从而激活NF-κB 信号通路和MAPK 信号通路，MAPK 信号通路下游包括细胞外信号调节激酶（ERK1/2）、C-Junn 端激酶（JNK）和P38，最终诱导促炎细胞因子TNF-α、IL-1 和IL-12 表达，提高炎症细胞趋化因子与IFN-α/β 诱导的基因产物表达，TLRs与TRIF结合则诱导IRF3/7表达从而引起IFN-β表达升高，说明Toll样受体信号通路是RA炎症发生的机制之一［23-24］。RA 病理发展特征为滑膜肿瘤样扩张，随后邻近关节软骨和骨被破坏，主要原因为活化的纤维细胞样滑膜细胞过度增殖、迁移和激活，侵袭邻近软骨和营养骨血管，因此抑制滑膜细胞增殖和迁移是治疗RA 的理想靶点［25-27］。研究表明PPAR-γ可能参与滑膜细胞表达调控，RA 成纤维样滑膜细胞中，PPAR-γ 表达下调，且PPAR-γ 抑制剂增加滑膜细胞增殖和迁移，而PPAR-γ激动剂表达抑制滑膜细胞增殖和迁移，为后续通过PPAR-γ 信号通路对RA 进行治疗提供了新思路［28-29］。GO 与KEGG 富集分析发现RA 发生与淋巴细胞、IL-17 信号通路、Toll样受体信号通路激活促进炎症因子释放及PPAR 信号通路调控纤维细胞样滑膜细胞密切相关。

机器学习是一门多领域交叉学科，涉及概率论、统计学、逼近论、凸分析、算法复杂度理论等多学科，其常见算法包括决策树算法、朴素贝叶斯算法、支持向量机算法、随机森林算法和人工神经网络算法等［30］。本研究运用生物信息学中强大的特征选择算法（SVM-RFE 和LASSO 回归法）分别处理芯片数据，利用SVM 在数据集训练得到的权重向量对特征进行排序，剔除无用特征，用剩余特征再次训练模型进行下一次迭代，选出需要的特征数，保留交叉验证后误差最小的疾病特征基因，运用LASSO回归法将处理的差异表达基因进行特征筛选，剔除无相关性或相关性不显著的疾病特征基因，得到相关性显著的疾病特征基因。将两种算法得出的结果取交集，并对所得结果进行检验，绘制ROC 曲线，用检测基因样本对得到的结果进行检验，最终得到AUC>80%的5 个基因为疾病特征基因，即SLAMF8、IGHM、CXCL10、AIM2、AKR1B10。

免疫细胞浸润是炎症反应的前提，免疫细胞与RA 基因表达可能存在密切联系，研究发现调节免疫细胞浸润情况可有效控制RA 发生。因此本研究通过CIBERSORT 算法分析免疫细胞在RA 滑膜中的表达，并通过将筛选的疾病特征基因表达与免疫浸润表达进行相关性分析，结果发现RA 滑膜中浆细胞、CD8 T 细胞、记忆静息T 细胞CD4、记忆激活T细胞CD4、T 细胞滤泡辅助器、激活的NK 细胞、单核细胞、M1 巨噬细胞、活化的肥大细胞9 种免疫细胞与正常滑膜存在明显差异，这种差异与疾病特征基因显著相关，发现CD8 T 细胞、浆细胞、滤泡辅助细胞、激活记忆T 细胞CD4 在SLAMF8、IGHM、AIM2、CXCL10 基因中表达较高，且该基因与这些免疫细胞显著相关，活化的肥大细胞、记忆静息CD4 T 细胞、单核细胞、激活的NK细胞在AKR1B10基因中表达较高，且该基因与这些免疫细胞显著相关。此外SLAMF8 还有M1 巨噬细胞和B 记忆细胞，IGHM 基因还有M1 巨噬细胞、M0 巨噬细胞和γδT 细胞，AIM2 还有M1 巨噬细胞与γδT 细胞，CXCL10 还有M1 巨噬细胞，不仅证实了5 个疾病特征基因对RA免疫应答的重要性，也可以看出不同基因影响不同免疫细胞在体内的表达，从而影响RA发生。

综上，本研究通过生信分析与机器学习法得到了RA 的主要疾病特征基因及免疫浸润情况，并阐明了不同特征基因与免疫细胞表达的相关性，说明可通过判断患者特征基因表达预测RA 免疫细胞表达，对不同患者进行针对性治疗，提高疾病治愈率，提供临床诊疗思路。