基于生物信息學探討骨關節炎患者LncRNA表達及CeRNA互作網絡

陳鋒 章曉云,2 陳躍平 廖建釗 周毅 李嘉駿 吳煒鋒

(1廣西中醫藥大學附屬瑞康醫院，廣西 南寧 530011；2江西中醫藥大學)

骨關節炎(OA)是一種以關節軟骨退變鈣化、軟骨下骨重塑、骨贅形成、滑膜炎癥和血管侵襲為特征的慢性退行性疾病〔1～3〕。伴隨全球人口老齡化的不斷發展，OA發病率逐年遞增，已成為全球范圍內的公共健康問題〔4〕。研究表明，OA的發生與創傷、遺傳、衰老、肥胖等多種因素有關，但具體發病機制尚不明確。目前臨床上非手術治療方法包括口服硫酸氨基葡萄糖、皮質類固醇、透明質酸鈉、多西環素和基質金屬蛋白酶抑制劑等，雖然這些藥物均有不同程度的療效，但均無法完全阻斷OA的進展〔5,6〕；而手術治療雖能矯正關節畸形、改善關節功能，但存在術后感染、血栓和二次手術等潛在風險〔7,8〕。因此，全面挖掘OA的潛在發病機制，尋找新的治療靶點，具有非常重要的意義。

長鏈非編碼RNA(LncRNA)是一類長度大于200nt且不能編碼蛋白的RNA。首次發現時認為其沒有編碼功能，被定義成遺傳暗物質。但近年研究證實LncRNA在生命活動過程中發揮著重要的調控作用，與疾病的發生發展密切相關〔9〕。研究表明LncRNA在OA發展過程中可作為競爭性內源RNA(CeRNA)調節微小RNA(miRNA)的表達模式和生物學特征〔10～12〕。即LncRNA可以海綿吸附miRNA來影響信使RNA(mRNA)的表達，從而參與OA的免疫反應、炎癥反應、細胞代謝等生物學過程。因此本研究通過GEO數據庫挖掘與OA相關的芯片,借助R語言分析差異基因,并構建LncRNA-miRNA-mRNA網絡來深入研究LncRNA在CeRNA調控網絡中的作用機制。

1 資料和方法

1.1篩選OA的相關芯片 以“Osteoarthritis”、“Genome-wide Expression Profiles”為關鍵詞，在GEO數據庫中檢索與OA相關的芯片,獲得編號為GSE51588的芯片矩陣文件和編號為GPL13497的基因注釋文件,該芯片中包含50個樣本,包括40例OA患者的關節組織和10例健康對照。

1.2芯片數據分析 利用Perl對數據進行基因重注釋，并利用R語言校正、分類，limma包對編碼蛋白質的mRNA及LncRNA進行差異分析，設置差異倍數(FC)>0.5，P<0.05，最后運用pheatmap包對存在差異的LncRNA繪制差異熱圖。

1.3LncRNA、miRNA、mRNA互作預測 利用mircode數據庫篩選與差異LncRNA互作的miRNA。再從公共數據庫(miRTarBase、TargetScan、miRDB)預測與miRNA互作的靶基因，取3個數據庫均能預測到的靶基因，與步驟1.2所得的差異mRNA取交集。整合得到LncRNA-miRNA-mRNA互作網絡，并導入Cytoscape繪制CeRNA網絡圖。

1.4蛋白互作網絡構建 為進一步探究OA的Hub基因，將差異基因導入String數據庫中，限定研究物種為“Homo Sapiens”獲得蛋白互作關系，設置連接評分>0.97，再導出蛋白互作關系的數據文件，選取文件中的node1、node2和連接評分三列導入到Cytoscape中，最后利用Cytoscape軟件的“Network Analyzer”工具進行可視化處理，構建PPI網絡圖，并根據度值篩選出Hub基因。

1.5Hub基因與核心LncRNA的相關性分析 選取CeRNA調控網絡中度值最高的LncRNA與其互作的Hub基因，保留各自對應的疾病組數據，利用R語言對數據進行校正后，設置相關系數的絕對值>0.4，P<0.01，最后運用limma包對存在統計學意義的基因繪圖。

