基于L1000数据集筛选骨质疏松症的潜在治疗化合物

梁鹏晨,周紫艳,常庆,易清清,孙苗苗,唐晔翎,曹励欧,杨洁

(1.上海大学微电子学院,上海 201800;2.上海健康医学院附属嘉定区中心医院临床科研中心,上海 201800;3.上海中医药大学研究生院,上海 201203;4.上海健康医学院附属嘉定区中心医院肾内科,上海 201800;5.上海交通大学医学院附属仁济医院肾内科,上海200127;6.上海健康医学院附属嘉定区中心医院GCP办公室,上海 201800)

骨质疏松症是人类常见的骨骼疾病,其特征是骨量低,骨组织恶化,骨结构破坏,骨强度降低,导致骨折风险增加。根据世界卫生组织的诊断分类,骨质疏松症的定义是髋部或腰椎的骨密度低于青年参考人群平均骨密度的2.5个标准差。骨重建是在去除旧骨(骨吸收)和形成新骨(骨形成)之间的平衡过程,它分别由破骨细胞和成骨细胞两种不同的细胞完成[1]。当这两个过程脱钩时,骨吸收超过骨形成,从而导致骨质疏松[2]。骨质疏松症的治疗药物主要分为两类,一类是减缓骨吸收的抗吸收药物,另一类是刺激骨形成的合成代谢药物。目前临床上用于治疗骨质疏松症的药物有双膦酸盐、核因子κB受体活化因子配体(RANKL)抑制性抗体、甲状旁腺素类似物和抗硬化素抗体等[3]。虽然这些药物是有效的,但大多数都有一些局限性和不良反应。考虑到这一点,本研究利用生物信息学的方法筛选骨质疏松症的潜在治疗化合物。

L1000是一种新的低成本和高通量的基因表达谱分析方法,用于评估由微阵列产生的全基因组基因表达特征[4]。基于L1000分析方法,集成网络细胞特征库项目团队开发了新一代的“连接图”(connectivity map,CMap)数据库。新一代的CMap数据库存储了百万多个、细胞被“扰动”处理后的基因表达谱,这些“扰动”包括小分子化合物、基因过表达(cDNA)、基因敲低(shRNA)等干扰[4]。这一庞大的CMap数据库为理解未知化合物的作用机制和寻找现有药物的新适应证提供了巨大的机会[5]。

本研究利用上传的骨质疏松症差异基因表达谱,在CMap数据库的L1000数据集中,筛选相应潜在的治疗骨质疏松症的化合物。并进一步通过分子对接,ADMET(吸收、分布、代谢、排泄、毒性)计算分析,细胞增殖活性实验,评价这些化合物的作用靶点和药理学特性,以期为临床提供安全有效的骨质疏松症的潜在治疗药物。

1 材料与方法

1.1 材料

