类风湿关节炎差异表达基因筛选、生物信息学分析及其潜在治疗药物的预测

孙胜南，廖苑君，蓝树金，李让，赵小蕾，覃继恒，饶绍奇

1广东医科大学医学系统生物学研究所，广东东莞 523808；2广东医科大学公共卫生学院

类风湿关节炎(RA)是一种慢性炎症性关节自身免疫性疾病，病理特征为关节滑膜细胞增生，新生血管形成血管翳，炎症细胞浸润，软骨和骨结构破坏[1,2]。RA患者存在持续、反复的滑膜炎症，以及严重的关节损伤，关节炎严重者可致残，使患者劳动力下降，严重影响生活质量[3]。文献[4]报道，每年每10万人中有5～50人罹患RA，其中最常发生于妇女和老年人。文献[5]报道，RA患病率和发病率呈逐年上升趋势，RA与其他风湿病是一个日益严重的公共卫生问题。虽然抗风湿药物甲氨蝶呤等能减弱或缓解RA的病情[6]，但患者治愈率仍不高。因此，RA早期阶段的检测非常必要[7]，然而现有的RA生物标志物预测价值有限，迫切需要进一步深入研究RA发生发展的分子机制，从而实现RA的早发现、诊断、治疗，以提高治疗效果，改善患者预后。基因芯片是基于碱基互补配对原则，利用化学的合成等技术，出现的高通量测序方法。至今为止，已经建成了很多的基因序列、表达量公共数据库，但这些基因表达量数据在不同疾病模型及生物学发育过程中所代表的生物学意义却知道的少之又少。大量资料表明，疾病的发生是受特异基因表达量的变化影响，这些基因变化的表达水平与人类正常的表达水平有明显的不同[8]，利用统计学、生物学、计算机等不同领域的帮助，提取不同生理、病理条件下基因表达量所代表的意义，研究不同基因所在通路，从而发现关键性的致病基因，这对RA的早期发现和诊断有巨大的意义。本研究通过筛选RA差异表达基因，并对差异表达基因进行生物信息学分析，且预测潜在的RA治疗药物，为RA的分子发生机制、临床诊断和药物治疗提供参考。

1 资料与方法

