基于機(jī)器學(xué)習(xí)的成纖維細(xì)胞生長因子受體激酶抑制劑虛擬篩選模型

丁俊濤，劉博，吳建章，李物蘭

1.溫州醫(yī)科大學(xué) 第一臨床醫(yī)學(xué)院（信息與工程學(xué)院），浙江 溫州 325035；2.溫州醫(yī)科大學(xué)附屬眼視光醫(yī)院，浙江 溫州 325000

成纖維細(xì)胞生長因子受體（fibroblast growth factor receptors, FGFR）是一種在許多生物學(xué)過程中發(fā)揮重要作用的受體酪氨酸激酶，包括FGFR1、FGFR2、FGFR3、FGFR4四種亞型[1]。FGFR信號傳導(dǎo)的異常激活在不同類型腫瘤（如膽管癌、子宮內(nèi)膜癌、尿路上皮癌和肺癌等）的發(fā)生和進(jìn)展中起重要作用[2-4]，F(xiàn)GFR抑制劑具有治療這些疾病的潛力。目前FGFR激酶家族中已有3種小分子靶向抑制劑在近三年被批準(zhǔn)用于癌癥治療[5]，但是隨后發(fā)現(xiàn)它們均表現(xiàn)出較強(qiáng)的高血磷癥、腹瀉等不良反應(yīng)，因此急需尋找出新穎先導(dǎo)化合物骨架用于改構(gòu)成安全性更好的抑制劑。

在新藥發(fā)現(xiàn)領(lǐng)域，與實驗篩選相比，虛擬篩選方法在經(jīng)濟(jì)成本、時間效率上極具優(yōu)勢[6]。目前虛擬篩選可分為傳統(tǒng)經(jīng)典的計算機(jī)輔助藥物設(shè)計（computer aided drug design, CADD）[7]和現(xiàn)代新穎的人工智能藥物設(shè)計（artificial intelligence drug design, AIDD）[8]。在面對千萬級以上大數(shù)據(jù)庫的虛擬篩選方面，AIDD的效率遠(yuǎn)高于CADD。在已報道的FGFR抑制劑中，部分是基于CADD的新藥發(fā)現(xiàn)，未見基于AIDD的抑制劑研究。因此，本研究建立了基于AIDD的FGFR激酶抑制劑虛擬篩選模型，選擇FGFR四種受體中與腫瘤等疾病關(guān)系密切的FGFR1，將AIDD篩選得到的化合物用CADD進(jìn)一步進(jìn)行了FGFR1激酶抑制劑的虛擬篩選和分子動力學(xué)模擬研究，旨在為FGFR抑制劑的研究提供高效AIDD篩選模型和苗頭化合物。

1 材料和方法

1.1 FGFR激酶抑制劑數(shù)據(jù)整理由于FGFR1-4激酶結(jié)構(gòu)，尤其是FGFR1-3，高度相似，因此目前已報道的FGFR激酶抑制劑絕大部分是泛FGFR抑制劑。將BindingDB公共藥物實驗數(shù)據(jù)集[9]導(dǎo)入Mysql[10]軟件進(jìn)行整理查詢，獲得2196條FGFR激酶抑制劑實驗數(shù)據(jù)，并將半數(shù)抑制濃度小于100 nmol/L的記為活性數(shù)據(jù)，標(biāo)簽值為1，其余的記為非活性數(shù)據(jù)，標(biāo)簽值為0。在對數(shù)據(jù)進(jìn)行清洗、去重、標(biāo)簽等操作后，最終得到活性數(shù)據(jù)1275條，非活性數(shù)據(jù)921條，作為后續(xù)處理和訓(xùn)練的數(shù)據(jù)集合。

1.2 分子的特征工程采用兩種分子特征表示方法：①分子指紋MACCS[11]：包含166位分子指紋，指紋中的每一位都代表特定的子結(jié)構(gòu)，可使用RDKit 根據(jù)分子的簡化分子線性輸入規(guī)范（simplified molecular input line entry system, SMILES）[12]計算獲得；②分子描述符：使用RDKit根據(jù)分子SMILES計算獲得，內(nèi)容包括分子量、脂水分配系數(shù)、拓?fù)錁O性表面積和可旋轉(zhuǎn)鍵數(shù)等，從不同角度描述分子性質(zhì)，RDKit描述符合共有206維。

