基于生物信息学探究青娥丸治疗骨质疏松症的作用机制※

王 浩 谭国庆 薛海鹏 徐展望△

(1.山东中医药大学2019级硕士研究生，山东 济南 250035；2.山东中医药大学附属医院脊柱骨科，山东 济南 250035)

骨质疏松症(osteoporosis，OP)是一种骨量减少、骨密度减低、骨脆性增加，从而易发生骨质减少性骨折乃至威胁生命的全身性疾病[1-2]。在人口老龄化形势不断加剧的情况下，OP的预防与治疗成为人类健康管理的重要研究方向。现代医学治疗OP的药物一般包括骨矿化促进剂(钙剂、维生素D)、骨吸收抑制药(双膦酸盐、降钙素)、骨形成促进剂(氟制剂)，但其过高的费用、较多的毒副作用等问题难以解决[3]。随着中医学的发展，中医药治疗OP体现出不良反应发生率低、辨证论治特异性强、多种疗法相结合等独特优势，得到了广泛关注。

青娥丸最早见于《太平惠民和剂局方》，全方由杜仲、补骨脂、核桃仁、大蒜4味组成，方中杜仲温补肝肾，强筋壮骨；补骨脂补肾壮阳，固精缩尿，温脾止泻；核桃仁补肾温肺，润肠通便；大蒜温中行滞，解毒消肿。全方共奏补肾强骨生髓之效，是治疗OP的经典名方。现代医学研究认为，青娥丸有抑制破骨细胞的骨吸收活动、促进成骨细胞的增殖、增强骨的质量和数量等多方面作用[4-5]。本研究基于生物信息学和网络药理学方法，科学系统地探究青娥丸治疗OP的作用机制，为该方在临床应用中提供理论支持，并进一步探索中药组方治疗OP的新方向和新思路。

1 资料与方法

1.1 高通量基因表达数据库(Gene Expression Omnibus，GEO)基因芯片数据分析OP差异表达基因 以“osteoporosis”为关键词检索GEO中有关OP的实验数据集，获得满足条件的编号为GSE100609的数据集，其中包含了4例从OP女性患者单核细胞中分离的破骨细胞前体细胞为实验组，4例属于非OP女性单核细胞中分离的破骨细胞前体细胞为对照组。利用生物信息学之RMA(robust multiarray average)算法背景校正和矩阵数据归一化处理该数据集，利用R语言技术分析芯片数据，并结合P值和差异倍数(FC)进行筛选，其中显著差异基因的筛选条件为P< 0.05，logFC绝对值>0.5，最终获得OP的差异表达基因。

1.2 OP疾病靶基因预测 以“osteoporosis”作为关键词，采用疾病基因数据库(DisGeNET)、人类基因数据库(GeneCards)、人类孟德尔遗传综合数据库(OMIM)、人类疾病数据库(MalaCards)、遗传药理学与药物基因组学数据库(PharmGKB)等多个数据库对OP已知靶基因进行检索与筛选，合并各数据库的疾病靶基因，将合并后的疾病靶基因与R语言技术分析芯片获得的差异表达基因取并集，删除重复的靶基因后，即获得已明确的OP发病过程的疾病靶点。

1.3 青娥丸有效成分筛选及作用靶点 通过中药系统药理学数据库与分析平台(TCMSP)和中药分子机制的生物信息学数据库(BATMAN-TCM)，检索杜仲、补骨脂、核桃仁、大蒜4味中药化学成分，TCMSP数据库中以口服生物利用度(oral bioavailability，OB)≥30%和类药性(drug likeness，DL)≥0.18作为筛选条件，BATMAN-TCM数据库中以Score off≥20和P≤0.05作为筛选条件，筛选出青娥丸有效成分并将结果导入蛋白质数据库(Unitprot)进行转换处理，得到青娥丸有效成分的作用靶点。

1.4 映射靶基因获取 将所筛选的青娥丸有效成分的作用靶点与OP疾病靶基因交叉生成Venny图，得出青娥丸与OP的映射靶基因。

