類風濕關節炎外周血單核細胞circRNA-miRNA-mRNA網絡構建及關鍵基因預測

郭錦晨，劉健，王莖，周巧，黃旦

（1.安徽中醫藥大學新安醫學教育部重點實驗室，安徽 合肥 230038；2.安徽中醫藥大學第一附屬醫院風濕病科，安徽 合肥 230031）

類風濕關節炎（rheumatoid arthritis，RA）是一種病因尚未完全闡明以滑膜炎和血管炎為基本病理改變的慢性全身性炎癥性自身免疫性疾病[1]，其病理機制尚未完全闡明，暫無根治之藥，因此探索可靠、準確的生物學標志物對其早期診斷及療效評估有重大意義。研究證實[2]，單核細胞/巨噬細胞可通過分泌腫瘤壞死因子-α（TNF-α）、白細胞介素1β（IL-1β）、白細胞介素6（IL-6）等炎性細胞因子來維持RA 患者的炎性環境，還可引起大量的免疫細胞遷移至關節炎，最終導致關節發生不可逆的損傷。研究發現[3]，環狀RNA（circular RNA，circRNAs）在RA 患者外周血單核細胞中存在差異表達，部分學者認為可能成為診斷RA 的新型生物標志物，微小核糖核酸（microRNA，miRNA）在免疫細胞分化、成熟、穩態調節方面發揮著重要的作用，部分miRNAs 已證實與T、B淋巴細胞的激活、炎癥因子和趨化因子的產生等自身免疫過程相關。circRNA 富含miRNA 的結合位點，大量研究表明可以起到海綿作用作為競爭內源ceRNAs，通過其miRNA 反應元件與miRNA 競爭性結合，然后間接調節miRNA 靶向信使RNA（mRNA）的表達，與RA 的發生發展密切相關[4]。本研究通過對基因表達綜合數據庫（Gene Expression Omnibus，GEO）中RA 患者和正常對照者的外周血單核細胞測序數據進行分析，研究circRNA、miRNA 及mRNA 差異表達譜，并構建內源性競爭RNA（ceRNA）網絡，為探索RA 診斷和治療的新靶點提供理論依據。

1 資料與方法

1.1 資料來源 在GEO 據庫（http://www.ncbi.nlm.nih.gov/geo）下載RA 患者外周血單核細胞相關數據集：circRNA數據集（GSE189338）、miRNA數據集（GSE124373）和mRNA 數據集（GSE55457）。采用“Limma”軟件包校正微陣列數據表達水平，采用“edgeR”軟件包對RA 組與正常對照組樣本進行差異表達分析，采用“pheatmap”軟件包進行熱圖繪制，LogFC 過濾閾值為1，校正后P（false discoveryrate，FDR）＜0.05 設定為篩選標準，如果表達數值很大，需要對數據取Log2。

1.2 circRNA-miRNA-mRNA 網絡的構建 根據組織樣本信息，將circRNA 數據集GSE189338、miRNA數據集GSE124373 和mRNA 數據集GSE55457 分別分為RA 組和正常對照組，在3 個數據集中分別提取差異表達的circRNA、miRNA 和mRNA。將RA患者差異表達的circRNA 通過CSCD 數據庫（https://gb.whu.edu.cn/cscd）數據庫預測相應的靶標miRNA，結合GSE124373 中篩選得到的差異表達miRNA 進一步篩選目標miRNA。然后，利用miRDB（http://www.mirdb.org）和TargetScan（http://www.tar getscan.org/vert_72）數據庫預測RA 患者目標miRNA對應的mRNA，再與GSE55457 中RA 患者差異表達的mRNA 進一步取交集以識別目標mRNA。采用“venn”軟件包對各差異表達的miRNA、mRNA 與靶標取交集，基于ceRNA 理論通過Cytoscape 3.7.1 軟件構建ceRNA 調控網絡。

