林適 林燕平 黃佳純 袁嘉堯 連曉航 林賢燦 萬雷 黃宏興*

1. 廣州中醫(yī)藥大學(xué)第三臨床醫(yī)學(xué)院，廣東 廣州 510000

2. 廣州中醫(yī)藥大學(xué)第三附屬醫(yī)院，廣東 廣州 510240

骨質(zhì)疏松癥(osteoporosis, OP)是一種表現(xiàn)為骨強(qiáng)度下降的全身性骨骼疾病，與代謝相關(guān)，其特點是骨量降低、骨組織微結(jié)構(gòu)退化，進(jìn)而導(dǎo)致骨脆性增加和骨折風(fēng)險升高[1-2]。骨質(zhì)疏松癥是老年人骨折最常見的原因[3]。骨質(zhì)疏松癥和骨質(zhì)疏松性骨折在世界范圍內(nèi)構(gòu)成巨大的醫(yī)療、公共衛(wèi)生和經(jīng)濟(jì)負(fù)擔(dān)[4]。由于人口老齡化的加劇和過去二三十年來我國生活方式發(fā)生的巨大變化，預(yù)計在不久的將來，骨質(zhì)疏松癥患者數(shù)量會激烈增長[5]，我國公共衛(wèi)生將面臨嚴(yán)峻挑戰(zhàn)。高通量技術(shù)的最新進(jìn)展使得對許多生物化學(xué)骨生物標(biāo)志物的識別成為可能，有助于對骨質(zhì)疏松癥患者的診斷和管理[6]。例如血清PINP、CTX水平是臨床常用的骨形成和骨吸收標(biāo)志物[7]。然而，確定其他有用的生物標(biāo)志物或靶標(biāo)，以改善骨質(zhì)疏松癥的診斷和管理仍是當(dāng)務(wù)之急。

非編碼RNA包括微小RNA(miRNA)，占人類基因轉(zhuǎn)錄組的80 %～90 %，在各種生物過程中發(fā)揮著重要作用[8]。研究[9-10]發(fā)現(xiàn)，miRNA與骨質(zhì)疏松癥相關(guān)，其主要參與調(diào)控成骨細(xì)胞分化和破骨細(xì)胞分化等進(jìn)程。一系列miRNA被證實影響骨質(zhì)疏松癥的發(fā)病機(jī)制，如miR-415、miR-93、miR-125b等，有望成為骨質(zhì)疏松癥診斷及治療的新靶點[11]。

隨著生物信息學(xué)的發(fā)展，不同算法和研究策略應(yīng)用于識別基因網(wǎng)絡(luò)的潛在調(diào)控機(jī)制，使得許多疾病得到全面深入的了解。加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析 (weighted gene co-expression network analysis, WGCNA) 就是這樣一種算法，可以將生物學(xué)上高度相關(guān)的基因聚類到模塊中，描述模塊與臨床特征之間的相關(guān)性，已被廣泛應(yīng)用于疾病相關(guān)標(biāo)志物的篩選[12]。在本研究中，基于GEO數(shù)據(jù)庫中的基因表達(dá)數(shù)據(jù)，構(gòu)建骨質(zhì)疏松miRNA共表達(dá)模塊，將模塊與臨床性狀結(jié)合，得到骨質(zhì)疏松相關(guān)度最高的模塊并進(jìn)一步與獲得的差異miRNA進(jìn)行交集處理，然后獲取其靶標(biāo)mRNA，再進(jìn)一步進(jìn)行miRNA-mRNA調(diào)控網(wǎng)絡(luò)構(gòu)建，并對靶標(biāo)mRNA進(jìn)行基因功能及通路富集分析，以此為骨質(zhì)疏松癥的診斷及治療提供新的見解。

1 資料和方法

1.1 數(shù)據(jù)下載