1.5 构建“青娥丸药物-有效成分-靶点”网络图 将映射靶基因导入Cytoscape3.7.2软件，构建“青娥丸药物-有效成分-靶点”网络图，并进行网络拓扑分析。

1.6 构建蛋白质-蛋白质相互作用网络图(PPI网络图) 抽取映射靶基因，借助Cytoscape3.7.2软件构建青娥丸治疗OP的PPI网络图，根据节点连接度(Degree)值的大小排名顺序进行筛选，得到青娥丸治疗OP的关键靶基因，所获取的靶基因则是直接或间接调控发挥治疗OP作用的关键靶点，揭示出青娥丸治疗OP的关键机制。

1.7 基因本体(GO)生物过程及京都基因与基因组百科全书(KEGG)通路富集分析 将映射靶基因导入基因功能注释数据库(DAVID)，进行KEGG通路及GO生物过程的富集分析。以P≤0.05，FDR≤0.05为筛选标准，将筛选出的排名前20个GO条目借助R语言ggplot2包绘制富集分析气泡图，借助Cytoscape3.7.2软件构建映射靶基因在KEGG通路上的富集网络图。

2 结果

2.1 OP相关差异表达基因检索 通过OP差异表达基因分析结果显示，发现542个mRNA存在差异表达，其中上调188个，下调354个，用R语言中heatmap与limma程序包绘制芯片原始的差异基因热图与火山图(见图1、2)。

图1 OP相关芯片差异基因聚类图

2.2 OP疾病靶基因筛选 通过DisGeNET、GeneCards、OMIM、MalaCards、PharmGKB数据库对OP已知靶基因进行检索与筛选，合并后共获得3576个疾病靶基因，将合并后的疾病靶基因与542个差异表达基因取并集并删除重复值，获得已明确的OP发病过程的疾病靶点共3811个。

2.3 青娥丸有效成分筛选及作用靶点 在TCMSP及BATMAN-TCM检索到4味中药活性成分共计48个，相对应靶点共计1757个，删除重复值后得到青娥丸最终靶点973个。

