面向虛擬篩選的GPU加速的分子對接方法

胡海峰,王領悅,唐詩迪,胡鳴珂,吳建盛

(1.南京郵電大學 通信與信息工程學院,南京 210023;2.南京郵電大學 地理與生物信息學院,南京 210023;3.南京林業大學 經濟管理學院,南京 210037)

虛擬篩選作為計算機輔助藥物設計的實用技術,在現代藥物研發中起著重要作用[1]。虛擬篩選在計算平臺上模擬藥物篩選的過程,可快速從幾十至上百萬的小分子化合物(配體)中,篩選出有可能和受體結合的活性化合物,大大降低實際篩選小分子化合物數目,縮短藥物研發周期,降低藥物研發的成本[2-3]。虛擬篩選的方法可分為基于結構的虛擬篩選和基于配體的虛擬篩選[4]。基于結構的虛擬篩選利用分子對接技術,在已知蛋白質受體的三維結構的基礎上,在結合位點處自動匹配化合物數據庫中的小分子配體,用打分函數對可能的結合模式進行能量計算,最終得到小分子化合物結合活性的排名[5-7]。而基于配體的虛擬篩選方法不需要知道目標蛋白質受體的結構,代表性的方法包括藥效團模型、定量構效關系和結構相似性等方法[8-9]。相對于基于配體的虛擬篩選,基于結構的虛擬篩選能避免活性化合物結構微小變化引起的活性改變[10],對小分子配體的活性值預測更加準確。

常用的分子對接軟件包括AutoDock4[11]、AutoDock Vina[12]、AutoDock Vina 1.2.0[13]和Smina等[14]可以進行基于結構的虛擬篩選,其中,AutoDock Vina提高了結合模式預測的平均準確度,通過使用更簡單的打分函數加快了搜索速度,是應用最廣泛的分子對接軟件。這些分子對接軟件分別采用不同的搜索算法和打分函數進行分子結合活性預測,需要耗費大量時間和計算資源來完成分子對接計算,然而在實際虛擬篩選時,篩選的候選化合物越多,可以更大可能的搜索到合適小分子配體,提高先導化合物質量[15]。因而面對大規模分子對接時,過長的篩選時間制約了分子對接軟件的應用。例如AutoDock Vina在常規計算平臺上篩選10億個化合物的數據庫時,每個配體的平均對接時間為15 s,篩選所有化合物的時間大概需要476年,目前的分子對接軟件無法滿足現代藥物研發的需求[15]。

近年來研究人員針對分子對接軟件對接時間過長的缺點進行了兩方面改進。一方面,提出基于AutoDock Vina改進版本的Qvina1和Qvina2分子對接軟件[16-17],Qvina1提出了啟發式算法減少不必要的搜索次數,改進了局部搜索算法,與AutoDock Vina相比提高了運行速度,但對接精度方面有所下降。在此基礎上,Qvina2優化了啟發式算法,提高了對接精度,與AutoDock Vina相比實現了20.49倍的最大加速,是已知AutoDock Vina系列中最為高效的對接軟件。另一方面,虛擬篩選中的分子對接計算特別適合使用并行計算加速,GPU(Graphic processing unit)因其強大的并行處理硬件架構特別適合應用于分子對接軟件[18-19]。GPU具有數千個計算核心,可提供強大的計算性能,性價比較高且易于開發,并且有完善的開發標準,例如并行計算架構(CUDA)和開放運算語言(OpenCL)[20]。近年來,考慮到分子對接軟件的特點,基于GPU的異構體系也被應用于改進分子對接軟件[21],例如實現了AutoDock4的加速版本AutoDock4-GPU[22],在細粒度上并行了拉馬克遺傳算法,取得了很好的加速比。但AutoDock4的搜索效率很低[13],制約了AutoDock4-GPU的加速性能。

本文提出一種基于GPU的QVina 2并行化方法Qvina2-GPU,在目前AutoDock Vina系列中最為高效的Qvina2基礎上,利用GPU硬件高度并行體系加速分子對接過程。具體包括增加初始化分子構象數量,大量線程分別對應不同的初始構象,使得每個線程獨立承擔蒙特卡羅搜索算法子任務,獨立搜索最佳的分子構象,即通過增加蒙特卡羅的迭代搜索的廣度以減少每次蒙特卡羅迭代搜索深度,并利用Wolfe-Powell準則改進局部搜索算法,相對于QVina 2中基于Armijo準則搜索算法,提高了對接精度,進一步減少蒙特卡羅迭代搜索深度,從而顯著減少了蒙特卡羅構象搜索時間,減少了分子對接軟件的整體運行時間。 本文的貢獻歸納如下:

1)提出基于GPU的QVina 2并行化方法Qvina2-GPU,設計QVina2-GPU異構并行架構,增加初始化分子構象數量,擴展蒙特卡羅的迭代局部搜索中線程的并行規模,每個線程獨立承擔蒙特卡羅搜索算法子任務,通過增加蒙特卡羅迭代搜索廣度以減少每次蒙特卡羅迭代搜索深度。

2)提出基于Wolfe-Powell準則的局部搜索算法,通過迭代步長的優化,改進了QVina 2中基于Armijo準則搜索算法,提高了分子對接精度,進一步減少蒙特卡羅迭代搜索深度。

3)使用公開數據庫中140個小分子復合物[22]作為測試數據集,實驗結果驗證了在保證對接精度的條件下,我們提出的QVina2-GPU對Qvina2的平均加速比達到5.18倍,最大加速比達到12.28倍,具有很大的實際應用前景。

1 QVina-GPU異構并行架構

QVina2-GPU異構并行架構如圖1所示,從圖中可以看出,QVina2-GPU整體架構由主機端和設備端組成,主機端由數據準備模塊和優化輸出模塊組成,這部分軟件運行在CPU平臺上,主要完成環境配置、數據輸入以及篩選最優分子構象的功能;設備端運行蒙特卡羅并行算法模塊,利用GPU并行化處理能力進行算法加速,主要思想是通過大量增加蒙特卡羅迭代搜索的初始態數量,以減少每次搜索深度,這樣在保證全局優化搜索結果一致的情況下,大量減少了分子對接軟件運行時間。

圖1 QVina2-GPU的異構并行架構Fig.1 Heterogeneous parallel architecture of QVina2-GPU

數據準備模塊包括讀取文件、配置環境、初始化數據和分配設備內存四個操作,為設備端基于GPU的蒙特卡羅并行算法的輸入做準備。

1)讀取文件模塊用于讀取配體和受體的格式文件和配置文件。格式文件存儲受體和配體原子坐標、部分電荷和原子種類信息;配置文件主要包含對接中心、對接盒子的體積、線程數和搜索深度,其中線程數N是設備端子任務的數量,搜索深度D決定蒙特卡羅并行算法迭代次數的大小。

2)設置Opencl模塊用于配置QVina2-GPU的異構并行環境,配置環境主要包括識別并選擇平臺和設備,創建OpenCL上下文、命令隊列、程序對象和內核等OpenCL環境。

3)初始化數據包括生成網格緩存和隨機數生成表,分別用于計算分子構象能量以及生成初始化分子構象,作為蒙特卡羅迭代搜索的初始態輸入。

4)分配設備內存模塊根據設定的子任務數量N,將生成的網格緩存、隨機數生成表、初始化分子構象分配在只讀內存的設備存儲器中,便于設備端讀取數據。

優化輸出模塊對設備端基于GPU的蒙特卡羅并行算法的輸出數據進行處理。設備端每個線程運行蒙特卡羅并行算法得到的最優分子構象分配在全局存儲器中,從N個分子構象中選取前k個最佳分子構象返回給主機端的優化輸出模塊,該模塊對k個最佳構象進行精細的局部搜索,并根據搜索結果輸出包含全局最優分子構象的配體文件。

2 基于GPU的蒙特卡羅并行算法

從圖1可以看出設備端的基于GPU的蒙特卡羅并行算法主要兩個部分:1)分配N個蒙特卡羅迭代搜索初始態;2)GPU子任務執行的基于Wolfe-Powell準則的蒙特卡羅搜索算法。

2.1 初始化蒙特卡羅迭代搜索

分子構象是指配體小分子和靶標蛋白質結合時配體的位置和姿態,蒙特卡羅搜索算法目的是在分子構象空間中找到一個最佳配體分子構象,并在此基礎上通過打分函數來計算配體和靶標結合的能量值,分子構象Ci∈iNC由一系列變量構成,具體表示為

(1)