1.3 ceRNA 網絡節點熱圖和箱線圖分析 采用“pheatmap”軟件包對RA 患者外周血單核細胞ceRNA 網絡節點中circRNA、miRNA 和mRNA 進行熱圖繪制?！皉eshape2”包是由Hadley Wickham 開發的一個R 包，利用reshape2 包可以融化數據進行重新整合，并對長寬數據格式進行相互轉換，最后使用“ggpubr”包繪制circRNA、miRNA 和mRNA 箱線圖。

1.4 網絡中mRNA 的GO 和KEGG 功能富集分析為進一步分析mRNA 的功能和主要作用通路，將上述篩選得到的mRNA 輸入DAVID 數據庫（https://david.ncifcrf.gov/），選定物種為“homo sapiens”，設定閾值P＜0.05，將RA 關鍵基因從生物過程（biological process，BP）、細胞組成（cellular component，CC）、分子功能（molecular function，MF）3 方面進行基因本體功能（GO）富集分析。利用R 語言“clusterProfiler”包對RA 患者關鍵靶點進行KEGG（https://www.kegg.jp/）富集分析，并使用“ggplot2”包繪制氣泡圖。

2 結果

2.1 差異表達的circRNA、miRNA 和mRNA 在GSE189338 數據集發現13 個上調circRNA 和45個下調circRNA（圖1A）；在GSE124373 數據集發現24 個上調miRNA 和11 個下調miRNA（圖1B）；在GSE55457 數據集發現112 個上調mRNA 和38 個下調mRNA（圖1C）。從CSCD 數據庫中檢索相關circRNA 數據，24 個circRNA 的基本結構模式見圖2，24 個circRNA 的基本信息，包括位置、基因組長度、靶向基因名等見表1。

圖1 RA 相關差異表達circRNA、miRNA、mRNA 熱圖

圖1 RA 相關差異表達circRNA、miRNA、mRNA 熱圖（續）

表1 24 個差異表達circRNA 的基本特征

圖2 篩選出的24 個circRNA 基本模式圖

2.2 circRNA-miRNA-mRNA 網絡的構建 通過CSCD 數據庫對circRNA 進行靶向預測，篩選出1378 個靶向miRNA。將這1378 個靶標miRNA 與miRNA 數據集GSE124373 得到的35 個差異miRNA取交集，篩選出15 個目標miRNA（圖3A）。再利用TargetScan 和miRDB 數據庫對這15 個目標miRNA進行靶向預測，篩選出7972 個靶標mRNA，將這7972 個靶標mRNA 與mRNA 數據集GSE55457 得到的150 個差異mRNA 取交集，篩選出60 個目標mRNA（圖3B）。根據ceRNA 網絡關聯數據并運用Cytoscape 3.7.1 軟件，最終篩選出10 個circRNA、8個miRNA 和38 個mRNA 構建circRNA-miRNA-mRNA可視化網絡，見圖4。

圖3 差異表達miRNA 及mRNA 篩選

圖4 可視化ceRNA 網絡

2.3 ceRNA 網絡節點熱圖和箱線圖 通過對ceRNA 網絡節點進行熱圖和箱線圖分析，發現6 個上調circRNA，4 個下調circRNA，其中，hsa_circ_0009012、hsa_circ_0003695、hsa_circ_0001589、hsa_circ_0092-125、hsa_circ_0074598 在RA 患者中高表達，hsa_circ_0058230、hsa_circ_0009127、hsa_circ_0067-492、hsa_circ_0001543 在RA 患者中低表達；發現4個上調miRNA，4 個下調miRNA，其中hsa-miR-6165、hsa-miR-98-5p、hsa-miR-619-5p、hsa-miR-186-5p 在RA 患者中高表達，hsa-miR-4516、hsamiR-3613-3p、hsa-miR-766-3p、hsa-miR-6737-3p 在RA 患者中低表達；發現27 個上調mRNA，11 個下調mRNA，其中，ADAM28、TREML2、PDCD1LG2、CD38、CD79A、RASGRP1、DAZL、ST8SIA1、E2F8、KIF23、TMEM156、BARD1、MS4A1、PNOC、LAMP3、PIM2、SEMA4D 在RA 患者中高表達，LGALSL、CYTH3、ADM2、KCNK3、TRHDE、MXD1、NR4A3 在RA 患者中低表達，見圖5、圖6。

