基于加权基因共表达网络分析和分子对接探讨青娥丸治疗腰椎间盘突出症的作用机制*

金晖,葛殊玮,彭广,吴洁,吴睿哲,曾景奇,王凡

湖南中医药大学第二附属医院,湖南 长沙 410005

腰椎间盘突出症(lumbar disc herniation,LDH)是一种常见的脊柱科疾病,主要表现为腰背部疼痛、下肢放射性疼痛、麻木以及腰部活动障碍等症状。该病的发病机制主要是腰椎间盘发生退变,在其他因素的共同作用下,导致纤维环破裂,髓核突出,从而压迫或刺激神经根[1]。目前,LDH的治疗手段主要分为非手术治疗和手术治疗两大类。尽管手术治疗能够迅速缓解疼痛、恢复腰椎功能,但手术风险较大,且易于复发,这使得非手术治疗成为一种更为安全有效的治疗方式[2]。在非手术治疗中,中医药疗法占据着重要地位,不仅具有方便、经济、高效等优势,而且在延缓LDH进程、改善症状和提升生活质量方面发挥着不可替代的作用[3]。

青娥丸作为补肾强腰的代表方剂,具有温补肝肾、活血化瘀、消肿止痛的功效[4]。多项研究已证实,青娥丸治疗腰椎间盘突出症的疗效显著,能够有效缓解临床症状,并改善患者的生活质量[5]。目前,青娥丸治疗LDH的作用机制尚未完全阐明。近年来,随着加权基因共表达网络分析(weighted gene co-expression network analysis,WGCNA)和网络药理学技术的不断发展,这些工具已广泛应用于药物作用机制的研究。因此,本研究拟采用WGCNA和分子对接技术探索青娥丸治疗LDH的作用机制,以期为该病的临床治疗提供新思路。

1 资料与方法

1.1 基因芯片数据获取以“lumbar disc herniation”为关键词,在GEO数据库(www.ncbi.nlm.nih.gov/geo)[6]中检索并下载数据集“GSE124272”。

1.2 加权基因共表达网络分析首先,根据绝对中位差(median absolute deviation,MAD)剔除数据集“GSE124272”的离群值。然后,使用R 4.2.2软件的“WGCNA”包[7]分析数据。通过pickSoft Threshold函数确定合适的软阈值(β),采用Dynamic Tree Cut算法确定基因模块。为保证结果的准确性,设定每个模块至少包含300个基因。通过模块基因与性状的关联分析,筛选与LDH发病高度相关的模块。通过模块基因与表型数据(如疾病状态、身高、年龄、性别等)的关联分析,以Module Membership(MM)≥0.8为条件筛选重要基因。

