潘俊曦 田天照 劉保新 張勝 馮慶輝 彭志華 蔡迎峰

廣州醫(yī)科大學(xué)附屬中醫(yī)醫(yī)院骨傷科，廣東 廣州 510000

絕經(jīng)后骨質(zhì)疏松癥(postmenopausal osteoporosis，PMOP)是常見的代謝性骨病，主要發(fā)生于絕經(jīng)后女性，由于絕經(jīng)后雌激素水平下降，骨吸收活動(dòng)增強(qiáng)，導(dǎo)致全身骨量下降、骨微結(jié)構(gòu)破壞、骨脆性增加，繼發(fā)骨質(zhì)疏松性骨折具有較高的病殘率、病死率，增加社會(huì)、家庭的經(jīng)濟(jì)及心理負(fù)擔(dān)，PMOP已成為我國重要的公共衛(wèi)生問題之一[1]。糖尿病是一種以高血糖為特征的全身代謝性疾病，可導(dǎo)致外周血管、神經(jīng)末梢等發(fā)生病變，引發(fā)多系統(tǒng)并發(fā)癥，主要包括1型與2型糖尿病[2-3]。中國是世界上2型糖尿病(type 2 diabetes mellitus，T2DM)患者最多的國家，2013年統(tǒng)計(jì)約有1.36億患者，占全球患病總數(shù)的1/3[4]。近年來多項(xiàng)研究表明T2DM患者有較低的骨密度值和較高的髖部及椎體骨折風(fēng)險(xiǎn)，特別是患有T2DM的絕經(jīng)后女性患者[5-6]。T2DM被認(rèn)為與PMOP發(fā)病具有密切聯(lián)系，但兩者間的病理機(jī)制聯(lián)系目前尚不明確，主要認(rèn)為是一些促進(jìn)骨形成因素(如胰島素、肥胖)、抑制骨形成因素(高血糖、微血管病變)等相互綜合，導(dǎo)致骨代謝平衡被打破[7]。微小核糖核苷酸(microRNA，miRNA)是一種單鏈內(nèi)源性非編碼RNA分子，通過與信使RNA的3′非翻譯區(qū)互補(bǔ)序列相結(jié)合參與調(diào)控轉(zhuǎn)錄后翻譯[8-9]。研究表明，miRNA參與調(diào)控PMOP與T2DM的病理過程，Xu等[10]發(fā)現(xiàn)miR-125a-5p 通過靶向STAT3 抑制肝臟糖異生減輕T2DM進(jìn)展；Suarjana等[11]發(fā)現(xiàn)miR-21參與抑制破骨細(xì)胞增殖，改善PMOP。深入挖掘miRNA調(diào)控疾病病理機(jī)制是目前的重要研究熱點(diǎn)，探索PMOP與T2DM之間的關(guān)聯(lián)miRNA及其調(diào)控作用將有助于深化對(duì)疾病聯(lián)系的認(rèn)知，發(fā)現(xiàn)新的治療靶點(diǎn)與方向。

機(jī)器學(xué)習(xí)是涉及概率學(xué)、統(tǒng)計(jì)學(xué)、近似理論、人工智能的多領(lǐng)域交叉學(xué)科，通過計(jì)算機(jī)實(shí)現(xiàn)實(shí)時(shí)模擬人類學(xué)習(xí)過程[12-13]。隨機(jī)森林是一種經(jīng)典機(jī)器學(xué)習(xí)算法，通過建立決策樹分類器模型，對(duì)分類變量進(jìn)行反復(fù)迭代評(píng)分，產(chǎn)生高精確性分類器，幫助發(fā)現(xiàn)分類變量中的重要聚類及獨(dú)立因子，已成為生物醫(yī)學(xué)領(lǐng)域挖掘生物標(biāo)志物的重要算法[14]。本研究借助隨機(jī)森林算法，通過挖掘PMOP與T2DM患者miRNA芯片數(shù)據(jù)，篩選聯(lián)系兩疾病的關(guān)聯(lián)miRNA，并對(duì)其調(diào)控靶向基因的分子機(jī)制進(jìn)行探索，以期為防治PMOP與T2DM提供新的思路方向。

1 材料與方法

1.1 芯片數(shù)據(jù)獲取與預(yù)處理