圖5 ceRNA 網絡節點circRNA、miRNA、mRNA 熱圖分析

圖5 ceRNA 網絡節點circRNA、miRNA、mRNA 熱圖分析（續）

圖6 ceRNA 網絡節點circRNA、miRNA、mRNA 箱線圖分析

2.4 靶點通路分析 將上述篩選出的38 個差異mRNA通過DAVID 數據庫進行GO 功能富集分析，共得到GO 條目289 個，在BP 中，基因主要富集于B cell activation、B cell differentiation、lymphocyte differentiation、mononuclear cell differentiation、B cell prolifera tion 等252 個條目中；在CC 中，基因主要富集在external side of plasma membrane、RNA polymerase II transcription repressor complex、lamellar body、pericentriolar material、ruffle 等5 個條目中；在MF 中，基因主要富集在fibronectin binding、cytokine binding、receptor ligand activity、signaling receptor activator activity 等32 個條目中；各類別排名前10 的GO 條目繪制氣泡圖，見圖7。KEGG 通路富集篩選得到11 條信號通路（P＜0.05），見圖8。

圖7 差異表達mRNA GO 富集分析

圖8 差異表達mRNA KEGG 富集分析

3 討論

RA 是一種病因和發病機制尚未明確的慢性炎癥性自身免疫性疾病，免疫炎癥反應可破壞全身多關節軟骨和骨骼，氧化應激反應可進一步加重炎癥組織損傷，最終導致患者喪失勞動力、甚至致殘[5]。CircRNA 是一類有較高的穩定性和進化保守性的內源性表達的具有閉合環狀結構非編碼RNA，越來越多的證據表明[6]，包括circRNAs、miRNA 和mRNA在內的調控網絡在RA 的發病過程中起著重要的調節作用，參與自身免疫和炎癥反應的調節。如Li B等[7]研究發現，hsa_circ_0001859 可以抑制miR-204/211 活性而上調激活轉錄因子2（activating transcription factor 2，ATF2），從而促進RA 的炎癥發展。有專家對RA 患者PBMCs 進行高通量測序發現，在RA中很多差異表達的circRNA 可能是RA 的潛在診斷生物標志物，并且在RA 的發病過程中可能起調控作用[8]。因此，本研究基于ceRNA 理論構建circRNA-miRNA-mRNA 網絡篩選參與調控RA 的潛在mRNA，為RA 臨床治療及診斷提供新靶標。

本研究中首先從GEO 數據庫獲得RA 相關circRNA、miRNA 和mRNA 測序數據集進行差異分析，篩選出差異表達circRNA、miRNA、mRNA 分別有58 個、35 個、150 個?；赾eRNA 理論進一步構建circRNA-miRNA-mRNA 網絡，篩選出10 個circRNA、8 個miRNA 和38 個mRNA 組成的ceRNA 調控網絡。通過對ceRNA 網絡節點進行箱線圖分析，發現hsa_circ_0009012、hsa_circ_0003695 等circRNA在RA患者中高表達，hsa_circ_0058230、hsa_circ_00091273 等circRNA 在RA 患者中低表達，hsa-miR-6165、hsa-miR-98-5p 等miRNA 在RA患者中高表達，hsa-miR-4516、hsa-miR-3613-3p 等miRNA 在RA 患者中低表達，ADAM28、TREML2、PDCD1LG2、CD38、CD79A、RASGRP1、DAZL 等mRNA在RA 患者中高表達，LGALSL、CYTH3、ADM2、KC NK3、TRHDE 等mRNA 在RA 患者中低表達。通過對差異基因集進行GO 和KEGG 富集分析，發現分子功能主要富集在纖維蛋白結合、細胞因子結合、受體配體活性、信號受體激活劑活性等，生物過程主要富集在B 細胞活化、分化及增殖、淋巴細胞分化、單核細胞分化等。KEGG 通路主要富集在造血細胞譜系、原發性免疫缺陷、B 細胞受體信號通路、Toll 樣受體信號通路、T 細胞受體信號通路等。