2.2 子任務執行的基于Wolfe-Powell準則的蒙特卡羅搜索

每個線程運行的基于Wolfe-Powell準則蒙特卡羅搜索算法表示為GPU子任務,蒙特卡羅的算法流程如圖2所示,大體可以分為5步,包括隨機擾動、啟發式條件、局部搜索和模擬退火準則四個模塊。下面以線程Ti處理初始分子構象Ci∈￡為例說明子任務執行基于Wolfe-Powell準則的模特卡羅搜索過程。

圖2 基于Wolfe-Powell準則的蒙特卡羅搜索算法流程Fig. 2 Flow of Monte Carlo search algorithm based on Wolfe Powell criterion

Step 1首先對用隨機抖動函數R(·)對初始化分子構象Ci進行處理:

(2)

(3)

(4a)

(4b)

1)計算搜索方向pk:Bkpk=-f(xk),其中,f(·)為能量計算函數,xk和Bk分別是第k步搜索構象和海森矩陣。

2) 確定基于Wolfe-Powell準則的步長αk:

(5)

上式表達的是精確線性搜索中可接受的搜索步長αk,實際執行中,我們采用非精確線性搜索中滿足Wolfe-Powell準則的搜索步長αk,Wolfe-Powell準則描述具體見2.3節。

3) 更新搜索構象xk+1:xk+1=xk+αkpk,其中,pk和αk分別是上面計算的搜索方向和步長。

4) 更新海森矩陣Bk+1:

(6)

其中,sk=αkpk,yk=f(xk+1)-f(xk), 海森矩陣更新需要耗費計算資源,減少BFGS迭代搜索的深度可以有效的減少海森矩陣的更新,有效提高計算效率。如果‖f(xk+1)‖≤ε或迭代次數k=T,ε為預定值,則BFGS迭代結束,并設局部搜索的新構象

(7)

2.3 基于Wolfe-Powell準則的局部搜索算法

蒙特卡羅搜索算法中的局部搜索算法是找到最佳配體分子構象的重要模塊,相對于其他模塊耗費時間較長。QVina2利用基于不精確線性搜索Armijo準則局部搜索對分子構象進行優化,Armijo準則的主要思想是確定搜索方向pk,找到合適的步長αk,保證目標函數值下降,且下降量滿足不等式關系:

(8)

其中,c1為常數且滿足c1∈(0,1)。在步長更新的過程中,公式(8)主要目的是限制步長,確保下一個搜索構象的能量有足夠的下降,但是該條件不能確保步長值接近最佳步長,因為只要步長足夠小就可以滿足公式(8),如果步長過小會增加BFGS搜索迭代的次數。為了克服Armijo準則這一缺陷,我們提出使用基于Wolfe-Powell準則的局部搜索算法。Wolfe-Powell屬于不精確線性搜索準則,相比于Armijo準則,Wolfe-Powell準則更適合BFGS算法,可以保證海森矩陣迭代的正定性,增加了對步長的限制條件:

(9)

其中,c2為常數且滿足c2∈(0,1)。公式(9)的作用是限制過小的步長,保證能量函數的方向導數充分下降,找到更合適的步長,從而以更小的迭代次數找到更佳的分子構象,提高對接精度。因而,QVina2-GPU通過基于Wolfe-Powell準則局部搜索可以進一步減少蒙特卡羅迭代算法的搜索深度,在保證全局優化搜索結果一致的情況下,提高軟件對接速度。Wolfe-Powell準則確定步長的迭代過程如下:

Step 1給定c1∈(0,0.5),c2∈(c1,1),令α0=1,n=0,設定步長迭代次數L。

Step 2f(xk+1)=f(xk+αnpk),其中pk和xk分別是已知的搜索方向和搜索構象,αn為第n次迭代的步長。

Step 3若αn滿足公式(8)和(9)或步長迭代次數n=L,則αk=αn并退出步長更新迭代。若αn不滿足公式(8),令b=αn,αn=αn/2,n=n+1,返回Step 3。

Step 4若αn不滿足公式(9),令αn=min(2αn, 0.5*(αn+b)),n=n+1,返回Step 3。