1.3 基于機(jī)器學(xué)習(xí)虛擬篩選模型的構(gòu)建采用隨機(jī)森林和支持向量機(jī)兩種機(jī)器學(xué)習(xí)方法建立虛擬篩選模型：①隨機(jī)森林：建立由決策樹組成的“森林”，利用多棵決策樹的結(jié)果對樣本進(jìn)行訓(xùn)練并預(yù)測的分類器，通過Sklearn[13]完成模型建立并對隨機(jī)森林超參數(shù)進(jìn)行調(diào)試和訓(xùn)練。②支持向量機(jī)：使用監(jiān)督學(xué)習(xí)模式對輸入數(shù)據(jù)進(jìn)行二元分類的線性分類器，其算法核心思想是求解輸入數(shù)據(jù)的最大邊距超平面以完成分類。使用Sklearn完成模型構(gòu)建并完成訓(xùn)練測試。隨機(jī)森林和支持向量機(jī)兩種模型的構(gòu)建都通過Python和Sklearn完成，模型的整體結(jié)構(gòu)包括化合物特征提取、數(shù)據(jù)載入、劃分測試集和驗證集、模型擬合訓(xùn)練、參數(shù)迭代調(diào)整以及最終的活性預(yù)測。

1.4 基于機(jī)器學(xué)習(xí)虛擬篩選模型的評價指標(biāo)使用準(zhǔn)確率、精準(zhǔn)率、召回率、曲線下面積（area under curve, AUC）4個指標(biāo)驗證機(jī)器學(xué)習(xí)模型在數(shù)據(jù)集上的綜合性能。對于二分類問題，預(yù)測結(jié)果有4種劃分，分別是：真陽（true positive, TP）、真陰（true negative, TN）、假陽（false positive,FP）、假陰（false negative, FN）。本研究中所用的指標(biāo)計算方法如下：

準(zhǔn)確率：所有預(yù)測正確（TP與TN）占總的比率，用于判斷模型預(yù)測分類是否準(zhǔn)確。

精準(zhǔn)率：預(yù)測為正且實際為正占全部預(yù)測為正的比率。

召回率：預(yù)測為正且實際為正占全部實際為正的比率。

AUC：對受試者工作特征曲線（receiver operating characteristic curve, ROC）進(jìn)行積分計算，用來定量描述分類器的好壞，AUC的值越大則分類性能越強(qiáng)。

1.5 基于分子對接的虛擬篩選在用機(jī)器學(xué)習(xí)模型篩選獲得潛在的活性化合物后，選擇FGFR家族中的FGFR1，進(jìn)一步采用兩級遞進(jìn)的基于分子對接的虛擬篩選，即依次用Autodock Vina軟件[14]和Glide 方法的XP（Extra Precision）精度進(jìn)行FGFR1激酶抑制劑的虛擬篩選。使用的模板蛋白結(jié)構(gòu)從PDB網(wǎng)站獲取，蛋白結(jié)構(gòu)的PDBID為4V05[15]。對蛋白結(jié)構(gòu)模型進(jìn)行加氫、能量優(yōu)化等操作后，以模型內(nèi)包含的AZD4547的結(jié)合位置為中心，選取60 ?×60 ?× 60 ?大小的格子作為對接的區(qū)域進(jìn)行虛擬篩選。

1.6 分子動力學(xué)模擬使用Gromacs2020.4軟件進(jìn)行分子動力學(xué)模擬，使用Charmm36力場[16]為復(fù)合體系的每個蛋白質(zhì)生成拓?fù)浜蛥?shù)文件，小分子配體使用CGenff力場參數(shù)。然后將體系溶解于TIP3P十二面體水盒中，距離體系表面10 ?處，加入適量的反離子使體系中和，接著采用最陡下降法進(jìn)行能量最小化50000步，然后進(jìn)行100 ps的NVT （Canonical Ensemble，恒定粒子數(shù)、體積和溫度）模擬和100 ps的NPT（Constant-pressure Constant-temperature，恒定粒子數(shù)、壓力和溫度）模擬，溫度限制在300 k，最后用MD模擬每個系統(tǒng)的NPT，持續(xù)時間為500 ns。所有仿真步驟均為 2 fs，軌跡每2 ps保存1次，以供后續(xù)分析。

2 結(jié)果

2.1 BindingDB中活性化合物理化性質(zhì)分析對BindingDB庫內(nèi)實驗數(shù)據(jù)進(jìn)行分析，庫內(nèi)化合物分子質(zhì)量的分布有一定的差異（見圖1A），活性化合物的分子量分布大多集中在442～557之間，且其中位數(shù)略高于非活性化合物，符合類藥性五原則對于分子量的要求。化合物的拓?fù)錁O性表面積 （topological polar surface area, tPSA）分布也有一定差異（見圖1B），活性化合物的tPSA較高，集中在89～117之間。另外，活性化合物的定量評估類藥性（quantitative estimate of drug-likeness,QED）多集中于0.5，而非活性化合物集中于0.4（見圖1C），說明活性化合物有更好的成藥性。最后，活性化合物結(jié)構(gòu)中成環(huán)多數(shù)為5個，而非活性化合物則為4個（見圖1D）。環(huán)的數(shù)量也可以導(dǎo)致化合物親疏水性質(zhì)和與靶點親和力的變化。