CD38 一種廣泛分布于不成熟的造血細胞和活化的T、B 淋巴細胞表面的Ⅱ型穿膜糖蛋白，參與T細胞的共刺激，調控B 細胞的成熟、分化與凋亡，在細胞活化、信號轉導、介導細胞因子產生等過程中發揮重要的作用[9]。研究發現[10]，RA 患者外周血淋巴細胞中CD38+細胞表達顯著升高，CD38 能夠調節IL-1α 和IL-1β 的分泌水平介導RA 免疫紊亂。CD79A 是參與B 淋巴細胞信號轉導的基因，富集于T 細胞受體復合物，可與B 細胞受體結合并在識別抗原后產生信號[11]。RAS 鳥嘌呤釋放蛋白1（RAS guangreleasing protein 1，RASGRP1）是一種作用于MAPKs 和T 細胞受體等信號通路的RAS 鳥氨酸交換因子，可結合鈣離子參與細胞內外信號傳遞[12]。E2F8 通過介導CyclinD1 在調節細胞周期、增殖、分化及凋亡和DNA 損傷修復過程中發揮重要作用，參與Notch 與NF－κB 信號通路的調控[13]。溶酶體相關膜蛋白3（LAMP3）是自噬和促進脂質生成過程中不可缺少的一個分子，在特定的細胞核分化階段表達，在脂代謝、細胞降解途徑起重要作用[14]。PIM2 一種對RA 有抑炎作用的絲/蘇氨酸蛋白激酶，在RA 中發現外源性過表達PIM2 能夠活化哺乳動物雷帕霉素靶蛋白1（mam malian target of rapamycin 1，mTORC1），從而阻斷NF-κB 和COX-2 等信號通路介導的炎癥反應[15]。信號素4D（Semaphorin 4D，Sema4D）是一個跨膜型同源二聚體蛋白質，在神經細胞導向、免疫調節、骨代謝及血管再生等過程中發揮重要作用，免疫系統所有T 細胞、活化B 細胞均表達Sema4D，血清Sema4D 過表達可能與RA 患者肺間質病變（interstitial lung disease，ILD）的發生存在一定聯系[16]。NR4A3 屬于核激素受體（NHRs）超家族中一種具有轉錄激活與抑制的雙向轉錄調控因子，能夠參與調節細胞生物學行為、炎癥反應、應激反應、能量代謝、血內皮細胞損傷等[17]。TLRs 是一個對適應性免疫、固有免疫都具有重要作用的單跨膜受體家族，Toll 樣受體信號通路由髓樣分化因子（MyD88）依賴與非依賴性通路組成，均經NF-κB 誘導炎性細胞因子而產生[18]。B 細胞受體信號通路在RA 中起著重要作用，可分泌TNF-α、IL-6 等細胞因子、CXCL13 等趨化因子引發關節炎癥，還可充當抗原遞呈細胞活化抗原特異性T 細胞，導致關節組織炎性損傷。T 細胞受體信號通路屬于免疫調節的重要通路之一，活化的CD4+T 細胞與人類白細胞抗原（HLA）、組織相容性復合體（MHC-Ⅱ）等分子相互作用可誘導PI3K-Akt 信號通路，導致CD4+T 細胞的成熟，從而激活CD8+T 細胞抗原促進炎癥發生。此外，CD4+T 細胞的凋亡異常還會導致RA 患者大量分泌致炎因子，加重炎癥反應，促進疾病活動。

綜上所述，RA 是一種由多種因素導致的免疫系統疾病，發病機制復雜。本研究篩選出RA 患者外周血單核細胞中差異表達的lncRNA、circRNA和mRNA，構建了ceRNA 調控網絡，分析了差異表達mRNA 的生物學功能和調控通路情況。