通過基因表達(dá)數(shù)據(jù)庫(gene expression omnibus，GEO)檢索下載GSE70318基因芯片作為訓(xùn)練集數(shù)據(jù)，從中篩選PMOP與T2DM關(guān)聯(lián)miRNA。GSE70318芯片基于GPL20631平臺(tái)進(jìn)行檢測，包含57個(gè)絕經(jīng)后女性血清樣本miRNA測序數(shù)據(jù)，分別為19個(gè)T2DM患者、19個(gè)PMOP患者以及19個(gè)同時(shí)診斷為PMOP與T2DM的患者。下載GSE74209芯片作為測試集數(shù)據(jù)驗(yàn)證算法結(jié)果，該芯片包含12個(gè)PMOP患者測序數(shù)據(jù)。使用R語言preprocess Core軟件包中對(duì)GSE70318原始數(shù)據(jù)進(jìn)行分位數(shù)歸一化，得到標(biāo)準(zhǔn)化表達(dá)矩陣并進(jìn)行基因名重注釋，當(dāng)出現(xiàn)重復(fù)的基因名時(shí)對(duì)其進(jìn)行合并取均值處理。

1.2 隨機(jī)森林算法篩選PMOP與T2DM關(guān)聯(lián)miRNA

使用Randomforest軟件包構(gòu)建隨機(jī)森林模型，通過隨機(jī)生成大量的分類樹并對(duì)每棵樹的miRNA分類結(jié)果進(jìn)行迭代評(píng)分獲得分類結(jié)局，最終對(duì)所有單棵樹的分類結(jié)果進(jìn)行綜合判定。使用Caret軟件包對(duì)模型中miRNA重要性進(jìn)行排序，選取排名前10的miRNA，并根據(jù)miRNA的表達(dá)量繪制差異miRNA散點(diǎn)圖。

1.3 評(píng)價(jià)納入miRNA的預(yù)測性能

接受者工作特征(receiver operator characteristic，ROC)曲線是用于評(píng)價(jià)分類器預(yù)測性能的經(jīng)典工具。本研究中以驗(yàn)證集數(shù)據(jù)對(duì)表達(dá)量顯著差異的miRNA進(jìn)行ROC曲線繪制，計(jì)算其曲線下面積(area under curve，AUC)，評(píng)價(jià)miRNA的預(yù)測性能。AUC值范圍在0.5～1，其值越大代表分類器預(yù)測性能越好，設(shè)定AUC>0.7為條件，認(rèn)為miRNA在模型中的分類結(jié)果準(zhǔn)確。

1.4 miRNA靶向基因預(yù)測篩選

使用Targetscan[15]、miRWalk2.0[16]和DIANA TOOLS[17]3個(gè)公共數(shù)據(jù)庫預(yù)測PMOP與T2DM關(guān)聯(lián)miRNA靶向基因，并對(duì)多數(shù)據(jù)庫預(yù)測結(jié)果取共交集，獲得miRNA靶向基因[18]。

1.5 構(gòu)建miRNA-靶向基因調(diào)控互作網(wǎng)絡(luò)

將靶向基因上傳至蛋白互作分析數(shù)據(jù)庫STRING[19]進(jìn)行蛋白互作分析，預(yù)測靶向基因的蛋白-蛋白互作(protein-protein interaction，PPI)關(guān)系，使用網(wǎng)絡(luò)構(gòu)建工具Gephi軟件構(gòu)建miRNA-靶向基因調(diào)控互作網(wǎng)絡(luò)。

1.6 基因本體論與信號(hào)通路富集分析

基因本體論(gene oncology，GO)[20]和京都基因與基因組百科全書(Kyoto Encyclopedia of Genes and Genomes，KEGG)[21]信號(hào)通路富集分析是分別從基因產(chǎn)物功能與分子信號(hào)通路水平探索基因參與生物學(xué)功能調(diào)控作用的重要方法。本研究通過將靶向基因上傳至DAVID6.8數(shù)據(jù)庫進(jìn)行GO與KEGG富集分析，并使用R語言對(duì)富集分析結(jié)果進(jìn)行可視化。

2 結(jié)果

2.1 隨機(jī)森林算法篩選miRNA

