生物信息学预测类风湿关节炎核心基因与互作miRNA

丁 晓，郝 颖，王维山，孟德峰，王梦雨，王 超

类风湿性关节炎(rheumatoid arthritis，RA)是一种未知病因的自身免疫性疾病，导致关节疼痛、肿胀、僵硬和功能丧失，严重影响患者的生活质量[1-2]。在成年人中，RA患者约占0.1%～0.5%[3]。肿瘤坏死因子在RA发病过程中具有关键作用，。目前临床上肿瘤坏死因子抑制剂英夫利昔单抗已广泛用于RA治疗。然而，仍有30%的患者对肿瘤坏死因子不具有敏感性。因此对于RA 的治疗仍需寻找预测性的生物标志物或治疗靶点。有研究[4]表明滑膜炎是RA发病的重要机制。滑膜细胞中的成纤维细胞样滑膜细胞可能到导致RA患者关节的局部损伤，RA患者关节疼痛与功能降低与滑膜炎有重要的关联。然而，目前对与RA发病机制相关的滑膜中基因和蛋白质表达的机制仍然知之甚少，这对于改善RA的预防、诊断和治疗是必要的。因此，该研究探索RA发病相关的核心基因，并预测与核心基因相互作用的miRNA，为RA的诊断及治疗药物的研制提供较为可靠的路径及作用靶点。

1 材料与方法

1.1 基因芯片数据通过GEO数据库下载GSE55235滑膜组织基因芯片原始数据。该芯片为Affymetrix Human Genome U133A芯片，该芯片包含10个正常滑膜数据(GSM1332201，GSM1332202，GSM1332203，GSM1332204，GSM1332205，GSM1332206，GSM1332207，GSM1332208，GSM1332209和GSM 1332210)和10个RA患者滑膜数据(GSM1332221, GSM1332222,GSM1332223,GSM1332224,GSM1332225,GSM1332226,GSM1332227,GSM1332228,GSM1332229,and GSM1332230)

1.2 筛选差异基因利用R语言3.5.0对芯片进行数据标一化，通过R语言中的limma包筛选RA和正常滑膜组织的差异表达基因，采用t检验分析RA和正常组织之间的表达，以调整后P<0.05, |log2FC|>2作为筛选条件。对差异基因进行聚类分析，利用R语言gplots包绘制热图。