綜上所述,基于Wolfe-Powell準則的蒙特卡羅搜索算法如下面所述,算法流程大體可以分為5步,包括隨機擾動、啟發式條件、局部搜索和模擬退火準則四個模塊。基于Wolfe-Powell準則的蒙特卡羅搜索算法交給GPU的子任務執行,這樣大量的子任務可以并行執行搜索算法,數量眾多的搜索線程可以實現對分子構象空間的有效搜索,可以減少單獨線程蒙特卡羅搜索算法的搜索深度,在保證全局優化搜索結果一致的情況下,顯著減少了分子對接軟件運行時間。

3 實驗及性能分析

3.1 實驗環境及數據集

QVina2-GPU的異構并行架構通過開放運算語言OpenCL實現,實驗環境的CPU為4核的Intel(R) Xeon(R) Gold 6130 CPU @2.10GHz,GPU為NVIDIA GeForce RTX 3090,操作系統為Windows 10,編譯器為Visual Studio 2019,CUDA的版本為11.1,OpenCL的版本為3.0。

圖3顯示了配體MGT1484和藥物靶標受體PROTEIN分子對接結果的3D結構示意圖,圖中左側框圖為分子對接過程的對接口袋,右側框圖為具體的對接過程,受體和配體之間的虛線為分子作用力。分子對接軟件就是通過打分函數計算分子對接過程的結合能量,來預測藥物靶標受體和配體結合的可能性,選出具有可能成藥的活性化合物。

圖3 藥物靶標受體和配體分子對接示意圖Fig.3 Schematic diagram of drug target receptor and ligand molecular docking

算法: 基于Wolfe-Powell準則的蒙特卡羅搜索算法輸入:N個隨機初始構象{C1, …,CN}:搜索深度:D;set d=0;Parallel for Ci(i=0, ..., N) dowhile d

使用包含140個藥物靶標受體-配體復合物的數據集作為實驗數據集,其中85個來自Astex多樣性集[23]、35個來自CASF-2013[24]和20個來自蛋白質數據庫[25],該數據集包含了廣泛的配體復雜度和靶標性質。數據集根據原子數目Natom分為簡單復合物、中等復雜復合物和復雜復合物,其中,簡單復合物的數量為46,原子數目Natom范圍為[5,23],中等復雜復合物的數量為51,Natom范圍為[24,36],復雜復合物的數量為43,Natom范圍為[37,108],三種復合物對應配體復雜度越來越高。

3.2 實驗指標和對比算法

對接精度和對接時間是評估分子對接軟件性能最重要的兩個指標,對接精度主要由對接分數Score和和RMSD(Root mean square deviation)決定,1)對接分數Score表示配體和受體之間的結合能,Score越低則結合越穩定,2)RMSD表示分子對接軟件搜索得到的構象與生物實驗獲得的X-ray結構(標準)之間的相似性,RMSD越低則表示搜索到的構象精度越高,當RMSD小于2?(10-10m)時,該構象是可以被接受,表示成功預測。3)對接時間是整個分子對接軟件的運行時間,對接時間越短,對接精度越高,表示分子對接軟件性能越好。

在相同的實驗環境中比較QVina2與QVina2-GPU分子對接軟件的性能,兩種算法的介紹如下:

1)QVina2(Quick Vina2):是一種快速準確的分子對接軟件,并提出了啟發式算法改進了局部搜索算法,與傳統的分子對接軟件相比提高了運行速度。

2)QVina2-GPU:是我們在QVina2對接軟件基礎上,提出的一種異構并行架構并實現了基于Wolfe-Powell準則的蒙特卡羅搜索,利用GPU硬件高度并行體系加速分子對接過程。

在分子對接軟件運行前,需要先設置配置文件。QVina2的配置文件中的配置參數包括搜索空間的中心、對接盒子的體積,另外還包括兩個重要的參數:初始化分子構象數量N和搜索深度D,二者反應的分子構象搜索的廣度和深度,其中分子對接軟件的運行時間主要由搜索深度D決定,

3)QVina2默認參數N為8,搜索深度D隨著配體分子的復雜度增加而增大,對于簡單分子復合物、中等復雜復合物和復雜復合物的平均搜索深度分別為16 170、29 663和41 212。

QVina2-GPU利用GPU并行化處理能力大量增加初始化分子構象數量N,使得大量線程分別對應不同的分子初始構象,使得每個線程獨立承擔蒙特卡羅搜索算法子任務,獨立搜索最佳的分子構象。以減少每個子任務的每次搜索深度D,在保證全局優化搜索結果一致的情況下,顯著減少了分子對接軟件運行時間。下面將重點討論初始化分子構象數量N和搜索深度D對性能的影響。