1.3 青娥丸活性成分及作用靶点筛选利用中药系统药理学数据库与分析平台(Traditional Chinese Medicine Systems Pharmacology Database and Analysis Platform,TCMSP)和中药分子机制生物信息学分析工具(Bioinformatics Analysis Tool for Molecular Mechanism of Traditional Chinese Medicine,BATMAN-TCM)数据库检索青娥丸(杜仲、核桃仁、补骨脂和大蒜)的活性成分[8-9]。其中,TCMSP数据库的筛选条件为口服生物利用度(oral bioavailability,OB)≥30%和类药性(drug likeness,DL)≥0.18,BATMAN-TCM数据库的筛选条件为Score≥20和P≤0.05。然后,基于Uniprot数据库(https://beta.uniprot.org/)获取青娥丸活性成分的靶点信息并进行标准化处理。

1.4 “药物-活性成分-靶点”网络构建将青娥丸活性成分和作用靶点信息导入Cytoscape 3.9.1软件[10],去除相关度不高的成分和靶点数据,得到“药物-活性成分-靶点”网络;对该网络进行拓扑分析,并按照度值(Degree)从大到小的顺序进行排序。

1.5 LDH疾病靶点获取以“lumbar disc herniation”为检索词,通过DisGeNET数据库(https://www.disgenet.org/)和GeneCards数据库(https://www.genecards.org/)获取LDH相关疾病靶点。Score值用于衡量靶点与疾病的相关性,在DisGeNET数据库中,选取Score值大于0的靶点;在GeneCards数据库中,选取Score值大于5的靶点。最后,将两个数据库得到的靶点进行合并去重整理。

1.6 基因本体(Gene Ontology,GO)和京都基因与基因组百科全书(Kyoto Encyclopedia of Genes and Genomes,KEGG)富集分析将“1.2”项的模块基因、“1.3”项的青娥丸作用靶点及“1.5”项的LDH疾病靶点取交集,作为青娥丸治疗LDH的关键靶点,并采用Venny 2.1.0(https://bioinfogp.cnb.csic.es/tools/venny/index.html)绘制韦恩图。将“1.3”项的青娥丸作用靶点与“1.5”项的LDH疾病靶点取交集,随后将交集基因与“1.2”项的模块基因进行合并去重整理,得到的靶点通过R软件中的ClusterProfiler包[11]进行GO功能富集和KEGG通路富集分析,根据P值进行升序排序,结果采用ggplot2包进行可视化展示。

1.7 分子对接将“1.4”项度值排序前5位的活性成分作为配体,“1.6”项的关键靶点作为受体。从TCMSP网站下载活性成分的mol2文件,并从PDB数据库(https://www.rcsb.org/)下载关键靶点的PDB文件。采用AutoDockTools 1.5.7软件[12]进行分子对接,PyMOL 2.4.1软件进行可视化展示。

2 结果

2.1 WGCNAGSE124272数据集共涉及5 495个基因。加权基因共表达网络见图1,图中尺度独立性(scale independence)越高代表该网络越接近于无尺度网络分布,通常选择纵轴高于0.9的数值[13];平均连接度(mean connectivity)越趋于平稳代表模块网络之间越紧密,因此选择最佳软阈值β=8。然后,对基因相关性大于0.9的模块进行合并,共确定9个模块,见图2。通过模块基因与性状的关联分析,发现blue模块(R=0.63,P=0.008)可能是LDH发生发展的关键模块,见图3。根据模块基因与表型数据的相关系数(MM)筛选模块内的重要基因,得到310个基因,见图4。

图1 软阈值确定图

图2 聚类分析模块图

图3 模块基因与性状的关联分析图

图4 模块基因与表型数据的关联分析图

2.2 青娥丸活性成分及作用靶点基于TCMSP和BATMAN-TCM数据库得到65个青娥丸活性成分,其中杜仲27个,补骨脂13个,核桃仁4个,大蒜21个。通过Uniprot数据库获取997个青娥丸活性成分的作用靶点。采用Cytoscape 3.9.1软件构建“药物-活性成分-靶点”网络,通过Network Analyzer插件对该网络进行拓扑分析,Degree值排名前5位的活性成分是:补骨脂素、鞣花酸、原人参二醇、异补骨脂素、槲皮素,这些成分可能是青娥丸发挥药理作用的关键成分,见图5。

注:DZ:杜仲;HTR:核桃仁;BGZ:补骨脂;DS:大蒜。

2.3 LDH疾病靶点通过DisGeNET和GeneCards数据库得到1 127个LDH相关疾病靶点。将310个blue模块基因、997个青娥丸作用靶点及1 127个LDH疾病靶点取交集,得到5个共同靶点,即丝裂原活化蛋白激酶14(mitogen activated protein kinase 14,MAPK14)、基质金属蛋白酶-9(matrix metalloproteinase-9,MMP-9)、白细胞介素-1β(interleukin-1β,IL-1β)、碱性磷酸酶(alkaline phosphatase,ALP)、活性素受体2A(activin receptor 2A,ACVR2A),这5个靶点可能是青娥丸治疗LDH的关键靶点,见图6。

图6 blue模块基因、青娥丸作用靶点及LDH疾病靶点韦恩图

2.4 GO功能富集分析和KEGG通路富集分析将997个青娥丸作用靶点与1 127个LDH疾病靶点取交集,得到123个共同靶点;将这些共同靶点与310个blue模块基因合并去重,最终筛选得到380个基因,将其导入R软件中,使用ClusterProfiler包进行GO和KEGG富集分析。GO富集分析得到 1 971 条生物过程(biological process,BP)、65条细胞成分(cellular component,CC)和120条分子功能(molecular function,MF),排名前8位的条目如图7所示。其中,BP主要涉及细胞对化学应激的反应、骨化、对营养水平的反应、神经元死亡调节、肌肉细胞增殖等;CC主要涉及分泌颗粒膜、分泌颗粒腔、细胞质囊泡腔、囊泡腔等;MF主要涉及受体配体活性、信号受体激活剂活性、细胞因子受体结合、生长因子结合等。KEGG富集分析共得到125条通路,主要涉及MAPK信号通路、破骨细胞分化、TNF信号通路、细胞凋亡、Toll样受体信号通路等,排名前20位的通路如图8所示。

图7 GO功能富集分析图

图8 KEGG通路富集分析图

2.5 分子对接采用AutoDock软件将青娥丸关键成分(补骨脂素、鞣花酸、原人参二醇、异补骨脂素、槲皮素)与筛选出的关键基因(MAPK14、MMP-9、IL-1β、ALP、ACVR2A)进行分子对接,对接结合能见图9。在分子对接过程中,结合能是判断分子间结合能力的重要指标之一,其数值越小表示分子之间结合越强。研究表明,当结合能小于0时,受体与配体可以自发结合,且数值越小自发结合的倾向越明显[14]。结果显示,青娥丸关键成分与关键基因的结合能均小于-4.0 kcal·mol-1,这表明他们具有较好的结合活性,其中MMP-9-补骨脂素结合能为-7.01 kcal·mol-1、MMP-9-异补骨脂素结合能为-6.62 kcal·mol-1、ACVR2A-原人参二醇结合能为-6.85 kcal·mol-1、ALP-鞣花酸结合能为-6.18 kcal·mol-1。选择5个关键成分与某个关键基因结合能最低的对接结果进行可视化展示,其中MMP-9-补骨脂素结合能力最强,并在ALA-417、TYR-420、LEU-418、MET-422处形成氢键,见图10。

图9 分子对接结合能热图

图10 分子对接结果可视化展示图

3 讨论

LDH引起的腰背部疼痛严重影响了患者的日常生活,而LDH的发病过程通常与椎间盘退变(intervertebral disc degeneration,IDD)密切相关。研究表明,IDD的患病率随着年龄增长而逐渐上升,40岁以上人群的患病率超过90%[15]。椎间盘组织的稳定性依赖于中央的髓核细胞(nucleus pulposus cells,NPCs)及其产生的细胞外基质(extracellular matrix,ECM)。IDD发生发展过程中,NPCs凋亡增加,炎症因子和基质金属蛋白酶等物质分泌增加,从而加速邻近细胞凋亡和ECM降解,导致髓核脱水,无法实现流体增压,从而引起椎间盘受力不均,周围纤维环的轴向压力增加,导致纤维环破裂引发LDH[16]。

目前,LDH的手术治疗存在感染、出血等并发症风险,而中医药凭借其便捷、经济、高效等优点,在延缓IDD进程方面具有独特优势。LDH归属于中医学“腰腿痛”“腰痹”“骨痹”等范畴,常见证型有气滞血瘀、肾阴不足等,治疗以活血化瘀、补肾强腰为主[17-18]。青娥丸是治疗肾虚腰痛的经典方剂,首载于宋代《太平惠民和剂局方》,方中补骨脂、杜仲、核桃仁补肾温阳,大蒜补肾止痛,四药配伍,共奏补肾益腰之效[4]。临床上,青娥丸治疗LDH疗效显著。

本研究基于WGCNA和分子对接探索青娥丸治疗LDH的分子机制。通过整合中药作用网络与疾病生物网络,旨在深入理解药物与机体的相互作用,并科学预测中药治疗疾病的有效机制[19]。通过数据挖掘得到青娥丸治疗LDH的5个关键靶点,包括MAPK14、MMP-9、IL-1β、ALP、ACVR2A。其中,MAPK14作为MAPK家族的成员,在细胞应激、炎症反应和细胞凋亡等生理病理过程中发挥重要作用,并参与IDD的发展[20]。MMP-9属于基质金属蛋白酶家族,参与细胞增殖、细胞迁移、血管生成、肿瘤侵袭和炎症反应等多种生理病理过程[21]。研究显示,MMP-9可以降解椎间盘基质的胶原蛋白和其他成分,从而促进IDD和LDH的进展[22]。此外,MMP-9可调节炎症反应和神经元损伤,从而加剧疼痛症状[23]。IL-1β作为一种促炎细胞因子,可通过激活MAPK14信号通路诱导MMP-9表达[24]。ACVR2A是一种重要的受体蛋白,能够抑制NF-κB信号通路激活,减少炎症因子产生,进而延缓LDH的发展[25]。ALP是一种膜结合型酶,参与骨代谢和骨发育,同时调节其他信号通路和细胞因子的产生,从而延缓LDH的进展[26]。因此,本研究初步揭示了青娥丸治疗LDH的分子机制。

本研究发现青娥丸的5个关键成分,其中补骨脂素和异补骨脂素是从中药补骨脂中提取的植物雌激素,他们与骨代谢、胰岛素分泌、葡萄糖代谢、体质量控制等生理过程密切相关。研究表明,补骨脂素能在一定程度上缓解IL-1β诱导的LDH;异补骨脂素可以调控p38MAPK信号通路,从而增强抗氧化和抗炎活性,达到改善LDH的作用[27-28]。鞣花酸是一种天然多酚化合物,具有抗氧化、抗炎、抗癌和抗菌等多种生物活性[29]。原人参二醇是一种类固醇化合物,具有抗肿瘤、抗炎和抗氧化等作用,可以抑制MAPK和NF-κB信号通路,从而减轻椎间盘退化程度,促进椎间盘细胞的生存与增殖[30]。槲皮素是一种天然黄酮类化合物,具有抗氧化、抗炎、抗癌、抗衰老等活性,可以延缓IDD进程,但其量效关系和安全性仍需进一步研究[31]。

综上所述,青娥丸可能通过补骨脂素、鞣花酸、原人参二醇等活性成分作用于MAPK14、MMP-9、IL-1β等关键靶点,调控MAPK信号通路、破骨细胞分化等通路从而发挥治疗LDH的作用。这凸显了青娥丸的多成分、多靶点、多途径等特点,为进一步开发复方或化合物治疗LDH提供理论依据。本研究结果尚需通过体内外实验进一步验证。