用生物信息学分析预测绝经后骨质疏松症核心基因与互作miRNA的研究

柴毅 谭峰 樊巧玲

南京中医药大学，江苏 南京 210046

骨质疏松症(osteoporosis，OP)是最常见的骨骼疾病，以骨量低，骨组织结构破坏，最终导致骨脆性增加，骨强度下降及骨折风险增加，易发生骨折为特征的全身性骨病[1]。在美国，每年约有150万人次发生骨折，绝大多数都发生于绝经后妇女[2]。绝经后骨质疏松症(postmenopausal osteoporosis，PMOP)是与年龄相关的衰老性疾病，多发生于绝经2年以上，70岁以下的妇女。PMOP的病变是一个隐性的过程，全世界约有50%的绝经后妇女受到影响，被认为是老年人发病率最高的疾病之一[3]。PMOP已成为全球经济的负担之一，积极开展对PMOP的预防和治疗是公共卫生的重要任务。

OP的诊断基于全面的病史记录、体格检查、骨密度测定、影像学检查和必要的生化测定，OP的诊断主要基于双能X线吸收检测法(dual-energy X-ray absorptiometry，DXA)骨密度测量结果与脆性骨折[4]。然而DXA检测也存在一些问题，例如DXA检测成本较高，不同设备的DXA检测结果存在差异[5-6]。挖掘与OP或PMOP病理相关的核心基因，特别是高敏感性与高特异性的基因是预防和治疗该疾病性价比较高的途径之一。早期的研究发现了一些与OP有关的基因。例如，细胞周期蛋白E1(Cyclin E1，CCNE1)是细胞周期的调控因子。CCNE1参与了骨代谢过程，CCNE1在PMOP B细胞中的表达呈下调趋势[7]。又如细丝蛋白α(filamin A alpha，FLNA)是参与破骨细胞生成过程的关键因子，高表达的FLNA可促进破骨细胞生成[8]。微小核糖核酸(micro-ribonucleic acid，miRNA)是一类非编码蛋白的小RNA，miRNA通过抑制特定靶点的mRNA调控基因表达[9]。一些miRNA亦可作为OP或PMOP的敏感标志物或治疗药物的靶点，因而针对miRNA在OP或PMOP中的机制挖掘是很有必要的[10-11]。目前，治疗OP或PMOP的相关机制的探索很少深入到非编码RNA的作用层面。此外，利用生物信息学研究方法在OP或PMOP方面有关核心基因挖掘的研究数量相对较少，且现有关于鉴定OP或PMOP核心基因的生物信息学研究对与核心基因相关互作的miRNA预测的报道更为稀缺。

本研究通过对基因芯片GSE57273进行生物信息学分析筛选核心基因，并预测与这些核心基因相互作用的miRNA，为PMOP建立新的科学假说以及后续更深入的研究提供依据，并为PMOP诊断以及治疗药物的研发提供较为可靠的路径和作用靶点。

1 材料与方法

1.1 数据收集