3.3 兩種準則對接精度比較

比較對象QVina2使用的是基于Armijo準則的蒙特卡羅搜索算法,我們提出的QVina2-GPU使用的是基于Wolfe-Powell準則的蒙特卡羅搜索算法。因而,本次實驗討論基于Wolfe-Powell準則的蒙特卡羅搜索算法和基于Armijo準則的蒙特卡羅搜索算法的對接精度,即對接分數Score和和RMSD指標。

在140個配體小分子復合物上進行對接實驗,對比兩種準則下對對接分數Score和RMSD的影響,實驗結果如圖4所示,其中圖4a表示140個配體小分子在兩種準則條件下的Score,其中Score的單位是kcal/mol,表示每摩爾分子千卡能量,圖4b表示兩種準則條件下的RMSD指標。兩幅圖的縱坐標表示基于Armijo準則的指標,橫坐標表示基于Wolfe-Powell準則的指標。顏色條表示配體中的原子數,顏色條從深到淺對應原子數目由少至多,即配體分子復雜度越來越高。

圖4 兩種準則的Score和RMSD對比圖 (Score和RMSD越小越好)Fig.4 Comparison diagram of score and RMSD of the two criteria (the smaller the score and RMSD, the better)

結果顯示數據集中60.71%的復合物位于圖4a的對角線上方,表明使用Wolfe-Powell準則得到的對接分數Score小于Armijo準則的Score,即Wolfe-Powell準則性能優于Armijo準則;有39.29%的復合物位于圖4a的對角線上,表明用Wolfe-Powell準則得到的Score與Armijo準則相同;基于Wolfe-Powell準則的蒙特卡羅搜索算法有63.57%構象位于圖4b中垂直紅色虛線左邊區域內,表示RMSD指標小于2?(10-10m)可被接受,而基于Armijo準則蒙特卡羅搜索算法只有59.28% 構象位于圖4b 水平紅色虛線下方區域內,可被接受。結果表明,相對于QVina2中基于Armijo準則蒙特卡羅搜索算法,我們提出的基于Wolfe-Powell準則的蒙特卡羅搜索算法具有更高的對接精度,在Score和RMSD指標中均明顯優于基于Armijo準則蒙特卡羅搜索算法。

3.4 QVina2-GPU中初始化分子構象數和搜索深度對性能的影響

不失一般性,從實驗數據集中隨機選取8個復合物,QVina2的初始化分子構象數N和搜索深度D按照配置文件設置并固定不變,QVina2-GPU的搜索深度D設置為1,初始化分子構象數N從100增大到8 000,討論N變化對QVina2-GPU的對接分數Score、RMSD和對接時間的影響。

如表1-表3所示,復合物從上到下對應配體復雜度越來越高,表1-表2可以看出初始化分子構象數N為800可以保證QVina2-GPU的Score和RMSD指標與QVina2相當,繼續增加初始化分子構象數N,QVina2-GPU的Score和RMSD指標將優于QVina2。表3的實驗結果表明在所有情況下QVina2-GPU對接時間相對于QVina2有了較大的降低,而且QVina2-GPU對接時間隨著初始化分子構象數N的增加而增加,對于復雜復合物這種增加的趨勢更加明顯,因此,可以得出:增加初始化分子構象數N,可以有效的提高分子對接精度,但同時也增加了對接時間開銷。綜上所述,為了在保證精度前提下盡量減少分子對接時間,下面的實驗中,將QVina2-GPU的初始化分子構象數N設置為800。

表1 初始化分子構象數N 對Score的影響(Score越小越好,標粗表示性能達到或優于QVina2)Table 1 Effect of initial molecular conformation number N on score(The smaller the score, the better. The bold numbers mean that the performance reaches or exceeds QVina2)

表2 初始化分子構象數N 對于RMSD的影響(RMSD越小越好,標粗表示性能達到或優于QVina2)Table 2 Effect of initial molecular conformation number N on RMSD(The smaller the RMSD, the better. The bold numbers indicate that the performance reaches or exceeds QVina2)

