基于生物信息學分析探討類風濕關節炎與骨關節炎之間的關系①

章曉云 李華南 陳 鋒 柴 源 甘 斌 陳丁鵬 （江西中醫藥大學，南昌 330004）

類風濕關節炎（rheumatoid arthritis，RA）是一種自身免疫性疾病，可引起慢性炎癥表現、肥胖、心血管疾病等，影響全球0.3%～1.0%的人口，在中國每年至少500 萬人被診斷出RA［1-2］。骨關節炎（osteoarthritis，OA）是一系列能引起軟骨損傷和變性的退行性關節炎的總稱，影響了全球約2.4億人口，無論是從生活質量還是疾病負擔方面，都已成為全球范圍內嚴重的公共衛生問題［3-5］。兩者是目前最常見的關節炎類型，病因均未明確，暫無有效的根治方法，給個人和社會帶來了巨大的經濟壓力。

盡管兩種疾病的病因學不同，但均以滑膜炎癥反應、關節軟骨破壞、骨質流失等為特征，最終導致關節進行性退變，臨床上許多患者同時患有這兩種疾病。研究表明RA 是OA 發病的重要危險因素之一，RA 患者罹患OA 的風險比正常人高2.75 倍［4］。隨著基礎研究的不斷深入，已經發現兩種疾病之間存在諸多共同發病機制和受調控因子，在蛋白質水平和組織炎癥評分水平上具有很強的連續性［6-8］。因此，本研究通過GEO數據庫查找RA和OA血清基因芯片表達譜進行生物信息學分析，篩選出差異miRNA 并進行靶基因預測，最后交集篩選出核心mRNA，為確定RA 和OA 之間的共同發病機理創造條件，并為RA 和OA 治療藥物的研發提供有效的作用靶點。

1 資料與方法

1.1 數據收集 以“Rheumatoid arthritis”“Osteoar?thritis”“miRNA”為關鍵詞，在GEO 公共數據庫（https：//www.ncbi.nlm.nih.gov/geo/）［9］檢索與RA 及OA 相關的芯片，獲得編號為GSE115885 的RA 血清基因表達芯片，共包含15 例RA 患者和5 例健康對照，其所處平臺為GPL25134；獲得編號為GSE105027 的OA 血清基因表達芯片，共含12 例OA患者和12例健康對照，其所處平臺為GPL21575。

1.2 差異miRNA 的數據分析 利用Perl（V5.30.0.1）對數據進行重注釋，并利用R 語言校正、分類，limma 包對miRNA 進行差異分析，以P<0.05 和|logFC|≥0.5 作為過濾條件，篩選差異miRNA（differentially expressed miRNAs，DEmiRNAs）。最后運用Venny 平臺（http：//bioinfogp. cnb. csic. es/tools/venny/）將所得RA 和OA 的DEmiRNAs 進行映射取交集，得到交集DEmiRNAs。

1.3 miRNA 轉錄因子預測及功能富集分析 Fun?Rich（V3.1.3）［10］是一個用于基因和蛋白質功能富集和相互作用網絡分析的軟件。利用FunRich 軟件將篩選出的交集miRNA 進行轉錄因子預測及功能富集分析，并根據P值顯示前10位進行可視化。

1.4 靶基因預測 通過miRTarBase 數據庫（http：//mirtarbase. mbc. nctu. edu. tw/php/index. php）、Target Scan 數據庫（http：//www. targetscan. org/vert_72/）和miRDB數據庫（http：//mirdb.org/）分別預測DEmiRNA的靶基因，使用Venny 平臺取3 個數據庫均能預測到的靶基因，生成交集基因，并與步驟1.2所得的交集DEmiRNAs 構成關系網絡。整合得到miRNAmRNA 調控網絡，并導入Cytoscape（V3.7.2）繪制網絡圖。

