基于生物信息學(xué)挖掘骨關(guān)節(jié)炎潛在的關(guān)鍵生物標(biāo)志物*

應(yīng) 璞，路 通，許 岳，周 燁，戈昱凡，薛 燚，繆逸鳴

南京中醫(yī)藥大學(xué)常熟附屬醫(yī)院骨科,江蘇常熟 215500

骨關(guān)節(jié)炎(OA)是一種臨床上最常見的以關(guān)節(jié)軟骨退變和破壞為特征的慢性關(guān)節(jié)疾病，其臨床表現(xiàn)為關(guān)節(jié)疼痛、僵硬、腫脹伴關(guān)節(jié)活動受限等[1-2]。OA發(fā)病機制復(fù)雜，涉及關(guān)節(jié)軟骨的機械刺激及退變、滑膜組織的病變、軟骨下骨微循環(huán)改變、炎癥反應(yīng)等[3-5]，但目前OA的具體機制仍未闡明。近年來，隨著微陣列技術(shù)和生物信息學(xué)的發(fā)展，通過對OA的高通量微陣列和基因芯片數(shù)據(jù)進行研究，越來越多的OA相關(guān)差異表達基因(DEGs)被發(fā)現(xiàn)[6-9]。這些基因的發(fā)現(xiàn)，有望為闡明OA發(fā)病的分子機制和治療靶點提供新的理論指導(dǎo)。本研究從基因表達綜合數(shù)據(jù)庫(GEO)下載了mRNA數(shù)據(jù)集和微小RNA(miRNA)數(shù)據(jù)集，選取OA滑膜組和正常滑膜組進行比較與分析，分析OA中差異表達的mRNA和miRNA。本研究旨在通過生物信息學(xué)方法篩選出與OA發(fā)生、發(fā)展相關(guān)的關(guān)鍵mRNA和miRNA，建立mRNA-miRNA共表達網(wǎng)絡(luò)，從而進一步探索OA的發(fā)病機制和潛在生物標(biāo)志物。

1 材料與方法

