类风湿关节炎基因表达谱芯片的生物信息学分析①

袁文华 解琪琪 刘 正 沈 伟 冯晓飞 周海宇

(兰州大学第二医院骨科，甘肃省骨与关节疾病研究重点实验室，兰州 730000)

类风湿关节炎(rheumatoid arthritis，RA)是关节慢性炎症性疾病，发病率为0.5%～1%，可引起软骨和骨损伤，甚至导致残疾[1，2]。目前对于大部分RA患者均采用以延缓病变进程为主要目的的治疗，然而并不能够逆转或者阻止病变的进展，RA除了降低患者生活质量外，还给个人、家庭及社会带来了巨大的负担[3，4]。关于RA发生发展的机制尚不完全明确，其中遗传和滑膜炎症扮演重要角色，进一步阐明其发病机制是当前RA研究领域的难题[5-7]。近年来生物信息学作为生命科学领域的新兴学科，为疾病分子机制的实验研究提供了可能的依据和可行的思路[8，9]。本研究通过生物信息学相关方法对GEO数据库中3个RA与正常滑膜组织基因芯片数据集进行整合分析，筛选与RA相关的差异表达基因(differentially expressed genes,DEGs)，然后进行GO功能注释、KEGG通路富集分析、蛋白-蛋白相互作用(PPI)网络构建，识别参与RA发生发展的关键基因，为RA发病机制的进一步研究提供生物信息学依据。

1 材料与方法

1.1材料 生物信息学分析数据集GSE1919、GSE55235、GSE55457均来源于NCBI的GEO(Gene Expression Omnibus)数据库，均为滑膜组织总RNA的mRNA表达谱芯片分析[10-12]；GSE1919是基于GPL91平台[HG_U95A]的Affymetrix Human Genome U95A Array，该芯片数据包括5例RA滑膜组织样本、5例骨性关节炎组织样本及5例正常人滑膜组织样本，本研究纳入5例RA及5例正常滑膜组织；GSE55235是基于GPL96平台[HG-U133A]的Affymetrix Human Genome U133A Array，该芯片数据包括10例RA滑膜组织样本、10例骨性关节炎组织样本及10例正常人滑膜组织样本，本研究纳入10例RA及10例正常滑膜组织；GSE55457同样是基于GPL96平台[HG-U133A]的Affymetrix Human Genome U133A Array，该芯片数据包括10例RA滑膜组织样本、13例骨性关节炎组织样本及10例正常人滑膜组织样本，本研究纳入10例RA及10例正常滑膜组织。

1.2方法

1.2.1数据处理及筛选DEGs 3个数据集的原始CEL文件均通过R语言affy软件包分析，进行标准化预处理，得到每个样本的探针表达矩阵，利用R语言limma工具包输出同时满足|log2 fold change(log2FC)|>1且P<0.05的DEGs，然后根据相应平台文件所对应的R语言hgu133a.db、hgu95av2.db包将探针名转化为基因名；将3个数据集所得DEGs取交集，然后对交集所得基因(以|log2FC|为标准)使用heatmap.2工具包绘制热图，将每个DEG的差异表达情况进行直观展示[13，14]。