在確定初始化分子構象數N后,我們將搜索深度D從1增大到30,探究D對于QVina2-GPU的對接精度和對接時間的影響,表4-表6可以看出對于簡單復合物如5tim等,增大搜索深度D對Score、RMSE和對接時間的影響很小;對于中等復雜復合物如1hvy等,隨著搜索深度D增加,對接精度指標Score和RMSE總體上有所減少,但對接時間緩慢增加;對于復雜復合物3er5等,增大搜索深度D,對接精度指標Score和RMSE同樣總體上有所減少,但對接時間增加很快,會導致很大的時間開銷。

表4 搜索深度D 對Score的影響(Score越小越好,標粗表示性能達到或優于QVina2)Table 4 Effect of search depth D on score(The smaller the score, the better. The bold numbers indicate that the performance reaches or exceeds QVina2)

表5 搜索深度D 對RMSD的影響(RMSD越小越好,標粗表示性能達到或優于QVina2)Table 5 Effect of search depth D on RMSD(The smaller the RMSD, the better. The bold numbers indicate that the performance reaches or exceeds QVina2)

表6 搜索深度D 對于對接時間的影響(對接時間越小越好,標粗表示性能達到或優于QVina2)Table 6 Effect of Search Depth D on docking time(The smaller the docking time, the better. The bold numbers indicate that the performance reaches or exceeds QVina2)

綜上所述,為了在保證精度前提下盡量減少分子對接時間,特別是考慮到搜索深度D對于復雜復合物對接時間的影響特別明顯。因此,接下來的實驗,我們將數據集140個復合物的初始化分子構象數N設為800,搜索深度D設為1,以最少的對接時間保證了分子對接精度。

3.5 QVina2-GPU與QVina2的性能對比

首先對比了QVina2-GPU和QVina2在140個復合物上的對接分數Score和RMSD對接精度指標,兩幅圖的縱坐標表示QVina2的性能指標,橫坐標表示QVina2-GPU的性能指標。顏色條表示配體中的原子數,顏色條從深到淺對應原子數目由少至多,即配體分子復雜度越來越高。

對接分數Score結果如圖5a所示,大多數復合物分布在對角線周圍,其對接分數的皮爾森相關系數為0.914,表示它們之間存在極強相關性,即QVina2-GPU和QVina2對接分子Score的指標基本相同;圖5b顯示大多數復合物落入左下角區域,QVina2(水平紅色虛線下方區域內)和QVina2-GPU(垂直紅色虛線下方區域內)分別有58.28%和55%的復合物的預測結果可被接受,結果表明本文提出的QVina2-GPU在初始化分子構象數N為800和搜索深度D為1的條件下,達到了QVina2的對接精度。

圖5 QVina2-GPU與QVina2的Score和RMSD對接精度對比圖 (N=800,D=1)Fig. 5 Comparison chart of score and RMSD docking accuracy between QVina2-GPU and QVina2(N=800,D=1)

下面重點討論一下QVina2-GPU和QVina2在140個復合物上的對接時間指標,這里我們使用兩個對接時間指標,對接時間加速比Acc和蒙特卡羅時間加速比Accd,定義如下:

(10)

式中,Tqvina2是QVina2分子對接計算的對接時間,Tqvina2-gpu是QVina2-GPU分子對接計算的對接時間。QVina2與QVina2-GPU軟件運行中,蒙特卡羅迭代搜索算法都是最耗時的部分,因而,Tmc為QVina2中基于CPU的蒙特卡羅算法的運行時間,Tmc-gpu為QVina2-GPU的蒙特卡羅運行時間。對接時間加速比Acc和蒙特卡羅時間加速比Accd數值越大,說明我們提出的QVina2-GPU相對于QVina2有較高的加速比,即運行時間相對于QVina2有顯著的減少。