1.1材料 本研究采用GE數(shù)據(jù)庫(https://www.ncbi.nlm.nih.gov/geo/)下載mRNA數(shù)據(jù)集GSE55457、GSE82107、GSE55235和miRNA數(shù)據(jù)集GSE143514。其中GSE55457數(shù)據(jù)集包括10例OA滑膜組織和10例正常滑膜組織，GSE82107數(shù)據(jù)集包含10例OA滑膜組織和7 例正常滑膜組織，GSE55235數(shù)據(jù)集包含26例OA滑膜組織和20例正常滑膜組織，GSE143514數(shù)據(jù)集包含5例OA滑膜組織和3例正常滑膜組織。

1.2方法

1.2.1差異基因篩選 所有芯片數(shù)據(jù)由R語言軟件的affy程序包預(yù)處理，隨后通過limma程序進行篩選，OA與正常滑膜之間的差異表達基因(DEGs)，按log2FC>1且P<0.05的標(biāo)準(zhǔn)進行篩選(FC為差異倍數(shù))。重疊差異表達基因數(shù)據(jù)集采用R語言軟件的VennDiagram包Venn圖顯示。

1.2.2差異基因GO和KEGG富集分析 利用在線分析工具Metascape(http://metascape.org/)對差異表達基因進行GO功能富集分析和KEGG通路富集分析，以P<0.05且基因計數(shù)≥5為差異有統(tǒng)計學(xué)意義。

1.2.3蛋白質(zhì)-蛋白質(zhì)相互作用(PPI)網(wǎng)絡(luò)分析及miRNA-mRNA共表達網(wǎng)絡(luò)分析 利用在線分析數(shù)據(jù)庫STRING ( http://string-db.org/)將篩選出的差異表達基因?qū)霐?shù)據(jù)庫，進行差異基因PPI網(wǎng)絡(luò)分析,并利用Cytoscape軟件(https://cytoscape.org)構(gòu)建PPI網(wǎng)絡(luò)。

2 結(jié) 果

2.1篩選差異表達的mRNA和miRNA 對mRNA數(shù)據(jù)集GSE55457、GSE82107和miRNA數(shù)據(jù)集GSE143514原始數(shù)據(jù)進行處理后，按篩選標(biāo)準(zhǔn)進行分析，最后從GSE55457中篩選出626個差異表達mRNA，其中313個上調(diào)、313個下調(diào)。按同樣的篩選標(biāo)準(zhǔn)從GSE82107中篩選出1 647個差異表達mRNA，其中498個下調(diào)、1 147個上調(diào)，從GSE143514中篩選出142個差異表達miRNA，其中22個下調(diào)、120個上調(diào),見圖1。隨后，利用Venn圖比較GSE55457和GSE82107兩個數(shù)據(jù)集的差異表達mRNA，得到99個相同的mRNA，見圖2。

注：A為差異表達mRNA數(shù)據(jù)集GSE55457；B為差異表達mRNA數(shù)據(jù)集GSE82107；C為差異表達miRNA數(shù)據(jù)集GSE143514。

2.2鑒定相關(guān)的潛在mRNA 首先對數(shù)據(jù)集GSE55235的數(shù)據(jù)進行標(biāo)準(zhǔn)化處理，然后通過limma軟件包從GSE55235中篩選出334個差異表達的mRNA，見圖3。為了使網(wǎng)絡(luò)符合無尺度網(wǎng)絡(luò)分布,利用權(quán)重基因共表達網(wǎng)絡(luò)分析(WGCNA)軟件包中的函數(shù)pickSoftThreshold計算權(quán)重值，并選擇軟閾值β=16用來構(gòu)建共表達網(wǎng)絡(luò)，見圖4A，并采用平均連鎖層次聚類方法對基因進行聚類，見圖4B。同時，對模塊進行聚類分析，然后根據(jù)各模塊特征基因與樣本類型的相關(guān)性，發(fā)現(xiàn)綠色模塊與樣本類型的相關(guān)性最大，見圖4C。這些結(jié)果提示綠色模塊中的27個基因可能與OA的發(fā)生密切相關(guān)。此外，比較綠色模塊中篩選得到的27個mRNA與圖2中篩選得到的99個mRNA，最終獲得OA的19個交叉的關(guān)鍵mRNAs和27個綠色模塊mRNAs，見圖4D。19個交叉的關(guān)鍵mRNA為APOLD1、CD52、趨化因子配體(CXCL)13、CXCL2、CXCL3、DUSP1、FKBP5、FOSB、IGJ、IGLL5、JUN、基質(zhì)金屬蛋白酶(MMP)1、ZB1、NFIL3、NKG7、NR4A2、磷脂酶Cδ3(PLCD3)、POU2AF1、TNFRSF17。和27個綠色模塊mRNAs為BCL6、RHOB、NFIL3、KLF4、PLCD3、CEBPD、GADD45A、BTG1、GADD45B、EIF1、APOLD1、CD52、NKG7、CXCL13、CXCL3、CXCL2、DUSP1、MMP1、JUN、FKBP5、NR4A2、FOSB、IGJ、TNFRSF17、POU2AF1、MZB1、IGLL5。

圖2 GSE55457和GSE82107共同差異表達mRNA的Venn圖

2.3mRNA的GO分析和KEGG分析 GO功能富集分析顯示，19個關(guān)鍵mRNAs主要參與免疫相關(guān)過程。例如體液免疫反應(yīng)、細(xì)胞對趨化因子的反應(yīng)、白細(xì)胞趨化、粒細(xì)胞趨化、粒細(xì)胞遷移、細(xì)胞對脂多糖的反應(yīng)、B細(xì)胞活化和白細(xì)胞遷移、抗菌體液反應(yīng)以及對細(xì)菌的防御反應(yīng)等，見圖5A。KEGG通路富集分析顯示，19個關(guān)鍵mRNAs主要參與白細(xì)胞介素(IL)-17信號通路、腫瘤壞死因子(TNF)信號通路、轉(zhuǎn)化生長因子-β(TGF-β)信號通路、NOD樣受體信號通路、趨化因子信號通路、白細(xì)胞介素信號通路和PI3K/AKT信號通路等，見圖5B。

圖3 OA中差異表達mRNA數(shù)據(jù)集GSE55235火山圖

注：A為各種軟閾值功率的網(wǎng)絡(luò)拓?fù)浞治觯籅為基因模塊的聚類樹狀圖；C為OA和正常樣本之間基因模塊的相關(guān)性分析；D為圖2中99個篩選得到的mRNA和27個綠色模塊篩選得到的共同差異表達mRNA的Venn圖。

2.4PPI網(wǎng)絡(luò)分析及miRNA-mRNA共表達網(wǎng)絡(luò)分析 利用STRING數(shù)據(jù)庫和Cytoscape軟件構(gòu)建一個由11個節(jié)點和24條邊組成的PPI網(wǎng)絡(luò)，見圖6A。隨后通過MCODE插件獲得PPI評分最高的子網(wǎng)絡(luò)，該子網(wǎng)絡(luò)由CXCL2、JUN、CXCL3、MMP1、PLCD3共5個mRNA組成，見圖6B，說明CXCL2、JUN、CXCL3、MMP1、PLCD3這5個基因可能是OA最關(guān)鍵的基因。同時利用DIANA Tools在線數(shù)據(jù)庫(http://diana.imis.athena-innovation.gr/DianaTools/index.php)分析預(yù)測5個關(guān)鍵基因的靶向miRNA。然后，將獲得的所有靶向miRNA與從GSE143514中篩選出142個差異表達miRNA，見圖1C，進行比較，獲得5個關(guān)鍵基因的24個共同miRNA，見圖6C。此外，還預(yù)測了39個miRNA-mRNA基因?qū)Α８鶕?jù)預(yù)測結(jié)果，利用Cytoscape軟件構(gòu)建了一個由29個節(jié)點和39條邊組成的mRNA-miRNA共表達網(wǎng)絡(luò)，見圖6D、表1。

注：A為GO功能富集分析；B為KEGG功能富集分析。

表1 OA的mRNA-miRNA共表達網(wǎng)絡(luò)中的miRNAs和mRNAs

注：A為PPI網(wǎng)絡(luò)；B為PPI網(wǎng)絡(luò)5個mRNA組成；C為靶向miRNA的差異表達篩選；D為mRNA-miRNA共表達網(wǎng)絡(luò)。

3 討 論

OA發(fā)病機制復(fù)雜，目前已有多個基因和miRNA發(fā)現(xiàn)參與OA的發(fā)生、發(fā)展，但其發(fā)病機制仍未被闡明[10-14]。基于高通量測序基因芯片技術(shù)和生物信息學(xué)技術(shù)的發(fā)展，能夠篩選出OA的潛在的關(guān)鍵基因和關(guān)鍵miRNA，從而深入研究其分子機制并提供潛在的治療靶點。

本研究整合了4個GEO數(shù)據(jù)集，并通過篩選OA和健康對照組的滑膜組織中差異表達的mRNA和miRNA，確定了19個關(guān)鍵mRNA。通過GO功能富集分析和KEGG通路分析，這些關(guān)鍵mRNA主要參與免疫相關(guān)過程，例如體液免疫反應(yīng)、細(xì)胞對趨化因子的反應(yīng)、白細(xì)胞趨化、粒細(xì)胞趨化等，主要參與IL-17信號通路、TNF信號通路、TGF-β信號通路等。HAN等[15]通過研究5個GEO數(shù)據(jù)庫獲得了5個OA的關(guān)鍵差異表達基因(TLR7、CSF1R、APOE、C1QA、CCL5)，這些基因廣泛參與炎癥反應(yīng)調(diào)控、血管形態(tài)發(fā)生、內(nèi)皮細(xì)胞遷移、體液免疫應(yīng)答等生物學(xué)過程。MEEHAN等[16]證實OA患者中關(guān)節(jié)滑液中的IL-4、IL-6、IL-8、IL-15水平顯著高于血清中的水平，證實OA患者中的滑液中存在大量趨化因子和促炎細(xì)胞因子。目前的研究認(rèn)為，與OA相關(guān)的促炎因子，包括IL-1β、TNF-α、IL-6、IL-15、IL-17和IL-18，以及減少OA相關(guān)的關(guān)鍵抗炎因子，包括IL-4、胰島素樣生長因子和TGF-β等[3]。這些研究的結(jié)果與本研究結(jié)果是一致的，證實了炎癥和相關(guān)免疫反應(yīng)在OA的發(fā)生、發(fā)展中有著重要的作用。

本研究通過構(gòu)建PPI網(wǎng)絡(luò)分析獲得的5個關(guān)鍵基因分別為CXCL2、CXCL3、JUN、MMP1、PLCD3。其中CXCL2和CXCL3都是趨化因子CXC亞家族成員，均為炎性趨化因子，對多形核白細(xì)胞具有趨化作用。STONE等[17]發(fā)現(xiàn)，OA半月板細(xì)胞中趨化因子(CXCL1、CXCL2)的表達顯著增加，半月板細(xì)胞的促炎刺激增加了基質(zhì)金屬蛋白酶的產(chǎn)生和分解代謝基因的表達，在關(guān)節(jié)損傷后OA的發(fā)展中發(fā)揮重要作用。LYNCH等[18]評估了43種生物標(biāo)志物，結(jié)果發(fā)現(xiàn)OA軟骨組織中趨化因子CXCL3顯著高于正常軟骨組織。于雅麗等[19]發(fā)現(xiàn)JUN基因被確定為有價值的OA生物標(biāo)記物，這與本研究結(jié)果一致。c-JUN蛋白是一種轉(zhuǎn)錄因子，各種物理、化學(xué)、生物學(xué)刺激及細(xì)胞應(yīng)激反應(yīng)都可以通過促進c-JUN蛋白表達活化，進而對細(xì)胞增殖、凋亡等生物學(xué)過程進行調(diào)控。XU等[20]發(fā)現(xiàn)膝OA患者的軟骨下骨細(xì)胞中ERK/c-JUN N-末端激酶(JNK)信號通路被激活，并參與調(diào)控成骨，促進膝OA的發(fā)生、發(fā)展。目前已有多項研究表明軟骨組織或滑液中MMP-1的高水平可能與膝關(guān)節(jié)OA的發(fā)病密切相關(guān)[16,21-22]。PLCD3其被發(fā)現(xiàn)能夠促進腫瘤細(xì)胞的增殖、遷移和侵襲[23]，然而在OA的發(fā)病機制中的作用較少報道，需要后續(xù)更多研究探討PLCD3在OA致病機制中的作用。

本研究共發(fā)現(xiàn)39對mRNA-miRNA共表達網(wǎng)絡(luò)，其中miR-34a-5p在5組mRNA-miRNA均出現(xiàn)，提示其可能在OA的發(fā)病機制中存在重要作用。SONG等[24]指出，下調(diào)miR-34a-5p可促進OA細(xì)胞增殖，抑制細(xì)胞凋亡和自噬；ENDISHA等[25]也指出，miR-34a-5p可以促進OA期間的關(guān)節(jié)破壞。這些研究表明，miR-34a-5p是OA細(xì)胞表達的潛在調(diào)控因子[24-26]。另外，miR-20a-5p在4組mRNA-miRNA均出現(xiàn)，提示其可能在OA的致病機制中也存在重要作用[27]。ROSS等[28]通過整合mRNA和miRNA測序數(shù)據(jù)，創(chuàng)建了miRNA-mRNA相互作用網(wǎng)絡(luò)，發(fā)現(xiàn)miR-20a-5p可以顯著降低降解酶的產(chǎn)生以及炎癥相關(guān)基因的表達。這些miRNA的發(fā)現(xiàn)，為研究OA的分子機制提供了方向，也為OA的治療提供了潛在的靶點。

綜上所述，本研究通過生物信息學(xué)的方法發(fā)現(xiàn)了OA與炎癥及相關(guān)免疫反應(yīng)存在密切聯(lián)系。同時，也發(fā)現(xiàn)了5個可能參與OA發(fā)病機制的基因：CXCL2、CXCL3、JUN、MMP1、PLCD3。另外，還發(fā)現(xiàn)了一些OA相關(guān)的mRNA-miRNA共表達網(wǎng)絡(luò)，并提出miR-34a-5p、miR-20a-5p可能會成為未來研究OA分子機制和治療OA的潛在靶點。