通過對(duì)GSE70318芯片原始數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化及基因注釋，共獲得153個(gè)miRNA，構(gòu)建隨進(jìn)森林算法模型，對(duì)分類樹分類結(jié)果評(píng)分，篩選PMOP與T2DM關(guān)聯(lián)重要性排名前10的miRNA，分別為hsa-miR-188-3p、hsa-miR-181a-3p、hsa-miR-155-5p、hsa-miR-135a-5p、hsa-miR-382-3p、hsa-miR-32-3p、hsa-miR-576-3p、hsa-miR-942、hsa-miR-330-3p和hsa-miR-369-3p，結(jié)果如圖1所示。將上述miRNA分為PMOP-T2DM關(guān)聯(lián)組(19例)與對(duì)照組(38例，單純PMOP或T2DM)進(jìn)行對(duì)比，繪制表達(dá)量差異散點(diǎn)圖(見圖2)。結(jié)果表明，hsa-miR-188-3p、hsa-miR-181a-3p、hsa-miR-135a-5p和hsa-miR-369-3p表達(dá)量在關(guān)聯(lián)組與對(duì)照組間具有顯著統(tǒng)計(jì)學(xué)差異(P<0.01)。

圖1 miRNA排名前10位Fig.1 Ranking of top 10 miRNAs

圖2 miRNA表達(dá)差異散點(diǎn)圖Fig.2 Scatter diagram of miRNAs注：Group1：對(duì)照組；Group2：關(guān)聯(lián)組。

2.2 ROC曲線驗(yàn)證miRNA重要性

使用標(biāo)準(zhǔn)化的GSE74209芯片作為驗(yàn)證集數(shù)據(jù)，對(duì)顯著差異的miRNA繪制ROC曲線進(jìn)行預(yù)測性能驗(yàn)證。分別對(duì)hsa-miR-188-3p、hsa-miR-181a-3p、hsa-miR-135a-5p、hsa-miR-369-3p進(jìn)行驗(yàn)證，同時(shí)將上述4個(gè)miRNA作為一個(gè)集合模型共同驗(yàn)證以增加驗(yàn)證準(zhǔn)確性。ROC曲線結(jié)果如圖3所示，hsa-miR-369-3p獲得最高AUC值(0.757)，其余miRNA的AUC值為：hsa-miR-181a-3p(0.586)、hsa-miR-188-3p(0.688)、hsa-miR-135a-5p(0.500)、miRNA集合模型(0.667)，均低于0.7。

hsa-miR-369-3p的AUC值高于其余miRNA，認(rèn)為hsa-miR-369-3p在PMOP與T2DM疾病關(guān)聯(lián)性中預(yù)測性能更好；同時(shí)，將上述4個(gè)miRNA構(gòu)建集合模型進(jìn)行ROC驗(yàn)證，其AUC為0.667，認(rèn)為該集合模型未能提高預(yù)測準(zhǔn)確性。據(jù)此，選擇hsa-miR-369-3p作為關(guān)聯(lián)PMOP-T2DM的關(guān)鍵miRNA。

圖3 ROC曲線圖Fig.3 ROC diagram

2.3 關(guān)聯(lián)miRNA靶向基因預(yù)測及共表達(dá)網(wǎng)絡(luò)構(gòu)建

圖4 韋恩圖與調(diào)控互作網(wǎng)絡(luò) A：靶向基因篩選共交集；B：關(guān)聯(lián)miRNA-靶向基因調(diào)控網(wǎng)絡(luò)。Fig.4 Venn diagram and network of miRNA targetsA: Screening of cointersection of the targeted genes; B: Regulatory network of associated miRNA-targeted gene.

將hsa-miR-369-3p分別上傳至Targetscan、miRWalk2.0和DIANA TOOLS數(shù)據(jù)庫，設(shè)置物種為“Homo sapiens”(人類)，進(jìn)行靶向基因預(yù)測。Targetscan預(yù)測獲得2 324個(gè)靶向基因，miRWalk2.0獲得151個(gè)靶向基因，DIANA TOOLS獲得89個(gè)靶向基因，對(duì)3個(gè)數(shù)據(jù)庫獲得的靶向基因取共交集后最終獲得44個(gè)目標(biāo)靶向基因(圖4A)。