圖1 BindingDB數(shù)據(jù)庫FGFR1實驗數(shù)據(jù)理化性質(zhì)統(tǒng)計分析

2.2 模型評價通常認(rèn)為，準(zhǔn)確率和AUC值大于0.75，分類器有較好的準(zhǔn)確性和分類效果。在FGFR激酶抑制劑的模型中，隨機(jī)森林相比于支持向量機(jī)在準(zhǔn)確率、精準(zhǔn)率、召回率有更好的表現(xiàn)（見表1）。從ROC曲線來看，隨機(jī)森林拐點更靠近左上角（見圖2），說明隨機(jī)森林算法有更好的預(yù)測性能。

表1 機(jī)器學(xué)習(xí)模型評價指標(biāo)

圖2 機(jī)器學(xué)習(xí)模型的ROC曲線

2.3 虛擬篩選結(jié)果虛擬篩選的化合物庫來自陶素公司提供的包含1300萬個小分子的虛擬化合物庫。針對該庫使用隨機(jī)森林方法進(jìn)行第一級的虛擬篩選，將預(yù)測陽性概率為0.999以上的化合物認(rèn)為是潛在的活性化合物，共計篩選獲得10340個潛在的活性化合物，約占總數(shù)的0.7%，耗時約28.3 h。 隨后對10340 個化合物用Autodock Vina軟件和Glide方法逐級篩選FGFR1激酶抑制劑。用Autodock Vina軟件得到打分值小于-9.00 kcal/mol的分子共395個，占比3.8%。用Glide法從395個化合物中篩選得到3個打分小于-9.80 kcal/mol的化合物，占比0.76%。其對接打分值如表2所示，最優(yōu)化合物的打分結(jié)果為-11.12 kcal/mol。AZD4547（見圖3A）、化合物a（見圖3B）、化合物b（見圖3C）均與FGFR1激酶ALA-101形成氫鍵作用，分子都伸入口袋內(nèi)部疏水口袋，化合物a、化合物b以及AZD4547與FGFR1結(jié)合模式相同[17]。化合物a和b與陽性對照AZD4547一樣與FGFR1激酶蛋白有著較多的氫鍵相互作用，不同的是，化合物c僅由疏水作用維持結(jié)合（見圖3D），沒有與FGFR1形成氫鍵相互作用，也沒有鹵鍵等相互作用，導(dǎo)致其與FGFR1的結(jié)合效果較差。此 外，化合物a上Cl和F分別和殘基LYS-19、VAL-29形成鹵鍵，推測與FGFR1的結(jié)合力較AZD4547更好。

表2 隨機(jī)森林、Vina和Glide XP虛擬篩選結(jié)果

圖3 FGFR1激酶（PDBID為4V05）與AZD4547（A）、化合物a（B）、化合物b（C）、化合物c（D）的結(jié)合模式

2.4 分子動力學(xué)模擬結(jié)果將模板蛋白的配體AZD4547 與打分最高的3 個潛在活性化合物進(jìn)行 100 ns的分子動力學(xué)模擬分析。通過MM/GBSA方法對4個體系的較為穩(wěn)定的時間段70～80 ns進(jìn)行結(jié)合自由能計算（見圖4、表3），其中AZD4547結(jié)合自由能最優(yōu)，總的結(jié)合自由能為-46.72 kcal/mol，化合物a總的結(jié)合自由能為-38.13 kcal/mol，化合物b總的結(jié)合自由能為-42.87 kcal/mol，化合物c總的結(jié)合自由能為-13.68 kcal/mol，結(jié)合圖3中4個化合物的結(jié)合模式，說明氫鍵和鹵鍵是提升化合物與FGFR1親和力的關(guān)鍵，化合物d與FGFR1沒有氫鍵作用，所以結(jié)合自由能較高。對比結(jié)合口袋位置氨基酸殘基能量貢獻(xiàn)情況（見圖5），其中LEU21、VAL29、ALA49這3個殘基在4個體系中都貢獻(xiàn)了較多的結(jié)合自由能，從結(jié)合自由能角度解釋了化合物a、b、c的結(jié)合情況與AZD4547較為相似。但是化合物c與其他化合物不同的是，ASP178削弱了結(jié)合自由能，這也是化合物c結(jié)合自由能較差的一個關(guān)鍵因素。對比4組體系的均方根偏差（root mean square deviation,RMSD）（見圖6），4個小分子化合物在口袋內(nèi)較為穩(wěn)定，其中化合物a較AZD4547有更好的穩(wěn)定性，其RMSD波動范圍較小；化合物b與AZD4547相比波動性較大，但波動范圍均未超過0.4 ?；化合物c穩(wěn)定性最差，它們的RMSD結(jié)果與前期分子對接的結(jié)果相吻合。