通過GEO(gene expression omnibus)數(shù)據(jù)庫(https://www.ncbi.nlm.gov/gds)下載數(shù)據(jù)。miRNA數(shù)據(jù)：GSE74209數(shù)據(jù)集，mRNA數(shù)據(jù)：GSE7158數(shù)據(jù)集。GSE74209數(shù)據(jù)集包含6個OP樣本和6個正常樣本的miRNA表達(dá)數(shù)據(jù)；GSE7158數(shù)據(jù)集收集了14個高峰值骨量受試者和12個低峰值骨量受試者樣本的mRNA表達(dá)數(shù)據(jù)。

1.2 數(shù)據(jù)預(yù)處理、差異表達(dá)miRNA和差異表達(dá)mRNA篩選

通過Perl(www.perl.org/)處理下載得到的兩個數(shù)據(jù)集數(shù)據(jù)，并獲得mRNA和miRNA表達(dá)矩陣文件。利用R軟件中的limma(linear models for microarray data)工具包進(jìn)行差異表達(dá)miRNA(DE-miRNAs)和差異表達(dá)mRNA(DEGs)篩選。DEGs篩選標(biāo)準(zhǔn)為：|logFC|>1并且P-value<0.05；DE-miRNAs篩選標(biāo)準(zhǔn)為：|logFC|>1并且 adjustedP-value <0.05。

1.3 構(gòu)建miRNA加權(quán)共表達(dá)網(wǎng)絡(luò)

利用R中的“WGCNA”工具包構(gòu)建miRNA無尺度加權(quán)基因共表達(dá)網(wǎng)絡(luò)，鑒定與臨床性狀相關(guān)的共表達(dá)基因和模塊[13]。為獲得無尺度網(wǎng)絡(luò)，對數(shù)據(jù)進(jìn)行聚類以檢測異常值，并設(shè)置合適軟閾值，使網(wǎng)絡(luò)近乎接近無尺度網(wǎng)絡(luò)。采用層次聚類樹對基因模塊進(jìn)行識別，采用基于拓?fù)渲丿B矩陣的層次聚類對基因模塊進(jìn)行檢測。計算Pearson相關(guān)系數(shù)確定各模塊與臨床特征的相關(guān)性，獲取與臨床相關(guān)性最大的模塊并導(dǎo)出模塊內(nèi)miRNA。

1.4 獲取miRNA的mRNA靶點

將與臨床相關(guān)性最大模塊內(nèi)的miRNA和DE-miRNAs取交集得到OP相關(guān)的關(guān)鍵miRNA，并通過miRNet(https:www.mirnet.ca/miRNet/upload/MirUploadView.xhtml)預(yù)測關(guān)鍵miRNA下游靶基因。MiRNet是一個集成11個數(shù)據(jù)庫的miRNA-靶標(biāo)相互作用數(shù)據(jù)的在線工具，可用于預(yù)測miRNA的下游靶基因[14]。為增加預(yù)測結(jié)果的可靠性，本研究將預(yù)測靶基因和差異mRNA(DGEs)進(jìn)行取交集，作為關(guān)鍵miRNA在OP中發(fā)揮作用的目標(biāo)靶基因。

1.5 關(guān)鍵miRNA上游轉(zhuǎn)錄因子分析