Gene Expression Omnibus(GEO)数据库(https://www.ncbi.nlm.nih.gov/geo/),R.3.6.1(https://www.r-project.org/),pkCSM数据库(http://biosig.unimelb.edu.au/pkcsm/),TTD治疗靶点数据库(http://db.idrblab.net/ttd/),R studio(https://rstudio.com/),R语言GEOquery包,R语言limma包,R语言Pheatmap包,MetaScape数据库(https://metascape.org),CMap数据库(https://clue.io/),CMap Query匹配工具(https://clue.io/query),PDB数据库(http://www.rcsb.org/),Pubchem数据库(https://pubchem.ncbi.nlm.nih.gov/),AutoDock Tools(version 1.5.6,http://autodock.scripps.edu/resources/adt),Open Babel 软件(http://openbabel.org/),AutoDock Vina软件(http://vina.scripps.edu/)。

1.2 筛选骨质疏松症的差异表达基因及可视化

在GEO数据库中,获得编号为GSE35958的骨质疏松症转录组数据集以及探针平台文件GPL570,该数据集包含9个人骨髓间充质干细胞(human mesenchymal stem cells,hMSCs)样本的转录组数据,其中骨质疏松症老年患者的hMSCs的转录组数据有5组;非骨质疏松老年供体的骨髓hMSCs的转录组数据有4组。

利用R语言的GEOquery包,对GSE35958矩阵数据文件进行下载及整理;利用hgu133plus2.db包,将数据集的探针名转化成基因名;通过比较数据集各组基因与管家基因(GAPDH、ACTB)的表达量平均值差异,判断数据集的合理性。利用limma包内置函数normalizeBetweenArrays,对数据进行批次矫正。利用limma包对数据集进行差异分析(其中非骨质疏松样本作为对照组,骨质疏松样本作为实验组)。差异表达基因的筛选标准为log2(差异倍数)的绝对值>2.5,P<0.05。利用R语言内置函数绘制差异表达基因的可视化火山图。

1.3 差异表达基因的通路富集分析

将差异表达基因上传到MetaScape数据库[5]中,分析参数设置为Min Overlap=3、PValue Cutoff=0.01、Min Enrichment=1.5;进行GO生物过程、KEGG 通路、Reactome 基因集定制化通路分析,研究基因参与的相应细胞通路。

1.4 CMap化合物匹配

CMap是一个化学基因组学数据库,它收集用小分子化合物或其他“扰动物质”处理培养的人类细胞的基因表达谱,可以用来筛选逆转特定基因表达谱的小分子化合物[4]。通过新一代CMap数据库中的Query工具,将骨质疏松症的差异表达基因上传到数据库中,查询类型选择Gene expression(L1000)。匹配化合物通过CMap得分数值([+100,-100])反映其与上传基因表达谱的相似或相反关系。CMap得分数值为负数的化合物,显示其对上传疾病基因表达谱的负相关性,提示其潜在的治疗骨质疏松症作用。本研究中,CMap 得分小于-90的匹配化合物,视为潜在治疗化合物。

1.5 潜在治疗化合物的分子对接研究

检索文献及TTD治疗靶点数据库得到与骨质疏松症相关的靶蛋白,利用PDB数据库查询得到蛋白晶体结构,将靶蛋白的三维结构信息及对应原配体信息,整理成骨代谢相关靶蛋白库(表1)[6-7]。通过Pubchem数据库下载化合物3D结构的sdf格式文件。利用AutoDockTools 1.5.6 软件对靶蛋白进行除水分子及加氢准备,并分离蛋白中的原配体结构,将原配体所在位置作为对接活性口袋。利用Open Babel软件,将化合物sdf格式文件转换成对接所需的pdbqt格式。利用AutoDock Vina软件进行蛋白-小分子对接,获得结合自由能(affinity)值,该值越低说明对接越稳定[8]。将蛋白与其相应原配体对接的affinity值,作为判断对接活性的阈值。蛋白与原配体对接affinity值减去蛋白与化合物对接affinity值,作为相对得分。相对得分越高说明化合物和该靶蛋白的对接活性越好,该靶蛋白可能是其发挥作用的潜在靶点。利用R语言的Pheatmap包,绘制化合物与各靶点对接的相对得分热图。

表1 骨代谢相关靶蛋白库

1.6 潜在治疗化合物的ADME及毒性预测

用pkCSM数据库[9]的ADME 机器学习预测模块计算所选化合物的吸收、分布、代谢、排泄特征,包括水溶性、Caco-2通透性、人体肠道吸收、人体分布体积。利用pkCSM数据库的毒性预测模块,计算化合物的致突变性、最大耐受剂量、半数致死量、大鼠口服慢性毒性、肝毒性。

1.7 MC3T3-E1细胞增殖活力评价

选取传代培养后处于对数生长期的小鼠颅骨前成骨细胞(MC3T3-E1)经0.25%的胰蛋白酶消化,将细胞密度调整为5×103/mL。将化合物美西林、己烯雌酚用培养液配置成0、10、20、40、80和160 μg/mL的药液备用。将200 μL含有细胞的培养液接种在96孔板中(9组,每组3个复孔),待细胞贴壁后吸除上清液,分别加入不同浓度的药液。置恒温箱内培养1 d后弃去上清液,PBS冲洗3次,按比例放入CCK8与无血清培养基的混合液。继续培养1.5 h后吸取上清液100 μL,放入96孔细胞培养板中,在酶联免疫测试仪上测定450 nm以下的光密度(D)值,将0 μg/mL组作为对照孔(n=3)。

1.8 统计学分析

应用Graphpad Prism 7.0软件进行统计分析,数据用均值±标准差表示,多组间均数比较采用单因素方差分析,进一步两两比较采用Tukey检验,P<0.05为差异有统计学意义。

2 结果

2.1 骨质疏松症的转录组数据分析结果

在GEO数据集中获得骨质疏松症相关的表达谱数据集(5个骨质疏松组和4个健康对照组的hMSCs样本),利用R语言的limma包分析数据集,共获得195个与骨质疏松症相关的差异表达基因(图1)。其中骨质疏松组中表达量显著上调的基因有127个,差异显著性前10个基因分别是GNAO1、FKBP8、AP2A1、CBX8、CD74、FSCN1、PCDHGA3、GOLGA6L2、IL32、ADAMTSL4;表达量显著下调的基因有68个,差异显著性前10个基因分别是TMEM154、PPP3R1、B4GALT1、C20orf197、PCDHB5、EPHA3、CTSZ、COL14A1、ANKRD1、ZIC1。

蓝色圆点代表显著下调基因,红色圆点代表显著上调基因,黑色圆点代表无差异变化的基因图1 骨质疏松症差异表达基因的火山图

2.2 差异表达基因的通路富集分析结果

将195个骨质疏松症相关差异表达基因,通过MetaScape数据平台进行通路富集分析,富集显著的部分结果见图2。差异表达基因的GO 生物过程共有146个,涉及骨化、骨矿化、生物矿化组织发育、胶原生物合成过程的正调控、骨骼重塑、调节胶原蛋白代谢过程、骨吸收、血管发育、血管形态发生、骨髓白细胞细胞因子的产生、凋亡过程的积极调控、细胞外基质组织等。差异表达基因的KEGG信号通路共有9条,涉及突触小泡循环通路、内分泌和其他因子调节的钙重吸收通路。差异表达基因的Reactome 基因集共有9个,涉及细胞外基质组织、细胞外基质的降解等。

图2 差异基因的通路富集结果

2.3 CMap数据库匹配骨质疏松症潜在治疗化合物

CMap数据库通过计算上传基因表达谱的相似性,将化合物与疾病或生理表型联系起来。通过CMap分析,本研究搜索了与骨质疏松症基因表达谱呈负相关的化合物,CMap 得分小于-90的潜在治疗化合物共有10个(表2),其中与骨质疏松症基因表达谱匹配负相关性最高的己烯雌酚是人工合成的非甾体雌激素类药物。

表2 CMap匹配的骨质疏松症潜在治疗化合物

2.4 潜在治疗化合物作用的骨质疏松症靶点

将潜在治疗化合物与骨质疏松症相关靶点进行分子对接研究,通过比较化合物与各靶蛋白对接的相对得分(图3),预测化合物可能的骨质疏松症靶点。己烯雌酚的对接活性较好的靶点有CAⅡ、PGR、ERβ等;PLX-4720的对接活性较好的靶点有CAⅡ、PGR、JNK1等;美西林的对接活性较好的靶点有CAⅡ、PGR、TAK1等;KU-C103428N的对接活性较好的靶点有ERβ、FGFR1、AKT2等;土大黄苷的对接活性较好的靶点有P38δ、BMP2、TAK1等;马拉韦罗的对接活性较好的靶点有P38δ、ERβ、AKT1等;奎宁的对接活性较好的靶点有CAⅡ、PGR、TAK1等;MST-312的对接活性较好的靶点有CK、P38δ、BMP2等；利培酮的对接活性较好的靶点有AKT2、AKT1、P38δ等;CO-102862的对接活性较好的靶点有CAⅡ、PGR、ERβ等。

2.5 潜在治疗化合物的ADMET预测结果

潜在治疗化合物的吸收、代谢、分布和排泄特征预测结果见表3。水溶性表示该分子在25℃水中的溶解度,以log(mol/L)表示,该预测值越大,水溶解性越好;预测结果表明美西林、MST-312水溶解性能最好。细胞的Caco-2单层被广泛用作人肠黏膜的体外模型,以预测口服药物的吸收情况,结果表明,利培酮、奎宁、己烯雌酚、CO-102862具有较好的Caco-2通透性。人肠道吸收模型预测通过人体小肠吸收的化合物比例,预测值<30%的分子被认为吸收不好,预测结果表明所有化合物的人肠道吸收性较好,其中美西林、利培酮最为优异。人分布体积模型是药物总剂量需要均匀分布以获得与血浆中相同浓度的理论体积,分布体积[以log(L/kg)表示]越高,药物在组织中的分布越多,预测值<-0.15,则认为分布体积低,预测值>0.45,则认为分布体积高,预测结果表明奎宁、马拉韦罗、利培酮具有较好的分布体积,在组织中的分布多。总清除率[以log(mL/min/kg)表示]由比例常数CLtot衡量,主要是肝清除率和肾清除率的组合,它与生物利用度有关,对于确定达到稳态浓度的给药速率很重要。

表3 潜在治疗化合物的ADME性能

热图中方块的颜色越红,代表对应的分子和蛋白的对接自由能越低,对接越稳定,该蛋白可能是其潜在治疗靶点图3 潜在治疗化合物与骨质疏松症靶点的分子对接热图

通过pkCSM数据库的Toxicity模块对潜在治疗化合物的不同毒性指标进行预测,结果见表4。致突变性预测一种给定的化合物是否具有潜在致突变性。最大耐受剂量[以log(mg/kg/day)表示]提供了人体内化学品毒性剂量阈值的估计值。半数致死量(mol/kg)用于评估不同分子的相对毒性,是导致一组实验动物中50%死亡的给药剂量。大鼠口服慢性毒性[以log(mg/kg_bw/day)表示]旨在确定导致观察到不良作用时化合物的最低剂量。肝毒性预测给定化合物是否与肝功能受损相关。pkCSM数据库的毒性预测结果提示KU-C103428N 具有致突变性和肝毒性。最大耐受剂量较低的化合物有马拉韦罗、利培酮、奎宁。半数致死量较低的化合物有CO-102862、MST-312。大鼠口服慢性毒性较低的是马拉韦罗。根据以上预测结果,美西林具有较好的ADME性能和较低的毒性,综合评价较好。

表4 潜在治疗化合物的毒性预测

2.6 MC3T3-E细胞增殖活力测试结果

根据上述的综合预测表明,美西林有较好的ADME性能和较低的毒性,综合评价较好;因此进一步通过CCK-8实验,对美西林作用于MC3T3-E1细胞的活性进行研究。由于己烯雌酚有较多治疗骨质疏松症的实验报道,因此也对其进行了细胞活性研究。结果如图4所示,10、20、40、80和160 μg/mL美西林均能促进成骨细胞的增殖,其中20 μg/mL的美西林对成骨细胞增殖的促进最为显著;而己烯雌酚从40 μg/mL开始就显著抑制成骨细胞的增殖,提示该浓度开始的己烯雌酚对成骨细胞有一定的毒性。

与0 μg/mL组相比,*:P<0.05;Δ:P<0.01;#:P<0.001图4 美西林和己烯雌酚对MC3T3-E1细胞增殖的影响

3 讨论

本研究通过在CMap数据库中提交骨质疏松症的显著差异基因数据,在其中的L1000数据集匹配潜在的治疗骨质疏松症的化合物。其中CMap得分小于-90(负相关)的化合物共有10个:己烯雌酚、PLX-4720、美西林、KU-C103428N、土大黄苷、马拉韦罗、奎宁、MST-312、利培酮和CO-102862,这些化合物对骨质疏松症具有潜在的治疗作用。其中与骨质疏松症基因表达谱匹配负相关性最高的己烯雌酚是人工合成的非甾体雌激素类药物,而雌激素补充疗法是防治绝经后骨质疏松症的有效措施[10-11]。这说明本研究中利用骨质疏松症基因表达谱,在CMap数据库筛选骨质疏松症的潜在治疗药物的策略是合理的。李朝阳等[12]研究了己烯雌酚对去卵巢大鼠骨质疏松症的影响与作用机制,结果表明己烯雌酚能显著抑制去卵巢引起的骨高转换,使治疗组的骨量明显高于去卵巢组。许碧连等[13]研究了己烯雌酚对去卵巢大鼠不同部位骨骼的影响,结果表明己烯雌酚能有效预防去卵巢大鼠的骨丢失,但是对不同部位骨骼的作用不同。匹配负相关性最高的己烯雌酚能抑制骨转换,增加骨质疏松症模型大鼠的骨量以及预防骨丢失,这证实了我们筛选的可靠性。

Han等[14]从GEO数据库中收集了脓毒症相关基因微阵列数据,并通过L1000数据集进行了候选药物匹配,使用活化的巨噬细胞系分析所选候选药物的抗炎作用,发现CGP-60474可能是减弱内毒素血症过程的潜在治疗候选物。Sendama[15]通过筛选L1000数据库,发现鲁索替尼(Ruxolitinib)是治疗新冠肺炎的合理抗病毒候选药物。Zhang等[16]通过L1000数据库,筛选出可逆转糖尿病肾病中基因特征的候选药物,并对其中的候选药物PLK1抑制剂(BI-2536)进行了生物活性研究;实验结果表明BI-2536在体外抑制了系膜细胞增殖和细胞外基质,并改善了糖尿病肾病模型小鼠的蛋白尿和肾损伤。因此,在药物开发的早期阶段,基于L1000微干扰数据集的虚拟筛选策略是一种节约时间和成本的有效方式。

本研究预测所筛选化合物主要潜在的骨质疏松症相关靶点,发现各化合物具有多靶点的特性,与化合物对接活性较好的主要有CAⅡ、PGR、ERβ、BMP2等骨质疏松症的靶点。pkCSM数据库预测潜在治疗化合物的ADMET性质,结果显示美西林具有较好的ADME性能和较低的毒性,综合评价较好;成骨细胞增殖实验结果表明,美西林能促进成骨细胞的增殖,其中20 μg/mL时对成骨细胞增殖的促进最为显著。后续我们将进一步在去卵巢骨质疏松大鼠模型上研究美西林治疗骨质疏松症的疗效。