隨機森林算法預測絕經后骨質疏松癥與2型糖尿病關聯microRNA及其靶向基因調控作用

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

廣州醫科大學附屬中醫醫院骨傷科，廣東 廣州 510000

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

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

1 材料與方法

1.1 芯片數據獲取與預處理

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

1.2 隨機森林算法篩選PMOP與T2DM關聯miRNA

使用Randomforest軟件包構建隨機森林模型，通過隨機生成大量的分類樹并對每棵樹的miRNA分類結果進行迭代評分獲得分類結局，最終對所有單棵樹的分類結果進行綜合判定。使用Caret軟件包對模型中miRNA重要性進行排序，選取排名前10的miRNA，并根據miRNA的表達量繪制差異miRNA散點圖。

1.3 評價納入miRNA的預測性能

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

1.4 miRNA靶向基因預測篩選

使用Targetscan[15]、miRWalk2.0[16]和DIANA TOOLS[17]3個公共數據庫預測PMOP與T2DM關聯miRNA靶向基因，并對多數據庫預測結果取共交集，獲得miRNA靶向基因[18]。

1.5 構建miRNA-靶向基因調控互作網絡

將靶向基因上傳至蛋白互作分析數據庫STRING[19]進行蛋白互作分析，預測靶向基因的蛋白-蛋白互作(protein-protein interaction，PPI)關系，使用網絡構建工具Gephi軟件構建miRNA-靶向基因調控互作網絡。

1.6 基因本體論與信號通路富集分析

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

2 結果

2.1 隨機森林算法篩選miRNA

通過對GSE70318芯片原始數據進行標準化及基因注釋，共獲得153個miRNA，構建隨進森林算法模型，對分類樹分類結果評分，篩選PMOP與T2DM關聯重要性排名前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，結果如圖1所示。將上述miRNA分為PMOP-T2DM關聯組(19例)與對照組(38例，單純PMOP或T2DM)進行對比，繪制表達量差異散點圖(見圖2)。結果表明，hsa-miR-188-3p、hsa-miR-181a-3p、hsa-miR-135a-5p和hsa-miR-369-3p表達量在關聯組與對照組間具有顯著統計學差異(P<0.01)。

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

圖2 miRNA表達差異散點圖Fig.2 Scatter diagram of miRNAs注：Group1：對照組；Group2：關聯組。

2.2 ROC曲線驗證miRNA重要性

使用標準化的GSE74209芯片作為驗證集數據，對顯著差異的miRNA繪制ROC曲線進行預測性能驗證。分別對hsa-miR-188-3p、hsa-miR-181a-3p、hsa-miR-135a-5p、hsa-miR-369-3p進行驗證，同時將上述4個miRNA作為一個集合模型共同驗證以增加驗證準確性。ROC曲線結果如圖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，認為hsa-miR-369-3p在PMOP與T2DM疾病關聯性中預測性能更好；同時，將上述4個miRNA構建集合模型進行ROC驗證，其AUC為0.667，認為該集合模型未能提高預測準確性。據此，選擇hsa-miR-369-3p作為關聯PMOP-T2DM的關鍵miRNA。

圖3 ROC曲線圖Fig.3 ROC diagram

2.3 關聯miRNA靶向基因預測及共表達網絡構建

圖4 韋恩圖與調控互作網絡 A：靶向基因篩選共交集；B：關聯miRNA-靶向基因調控網絡。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數據庫，設置物種為“Homo sapiens”(人類)，進行靶向基因預測。Targetscan預測獲得2 324個靶向基因，miRWalk2.0獲得151個靶向基因，DIANA TOOLS獲得89個靶向基因，對3個數據庫獲得的靶向基因取共交集后最終獲得44個目標靶向基因(圖4A)。

將靶向基因導入STRING數據庫，設置關聯置信度為<0.40，物種為“Homo sapiens”，進行蛋白互作分析，將分析結果導入Gephi軟件中，構建miRNA-靶向基因調控互作網絡(圖4B)。對靶向基因在網絡中關聯信度進行計算，最終獲得網絡核心基因為：CYR61、CALD1、DDTT4和DUSP1。

2.4 GO與KEGG富集分析

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

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

3 討論

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

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

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

GO分析顯示，靶向基因主要富集于細胞外基質組織、運動行為、內質網組成、染色質結合等細胞基礎功能與結構環節，同時與鋅離子穩態、活性氧代謝等細胞循環活動密切相關，證明hsa-miR-369-3p主要參與調控細胞的基礎生命活動。KEGG結果中，礦物質吸收與骨代謝活動密切相關，鈣、磷等元素是維持骨代謝平衡的主要物質，研究發現鈣離子通道TRPV5/6是PMOP的治療靶點[34]；cAMP信號通路在細胞對胞外刺激反應中起重要作用，與鈣通道、鉀通道、鈉通道等蛋白磷酸化介導的興奮-收縮耦合反應相聯系[35]；代謝途徑是機體生命活動的基礎，研究發現腸道菌群對代謝物質的調控是防治PMOP與T2DM的潛在新靶點[36]。hsa-miR-369-3p可能通過調控上述信號通路中關鍵基因的信號轉導，介導PMOP骨代謝與T2DM內分泌代謝及相關細胞功能活動。

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