1.5 蛋白互作網絡構建 為進一步探究miRNA 所調控靶基因的作用機制，將靶基因導入String 數據庫（http：//string-db. org/cgi/input. pl）［11］，限定研究物種為“Homo Sapiens”獲得蛋白互作關系，設置連接評分（Combined score）>0.95，再導出蛋白互作關系的數據文件，選取文件中的node1、node2 和Com?bined score 三列導入Cytoscape 中，最后利用Cyto?scape 軟件的“NetworkAnalyzer”工具進行可視化處理，構建蛋白互作（PPI）網絡圖。

1.6 基因富集分析 利用DAVID 數據庫（https：//david. ncifcrf. gov/）［12］對mRNA 進行GO 功能富集分析，研究靶基因的主要生物功能；對潛在靶點進行KEGG 信號通路富集分析，研究靶基因的主要信號通路，P<0.05 代表富集結果顯著。最后利用R 語言繪制GO富集分析圖及KEGG富集分析氣泡圖。

2 結果

2.1 DEmiRNA 的鑒定 利用R語言對芯片進行重注釋和數據校正后，再對其進行差異分析。結果顯示，與健康對照相比，RA 患者的血清中共存在41個明顯改變的miRNAs，其中上調31個，下調10個；OA患者的血清中共存在265個明顯改變的miRNAs，其中上調124 個，下調141 個；對二者取交集后共得到4個交集DEmiRNAs，詳見表1。

表1 交集DEmiRNAsTab.1 Intersection DEmiRNAs

2.2 miRNA 調控的轉錄因子及富集分析結果 通過FunRich 軟件對4 個交集DEmiRNA 進行轉錄因子預測，共得到197個轉錄因子，根據其所調控的基因數與P值，選取富集程度顯著的10 個轉錄因子進行展示，見圖1。差異最顯著的10 個轉錄因子分別為LHX3、SP4、NFIC、VSX2、HOXA7、TCF3、MYC、HOXB4、ETS1、SP1，其基本信息見表2；對4 個交集DEmiRNAs 進行功能富集分析，其生物過程主要涉及肽代謝、轉錄、翻譯等；細胞成分主要涉及受體復合體、Ⅰ型膠原、Ⅸ型膠原等；分子功能主要涉及細胞外基質結構成分、蛋白質酪氨酸/絲氨酸/蘇氨酸磷酸酶活性、輔助轉運蛋白活性等，見圖2。

圖2 交集DEmiRNAs的功能富集分析Fig.2 Functional enrichment analysis of intersection DEmiRNAs

表2 關鍵轉錄因子的基本信息Tab.2 Basic information on key transcription factors

圖1 交集DEmiRNAs的轉錄因子預測Fig.1 Transcription factor prediction of intersection DEmiRNAs

2.3 “miRNA-基因”調控網絡預測 本研究運用miRTarBase、TargetScan、miRDB數據庫預測上述4個DEmiRNAs 的靶基因，通過Venny 平臺篩選靶基因在3 個數據庫的交集，共得到433 個靶基因，見圖3。最后將得到的靶基因與交集DEmiRNAs構成關系網絡導入Cytoscape 軟件構建miRNA-mRNA 調控網絡圖，見圖4。

圖3 交集DEmiRNAs的靶基因韋恩圖Fig.3 Target gene Venn diagram of intersection DEmiRNAs

圖4 交集miRNA-mRNA調控網絡Fig.4 Intersection miRNA-mRNA regulatory network

2.4 PPI網絡 借助String數據庫和Cytoscape軟件構建PPI網絡，見圖5。圖中共涉及40 個節點、43 條邊，其中節點代表蛋白，邊則代表蛋白間的互作關系，節點越大、顏色越深則度值越大；邊越粗，則蛋白間的關系越緊密。度值排名前5 的蛋白質為CDC5L、HNRNPU、GATAD2A、CHD4、ACTB。這些度值較大的蛋白質在整個網絡中發揮關鍵作用，可能是RA 跟OA 發生發展的關鍵基因，其基本信息見表3。

表3 關鍵靶點的基本信息Tab.3 Basic information on key targets