1.6GO與KEGG富集分析 利用DAVID數據庫對差異基因進行GO功能富集分析，研究OA的主要生物功能；對差異基因進行KEGG通路富集分析，研究OA的主要信號通路，P<0.05代表富集結果顯著。最后利用R語言繪制GO和KEGG富集分析氣泡圖。

2 結 果

2.1差異mRNA及LncRNA 利用R語言對芯片進行基因重注釋和數據校正后，再對其進行差異分析。結果顯示，與健康對照者相比，OA患者的膝關節組織中共存在3 292個明顯改變的mRNA，其中上調1 476個，下調1 816個；共存在67個明顯改變的LncRNA，其中上調45個，下調22個。在上調和下調的LncRNA中分別選取最顯著的5個來繪制差異熱圖，見圖1。

圖1 差異LncRNA在不同樣本中的表達

2.2CeRNA網絡 利用mircode等數據庫共預測到208個miRNA和1 936個mRNA，將預測到的mRNA與芯片中的3 292個差異mRNA取交集，再刪除不在該互作關系中的miRNA、LncRNA，整合后得到LncRNA-miRNA-mRNA互作網絡。將該網絡導入Cytoscape中繪制CeRNA網絡，見圖2。網絡中共包括310個節點(54個LncRNA節點、38個miRNA節點、218個mRNA節點)和1 085條邊，每條邊表示LncRNA與miRNA、miRNA與mRNA之間的關系。該網絡中度值排名前5位的LncRNA為鉀電壓門控通道亞家族Q成員1反義鏈(KCNQ1OT)1、X 染色體失活特異性轉錄本基因(XIST)(度值36)、反義Opa反應蛋白5RNA1 (OIP5-AS)1(度值33)、淋巴細胞白血病缺失基因(DLEU)1(度值27)、慢性淋巴細胞白血病上調(CLLU)1(度值27)。

圖2 CeRNA調控網絡

2.3PPI網絡 借助String數據庫和Cytoscape軟件構建PPI網絡圖，見圖3。圖中共包含30個節點、41條邊，其中節點代表蛋白基因，邊則代表蛋白基因間的互作關系，節點大小及顏色深淺表示其度值的大小，節點越大相應的度值越大，顏色由黃變橙相應的度值則由小變大，邊的粗細反映連接評分，邊越粗，評分越高。度值排名前5的蛋白基因為信號轉導和轉錄激活劑(STAT)3(度值6)、生長因子受體結合蛋白(GRB)2(度值6)、絲裂原激活的蛋白激酶(MAPK)1(度值6)、叉頭家族(FOX)O3(度值5)、原癌基因MYC(度值5)。這些度值較大的蛋白基因在整個網絡中起著關鍵的作用，為該網絡的Hub基因。

圖3 蛋白互作網絡

2.4共表達結果 選取CeRNA中度值最高的LncRNA(KCNQ1OT1、XIST)與Hub基因(STAT3、GRB2、MAPK1、FOXO3、MYC)分析其表達的相關性。除XIST與MAPK1之間呈正相關，其余均呈負相關，選取負相關部分存在統計學意義的數據進行展示，見圖4。

圖4 共表達結果

2.5GO富集分析和KEGG富集分析 為探究更具體的功能模式,利用DAIVD數據庫對關系網絡中的218個差異基因進行GO和KEGG富集分析。GO富集分析共得到237個條目,KEGG富集分析共得到45條信號通路,根據P值排序,選取P值最顯著的前20個條目以氣泡圖的形式展現,見圖5。圖中橫坐標為該通路富集的基因占人體總基因的比例，氣泡越大代表該通路富集的基因數越多，顏色越紅則代表富集程度越顯著。分析結果顯示,骨關節炎的生物過程主要包括對DNA結合轉錄因子活性的調控、對營養水平的反應、肌細胞的增殖等；信號通路主要涉及PI3K-AKT信號通路、MAPK信號通路、mTOR信號通路等。

圖5 GO和KEGG信號通路可視化結果

3 討 論

在OA的研究過程中，針對骨代謝、基質降解、氧化應激和免疫炎癥反應等方面的治療均已應用于臨床，但其療效并不理想〔13〕。而LncRNA在介導OA發病機制中起著重要調控作用，通過基因測序技術發現在軟骨分化過程中表達量發生改變的LncRNA高達數千個〔14,15〕。因此，深入研究LncRNA在OA發生發展中的作用及其調控機制意義重大。