1.2.2DEGs的基因本体论和通路富集分析 基因本体论(gene ontology，GO)是注释基因及其产物的重要方法和工具，有利于生物数据的整合和利用[15]。注释、可视化数据库DAVID(http://david.abcc.ncifcrf.gov/)是1个在线数据合成工具，为成功的高通量基因功能分析奠定基础。利用DAVID对DEGs进行GO和KEGG通路富集分析，并用R软件包GOplot进行结果可视化处理[16，17]。显著性基因富集的临界值设定为P<0.05。

1.2.3共表达分析 通过将所得58个DEGs导入STRING(Version：10.5)在线软件(https://string-db.org/)中，得出PPI关系(score>0.4)，并以TSV格式导出后，再将得到的源文件导入Cytoscape软件，使用插件NetworkAnalyzer、cytoHubba进行关键基因(Hub基因)的筛选，选用MCC算法，得分前10的基因选取为Hub基因[18-20]。

2 结果

2.1数据处理和筛选DEGs 对GSE1919、GSE55235、GSE55457数据集的DEGs进行标准化处理，标准化后以箱式图呈现(图1)，然后采用spearman算法对两组样本进行聚类分析，计算样本之间的相关系数，最终以R语言的pheatmap包进行结果可视化[21]。结果表明样本来源可靠(图2)；GSE1919筛选出262个DEGs，GSE55235筛选出527个DEGs，GSE55457筛选174个DEGs(图3)；图4A～C为3个数据集表达的聚类热图，将GSE1919、GSE55235、GSE55457数据集所得差异基因取交集得到58个DEGs(图4D)。

图1 样本表达校正箱式图Fig.1 Sample expression correction box diagramNote:Box plot of number of genes among RA patient samples and control samples in GSE1919，GSE55235，and GSE55457.Red dots represent mean values of gene expression after sample normalization.

图2 样本聚类图Fig.2 Sample clustering diagramNote:Heatmap of number of RA patient samples and control samples in GSE1919，GSE55235，and GSE55457.Green represents control group and purple represents RA group.

图3 差异基因火山图Fig.3 Differential gene volcano diagramNote:No significantly changed genes are marked as black dots，red dots represent up-regulated differential genes and green dots represent down-regulated genes.

图4 差异基因交集Venn图及热图分析Fig.4 Differential gene intersection Venn diagram and heatmap analysisNote:A-C.Heatmaps using hierarchical clustering of DEGs in GSE1919，GSE55235，and GSE55457；D.Venn diagram DEGs were selected with a fold change > 1 and P-value <0.05 among GSE1919，GSE55235，and GSE55457.

2.2GO和KEGG通路富集分析 采用DAVID数据库对58个DEGs进行GO功能富集分析结果表明，DEGs主要涉及免疫应答、趋化因子介导的信号通路、炎症反应及细胞表面受体信号通路等生物过程，介导趋化因子活性、CXCR3趋化因子受体结合、抗原结合、免疫球蛋白受体结合等分子功能，主要富集于质膜外区域；KEGG通路富集分析结果显示DEGs主要涉及细胞因子-细胞因子受体相互作用、趋化因子信号通路、Toll样受体信号通路等经典信号通路(图5)。

图5 差异基因GO分析及通路富集分析Fig.5 GO analysis and pathway enrichment analysis of differential genesNote:A.GO functional enrichment analysis of DEGs;B.Kyoto Encyclopedia of Genes and Genomes enrichment analysis of DEGs.

2.3差异表达编码蛋白质间相互网络分析及Hub基因表达验证 STRING工具和Cytoscape软件对58个显著DEGs进行分析结果表明，CCR5、CCL5、CCL19、CXCL9、CXCL10、CCR2、CXCL13、CD2、CD27、PTPRC为所得Hub基因(图6)。

图6 差异基因编码蛋白质的PPI分析和Hub基因Fig.6 PPI analysis of differentially encoded proteins and Hub geneNote:A total of 58 DEGs were imported into STRING online software (version：10.5) to identify interactions among nodes (the parameters were default in STRING database) and exported in TSV format.The obtained source file was imported into Cytoscape software.A.Data were analyzed using NetworkAnalyzer plugin;B.CytoHubba plugin was used to analyze hub genes with maximum correlation criterion (MCC).

3 讨论

RA是最常见的慢性关节炎之一，可导致软骨和骨损伤甚至残疾[1，2]。当前的治疗主要以延缓病变进程为主，并不能逆转或者阻止病变的进展[3，4]。且RA的病因尚不清楚，缺乏可靠的生物标志物评估其预后、治疗反应和毒性，导致治疗策略的制定在临床实践中十分困难[22]。因此，阐明RA的发病机制将有助于制定更有效的治疗策略，并可能找到潜在的治疗靶点。本研究利用生物信息学方法对3个数据集的RA和正常滑膜组织基因芯片进行整合分析，共筛选出58个DEGs。随后对DEGs进行GO功能注释及KEGG通路富集分析，寻找Hub基因，旨在为进一步研究RA的发病机制提供参考依据。

为了解DEGs在RA中所涉及的通路及功能，课题组进行了GO功能注释及KEGG通路富集分析。结果显示，DEGs主要涉及免疫应答、趋化因子介导的信号通路、炎症反应及细胞表面受体信号通路等生物过程，介导趋化因子活性、CXCR3趋化因子受体结合、抗原结合、免疫球蛋白受体结合等分子功能，主要富集于质膜外区域；上述生物过程及分子功能与免疫系统功能密切相关，而免疫系统功能的异常导致了RA的发生发展[23，24]。KEGG通路富集分析显示，DEGs主要富集到细胞因子-细胞因子受体信号通路、趋化因子信号通路和Toll样受体信号通路；细胞因子与细胞因子受体信号通路参与了RA的发病机制，并在RA中介导中性粒细胞自噬[25，26]；Toll样受体信号通路与炎症和自身免疫性疾病的发病机制密切相关[27]。大量研究表明，该信号通路在RA发生发展中起重要作用，且趋化因子信号通路同样在RA中起一定作用[28-33]。由此可见，课题组的分析结果与既往研究结论一致，可信度较高，提示这些DEGs和信号通路可能参与RA的发生发展，通过对这些DEGs的监测可能有助于评估RA的发生与进展，而这些信号通路有望成为潜在的药物治疗靶点。

CCR2、CCR5、CCL5、CCL19、CXCL9、CXCL10、CXCL13、CD2、CD27、PTPRC为所得Hub基因。CCR2和CCR5同属于趋化因子受体家族且互为同源物，CCR2位于3号染色体的趋化因子受体基因簇区，其编码的蛋白是单核细胞趋化蛋白-1(monocyte chemoattractant protein-1，MCP-1)的受体，是一种特异性介导单核细胞趋化性的趋化因子，而MCP-1主要参与RA等炎症性疾病的单核细胞浸润[34]。CCR5功能与CCR2类似，为趋化因子CCL5的受体[35，36]。早在1999年Suzuki等[37]就发现了CCR5通过选择性的募集促炎性T细胞到RA关节处从而发挥作用，之后大量研究也验证了两者在RA的进展中起重要作用[38-40]。CCL19是一种B 细胞特异性趋化因子，通过B 细胞上表达的相应受体(CXCR5、CXCR4和 CCR7)，促进 B 细胞在组织中的聚集和滑膜内生发中心的形成[41]；Sellam等[42]研究表明，RA患者血清中B细胞相关趋化因子CCL19表达水平升高，是B细胞介导的RA亚型新的生物标志物，可用于预测利妥昔单抗(RTX)治疗RA的临床疗效;有研究报道，RA患者的血液和滑液中趋化因子CCL19的水平升高，其通过CCR7/Rho信号通路刺激破骨细胞迁移和破骨活性，促进RA的发生发展[43，44]。CXC趋化因子家族，又称为γ-干扰素诱导的单核细胞趋化因子家族，为趋化因子大家族中的一个亚族，在炎症性疾病及肿瘤中发挥重要作用[45-47]。Tsubaki等[48]研究表明在RA的早期阶段，表达CXCR3的浆细胞大量被CXCL9阳性滑膜成纤维细胞招募至滑膜处；Kotrych等[49]通过检测422例RA患者和338例正常人的CXCL9的基因多态性，发现其与RA的关系密切。以上研究均说明CXCL9在RA的发生发展中发挥重要作用，但其具体作用机制不明，值得进一步深入探讨。CXCL10是一种蛋白质编码基因，其编码的促炎细胞因子参与多种生物过程，如外周免疫细胞的趋化、分化和活化等[50]；Laragione等[51]研究表明，CXCL10在RA的高侵袭性成纤维细胞样滑膜细胞(fibroblast-like synoviocyte,FLS)中的表达水平升高，发现CXCL10/CXCR3在调节RA患者FLS侵袭中具有一种新的自分泌/旁分泌作用，并提出CXCL10-CXCR3轴可能是一个新的潜在治疗靶点，可用于减少FLS侵袭及其相关的关节损伤和血管翳在RA中的侵袭和破坏；有研究报道，CXCL10在RA早期炎症中起关键作用，可作为RA早期疾病活动的生物标志物[52]。此外，CD2、CD27在RA中的作用研究较为深入，机制相对明确[53,54]。PTPRC又称受体型酪氨酸蛋白磷酸酶C(receptor-type tyrosine-protein phosphatase C，PTPRC)，为蛋白酪氨酸磷酸酶(protein tyrosine phosphatase，PTP)家族的成员，PTP参与调节细胞生长、分化、有丝分裂等多种细胞过程，已被证实是T细胞和B细胞抗原受体信号的重要调节因子，还可抑制JAK激酶，而这些功能均与RA的发生发展密切相关[1，4，22]。最新研究表明，PTPRC在RA中发挥作用，但具体作用机制尚不明确，值得深入研究[55，56]。CXCL13是一种通过与CXC趋化因子受体5(CXCR5)结合而吸引B淋巴细胞和滤泡辅助T淋巴细胞的趋化因子，涉及B淋巴细胞和滤泡辅助T淋巴细胞生发中心的形成，并可能将其募集到RA的滑膜组织中，产生抗体并引起局部免疫反应从而促进RA的发生发展，而血清CXCL13可作为RA早期炎症的生物标志物，通过监测CXCL13，可以让RA患者获得早期诊断及治疗，从而取得更佳的治疗效果[57-60]。

综上所述，本研究通过生物信息学方法整合RA表达谱芯片，分析其DEGs之间的相互关系，确定了可能导致RA发生发展的信号通路和Hub基因，其可能是RA潜在的治疗靶点。这些发现可能有助于提高研究人员对RA潜在分子机制的认识，为RA发病机制的进一步研究提供可能的参考依据。