从GEO公共数据库(https://www.ncbi.nlm.nih.gov/geo/)下载基因表达芯片GSE57273。该芯片包含3组PMOP药物干预前后样本。其所处的平台(Platforms)为GPL4133Agilent-014850 Whole Human Genome Microarray 4x44 K G4112F(Feature Number version)。另下载Series Matrix File以便后续使用。

1.2 差异基因的数据分析

使用GEO2R(https://www.ncbi.nlm.nih.gov/geo/geo2r/)和Morpheus(https://software.broadinstitute.org/morpheus/)在线分析软件进行差异基因(differentially expressed genes，DEGs)的甄别和筛选。GEO2R是GEO数据库自带的公共在线分析工具，它可将GEO数据进行复杂的R语言分析，从而呈现出每个基因的计算结果。经原始数据进行成组t检验统计学分析，以adj.P<0.01和|logFC|≥3作为DEGs的筛选条件[12]。

1.3 基因功能富集分析

GO(gene ontology)是常用的分析方法，其主要功能是注释基因或其产物并识别高通量基因组或转录组数据的特征生物学特性。GO按照生物途径(biological process，BP)、分子功能(molecular function，MF)、细胞定位(cellular component，CC)对基因进行注释和分类。此外，KEGG(Kyoto Encyclopedia of Genes and Genomes)数据库可供查询通路信息和信号通路检索等。KEGG通路分析是另一种常用的基因功能富集分析方法。本研究应用DAVID(Database for Annotation, Visualization and Integrated Discovery)进行在线分析提供所需的GO和KEGG生物功能富集数据。本研究使用的DAVID数据库版本为6.8，地址为https://david.ncifcrf.gov/，由美国国立变态反应与传染病研究所提供研究服务。使用Fisher Exact或EASE Score统计方法，GO各项以P<0.05且FDR<0.05为筛选条件，KEGG各项以P<0.05为筛选条件[13]。

1.4 蛋白互作网络与模型分析

STRING(Search Tool for the Retrieval of Interacting Genes)是一款可以用来呈现与评估蛋白互作(protein-protein interaction，PPI)的在线分析工具。STRING中的所有数据和下载文件都可以在“Creative Commons BY 4.0”许可下免费获取。本研究将所筛选的所有DEGs植入STRING(版本10.5，https://string-db.org/)分析工具试探它们之间潜在的联系。置信度(confidence score)≥0.4，互作最大值(maximum number of interactors)=0设为筛选条件[14]。此后，把STRING的计算结果导入Cytoscape(版本3.6.0)进行MCODE(Molecular Complex Detection)分析以挖掘PPI中连接最为紧密的集簇。本研究使用的MCODE版本为1.5.1，设置参数为degree=2，node score=0.2，k-core=2，max. depth=100[15]。

1.5 miRNA-基因调控网络的预测

CyTargetLinker可以扩展生物调控互作网络(regulatory interaction networks，RegINs)，由荷兰系统生物学联合会提供支持。它涵盖了miRNA—靶点、转录因子—靶点和药物—靶点之间的互作关系。本研究下载了人类物种基因数据集(https://projects.bigcat.unimaas.nl/cytargetlinker/regins/)。选用该数据集中基于实验验证的miRTarBase 4.4数据库(含20 942个RegINs)，基于预测功能的TargetScan 6.2数据库(含511 040个RegINs)和MicroCosm 5数据库(含541 039个RegINs)预测核心基因与miRNA之间的调控关系。

2 结果

2.1 DEGs的鉴定

本研究选用基因表达芯片GSE57273，经GEO2R初步分析共获得32 996个DEGs，随后由Morpheus分析并经条件筛选，最终获得841个DEGs，其中包含826个下调基因和15个上调基因。

2.2 GO功能与KEGG信号通路富集分析

根据GO的分析结果，本研究以P<0.05，FDR<0.05为筛选条件，并按照计数值从大到小排列，在BP、CC与MF类别中各选取前3项列为表1。可以看出，在生物学过程中，这些DEGs主要参与了基因表达，细胞大分子生物合成和RNA代谢过程；在细胞定位中，这些DEGs富集于核浆、细胞质基质以及粘附连接；从分子功能上看，这些DEGs具有使有机环状化合物结合、杂环化合物结合和核酸结合的作用。根据KEGG分析结果，本研究以P<0.05为筛选条件，按照计数值从大到小将DEGs富集的信号通路列为表2。结果显示，这些DEGs主要富集于癌症信号通路、病毒致癌通路、粘附斑、rap1信号通路和内质网蛋白加工通路。

表1 与PMOP相关的DEGs的GO富集分析

表2 与PMOP相关的DEGs的KEGG富集分析

图2 蛋白质互作网络的前3个集簇模块

2.3 核心基因的鉴定与筛选的集簇模块

通过STRING的PPI构建，经由Cytoscape对网络的计算工具得出所有DEGs的连接度(degree)(图1)。degree值表示网络中某一基因与周围基因的关系数量，因此degree越大代表与它相互作用关系的基因数量就越多。本研究按degree从高至低进行排序，以排名前10位的高degree基因定为核心基因，它们分别是HSP90AA1(degree=75)、EP300(degree=55)、SMARCA2(degree=44)、RANBP2(degree=41)、ASH1L(degree=36)、EIF4E(degree=35)、PTEN(degree=31)、CNOT6L(degree=30)、RPL7(degree=29)、KRAS(degree=29)。此外，MCODE分析发现17组集簇，共包含523个节点(node)与2 026条连线(edge)。本研究以score为依据展示前3组集簇并分析各个集簇所富集的通路(图2)。其中集簇模块A富集核糖体(hsa03010：Ribosome)与mRNA监视通路(hsa03015：mRNA surveillance pathway)(P<0.05)，集簇模块B富集泛素介导的蛋白质水解4通路(hsa04120：Ubiquitin mediated proteolysis 4)(P<0.05)，而集簇模块C未鉴定出具有统计学意义的通路(P>0.05)(表3)。

图1 DEGs的PPI网络

2.4 miRNA—核心基因调控网络模块的鉴定

本研究应用CyTargetLinker预测可以与上述筛选出的10个核心基因相互作用的miRNA。结果显示，在MicroCosm数据库中有258个预测的miRNA靶点互作关系，在TargetScan数据库中有1 171个预测的miRNA靶点互作关系，总共有875个节点与1 429条连线。另外，阈值(threshold)可对结果显示的可视化调控网络进行支持数据库的叠加筛选，通过调设阈值可控制调控网络的显示结果，其设置范围一般为1～3[16]。本研究将阈值设为2，结果显示共有37个miRNA与7个靶基因存在互作关系。这些基因与预测的miRNA如表3所示。

3 讨论

本研究通过对GEO数据库中的基因芯片GSE57273进行生物信息学分析获得DEGs，并分析了有关这些基因富集的生物过程、细胞定位、分子功能和信号通路为机制研究提供了理论依据与研究方向，通过挖掘与核心基因互作的miRNA为PMOP的研究提供新思路。

GO与KEGG分析有助于更深入地认识并筛选出DEGs的功能和作用。由于受到数据呈现空间的限制，本研究未能把这些DEGs所富集的生物学过程与信号通路全部列出。本研究由KEGG筛选出的富集通路以癌症通路为首。从参与的DEGs来看，该通路仅包含了少部分核心基因，由于该通路涉及了大部分非核心基因，这可能导致它们富集的通路与PMOP不同。另外，就仅涉及的核心基因而言，现有关于这些基因功能的研究也存在局限性，大部分集中于肿瘤领域的研究，而本研究结果为这些核心基因参与PMOP提供了一定依据，有助于拓展核心基因的功能。

表3 由CyTargetLinker扩展网络分析预测的与7个核心基因互作的miRNA

本研究罗列的核心基因主要富集于RNA代谢过程、基因表达过程、PI3K-Akt信号通路和癌症相关的信号通路等。HSP90AA1编码的蛋白质是一种功能类似于同型二聚体的诱导型分子。HSP90AA1参与细胞的生长发育过程，有研究证实使用HSP90AA1抑制剂可迅速导致细胞死亡，表明HSP90AA1在细胞活动中发挥重要作用[17-18]。因而，验证HSP90AA1是否参与PMOP的骨偶联的调控过程值得深入探究。EP300编码与腺病毒E1A相关的细胞p300转录共激活蛋白。它与组蛋白乙酰转移酶的功能类似，可以通过染色质重塑调节转录并且在细胞增殖和分化过程发挥重要作用。EP300在骨髓中高表达，并参与了细胞成骨分化和骨量减少的调控过程。EP300在髓核细胞中可受到骨形态发生蛋白2和骨形态发生蛋白7的调控，为EP300与PMOP的关联提供理论支撑[19-20]。由SMARCA2基因编码的蛋白属于SWI/SNF家族蛋白，并且这种蛋白与果蝇的brahma蛋白高度相似。该家族蛋白具有解旋酶和ATP酶活性，其通过改变基因周围的染色质结构来发挥调节基因转录的功能。SMARCA2在卵巢、大脑等处高表达，而这类蛋白的减少会影响间充质干细胞的成骨分化和成脂分化平衡，进而导致OP[21]。RANBP2编码与核孔复合物免疫定位有关的RAN结合蛋白，主要在睾丸、甲状腺、骨髓和大脑等处高表达。RAN是与核膜相关的RAS超家族的小GTP结合蛋白，它通过与蛋白质的相互作用来调控多种细胞功能。RANBP2细胞定位为核孔复合体，并且在神经视网膜中的表达非常丰富[22]。ASH1L编码转录激活因子的三空腔结构蛋白质，它同样在睾丸、甲状腺、骨髓和大脑等处高表达。现有研究证实ASH1L可以参与细胞分化以及骨髓造血的调控[23-24]。EIF4E基因编码的蛋白质是真核生物翻译起始因子4F复合物的组成部分。研究发现，EIF4E在病理状态下的骨髓间充质干细胞(bone marrow stromal cells，BMSCs)中表达水平会发生变化[25]。当然，EIF4E作为一种原癌基因，其表达和激活与转化和肿瘤发生密切相关。PTEN编码的蛋白质是磷脂酰肌醇-3，4，5-三磷酸3-磷酸酶，该基因被认为是在大量癌症中高频率突变的肿瘤抑制因子。近年来有研究已经发现PTEN/PI3K/AKT信号通路可能是调控BMSCs增殖和分化的途径之一[26]。CNOT6L可能参与包括细胞增殖的各种细胞活动，有研究证实CNOT6L可通过p53介导的信号通路途径来调节细胞周期阻滞和衰老[27-28]。RPL7编码一种由60S亚基组成的核糖体蛋白，该蛋白属于核糖体蛋白L30P家族，主要在卵巢与骨髓高表达。KRAS编码一种属于小GTP酶超家族成员的蛋白。有研究证实KRAS可以与骨形态发生蛋白4靶向性结合[29]。

此外，本研究还扩展了与核心基因相互作用的miRNA。与HSP90AA1相互调控的miRNA参与了类固醇生物合成过程。与EP300相互作用的miRNA参与D-谷氨酰胺和D-谷氨酸代谢和细胞周期调控过程。与SMARCA2靶向调控的miRNA富集于调节干细胞多能性的信号通路、肌动蛋白细胞骨架的调控与甲状腺激素信号通路等。与RANBP2互作的miRNA功能富集于蛋白的内质网加工过程和雌激素信号通路等。hsa-miR-137、hsa-miR-219-5p、hsa-miR-548 d-3p可调控ASH1L，参与了细胞氮化合物代谢过程。hsa-miR-134、hsa-miR-935可调控KRAS，目前尚未发现它们所富集的信号通路，因此仍需进一步研究探索。与EIF4E互作miRNA的功能富集于细胞周期调控和糖胺聚糖的生物合成等。

本研究预测的37个miRNA与核心基因为研究PMOP提供了新的角度。然而，本研究所得到的结论尚需要更多的研究验证这些核心基因和miRNA在PMOP病理进展中的特异性与可靠性。随着未来针对这些核心基因及miRNA不断地深入研究，它们或可应用于PMOP病理的筛查，或为PMOP药物研发提供作用靶点，或有助于对PMOP合并肿瘤患者进行基因诊断。此外，上述核心基因及miRNA或可作为遗传致病因子的筛查点并为此定制个性化的防治措施。由于PMOP机制复杂，影响因素众多，探寻并梳理具有逻辑关系和实际意义的基因及互作集簇，找出后续研究工作的切入点是需要进一步思考的问题。