基于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時對成骨細胞增殖的促進最為顯著。后續我們將進一步在去卵巢骨質疏松大鼠模型上研究美西林治療骨質疏松癥的療效。