圖6和圖7給出了140個復合物的對接時間加速比Acc和蒙特卡羅時間加速比Accd。為了區分不同復雜度配體分子的加速比情況,三維坐標的橫縱坐標分別表示配體的原子數目Natom和可旋轉鍵數目Natom,配體的復雜度由原子數目Natom和可旋轉鍵數目Natom決定,一般來說Natom和Natom越大,配體的復雜度越高。圖6和圖7中每個條形柱體分別表示配體的復雜程度和對應的對接時間加速比Acc和蒙特卡羅時間加速比Accd,對接時間加速比的藍色柱狀條表示加速比達到10的復合物,蒙特卡羅時間加速比的藍色柱狀條為加速比達到60的復合物,越靠近藍色加速比越高,圖5可以看出Acc分布范圍在3.11倍～12.28倍之間,圖7可以看出Accd分布范圍為6.48倍～60.28倍,其中,2bm2復合物加速效果最好,它的蒙特卡羅時間加速比Acc達到了60.28倍,對接時間加速比Accd達到了12.28倍。上面數據充分說明了我們提出的QVina2-GPU在不同類型的復合物上都有明顯的加速效果,顯著的減少了分子對接的計算時間,這對于需要大量篩選配體分子的生物實驗特別有意義。

圖6 140個復合物的對接時間加速比Fig.6 Docking time acceleration ratio of 140 composites

圖7 140個復合的物蒙特卡羅時間加速比Fig.7 Monte Carlo time acceleration ratio of 140 compounds

為進一步顯示QVina2-GPU加速比性能,三種配體分子復雜度下的平均加速比結果(見表7),結果表明,最大平均蒙特卡羅時間加速比為22.91倍,最大平均對接時間加速比為5.71倍,三種不同復雜度復合物條件下,QVina2-GPU都取得了令人滿意的加速比,總體平均蒙特卡羅時間加速比和平均對接時間加速比分別到達了18.63和5.18。從上面的分析可知,QVina2-GPU利用GPU并行化處理能力大量增加初始化分子構象數量N,使得大量線程分別對應不同的分子初始構象,使得每個線程獨立承擔蒙特卡羅搜索算法子任務,增加初始化分子構象數量N以減少每個子任務的每次搜索深度D,從而顯著減少了蒙特卡羅構象搜索時間,減少了分子對接軟件的整體運行時間。平均對接時間加速比隨著復雜度的增大而提高,這是由于隨著復雜度的增加,蒙特卡羅搜索算法運行時間在整個軟件運行時間比例加重,平均對接時間加速比隨著增加。

表7 三種配體復雜度的平均時間和平均時間加速比Table 7 Average time and average time acceleration ratio of three ligand complexities

3.6 QVina2-GPU與Vina-GPU性能對比

為了區分Vina-GPU[26]和QVina2-GPU兩種方法在性能上的不同,我們從Vina-GPU論文里獲取可執行程序,在相同的硬件平臺實驗環境下,對上文提到的隨機8個復合物對比QVina2-GPU與Vina-GPU的對接分數Score和RMSD對接精度指標和對接時間,如表8所示,實驗表明,對于簡單分子和中等復雜分子QVina2-GPU可以在對接精度基本相同的情況下加速Vina-GPU,對于復雜分子,我們提出的QVina2-GPU在對接精度略微下降的條件下,運行時間上有明顯的減少,相對于Vina-GPU的加速比最大可達32.92,這對于實際藥物篩選應用具有很大的意義。

表8 Vina-GPU與QVina2-GPU的性能對比Table 8 Performance comparison between Vina GPU and QVina2 GPU

4 結 論

本文提出基于GPU加速的分子對接軟件并行化方法QVina2-GPU,針對Qvina2在大型虛擬數據庫中篩選時間過長的情況,通過增加初始化分子構象數量來擴展蒙特卡羅的迭代局部搜索中線程的并行規模,減少蒙特卡羅迭代搜索算法的搜索步數,而且利用Wolfe-Powell準則改進局部搜索算法,提高了對接精度,進一步減少蒙特卡羅迭代算法的搜索深度,提高分子對接軟件的運行速度。在公開的藥物靶標受體-配體復合物的數據集上驗證了不同類型復合物上的平均模特卡羅時間加速比和平均對接時間加速比,結果表明QVina2-GPU在達到Qvina2分子對接精度的條件下,相對于QVina2有顯著的加速效果。QVina2-GPU的代碼可以在https://github.com/DeltaGroupNJUPT/QuickVina2-GPU上獲得,其中包含可執行程序和程序使用說明,便于科研人員使用。

由于QVina2-GPU計算構象的能量使用的依舊是QVina2的能量計算函數,是一種經驗公式,篩選出能量較低的復合物用于生物實驗中,存在較高的偏差[27],近年來已有學者研究基于深度學習的能量計算函數。未來,我們的工作將通過研究新的能量計算函數來提高QVina2-GPU的對接性能。