生物信息學預測類風濕關節炎核心基因與互作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相關實際意義的基因，仍是后續研究所思考的問題。