表3 FGFR1與4種化合物結(jié)合的自由能情況（kcal/mol）

圖5 FGFR1激酶化合物復(fù)合體系在70～80 ns時氨基酸貢獻(xiàn)自由能統(tǒng)計

圖6 FGFR1激酶化合物體系分子動力學(xué)模擬的RMSD

3 討論

近年由于人工智能的蓬勃發(fā)展，越來越多的機(jī)器學(xué)習(xí)和深度學(xué)習(xí)技術(shù)被應(yīng)用于藥物發(fā)現(xiàn)。在激酶領(lǐng)域中，受體酪氨酸激酶家族已有較為成功的人工智能結(jié)合發(fā)現(xiàn)藥物案例。在人工智能與新藥發(fā)現(xiàn) 中，主要困難在于藥物活性數(shù)據(jù)與陰性數(shù)據(jù)的缺乏，目前較為成功的案例如血管內(nèi)皮生長因子受體、表皮生長因子受體往往是具備較多實驗測試數(shù)據(jù)的靶點[18-20]。而本研究的FGFR抑制劑目前數(shù)據(jù)相對缺少，現(xiàn)有數(shù)據(jù)庫中活性數(shù)值模糊，對該領(lǐng)域機(jī)器學(xué)習(xí)發(fā)展有較大掣制，且尚未見有利用人工智能發(fā)現(xiàn)FGFR抑制劑的相關(guān)報道。因此本研究在數(shù)據(jù)準(zhǔn)備時，對已有藥物數(shù)據(jù)庫和文獻(xiàn)抑制劑數(shù)據(jù)進(jìn)行整理，并將其用于人工智能模型構(gòu)建和藥物篩選。本研究中使用的隨機(jī)森林和支持向量機(jī)兩種機(jī)器學(xué)習(xí)模型的準(zhǔn)確率分別為0.878和0.770，AUC分別是0.952和0.840，這表明本研究構(gòu)建的模型較好。與支持向量機(jī)算法構(gòu)建的模型相比，隨機(jī)森林構(gòu)建的模型在準(zhǔn)確率、精準(zhǔn)率、召回率方面總體更優(yōu)。同時機(jī)器學(xué)習(xí)模型從1300萬個小分子化合物庫篩選獲得10340個潛在的活性化合物，耗時約28.3 h，這說明機(jī)器學(xué)習(xí)模型在面對較大虛擬篩選藥物庫時可以快速準(zhǔn)確地從中篩選出可能的先導(dǎo)化合物，快速縮小篩選范圍，節(jié)約時間和經(jīng)濟(jì)成本，提高效率。

分子動力學(xué)模擬是一種有效分析小分子化合物與激酶結(jié)合模式和結(jié)合能力的方法。通過分子動力學(xué)模擬并計算激酶和化合物的結(jié)合自由能，a和b的圖6 FGFR1激酶化合物體系分子動力學(xué)模擬的RMSD結(jié)合自由能與AZD4547相近，在理論上說明化合物a、b可能有接近AZD4547的抑制活性。將結(jié)合自由能分解，最優(yōu)的3個化合物的結(jié)合自由能主要貢獻(xiàn)來源于范德瓦爾斯力，其次是靜電力，由于化合物a的靜電力貢獻(xiàn)較小，導(dǎo)致其最終的自由能情況較差。這說明對化合物的改造可以通過增強(qiáng)靜電力的貢獻(xiàn)和維持疏水作用來提升化合物抑制活性。對化合物周圍的氨基酸殘基能量貢獻(xiàn)進(jìn)行分析，在整個過程中，LEU21、ALA101兩個殘基對結(jié)合能的貢獻(xiàn)在幾個化合物體系中均有呈現(xiàn)，說明LEU21、ALA101兩個殘基是主要影響該系列先導(dǎo)化合物和FGFR1結(jié)合的關(guān)鍵殘基，在化合物b和FGFR1的復(fù)合體系中，GLU108和ASP178兩個殘基對結(jié)合能的貢獻(xiàn)起到反作用，這說明這兩個殘基可能阻礙化合物和FGFR1結(jié)合，因此化合物b的結(jié)合能較其他3個化合物相比更高。上述結(jié)果表示該模型的虛擬篩選性能較好。

總之，本研究提供了一種基于機(jī)器學(xué)習(xí)的FGFR激酶抑制劑虛擬篩選模型，其可以短時間內(nèi)高效篩選規(guī)模較大的小分子庫，快速高效地獲得潛在活性的FGFR激酶抑制劑。