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

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

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

1 材料與方法

1.1材料 本研究采用GE數(shù)據(jù)庫(kù)(https://www.ncbi.nlm.nih.gov/geo/)下載mRNA數(shù)據(jù)集GSE55457、GSE82107、GSE55235和miRNA數(shù)據(jù)集GSE143514。其中GSE55457數(shù)據(jù)集包括10例OA滑膜組織和10例正?；そM織，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語(yǔ)言軟件的affy程序包預(yù)處理，隨后通過(guò)limma程序進(jìn)行篩選，OA與正?；ぶg的差異表達(dá)基因(DEGs)，按log2FC>1且P<0.05的標(biāo)準(zhǔn)進(jìn)行篩選(FC為差異倍數(shù))。重疊差異表達(dá)基因數(shù)據(jù)集采用R語(yǔ)言軟件的VennDiagram包Venn圖顯示。

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

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

2 結(jié) 果

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

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

2.2鑒定相關(guān)的潛在mRNA 首先對(duì)數(shù)據(jù)集GSE55235的數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理，然后通過(guò)limma軟件包從GSE55235中篩選出334個(gè)差異表達(dá)的mRNA，見(jiàn)圖3。為了使網(wǎng)絡(luò)符合無(wú)尺度網(wǎng)絡(luò)分布,利用權(quán)重基因共表達(dá)網(wǎng)絡(luò)分析(WGCNA)軟件包中的函數(shù)pickSoftThreshold計(jì)算權(quán)重值，并選擇軟閾值β=16用來(lái)構(gòu)建共表達(dá)網(wǎng)絡(luò)，見(jiàn)圖4A，并采用平均連鎖層次聚類(lèi)方法對(duì)基因進(jìn)行聚類(lèi)，見(jiàn)圖4B。同時(shí)，對(duì)模塊進(jìn)行聚類(lèi)分析，然后根據(jù)各模塊特征基因與樣本類(lèi)型的相關(guān)性，發(fā)現(xiàn)綠色模塊與樣本類(lèi)型的相關(guān)性最大，見(jiàn)圖4C。這些結(jié)果提示綠色模塊中的27個(gè)基因可能與OA的發(fā)生密切相關(guān)。此外，比較綠色模塊中篩選得到的27個(gè)mRNA與圖2中篩選得到的99個(gè)mRNA，最終獲得OA的19個(gè)交叉的關(guān)鍵mRNAs和27個(gè)綠色模塊mRNAs，見(jiàn)圖4D。19個(gè)交叉的關(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個(gè)綠色模塊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共同差異表達(dá)mRNA的Venn圖

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

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

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

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

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

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

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

3 討 論

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

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

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

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

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