將靶向基因?qū)隨TRING數(shù)據(jù)庫，設(shè)置關(guān)聯(lián)置信度為<0.40，物種為“Homo sapiens”，進(jìn)行蛋白互作分析，將分析結(jié)果導(dǎo)入Gephi軟件中，構(gòu)建miRNA-靶向基因調(diào)控互作網(wǎng)絡(luò)(圖4B)。對(duì)靶向基因在網(wǎng)絡(luò)中關(guān)聯(lián)信度進(jìn)行計(jì)算，最終獲得網(wǎng)絡(luò)核心基因?yàn)椋篊YR61、CALD1、DDTT4和DUSP1。

2.4 GO與KEGG富集分析

將靶向基因上傳至DAVID6.8數(shù)據(jù)庫，設(shè)置命名類型為“Official Symbol”，物種為“Homo sapiens”，進(jìn)行GO生物學(xué)功能分析和KEGG信號(hào)通路富集分析，以P<0.05為差異具有統(tǒng)計(jì)學(xué)意義。GO富集分析(圖5A)顯示，靶向基因主要富集于信號(hào)轉(zhuǎn)導(dǎo)、內(nèi)質(zhì)網(wǎng)膜的組成部分、細(xì)胞外基質(zhì)組織、肌細(xì)胞細(xì)胞穩(wěn)態(tài)、蛋白質(zhì)加工、活性氧代謝過程、蛋白K48連接的泛素化、運(yùn)動(dòng)行為、細(xì)胞鋅離子穩(wěn)態(tài)、細(xì)胞生長的正調(diào)控、RNA剪接，著絲粒和染色質(zhì)結(jié)合等GO生物學(xué)過程。KEGG信號(hào)通路富集(圖5B)主要集中于內(nèi)質(zhì)網(wǎng)中的蛋白質(zhì)加工、礦物質(zhì)吸收、血管平滑肌收縮、溶酶體、cAMP信號(hào)通路、Ras信號(hào)通路、PI3K-Akt信號(hào)傳導(dǎo)途徑和代謝途徑。

圖5 GO與KEGG富集分析結(jié)果A：GO富集分析；B：KEGG信號(hào)通路富集分析。Fig.5 Results of GO and KEGG enrichment analysisA: Go enrichment analysis; B: Enrichment analysis of KEGG signal pathway.

3 討論

miRNA通過調(diào)控轉(zhuǎn)錄后翻譯，參與調(diào)控細(xì)胞的增殖、分化、凋亡等生命活動(dòng)，被認(rèn)為是診斷、治療疾病的新靶點(diǎn)[22]。人類基因組轉(zhuǎn)錄翻譯活動(dòng)約60%受miRNA調(diào)控，但目前所能識(shí)別并確認(rèn)功能活動(dòng)的miRNA僅有一小部分，miRNA在許多疾病病理過程中的作用目前仍未明確[23]。機(jī)器學(xué)習(xí)算法可以從大量變量中精確識(shí)別關(guān)鍵變量，這些變量具有合理的預(yù)測價(jià)值及臨床意義[24]，目前包括隨機(jī)森林、支持向量機(jī)、神經(jīng)網(wǎng)絡(luò)等多種不同的機(jī)器學(xué)習(xí)算法被不斷發(fā)展并應(yīng)用于生物醫(yī)學(xué)領(lǐng)域[25]。

本研究使用隨機(jī)森林算法建立分類樹模型，以同時(shí)關(guān)聯(lián)PMOP與T2DM為數(shù)據(jù)特征，對(duì)miRNA進(jìn)行多次分類并對(duì)分類結(jié)局進(jìn)行評(píng)分，篩選出10個(gè)關(guān)聯(lián)性最高的miRNA。通過對(duì)這10個(gè)miRNA進(jìn)行組間表達(dá)量差異分析，進(jìn)一步確定了4個(gè)顯著差異的miRNA(hsa-miR-188-3p、hsa-miR-181a-3p、hsa-miR-135a-5p、hsa-miR-369-3p)。對(duì)上述4個(gè)miRNA進(jìn)行ROC驗(yàn)證，結(jié)果提示hsa-miR-369-3p的AUC值最高，預(yù)測性能最好，可能是聯(lián)系PMOP與T2DM的關(guān)鍵miRNA。

