基于機器學習的成纖維細胞生長因子受體激酶抑制劑虛擬篩選模型

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

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

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

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

1 材料和方法

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

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

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

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

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

精準率：預測為正且實際為正占全部預測為正的比率。

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

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

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

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

2 結果

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

圖1 BindingDB數據庫FGFR1實驗數據理化性質統計分析

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

表1 機器學習模型評價指標

圖2 機器學習模型的ROC曲線

2.3 虛擬篩選結果虛擬篩選的化合物庫來自陶素公司提供的包含1300萬個小分子的虛擬化合物庫。針對該庫使用隨機森林方法進行第一級的虛擬篩選，將預測陽性概率為0.999以上的化合物認為是潛在的活性化合物，共計篩選獲得10340個潛在的活性化合物，約占總數的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所示，最優化合物的打分結果為-11.12 kcal/mol。AZD4547（見圖3A）、化合物a（見圖3B）、化合物b（見圖3C）均與FGFR1激酶ALA-101形成氫鍵作用，分子都伸入口袋內部疏水口袋，化合物a、化合物b以及AZD4547與FGFR1結合模式相同[17]?；衔颽和b與陽性對照AZD4547一樣與FGFR1激酶蛋白有著較多的氫鍵相互作用，不同的是，化合物c僅由疏水作用維持結合（見圖3D），沒有與FGFR1形成氫鍵相互作用，也沒有鹵鍵等相互作用，導致其與FGFR1的結合效果較差。此 外，化合物a上Cl和F分別和殘基LYS-19、VAL-29形成鹵鍵，推測與FGFR1的結合力較AZD4547更好。

表2 隨機森林、Vina和Glide XP虛擬篩選結果

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

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

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

圖5 FGFR1激酶化合物復合體系在70～80 ns時氨基酸貢獻自由能統計

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

3 討論

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

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

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