圖5 蛋白互作網絡Fig.5 Protein interaction network

2.5 GO 功能與KEGG 信號通路富集分析 將靶基因導入DAVID 數據庫，分析獲得疾病基因和蛋白信息所富集的通路，并將基因信息進行可視化處理，見圖6。結果顯示，靶基因的生物過程主要涉及mRNA代謝過程的調控、細胞周期過程的負調控、有絲分裂細胞周期等；細胞成分主要涉及轉錄因子復合物、轉移酶復合物、蛋白激酶復合物等；分子功能主要涉及泛素蛋白轉移酶活性、蛋白激酶調節劑活性、DNA 解旋酶活性等；信號通路主要涉及調控干細胞多能性的信號通路，Hippo信號通路、FoxO 信號通路、Apelin信號通路、p53信號通路等。

圖6 GO和KEGG信號通路可視化結果Fig.6 Visualization results of GO and KEGG signaling pathway

3 討論

近年來，miRNA 的研究引起了人們的極大興趣，已經發現了2 000 多種miRNAs，它們在真核基因表達調控中具有廣泛作用，參與細胞分化，增殖及凋亡等多種生物過程，人類基因組中約有1/3 的基因處于miRNA的控制之下［13-14］。研究表明miRNA通過調控自身抗體與炎癥細胞因子釋放、軟骨細胞自噬凋亡、調節免疫等機制在RA 與OA 疾病的發病過程中發揮至關重要的作用［13-15］。隨著以基因芯片為代表的高通量測序技術的發展，目前可在短時間獲取疾病的相關信息，進而從基因層面進行更全面的生物標志物篩選。本研究通過生物信息學分析了解RA 與OA 之間的相關性，了解兩病的相關發病機理，為藥物研發提供一定的參考資料。

本研究通過GEO 數據庫查找獲取RA 血清基因表達芯片GSE115885 和OA 血清基因表達芯片GSE105027，對兩種芯片進行生物學信息分析共獲得4 個DEmiRNAs，分別為hsa-miR-98-5p、hsa-miR-3613-3p、hsa-miR-6858-3p、hsa-miR-4281。研究結果顯示除hsa-miR-98-5p 對兩種疾病的作用變化趨勢相反，其他三種miRNA 的作用趨勢一致。相關研究表明RA 與OA 疾病發生和發展過程中涉及多種轉錄因子，主要涉及滑膜炎癥反應，軟骨與骨組織的破壞［16］。本研究對4 個交集DEmiRNA 轉錄因子預測結果顯示與其相關的轉錄因子共有197 個，差異最顯著的為LHX3、SP4、NFIC、VSX2、HOXA7、TCF3、MYC、HOXB4、ETS1、SP1。研究表明miRNAs參與了幾乎所有已知細胞過程的轉錄后調控機制，包括發育、分化、凋亡以及對病原體感染的先天性和適應性免疫反應，而轉錄因子與RNA 聚合酶Ⅱ形成的轉錄起始復合體在轉錄起始過程中具有重要的作用［17-18］。本研究結果對交集DEmiRNAs 進行功能富集分析顯示其生物過程主要涉及肽代謝、轉錄、翻譯等，與既往研究完全符合［17-18］。因此，筆者認為這4種交集DEmiRNAs與RA 和OA 的發病過程中有著緊密的關系，可以為闡明RA 與OA 之間的關系以及后續研究提供依據。