研究發(fā)現(xiàn)，hsa-miR-369-3p通過參與細(xì)胞的增殖、分化、遷移等生命活動(dòng)，參與調(diào)控人體生理病理活動(dòng)[26-27]；hsa-miR-369-3p的靶向基因參與的生物學(xué)進(jìn)程、信號(hào)通路是miR-369發(fā)揮作用的重要物質(zhì)基礎(chǔ)，可能是miR-369關(guān)聯(lián)PMOP與T2DM的分子基礎(chǔ)。CYR61是一種基質(zhì)細(xì)胞蛋白，可在礦化組織中發(fā)現(xiàn)，被認(rèn)為與骨組織再生、成骨分化有顯著關(guān)聯(lián)，Zhao等[28]對(duì)CYR61進(jìn)行敲除和過表達(dá)，證明CYR61通過Wnt信號(hào)通路調(diào)控成骨細(xì)胞血管形成，提高骨量，與PMOP的診斷預(yù)后密切相關(guān)；Feng等[29]發(fā)現(xiàn)CYR61在T2DM患者外周血中表達(dá)升高，可能是T2DM的生物標(biāo)志物。Diao等[30]通過對(duì)骨肉瘤患者測序芯片進(jìn)行差異分析，鑒別出其中的關(guān)鍵基因CALD1，認(rèn)為CALD1可能是骨肉瘤轉(zhuǎn)移的潛在靶點(diǎn)；Wang等[31]發(fā)現(xiàn)CALD1可能通過參與轉(zhuǎn)錄因子、miRNA的調(diào)控活動(dòng)，影響2型糖尿病微血管病變發(fā)展；DUSP1基因是MAPK信號(hào)通路、JNK信號(hào)通路等多個(gè)交互通路的明星分子，參與細(xì)胞自噬，細(xì)胞自噬被認(rèn)為是PMOP、T2DM病程進(jìn)展的一種重要調(diào)控機(jī)制[32-33]。

GO分析顯示，靶向基因主要富集于細(xì)胞外基質(zhì)組織、運(yùn)動(dòng)行為、內(nèi)質(zhì)網(wǎng)組成、染色質(zhì)結(jié)合等細(xì)胞基礎(chǔ)功能與結(jié)構(gòu)環(huán)節(jié)，同時(shí)與鋅離子穩(wěn)態(tài)、活性氧代謝等細(xì)胞循環(huán)活動(dòng)密切相關(guān)，證明hsa-miR-369-3p主要參與調(diào)控細(xì)胞的基礎(chǔ)生命活動(dòng)。KEGG結(jié)果中，礦物質(zhì)吸收與骨代謝活動(dòng)密切相關(guān)，鈣、磷等元素是維持骨代謝平衡的主要物質(zhì)，研究發(fā)現(xiàn)鈣離子通道TRPV5/6是PMOP的治療靶點(diǎn)[34]；cAMP信號(hào)通路在細(xì)胞對(duì)胞外刺激反應(yīng)中起重要作用，與鈣通道、鉀通道、鈉通道等蛋白磷酸化介導(dǎo)的興奮-收縮耦合反應(yīng)相聯(lián)系[35]；代謝途徑是機(jī)體生命活動(dòng)的基礎(chǔ)，研究發(fā)現(xiàn)腸道菌群對(duì)代謝物質(zhì)的調(diào)控是防治PMOP與T2DM的潛在新靶點(diǎn)[36]。hsa-miR-369-3p可能通過調(diào)控上述信號(hào)通路中關(guān)鍵基因的信號(hào)轉(zhuǎn)導(dǎo)，介導(dǎo)PMOP骨代謝與T2DM內(nèi)分泌代謝及相關(guān)細(xì)胞功能活動(dòng)。

本研究立足于探索miRNA在PMOP與T2DM中的聯(lián)系與調(diào)控作用，通過挖掘基因芯片數(shù)據(jù)，借助隨機(jī)森林算法篩選出關(guān)聯(lián)兩疾病的關(guān)鍵miRNA(hsa-miR-369-3p)，通過預(yù)測分析靶向基因參與的生物過程及信號(hào)通路，從分子機(jī)制水平理解hsa-miR-369-3p對(duì)PMOP與T2DM的潛在調(diào)控作用。但本研究僅從個(gè)別基因測序芯片數(shù)據(jù)出發(fā)，缺乏大樣本數(shù)據(jù)的驗(yàn)證支持，希望后續(xù)能繼續(xù)對(duì)hsa-miR-369-3p進(jìn)行實(shí)驗(yàn)驗(yàn)證以探究其具體調(diào)控機(jī)制，指導(dǎo)臨床。