本研究基于GEO數據庫聯合DAVID功能富集分析軟件，通過生物信息技術篩選調控差異表達miRNA的上游LncRNA及miRNA作用的下游靶基因，旨在探討OA發展中LncRNA表現出的調控網絡模式。CeRNA網絡分析得出有54個LncRNA、38個miRNA和218個mRNA與OA發生發展相關。對差異表達的mRNA進行GO和KEGG富集分析，得出這些基因所富集的功能主要包括DNA結合轉錄因子活性的調控、營養水平的反應、肌細胞的增殖等生物學過程，涉及PI3K-AKT、MAPK、mTOR等信號通路，這些通路在OA發生發展中均起著重要的作用〔16〕。

通過分析共得到5個關鍵LncRNA(KCNQ1OT1、XIST、OIP5-AS1、DLEU1、CLLU1)，與之互作的mRNA在GO和KEGG通路分析中均顯示與OA高度相關。這些結果表明，在OA的發生發展過程中，上述LncRNA可能發揮著重要作用。研究表明，KCNQ1OT1與miR-200a存在潛在結合位點，而在CeRNA網絡中Hub基因FOXO3又是miR-200a的潛在靶基因。KCNQ1OT1可通過作用miR-200a影響FOXO3的表達，進而調控細胞自噬、衰老、凋亡和炎癥反應等〔17〕；XIST在細胞增殖、分化和基因組維持中起著關鍵作用，研究表明降低XIST的表達，可顯著降低骨髓間充質干細胞(BMSC)的增殖和遷移能力，從而導致炎癥反應發生〔18～20〕。研究表明沉默OIP5-AS1能顯著抑制細胞增殖活性，通過上調OIP5-AS1可調節miR-410的表達，并進一步調控靶基因KLF10，從而介導PI3K/AKT信號通路，促進細胞增殖〔21,22〕；DLEU1和CLLU1可在轉錄中、轉錄后、表觀遺傳等方面調控基因的表達，參與膠質瘤、骨肉瘤、宮頸癌等疾病的發展進程，但其在OA中的潛在作用尚無相關研究，這就為今后通過作用DLEU1和CLLU1治療OA提供新的方向和線索〔23,24〕。

本研究結果說明OA的作用機制無法依靠某個特定靶點來調控，且各基因間存在錯綜復雜的互作關系：LncRNA與miRNA結合，通過隔離mRNA而下調miRNA的結合，最終促進mRNA的表達。STAT3是一類參與成骨細胞增殖、分化、炎癥和免疫反應的特異性基因，通過激活STAT3信號通路不僅能刺激MSC成骨分化，還參與Cdc42等生物因子介導的軟骨損傷過程。這表明STAT3在OA的軟骨下骨硬化和軟骨損傷過程中發揮了重要作用，因此抑制STAT3表達是阻止關節損傷進入不可逆階段的重要干預措施,對治療OA意義重大〔25,26〕；GRB2作為一種銜接蛋白，經研究證實與PI3K/AKT/NF-κB信號通路的轉導存在密切關聯，沉默GRB2能夠抑制上述通路的活化，進而減少組織的炎性浸潤發揮保護作用，GRB2有望成為OA治療的重要靶點〔27〕；MAPK1參與了細胞增殖分化、應激反應等活動，在細胞外信號轉導與細胞反應中發揮著關鍵作用〔28〕；MYC作為一種轉錄因子，可調控CyclinD1、CyclinE及P27等細胞周期相關基因的表達，與軟骨細胞的增殖、分化和凋亡密切相關〔29〕。本研究表明通過調節MAPK1、FOXO3及MYC水平減少細胞凋亡，并通過調節STAT3及GRB2表達，降低炎癥因子等物質的水平，從而改善OA的炎癥癥狀及軟骨退變情況。這些特異性指標蛋白的篩選可為進一步研究OA的臨床或動物實驗提供可靠的檢測指標，不僅能有效縮短實驗周期，還可降低實驗成本。但闡明LncRNA-miRNA-mRNA間的互作關系及LncRNA在OA中的作用機制是未來研究的重點。