随机森林算法预测绝经后骨质疏松症与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进行实验验证以探究其具体调控机制，指导临床。