通過FunRich(version 3.1.3)進(jìn)行關(guān)鍵miRNA上游轉(zhuǎn)錄因子分析。FunRich(http://www.funrich.org)是一個功能性的富集和相互作用網(wǎng)絡(luò)分析工具，可以用來分析基因集富集的轉(zhuǎn)錄因子[15]。

1.6 GO功能注釋和KEGG途徑富集分析

基于WebGestalt(WEB-based GEne SeT AnaLysis，http://www.webgestalt.org/)在線工具[16]進(jìn)行關(guān)鍵miRNA目標(biāo)靶基因GO功能注釋和KEGG通路富集分析，富集條目篩選標(biāo)準(zhǔn)為P-value<0.05。

1.7 構(gòu)建并分析miRNA-mRNA調(diào)控網(wǎng)絡(luò)

先從String在線數(shù)據(jù)庫(https://www.string-db.org/)中獲得關(guān)鍵miRNA目標(biāo)靶基因蛋白互作(PPI)關(guān)系文件，依次將關(guān)鍵miRNA和下游靶基因?qū)?yīng)關(guān)系文件、PPI互作關(guān)系文件并導(dǎo)入Cytoscape(version 3.7.1)中，利用Merge功能合并二者，構(gòu)建miRNA-mRNA調(diào)控網(wǎng)絡(luò)。最后，通過Cytoscape中MCODE方法提取核心調(diào)控網(wǎng)絡(luò)。

2 結(jié)果

2.1 DE-miRNAs和DEGs篩選

通過GSE74209數(shù)據(jù)集篩選得到142個差異miRNA,其中上調(diào)82個，下調(diào)60個；GSE7158數(shù)據(jù)集中篩選得到1 111個差異mRNA，其中上調(diào)625個，下調(diào)486個(基因熱圖及火山圖見圖1，綠色代表下調(diào)，紅色代表上調(diào)，黑色代表中間表達(dá)，熱圖中顏色越深顯著性越高)。

2.2 構(gòu)建WGCNA網(wǎng)絡(luò)及臨床相關(guān)性模塊獲取

利用WGCNA R工具包中的pickSoftThreshold函數(shù)計算軟閾值[17]。本研究中軟閾值設(shè)為6，尺度獨立性為0.9，平均連通性相對較高，以此構(gòu)建網(wǎng)絡(luò)更接近無尺度網(wǎng)絡(luò)(圖2A)?；诖藰?gòu)建WGCNA網(wǎng)絡(luò)分層聚類樹和共表達(dá)模塊(圖2B)。最終得到五個基因模塊，其中藍(lán)色模塊(90個miRNA)和骨質(zhì)疏松癥相關(guān)性最高為0.75(圖3A)，模塊內(nèi)基因顯著性較高(圖3B)。

2.3 關(guān)鍵miRNA獲取及下游靶基因預(yù)測

將藍(lán)色模塊內(nèi)的miRNA和DE-miRNA進(jìn)行取交集得到關(guān)鍵miRNA 6個，其中hsa-miR-3609、hsa-miR-664a-5p、hsa-miR-3182、hsa-miR-4297為下調(diào)miRNA，hsa-miR-3679-3p和hsa-miR-4500為上調(diào)miRNA(韋恩圖見圖4B)，通過miRNet預(yù)測得到hsa-miR-3609、hsa-miR-664a-5p、hsa-miR-3182、hsa-miR-4297、hsa-miR-3679-3p和hsa-miR-4500下游靶基因數(shù)分別為320、133、84、85、148和312(圖4A)。將所有預(yù)測靶基因和DEGs取交集得到52個關(guān)鍵miRNA在OP中發(fā)揮作用的目標(biāo)靶基因(韋恩圖見圖4C)。

2.4 關(guān)鍵miRNA上游轉(zhuǎn)錄因子預(yù)測

采用FunRich軟件預(yù)測關(guān)鍵miRNA上游轉(zhuǎn)錄因子，獲得差異較顯著的前5個轉(zhuǎn)錄因子(圖5)，分別為SP4、LHX3、NFIC、MYC和VSX2。

2.5 GO功能注釋和KEGG途徑富集分析

關(guān)鍵miRNA下游靶基因GO功能富集分析主要集中在對熱的反應(yīng)、DNA完整性檢查點、甘露糖基轉(zhuǎn)移酶活性等，通路富集分析主要富集在趨化因子信號通路、P53信號通路及炎癥、腫瘤等信號通路。見表1、表2。

表1 GO功能注釋Table 1 GO enrichment analysis

表2 KEGG通路富集分析Table 2 KEGG pathway enrichment analysis

2.6 構(gòu)建miRNA-mRNA調(diào)控網(wǎng)絡(luò)及核心網(wǎng)絡(luò)分析

通過Cytoscape構(gòu)建OP相關(guān)miRNA-mRNA調(diào)控網(wǎng)絡(luò)(圖6A)，MCODE方法提取核心調(diào)控網(wǎng)絡(luò)(圖6B)，發(fā)現(xiàn)hsa-mir-4500、CDKN1A、MAP2K7組成OP相關(guān)的核心調(diào)控網(wǎng)絡(luò)，其中hsa-mir-4500是上調(diào)差異miRNA, CDKN1A和MAP2K7均是下調(diào)的差異mRNA。

3 討論

目前治療骨質(zhì)疏松癥主要包括藥物防治、飲食管理和運(yùn)動療法，但是仍未有特效方法，這提示發(fā)掘骨質(zhì)疏松癥病因和分子機(jī)制對其防治極其重要。miRNA通過轉(zhuǎn)錄抑制和mRNA的失穩(wěn)來控制基因表達(dá)，與組織形態(tài)的發(fā)生、細(xì)胞凋亡等進(jìn)程和相關(guān)信號通路密切相關(guān)[18]，但是許多miRNA和mRNA的相互作用機(jī)制仍不清楚[19]。大量實驗研究[20-21]表明，miRNA與骨質(zhì)疏松癥的發(fā)生有關(guān)，主要是調(diào)節(jié)骨構(gòu)建、骨重吸收和成骨細(xì)胞分化之間的平衡。研究[22]證實miRNA表達(dá)差異與患者骨密度密切相關(guān)，miR-195、 miR-150是骨密度降低的高危因素，導(dǎo)致骨質(zhì)疏松癥的發(fā)生。本研究利用WGCNA分析和差異基因分析，篩選具有臨床相關(guān)性高同時差異表達(dá)的miRNA，充分結(jié)合兩種方法的優(yōu)勢。既往研究[23]中也有通過WGCNA方法篩選骨質(zhì)疏松癥臨床相關(guān)性高的lncRNA，但未結(jié)合差異表達(dá)分析。

本研究對GEO數(shù)據(jù)庫中的miRNA數(shù)據(jù)集進(jìn)行差異分析，并通過WGCNA分析獲得臨床相關(guān)性最高的miRNA模塊，最終篩選得到與OP發(fā)病相關(guān)性高且在OP中差異表達(dá)的6個關(guān)鍵miRNA，包括2個上調(diào)和4個下調(diào)差異表達(dá)miRNA，張等[24]發(fā)現(xiàn)過表達(dá)miR-664a-5p可促進(jìn)骨髓間充質(zhì)干細(xì)胞的成骨分化，而抑制miR-664a-5p則抑制骨髓間充質(zhì)干細(xì)胞的成骨分化。還未發(fā)現(xiàn)相關(guān)研究能夠證實其它幾個關(guān)鍵miRNA與OP有關(guān)，但是有相關(guān)研究證實其與腫瘤細(xì)胞遷移、增殖和分化相關(guān)，這些miRNA是需要關(guān)注的潛在研究方向。

探索miRNA上游轉(zhuǎn)錄因子有助于理解其作用機(jī)制，研究通過FunRich分析發(fā)現(xiàn)SP4、LHX3、NFIC、MYC和VSX2五個相關(guān)性最高的轉(zhuǎn)錄因子。特異性蛋白4(SP4)是具有鋅指結(jié)構(gòu)的轉(zhuǎn)錄因子，可以和具有G-C/T 序列結(jié)合，調(diào)控大量基因的表達(dá)，涉及多種生物學(xué)功能[25]。LHX3與神經(jīng)內(nèi)分泌密切相關(guān)，調(diào)控多種垂體激素，在生長發(fā)育中發(fā)揮重要作用[26]。NFIC屬于轉(zhuǎn)錄復(fù)制因子核因子Ⅰ家族中的一員[27]，可促進(jìn)大鼠牙濾泡細(xì)胞活力和成骨分化[28]。MYC作為上游信號，可以調(diào)控RANKL誘導(dǎo)的NFATc1，從而影響RANKL/RANK信號通路，敲減其表達(dá)可導(dǎo)致小鼠破骨細(xì)胞減少，而不影響成骨細(xì)胞數(shù)量，從而增加小鼠骨量[29]。

通過對關(guān)鍵miRNA進(jìn)行下游靶基因預(yù)測，并和GEO數(shù)據(jù)庫中差異表達(dá)mRNA整合得到52個關(guān)鍵miRNA目標(biāo)靶基因。GO和KEGG富集分析顯示這些基因涉及多個生物學(xué)進(jìn)程及通路，值得注意的顯著富集信號通路是趨化因子信號通路、P53信號通路。趨化因子通過控制炎癥過程中免疫細(xì)胞的遷移、定位和功能，在適應(yīng)性免疫系統(tǒng)中起著關(guān)鍵作用，目前已明確的趨化因子有27種[30]。由骨骼肌細(xì)胞或骨髓生態(tài)位的其他細(xì)胞釋放的趨化因子信號可以通過自分泌和旁分泌機(jī)制調(diào)節(jié)骨的形成和吸收[30]。趨化因子CCL3、CCL2均可以通過刺激RANKL的表達(dá)來促進(jìn)小鼠破骨細(xì)胞活性以及促進(jìn)骨吸收進(jìn)程[31-32]。趨化因子CX3CL1及其受體CX3CR1參與促進(jìn)破骨細(xì)胞生成，促進(jìn)骨吸收和BMD降低[33]。p53在骨質(zhì)疏松癥患者血清中表達(dá)豐富，p53下調(diào)可部分逆轉(zhuǎn)骨密度受損的結(jié)果[34]。Yu等[35]通過體外實驗證實p53可以抑制成骨細(xì)胞向成骨分化，而且白藜蘆醇可以逆轉(zhuǎn)p53誘導(dǎo)的成骨分化抑制。這些途徑有可能成為防治骨質(zhì)疏松癥的有效作用機(jī)制。

雖然本研究存在一些不足，如納入研究的數(shù)據(jù)集樣本量有限，可能導(dǎo)致篩選結(jié)果及生物學(xué)功能分析出現(xiàn)偏差，網(wǎng)絡(luò)構(gòu)建結(jié)果仍需進(jìn)一步實驗驗證等。但是本研究發(fā)揮了生物信息學(xué)的優(yōu)勢，通過mRNA差異表達(dá)譜及WGCNA分析，結(jié)合PPI互作關(guān)系，構(gòu)建了骨質(zhì)疏松癥miRNA-mRNA的潛在調(diào)控網(wǎng)絡(luò)，更加系統(tǒng)、全面地闡釋了骨質(zhì)疏松癥的分子機(jī)制，揭示了miRNA-mRNA在骨質(zhì)疏松癥中的表達(dá)是通過多途徑、多靶點進(jìn)行調(diào)控的，并且進(jìn)一步得到hsa-mir-4500、CDKN1A、MAP2K7組成的核心調(diào)控網(wǎng)絡(luò)，這些都將是未來的研究方向。同時，基于測序數(shù)據(jù)篩選miRNA，并對miRNA預(yù)測靶點和測序數(shù)據(jù)整合，具有一定的可靠性?？偠灾狙芯繛楣琴|(zhì)疏松癥的分子機(jī)制分析提供了新的思路與潛在靶點，對闡明其發(fā)病機(jī)制及防治藥物的開發(fā)具有重要的指導(dǎo)意義。