2.4 映射靶基因获取 将所筛选的青娥丸973个靶基因与2.2获取的3811个OP疾病靶基因交叉生成Venny图(http://bioinfogp.cnb.csic.es/tools/venny/index.html，见图3)，通过Venny图可以直观发现青娥丸对OP靶基因映射率为10%。该方共有973个有效靶点，其中映射靶点433个，作为青娥丸治疗OP的直接或间接作用靶基因，占该方作用靶点的44.5%左右，表明该复方在治疗OP方面具有显著特异性。

图2 OP相关芯片差异基因火山图

图3 “青娥丸-OP”映射靶基因Venny图

2.5 青娥丸药物-有效成分-靶点网络图构建及分析 借助Cytoscape3.7.2软件构建“青娥丸药物-有效成分-靶点”网络图(见图4)。按照Degree值排名顺序排名前4位的有效成分为：β-谷甾醇、豆甾醇、异补骨脂素、槲皮素，这4种有效成分与多个作用靶点相互关联，说明这些化合物是青娥丸治疗OP的主要有效成分。

图4 青娥丸药物-有效成分-靶点网络图

2.6 PPI网络图构建及分析 将433个映射靶基因导入Cytoscape3.7.2软件，根据Degre值≥45为筛选条件，得到排名前27位的映射靶基因为青娥丸治疗OP的关键靶基因，并构建青娥丸治疗OP的关键靶基因PPI网络图(见图5)，所获取的靶基因则是直接或间接调控发挥治疗OP作用的关键靶基因，揭示出青娥丸治疗OP的关键机制。

图5 青娥丸治疗OP关键靶基因PPI网络图

2.7 GO生物过程及KEGG通路富集分析 对2.4获得的映射靶基因采用DAVID数据库进行GO生物过程和KEGG信号通路分析，得到分子功能、生物过程和细胞成分条目188个，KEGG通路118条，以P≤0.05，FDR≤0.05为标准，并以P值升序排列，分别筛选出符合条件的20个生物过程和20条信号通路，并将GO生物过程以气泡图的形式展现出来(见图6)，将KEGG通路以表格形式展现出来(见表1)，并借助Cytoscape3.7.2软件构建映射靶基因在KEGG通路上的富集网络图(见图7)。

图6 青娥丸治疗OP映射靶基因的GO生物过程气泡图

表1 青娥丸治疗OP的KEGG富集通路

图7 青娥丸治疗OP映射靶基因在KEGG通路上的富集网络图

3 讨论

随着近年来现代医学对OP机制的认识，在中医理论指导下，中药复方青娥丸防治OP在临床和实验中均取得了一定成果[1]。经活性物质筛选及靶点预测发现，青娥丸具有诸多活性物质协同作用调控OP进程。植物甾醇作为一种天然活性成分，广泛存在于各种植物中，超过250种被学者研究报道，以β-谷甾醇、豆甾醇(stigmasterol)最常见[2]，其生理功能包括抗高脂血症、抗癌、抗炎、抗氧化和免疫调节等诸多方面[3]。植物甾醇通过抑制小肠对胆固醇的吸收来调节脂代谢，脂代谢异常引起肥胖可影响骨凋亡、骨形成、骨分化、骨微循环、骨吸收等多条信号通路的改变[4]。异补骨脂素(angelicin)是补肾中药补骨脂有效成分之一，通过激活骨矿化过程激活成骨细胞活性，调节骨代谢，预防骨质疏松。张军芳等[5]经体外细胞培养初步证明了异补骨脂素及锌制剂共同作用能促进成骨细胞增殖与分化。进一步试验表明，异补骨脂素在1×10-5mol/L浓度时能最大程度使大鼠成骨细胞增殖分化成熟，显著提高细胞碱性磷酸酶(ALP)，增加钙含量，促进骨钙素分泌[6]。补骨脂另一成分补骨脂酚(backuchiol)对于去势导致的雄性小鼠OP有较好治疗作用[7]，对雌激素受体(ER)α同样具有很强的结合亲和力，减少绝经后骨质流失[8]。研究表明，补骨脂酚干预OP的可能机制为抑制蛋白激酶B(Akt)和激活蛋白-1(AP-1)途径改善巨噬细胞刺激集落因子(M-CSF)和核因子κB受体活化因子配体(RANKL)诱导的破骨细胞分化和骨吸收[9]。槲皮素(quercetin)作为经典的抗衰老活性成分，能有效清除衰老的骨髓间充质干细胞，促进骨髓间质干细胞的增殖，刺激成骨细胞活性，改善骨质状态及结构[10]。槲皮素以其雌激素拮抗剂活性在绝经后OP治疗中应用广泛，临床OP常规治疗基础上，补充槲皮素明显改善骨密度、疼痛、骨代谢等指标[11]。青蛾丸具有肝肾同调、精血并补功效，组方中药富含多种活性成分，能同时影响骨代谢和骨吸收过程，对于改善慢性疼痛、骨丢失具有重要作用。

通过靶点映射发现，大量靶点共同存在于青娥丸及OP，且部分靶基因同时存在于诸多活性物质及通路，说明青娥丸通过调控复杂靶基因网络，协同各活性物质共同发挥作用。胰岛素(INS)为典型降血糖成分，INS缺乏易引起糖脂代谢紊乱，导致体内钙、磷元素代谢失调,使骨组织微结构破坏，加速OP进程[12]。INS一方面作用于成骨细胞表面受体，上调成骨细胞活性；另一方面调节血糖水平，抑制高糖对成骨细胞的毒性作用[13]。诸多研究发现，葡萄糖不仅对成骨细胞具有毒性作用，同时破骨细胞骨吸收活动对于葡萄糖具有依赖性，破骨细胞在7～25 mmol/L浓度时使骨吸收最活跃[14-15]，因此推测INS是青娥丸在糖尿病性OP临床治疗中具有良好疗效的关键基础靶基因。OP发生的根本原因是成骨与破骨之间平衡失调，而蛋白激酶B1(Akt1)正是控制这种平衡的传导中间体，调控成骨细胞与破骨细胞的成熟分化[16]。而且Akt1不同表达程度也会产生不同生理作用，Tu Q S等[17]使用脂联素下调Akt1活性来抑制破骨细胞作用，而当Akt1过度表达时又成功逆转了对破骨细胞分化的抑制作用。在成骨方面,地塞米松能间接活化Akt1诱导成骨细胞及骨细胞的凋亡,而使用Akt1抑制剂能明显降低成骨细胞凋亡[18]。白细胞介素6(IL-6)是生物学效应广泛的细胞炎症因子，广泛参与破骨细胞分化，最近数据显示IL-6尤其在绝经期后OP妇女中显著增加[19]。有关临床数据表明，雌激素较低时OP患者IL-6水平会高达210.0 ng/L[20]，进一步证明了IL-6在绝经后OP发病过程中的重要作用，其具体作用机制可能为绝经后妇女体内雌激素水平下降，导致雌激素抑制成骨细胞生成肿瘤坏死因子α(TNF-α)能力减弱，TNF-α的增加进一步刺激成骨细胞产生IL-6，进而发挥抑制成骨、促进破骨的作用[21]。血管内皮生长因子(VEGF)一方面通过血管构建对骨组织进行营养供给，增加骨量和治疗OP[22]；另一方面通过血管构建，维持骨组织稳定氧环境，抑制因缺氧分泌的炎症因子进一步加快骨流失。最重要的是VEGF能产生成骨性生长因子，调节成骨细胞的增殖分化及破骨细胞活性[23]。综上所述，众多关键靶点相互作用，产生复杂作用网络，调控炎性反应及糖脂代谢等诸多相关生理过程，这可能是青蛾丸干预OP的可能内在机制。

依托DAVID数据库进行基因可视化。KEGG通路富集分析发现诸多通路关于各种癌症，癌症相关通路普遍见于OP的诸多通路研究之中，癌症患者不仅因自身病变引起代谢异常而面临原发性OP的风险，而且还可能因特别是削弱性腺功能的相关治疗引起继发性OP[24]。研究显示，性激素剥夺治疗的前列腺及乳腺癌患者骨量丢失尤其明显[25]，临床数据显示前列腺癌去势治疗后OP患病率为37%。去势治疗后降低雄激素，经芳香化而来的雌激素也随之减少，进而破骨细胞增殖，活性上调，骨吸收增强[26]，因此在延长寿命治疗同时也应及时纠正骨量丢失。有研究显示，选择性抑制破骨细胞中缺氧诱导因子1α(HIF-1α)的表达，能减弱破骨细胞活性，骨吸收减少，改善骨量[27]。除改善骨代谢外，HIF-1一方面能调节组织、细胞适应低氧或缺氧环境，缺氧对大鼠成骨细胞的骨形成能力起到抑制作用[28]；另一方面，HIF-1也能调控血管生成构建，与VEGF协同作用促进骨血管为骨组织提供氧、营养、激素、细胞因子等物质[29]。越来越多证据表明，肝脏与骨之间有着复杂联系，非酒精性脂肪性肝病(NAFLD)更是能够通过多条通路影响骨骼发育成熟。临床研究表明，NAFLD患者腰椎、髋部等部位的骨密度(BMD)明显下降[30]。同时,NAFLD与绝经后女性BMD关系更密切，NAFLD被认为是OP的高危因素。NAFLD从肝纤维化、胰岛素抵抗(IR)、脂肪因子、胃饥饿素等多个方面对BMD产生影响。研究认为，脂毒性更易引起低骨量和OP, 胆固醇及其代谢产物能影响成骨细胞和破骨细胞活性，降低BMD。基于以上分析，青娥丸通过调控平衡多种细胞因子及代谢干预OP，尤其在性激素代谢方面调控作用显著，推测这是临床治疗绝经后OP疗效显著的原因。

本研究基于生物信息学相关技术手段，筛选青娥丸有效活性物质及作用靶点，并与OP相关疾病靶点映射分析获得关键靶基因，通过基因可视化寻找青娥丸可能作用通路，相对科学可靠地预测青娥丸干预OP的分子作用机制，为青娥丸的应用及OP的防治提供新思路。但这种信息整合的可靠性依托数据库的信息准确与否，具体作用机制仍需进一步试验研究的验证。