對自身靶基因的調控是miRNAs 在疾病發病過程中發揮作用的研究重點。本研究結果顯示4 個DEmiRNAs 共調控433個靶基因，其中有359個靶基因為共同調控，其中CDC5L、HNRNPU、GATAD2A、CHD4、ACTB 排名較靠前，在整個網絡中起著關鍵的作用。CDC5L 參與細胞周期調控，影響炎癥及細胞凋亡過程，因此其可通過調節炎癥細胞及軟骨細胞的凋亡從而影響RA 與OA 疾病的發生與發展［19-20］。研究表明HNRNPU 參與免疫基因轉錄調節過程，對RA 發病具有關鍵作用［21-22］。HNRNPU 可與功能性基因間重復RNA 元件相互作用，增強mRNA的穩定性，可促進炎癥基因的表達水平，從而可減輕RA 與OA 患者的炎癥反應［23］。CHD4是低氧誘導因子（hypoxia inducible factor，HIF）的共激活因子，CHD4 可增強HIF1α 的募集，從而促進HIF 目標基因轉錄［24］。研究表明HIF1α 介導了免疫炎癥細胞的激活，在RA 患者滑膜中表現出缺氧性并呈現HIF1α 上調［25］，HIF1α 也被認為是OA 的關鍵調節因子，因此筆者認為CHD4 可通過與HIF1 相互作用從而影響兩種疾病的發生。ACTB 參與軟骨細胞自噬與凋亡過程，從而影響RA 與OA 患者軟骨退變情況［26］。目前對于GATAD2A 的研究主要集中在癌癥疾病與精神分裂疾病，其可影響患者精神狀態［27］，尚未有關于其在RA 與OA 疾病中的研究，但RA 與OA 患者因為疾病的難治性可影響患者的心理及精神情緒，因此未來對于GATAD2A 與RA 與OA 疾病之間的關系有待進一步研究。

為充分了解受DEmiRNAs 調控的mRNA 在RA與OA 中涉及的共同通路及功能，本研究進行了GO及KEGG 分析。GO 結果顯示靶基因主要涉及mRNA代謝過程的調控、細胞周期過程的負調控、有絲分裂細胞周期等生物過程，主要涉及泛素蛋白轉移酶活性、蛋白激酶調節劑活性、DNA 解旋酶活性等分子功能，這些生物過程及分子功能均與炎癥細胞與軟骨細胞的增殖與分化，自噬與凋亡有關，從而影響RA 與OA 的炎癥變化與關節退變程度。KEGG 分析結果顯示信號通路主要涉及調控干細胞多能性的信號通路，Hippo 信號通路、FoxO 信號通路、Apelin 信號通路、p53 信號通路等。研究發現通過對誘導性多能干細胞進行編程，可使其具有對炎癥刺激作出反應的能力，并產生強大的自主調節的抗炎細胞因子，因此調控干細胞多能性的信號通路是未來研究RA 與OA 發病機制的主要途徑［28］。研究表明Hippo 信號通路及p53 信號通路均可影響細胞增殖和凋亡，抑制其表達可減少軟骨細胞凋亡以維持軟骨形態，FoxO 信號通路可調節軟骨細胞自噬和防御氧化應激，從而有效改善RA 與OA 患者的病情，減緩關節退變［29-31］。Apelin 可通過自分泌/旁分泌和內分泌途徑有效地調節炎癥，調節軟骨、滑膜、骨骼等組織和各種免疫細胞的變化，對于RA 和OA的發病機制至關重要［32］。

總而言之，RA 和OA 均存在龐大的疾病基因網絡，但這些網絡并非各自獨立，而是存在密不可分的聯系。本研究獲得RA 和OA 各自疾病基因網絡的交叉部分正是二者相關聯的部分，也是研究RA和OA 相互關系的突破點，同時也可為藥物同時干預兩種相互關聯的疾病提供較為可靠的路徑和作用靶點。且從側面反映出不同疾病可通過應用同一藥物治療的特色，與中醫辨證論治中的“異病同治”思想不謀而合，以此促進中藥現代化發展。此外本研究尚存在一些缺點：首先，由于篩選條件的限制，只能對主要轉錄因子、作用靶點、信號通路進行分析，在一定程度上使得研究結果具有局限性；其次，對差異基因信息進行整合處理依賴于生物信息學技術的發展，疾病數據庫的完整性、準確性直接決定著整合后信息的可靠性。故本文雖然能夠篩選出大量靶點和通路，但仍需后續結合體內外實驗進一步驗證與支持，使理論更加可靠。