1.1 芯片数据的筛选 截止至2019年8月，以“Rheumatoid arthritis(RA)，homo sapiens”为关键词从公共Gene Expression Omnibus(GEO)数据库(https://www.ncbi.nlm.nih.gov/GEO/)搜索目标基因芯片，经过筛选，最终选取GSE97779和GSE1919作为本次研究的目标数据基因芯片。GSE97779共有14个样本，由9例RA患者样本和5例健康样本组成；GSE1919共有15例样本，由5例RA患者样本、5例健康样本、5例骨关节患者样本(不纳入本研究)组成。

1.2 RA差异表达基因筛选 GEO2R(https://www.ncbi.nlm.nih.gov/geo/GEO2R/)是一个基于R语言的差异表达基因筛选工具，可以比较疾病和正常组之间的差异。通过输入序列号、定义组别、分配组别(疾病组和健康组)和单击“Top 250”，得到疾病组和健康组差异分析结果。在此结果中筛选差异表达基因，需同时满足以下两个条件:①P<0.05;②|log2FC|>1，即差异倍数(Fold change)FC>2或FC<-2。

1.3 RA差异表达基因蛋白质相互作用(PPI)网络构建及核心基因提取 将所得差异表达基因作为初始种子基因上传到String数据库(http://string-db.org/)，得到基因互作关系网络图，以置信分数(confidence)>0.4为筛选标准，将符合条件的互作基因关系导入到Cytoscape软件进行可视化，得到PPI网络。运用Cytoscape的cytoHubba插件，获取排名前十位的核心基因。

1.4 RA差异表达基因基因本体(GO)功能富集分析、京都基因与基因组百科全书(KEGG)通路富集分析 将筛选出的差异表达基因上传至DAVID数据库(https://david.ncifcrf.gov/)，进行GO功能富集分析和KEGG通路富集分析。GO分析包括生物学过程、细胞成分、分子功能三方面内容，主要对差异表达基因进行生物学过程、细胞成分、分子生物学功能进行分析；KEEG主要针对差异表达基因的信号通路进行富集分析。

1.5 RA潜在治疗药物筛选 药物-基因相互作用数据库(DGIdb，http：//dgidb.org/)是一个整合了来自组织和展示论文、数据库和网络资源的药物-基因相互作用和基因可药性信息的数据库。DGIdb将来自30个不同来源的内容标准化，包含关于查询药物、基因和药物-基因相互作用更详细的信息。将核心基因列表上传到DGIdb数据库基因清单窗口框，进行寻找基因-药物作用关系，得到PubMed id列表、相互作用类型和得分等结果数据，可直接将结果TSV格式下载到本地，根据得分筛选潜在治疗药物。

2 结果

2.1 RA差异表达基因筛选结果 按照筛选条件，得到GSE97779和GSE1919共有的100个差异表达基因，其中上调基因77个，下调基因23个(详见图1)。

2.2 RA差异表达基因的PPI网络构建和核心基因筛选结果 最终得到具有81个节点、383条边的PPI网络(图2A)，经过Cytoscape的cytoHubba插件筛选得到10个核心基因(图2B)，分别是C-C基序趋化因子受体5(CCR5)、选择素(SELL)、颗粒酶B(GZMB)、CD2分子(CD2)、白细胞介素2受体亚单位(IL2RB)、CD44分子(CD44)、白细胞介素受体7(IL7R)、白细胞介素2受体亚单位(IL2RG)、淋巴细胞胞质蛋白2(LCP2)、LCK原癌基因(LCK)。

注：Up_List1、down_List1表示GSE97779集合中的上调、下调差异表达基因，Up_List3、down_List4表示GSE1919集合中的上调、下调差异表达基因。

图1GSE97779和GSE1919共有的100个差异表达基因

注：A为RA差异表达基因的PPI网络，B为10个核心基因。

图2RA差异表达基因的PPI网络构建和关键基因筛选结果

2.3 RA差异表达基因GO功能富集分析、KEGG通路富集分析结果 RA差异表达基因GO功能富集分析结果以P<0.05为筛选条件，按照P值从大到小的顺序排列；在生物学过程中，这些差异表达基因主要参与了免疫反应、趋化因子介导信号通路、趋化作用、细胞趋化作用、炎症反应的正调控、炎症反应、T细胞受体信号通路、信号转导、白细胞迁移、细胞外基质、T细胞共刺激、调节细胞增殖、细胞间信号、细胞表面受体信号通路、细胞防御反应、调节免疫反应、T细胞分化、T细胞分化的调控、细胞黏附、适应性免疫反应和T细胞趋化性的正调控过程；由细胞组分可知，这些差异表达基因表达主要位于质膜外侧面、细胞外区域、原生质膜、α-βT细胞受体复合物、细胞外空间、细胞表面、膜、免疫突触和质膜部位；从分子功能上看，这些差异表达基因具有使趋化因子的活动、CXCR3趋化因子受体结合的作用。RA差异表达基因KEGG通路富集分析结果以P<0.05为筛选条件，按照P值从大到小的顺序排列，这些差异表达基因表达通路主要富集于细胞因子—细胞因子受体相互作用通路、趋化因子信号通路、类风湿关节炎、原发性免疫缺陷、自然杀伤细胞介导的细胞毒性、造血细胞系、Jak-STAT信号通路、T细胞受体信号通路。

2.4 RA潜在治疗药物预测筛选结果 根据PPI网络中筛选出的10个核心基因，分别为CCR5、SELL、GZMB、CD2、IL2RB、CD44、IL7R、IL2RG、LCP2、LCK。在DGIdb数据库中寻找相关潜在药物，最后搜索发现除LCP2外的9个基因都找到相关药物。由基因CCR5预测到的潜在治疗药物有胍基乙酸、依诺肝素、安克里维克、阿帕维洛克、塞尼克罗克、马拉维克、酚普罗考蒙、阿普洛克、马来酸维里韦洛克；由基因SELL5预测到的潜在治疗药物有双咪唑糖、里维潘塞尔、阿塞单抗；由基因GEZM5预测到的潜在治疗药物有(2S)-2-(甲胺基)-1-苯丙-1-醇；由基因CD25预测到的潜在治疗药物有阿利法西普、西普利单抗、奥曲肽、葡萄牛磺酸；由基因IL2RB5预测到的潜在治疗药物有艾地氯喹、地尼白介素、巴利昔单抗、赛尼哌、硫柳单抗；由基因CD445预测到的潜在治疗药物有是牛磺酸、干扰素GAMA-1B、甲酚酸、庆大霉素；由基因IL7R5预测到的潜在治疗药物有鲁索利替尼；由基因IL2RG5预测到的潜在治疗药物有赛尼哌、巴利昔单抗、地尼白介素；由基因LCK5预测到的潜在治疗药物有达沙替尼、波纳替尼、阿卡鲁替尼、十二醇、伊布替尼、萨拉卡蒂尼布、盐酸帕佐帕尼、伊洛拉塞提布、葡萄孢素。我们选取得分排名前三的药物马拉维克、阿利法西普、艾地氯喹作为治疗RA的候选治疗药物。

3 讨论

本研究通过对GEO数据库中的基因芯片GSE97779和GSE1919进行生物信息学分析获得差异表达基因，并分析了有关这些基因富集的生物学过程、细胞成分、分子功能和信号通路，为RA的机制研究提供了理论依据与研究方向。在我们的研究中，RA病例组和对照组之间共100个差异表达基因。功能分析显示，差异表达基因主要富集于细胞因子活性、G 蛋白偶联受体结合和细胞因子受体结合。KEGG通路分析显示，差异表达基因主要富集于细胞因子-细胞因子受体相互作用、趋化因子信号通路、病毒蛋白与细胞因子、细胞因子受体相互作用。

CD44抗原是一种参与细胞间相互作用、黏附和迁移的细胞表面糖蛋白。在人类中，CD44抗原由11号染色体上的CD44基因编码。CD44被认为粘连分子、吞噬糖蛋白-1、Hermes抗原、淋巴细胞归巢受体，一系列研究表明CD44可能是RA的标志物，与RA的发生发展相关[9～11]。这些研究表明，本研究的筛选结果与前人研究一致。CD2是在T细胞和NK细胞表面发现的一种细胞黏附分子，参与了CD8+T细胞和NK细胞介导的人类和其他几种哺乳动物细胞免疫的免疫调节[12]。CCR5是一种趋化因子受体，被多个参与炎症反应的免疫细胞广泛表达，我们推测CCR5可能通过这些炎症途径通路参与了RA发生。GZMB基因沉默通过调节炎症因子、凋亡相关因子、血管生成相关因子的表达，抑制MAPK信号通路，缓解RA引起的滑膜组织增生和关节软骨组织损伤[13]。LCP2是淋巴细胞胞质蛋白2，在PubMed数据库中未发现LCP2与RA有直接关系。但有研究证实其与肿瘤密切相关，包括结直肠癌[14,15]、多形性胶质母细胞瘤[16,17]、弥漫性大b细胞淋巴瘤[18]。癌症与RA可能存在共同的免疫通路。LCK是一种非受体src家族激酶，在细胞周期控制、细胞黏附、运动、增殖和分化等多种细胞过程中发挥重要作用。有研究表明，LCK是调节T细胞功能的关键分子之一，调节TCR信号的启动、T细胞的发育和T细胞内稳态[19]。我们推测LCK通过调节T细胞免疫通路影响RA的发生发展。到目前为止，还没有研究显示有SELL与RA的发生、发展有直接的联系。因此，对该基因的进一步研究可能会有新的发现，为RA的分子发病机制研究提供了新的方向。

功能富集分析中，100个差异表达基因主要富集于免疫、炎症、信号通路，如白细胞迁移、细胞黏附的正向调节、受体配体活性、细胞因子-细胞因子受体相互作用通路等。目前药物控制是RA患者的重要治疗手段，寻找有效、敏感的抗RA药物有助于改善RA患者的预后。10个核心基因共发现76种化学物质，其中选择得分前3的化合物为候选药物。马拉维克是一种一级选择性CCR5拮抗剂，仅被批准与其他抗逆转录病毒药物联合治疗艾滋病病毒感染的[20]，已有研究报道，马拉维克有助于免疫综合征的初始预防和免疫综合征的积极治疗[21]，减轻神经性疼痛症状的发展[22]。但其抗RA作用及作用机理至今尚未阐明。在本研究中，我们发现了马拉维克作为RA治疗剂的潜力。关于阿利法西普和艾地氯喹，已有研究报道它们与银屑病和新发型1型糖尿病有关[23,24]。我们的研究为RA的治疗提供了新的药物靶点，当然还需要更多的实验研究来证实这些发现。

总之，本文从GEO数据库获得RA相关数据集，基于GEO2R采用limma包分析筛选出100个差异表达基因、PPI网络构建及核心基因筛选，并对其进行GO和KEGG功能分析和潜在药物的预测。综上所述，我们采用大数据挖掘和计算生物学的综合生物信息学方法筛选出了10个与RA发生有关的核心基因(CCR5、SELL、GZMB、CD2、IL2RB、CD44、IL7R、IL2RG、LCP2、LCK)和3种潜在治疗药物。本研究为RA的发病机制、早期诊断和临床药物治疗提供有价值的参考。