1.3 差异基因GO功能富集和KEGG通路富集分析DAVID数据库(https://david.ncifcrf.gov/)是任何高通量基因功能分析的重要基础。利用DAVID对差异基因进行功能注释并对编码蛋白的功能和通路富集分析。本研究分析了RA芯片数据中显著上调和下调的差异基因，P<0.05为差异有统计学意义。将差异表达基因上传至DAVID6.8 数据库进行GO和KEGG分析。

1.4 差异基因蛋白互作网络分析及基因模块分析STRING数据库(http://string-db.org/)是一种软件系统，通常用于识别已知蛋白质和预测蛋白质之间的相互作用。每个节点是基因，蛋白质或分子，节点之间的连接代表这些生物分子的相互作用，其可用于鉴定RA中DEGs编码的蛋白质之间的相互作用和途径关系。中心节点中的相应蛋白质可以是具有重要生理调节功能的核心蛋白质或关键候选基因。将差异表达基因上传至String 10.5数据库进行分析，以综合分数(combined score)>0.4为筛选标准[5]。将String 10.5数据库结果导入Cytoscape v3.6.1插件MCOD进行模块分析分析以挖掘PPI中连接最为紧密的模块。筛选标准为degree cutoff=2，node score cutoff=0.2，K-score=2，max deapth=100，P<0.05差异有统计学意义。

1.5 miRNA-核心基因调控网络的预测使用Cytoscape v3.6.1自带的插件CyTargetLinker预测调控核心基因的miRNA，通过(https://projects.bigcat.unimaas.nl/cytargetlinker/linksets/)下载人类基因数据库，其中包括基于实验验证的miRTarBase 7.0数据库(包含502 652个互作网络)和基于预测功能的TargetScan 7.2数据库(包含264563个互作网络)，来预测核心基因与miRNA之间的调控关系。

2 结果

2.1 差异表达基因将原始数据进行处理，按照筛选条件发现RA患者滑膜和正常滑膜之间有605个差异表达基因，其中上调基因有314个，下调基因有291个。见图1。利用R语言gplots包对差异基因类聚分析，将差异基因可视化。见图2。

图1 RA差异表达基因的火山分析

2.2 差异基因GO分析和KEGG分析对上调基因进行GO分析，按照计数值从大到小排序，分别选取生物过程(biological process，BP)、分子功能(molecular function，MF)与细胞组分(cellular component，CC)中前5项进行分析(表1)。上调基因主要富集在免疫反应、信号转导、炎症反应、先天免疫反应和细胞黏附等生物学过程，参与质膜、膜的整体组成部分、细胞外外泌体、细胞质和细胞外区域的分析构成，影响蛋白质结合、蛋白质同源二聚化活性、同源蛋白结合、受体结合和钙离子结合等分子功能。下调基因未得到相关富集信息。对差异基因KEGG分析发现上调基因主要富集在结核、吞噬、趋化因子信号通路，细胞因子-细胞因子受体相互作用、破骨细胞分化、细胞黏附分子(cell adhesion molescules, CAMs)、利什曼病、金黄色葡萄球菌感染、甲型流感和癌症通路等相关信号通路。见图3。

图2 热图中的40个差异表达基因(包括20个上调基因和20个下调基因)红色：上调基因；绿色：下调基因

表1 上调基因的GO功能分析

分类条目差异基因数百分比(%)P值BPGO:0006955～免疫反应5618.065.49E-33BPGO:0007165～信号传导5016.132.95E-09BPGO:0006954～炎症反应4514.522.84E-24BPGO:0045087～先天免疫反应3410.974.59E-13BPGO:0007155～细胞粘附3110.003.00E-10CCGO:0005886～质膜12038.711.42E-11CCGO:0016021～膜的整体组成部分11737.745.34E-05CCGO:0070062～细胞外外泌体10032.261.44E-14CCGO:0005829～细胞质7423.874.01E-03CCGO:0005576～细胞外区域6019.353.11E-09MFGO:0005515～蛋白质结合17556.457.83E-04MFGO:0042803～蛋白质同源二聚化活性247.742.87E-03MFGO:0042802～相同蛋白结合216.772.75E-03MFGO:0005102～受体结合185.811.01E-04MFGO:0005509～钙离子结合185.819.60E-02

图3 差异表基因显著参与的KEGG信号通路

2.3 差异基因蛋白互作网络差异基因经String数据库构建蛋白互作网络分析，经Cytoscape v3.6.1计算出个基因的连接度(degree)，共包括552个节点和5163条相互作用边。见图4。Degree表示网络中基因与周围基因相连接数量，因此degree越大表示与它相关的基因数量就越多。因此筛选出排名前10的degree基因定位核心基因，分别为蛋白酪氨酸磷酸酶受体C(protein tyrosine phosphatase receptor type C, PTPRC, degree=129)、血管内皮生长因子A(vascular endothelial growth factor A, VEGFA, degree=119)、纤连蛋白1(fibronectin 1, FN1, degree=111)、整合素αM(integrin alpha M, ITGAM, degree=109)、表皮生长因子(epidermal growth factor receptor, EGFR, degree=102)、CD86(degree=102)、基质金属蛋白酶9(matrix metalloprotein，MMP9，degree=100)、整合素亚基β2(integrin subunit Beta 2, ITGB2, degree=199)、蛋白酪氨酸激酶结合蛋白(protein tyrosine kinase-binding protein,TYROBP, degree=97)、MYC(degree=89)。

通过Cytoscape v3.6.1插件MCODE进行子模块分析，共得到13个子模块，按照节点数量进行排序，选取前4个模块并对各模块参与基因进行通路分析。见图5。模块1主要富集在破骨细胞分化、PI3K-Akt信号通路、FoxO信号通路及类风湿关节炎等相关信号通路。

2.4 miRNA-核心基因网络模块的筛选通过CyTargetLinker预测上述10个核心基因的相互作用的miRNA，在miRTarBase 7.0数据库中有346个预测miRNA靶点互作关系，在TargetScan数据库中有127个预测的miRNA 靶点互作关系，总共有347个节点和473个边。通过调节数据库叠加阈值，控制调控网络显示结果，一般阈值设置为1～3。本文以2为阈值，共得到36个miRNA和3个核心基因存在互作关系。见图6。

图4 差异表达基因蛋白互作网络

图5 差异基因蛋白网络图的子网络分析(前4个模块)

3 讨论

本研究通过R语言对类风湿关节炎滑膜芯片进行分析获得相关差异基因，并通过对差异基因进行GO富集分析以及KEGG通路富集分析，为类风湿性关节炎的机制研究进一步提供研究方向及依据，通过筛选核心基因并预测互为作用的miRNA为RA的研究提供一种新的思路。

本研究共得到605个差异表达基因，其中有314个基因RA中表达上调，有291个基因在RA中表达下调。由于差异基因相对过多，只列举了前20个上调基因和前2个下调基因热图。对差异基因进行GO功能富集分析，在生物学过程中，差异基因主要参与免疫应答反应、炎症反应、信号传导等，这引起自身抗体的增强。在分子功能上，差异基因增强了受体活性，蛋白质结合和蛋白质同源二聚化活性。KEGG通路富集分析鉴定出吞噬，趋化因子信号通路，细胞因子-细胞因子受体相互作用，破骨细胞分化，CAMs等相关信号通路表明免疫反应在RA中具有重要的作用。

本研究还筛选了RA的10个关键基因，为RA 的诊断提供一些方向，degree最大的是PTPRC基因，与二肽基肽酶-4结合共同作为T细胞激活因子。蛋白酪氨酸蛋白磷酸酶是通过T细胞抗原受体激活T细胞所必需的。早期PTPase结构域具有酶活性，T细胞活化导致Src激酶相关磷蛋白1(src kinase-associated phosphoprotein 1,SKAP1)和非受体酪氨酸蛋白激酶的募集和去磷酸化[6-7]。有研究[8]表明PTPRC参与RA的发病进展，上调PTPRC可能有利于维持滑膜的正常功能。VEGFA在RA患者滑液和血清中大量表达，对滑膜血管翳的形成有特异作用[9]。FN1编码纤维蛋白，一种以血浆中可溶性二聚体形式存在的糖蛋白，有研究[10]表明在RA患者血液中FN1的表达增高，可为临床诊断提供依据。ITGAM和ITGB2是整合素家族的两种，ITGAM 编码整合素αMβ2的a链(即CD11b)，ITGB2码整联蛋白β链，广泛表达于多种免疫细胞亚群，在调节免疫耐受中具有重要的作用[11]。破骨细胞是导致RA关节破坏的关键之一，EGFR在破骨细胞形成及滑膜炎症过程中具有重要的作用[12]。有研究[13]表明EGFR在RA滑膜及关节液中表达。共刺激/抑制分子可以调控T细胞的活化，在免疫系统的动态平衡调节具有重要的作用。有研究[14]表明RA的发生与共刺激/抑制分子信号异常有关。CD86是最主要的共刺激分子之一，主要诱导Th2介导的体液免疫和Ig产生。CD86介导的共刺刺激是激发T细胞反应的关键。MMP9具有很强的降解软骨细胞外基质的作用，MMP9的表达增高对于RA中软骨的破坏具有重要的影响作用。c-myc是核内转录因子类原癌基因，其基因产物过度表达可以细胞增殖，c-myc的过度表达是刺激滑膜细胞和新生血管形成的重要原因。

图6 核心基因与miRNA互作网络

此外本研究还预测了与核心基因相互作用的miRNA。从miRTarBase 7.0数据库和TargetScan数据库仅预测到核心基因FN1、VEGF和EGFR相关的36个miRNA，为RA的研究提供新的角度。本研究通过疾病基因水平的差异，借助miRTarBase7.0数据库和TargetScan数据库，将基因水平与蛋白转录调控水平联系在一起，可以更深一层次的了解RA发病机制，提供多方面RA的治疗靶点，然而本研究得到的结论需要进一步的实验验证，随着研究的深入，PTPRC、VEGFA、FN1、ITGAM、EGFR、CD86、MMP9、ITGB2、MYC、TYROBP及核心基因FN1、VEGF和EGFR相关的36个miRNA或许可以应用于RA的病理筛查以及药物研发的靶点。由于RA的机制相对复杂，因素相关较多，探寻RA相关实际意义的基因，仍是后续研究所思考的问题。