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

胡海峰,王領(lǐng)悅,唐詩迪,胡鳴珂,吳建盛

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

虛擬篩選作為計(jì)算機(jī)輔助藥物設(shè)計(jì)的實(shí)用技術(shù),在現(xiàn)代藥物研發(fā)中起著重要作用[1]。虛擬篩選在計(jì)算平臺(tái)上模擬藥物篩選的過程,可快速?gòu)膸资辽习偃f的小分子化合物(配體)中,篩選出有可能和受體結(jié)合的活性化合物,大大降低實(shí)際篩選小分子化合物數(shù)目,縮短藥物研發(fā)周期,降低藥物研發(fā)的成本[2-3]。虛擬篩選的方法可分為基于結(jié)構(gòu)的虛擬篩選和基于配體的虛擬篩選[4]?；诮Y(jié)構(gòu)的虛擬篩選利用分子對(duì)接技術(shù),在已知蛋白質(zhì)受體的三維結(jié)構(gòu)的基礎(chǔ)上,在結(jié)合位點(diǎn)處自動(dòng)匹配化合物數(shù)據(jù)庫中的小分子配體,用打分函數(shù)對(duì)可能的結(jié)合模式進(jìn)行能量計(jì)算,最終得到小分子化合物結(jié)合活性的排名[5-7]。而基于配體的虛擬篩選方法不需要知道目標(biāo)蛋白質(zhì)受體的結(jié)構(gòu),代表性的方法包括藥效團(tuán)模型、定量構(gòu)效關(guān)系和結(jié)構(gòu)相似性等方法[8-9]。相對(duì)于基于配體的虛擬篩選,基于結(jié)構(gòu)的虛擬篩選能避免活性化合物結(jié)構(gòu)微小變化引起的活性改變[10],對(duì)小分子配體的活性值預(yù)測(cè)更加準(zhǔn)確。

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

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

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

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

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

3)使用公開數(shù)據(jù)庫中140個(gè)小分子復(fù)合物[22]作為測(cè)試數(shù)據(jù)集,實(shí)驗(yàn)結(jié)果驗(yàn)證了在保證對(duì)接精度的條件下,我們提出的QVina2-GPU對(duì)Qvina2的平均加速比達(dá)到5.18倍,最大加速比達(dá)到12.28倍,具有很大的實(shí)際應(yīng)用前景。

1 QVina-GPU異構(gòu)并行架構(gòu)

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

圖1 QVina2-GPU的異構(gòu)并行架構(gòu)Fig.1 Heterogeneous parallel architecture of QVina2-GPU

數(shù)據(jù)準(zhǔn)備模塊包括讀取文件、配置環(huán)境、初始化數(shù)據(jù)和分配設(shè)備內(nèi)存四個(gè)操作,為設(shè)備端基于GPU的蒙特卡羅并行算法的輸入做準(zhǔn)備。

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

2)設(shè)置Opencl模塊用于配置QVina2-GPU的異構(gòu)并行環(huán)境,配置環(huán)境主要包括識(shí)別并選擇平臺(tái)和設(shè)備,創(chuàng)建OpenCL上下文、命令隊(duì)列、程序?qū)ο蠛蛢?nèi)核等OpenCL環(huán)境。

3)初始化數(shù)據(jù)包括生成網(wǎng)格緩存和隨機(jī)數(shù)生成表,分別用于計(jì)算分子構(gòu)象能量以及生成初始化分子構(gòu)象,作為蒙特卡羅迭代搜索的初始態(tài)輸入。

4)分配設(shè)備內(nèi)存模塊根據(jù)設(shè)定的子任務(wù)數(shù)量N,將生成的網(wǎng)格緩存、隨機(jī)數(shù)生成表、初始化分子構(gòu)象分配在只讀內(nèi)存的設(shè)備存儲(chǔ)器中,便于設(shè)備端讀取數(shù)據(jù)。

優(yōu)化輸出模塊對(duì)設(shè)備端基于GPU的蒙特卡羅并行算法的輸出數(shù)據(jù)進(jìn)行處理。設(shè)備端每個(gè)線程運(yùn)行蒙特卡羅并行算法得到的最優(yōu)分子構(gòu)象分配在全局存儲(chǔ)器中,從N個(gè)分子構(gòu)象中選取前k個(gè)最佳分子構(gòu)象返回給主機(jī)端的優(yōu)化輸出模塊,該模塊對(duì)k個(gè)最佳構(gòu)象進(jìn)行精細(xì)的局部搜索,并根據(jù)搜索結(jié)果輸出包含全局最優(yōu)分子構(gòu)象的配體文件。

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

從圖1可以看出設(shè)備端的基于GPU的蒙特卡羅并行算法主要兩個(gè)部分:1)分配N個(gè)蒙特卡羅迭代搜索初始態(tài);2)GPU子任務(wù)執(zhí)行的基于Wolfe-Powell準(zhǔn)則的蒙特卡羅搜索算法。

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

分子構(gòu)象是指配體小分子和靶標(biāo)蛋白質(zhì)結(jié)合時(shí)配體的位置和姿態(tài),蒙特卡羅搜索算法目的是在分子構(gòu)象空間中找到一個(gè)最佳配體分子構(gòu)象,并在此基礎(chǔ)上通過打分函數(shù)來計(jì)算配體和靶標(biāo)結(jié)合的能量值,分子構(gòu)象Ci∈iNC由一系列變量構(gòu)成,具體表示為

(1)

2.2 子任務(wù)執(zhí)行的基于Wolfe-Powell準(zhǔn)則的蒙特卡羅搜索

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

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

Step 1首先對(duì)用隨機(jī)抖動(dòng)函數(shù)R(·)對(duì)初始化分子構(gòu)象Ci進(jìn)行處理:

(2)

(3)

(4a)

(4b)

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

2) 確定基于Wolfe-Powell準(zhǔn)則的步長(zhǎng)αk:

(5)

上式表達(dá)的是精確線性搜索中可接受的搜索步長(zhǎng)αk,實(shí)際執(zhí)行中,我們采用非精確線性搜索中滿足Wolfe-Powell準(zhǔn)則的搜索步長(zhǎng)αk,Wolfe-Powell準(zhǔn)則描述具體見2.3節(jié)。

3) 更新搜索構(gòu)象xk+1:xk+1=xk+αkpk,其中,pk和αk分別是上面計(jì)算的搜索方向和步長(zhǎng)。

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

(6)

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

(7)

2.3 基于Wolfe-Powell準(zhǔn)則的局部搜索算法

蒙特卡羅搜索算法中的局部搜索算法是找到最佳配體分子構(gòu)象的重要模塊,相對(duì)于其他模塊耗費(fèi)時(shí)間較長(zhǎng)。QVina2利用基于不精確線性搜索Armijo準(zhǔn)則局部搜索對(duì)分子構(gòu)象進(jìn)行優(yōu)化,Armijo準(zhǔn)則的主要思想是確定搜索方向pk,找到合適的步長(zhǎng)αk,保證目標(biāo)函數(shù)值下降,且下降量滿足不等式關(guān)系:

(8)

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

(9)

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

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

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

Step 3若αn滿足公式(8)和(9)或步長(zhǎng)迭代次數(shù)n=L,則αk=αn并退出步長(zhǎng)更新迭代。若α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準(zhǔn)則的蒙特卡羅搜索算法如下面所述,算法流程大體可以分為5步,包括隨機(jī)擾動(dòng)、啟發(fā)式條件、局部搜索和模擬退火準(zhǔn)則四個(gè)模塊?；赪olfe-Powell準(zhǔn)則的蒙特卡羅搜索算法交給GPU的子任務(wù)執(zhí)行,這樣大量的子任務(wù)可以并行執(zhí)行搜索算法,數(shù)量眾多的搜索線程可以實(shí)現(xiàn)對(duì)分子構(gòu)象空間的有效搜索,可以減少單獨(dú)線程蒙特卡羅搜索算法的搜索深度,在保證全局優(yōu)化搜索結(jié)果一致的情況下,顯著減少了分子對(duì)接軟件運(yùn)行時(shí)間。

3 實(shí)驗(yàn)及性能分析

3.1 實(shí)驗(yàn)環(huán)境及數(shù)據(jù)集

QVina2-GPU的異構(gòu)并行架構(gòu)通過開放運(yùn)算語言O(shè)penCL實(shí)現(xiàn),實(shí)驗(yàn)環(huán)境的CPU為4核的Intel(R) Xeon(R) Gold 6130 CPU @2.10GHz,GPU為NVIDIA GeForce RTX 3090,操作系統(tǒng)為Windows 10,編譯器為Visual Studio 2019,CUDA的版本為11.1,OpenCL的版本為3.0。

圖3顯示了配體MGT1484和藥物靶標(biāo)受體PROTEIN分子對(duì)接結(jié)果的3D結(jié)構(gòu)示意圖,圖中左側(cè)框圖為分子對(duì)接過程的對(duì)接口袋,右側(cè)框圖為具體的對(duì)接過程,受體和配體之間的虛線為分子作用力。分子對(duì)接軟件就是通過打分函數(shù)計(jì)算分子對(duì)接過程的結(jié)合能量,來預(yù)測(cè)藥物靶標(biāo)受體和配體結(jié)合的可能性,選出具有可能成藥的活性化合物。

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

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

使用包含140個(gè)藥物靶標(biāo)受體-配體復(fù)合物的數(shù)據(jù)集作為實(shí)驗(yàn)數(shù)據(jù)集,其中85個(gè)來自Astex多樣性集[23]、35個(gè)來自CASF-2013[24]和20個(gè)來自蛋白質(zhì)數(shù)據(jù)庫[25],該數(shù)據(jù)集包含了廣泛的配體復(fù)雜度和靶標(biāo)性質(zhì)。數(shù)據(jù)集根據(jù)原子數(shù)目Natom分為簡(jiǎn)單復(fù)合物、中等復(fù)雜復(fù)合物和復(fù)雜復(fù)合物,其中,簡(jiǎn)單復(fù)合物的數(shù)量為46,原子數(shù)目Natom范圍為[5,23],中等復(fù)雜復(fù)合物的數(shù)量為51,Natom范圍為[24,36],復(fù)雜復(fù)合物的數(shù)量為43,Natom范圍為[37,108],三種復(fù)合物對(duì)應(yīng)配體復(fù)雜度越來越高。

3.2 實(shí)驗(yàn)指標(biāo)和對(duì)比算法

對(duì)接精度和對(duì)接時(shí)間是評(píng)估分子對(duì)接軟件性能最重要的兩個(gè)指標(biāo),對(duì)接精度主要由對(duì)接分?jǐn)?shù)Score和和RMSD(Root mean square deviation)決定,1)對(duì)接分?jǐn)?shù)Score表示配體和受體之間的結(jié)合能,Score越低則結(jié)合越穩(wěn)定,2)RMSD表示分子對(duì)接軟件搜索得到的構(gòu)象與生物實(shí)驗(yàn)獲得的X-ray結(jié)構(gòu)(標(biāo)準(zhǔn))之間的相似性,RMSD越低則表示搜索到的構(gòu)象精度越高,當(dāng)RMSD小于2?(10-10m)時(shí),該構(gòu)象是可以被接受,表示成功預(yù)測(cè)。3)對(duì)接時(shí)間是整個(gè)分子對(duì)接軟件的運(yùn)行時(shí)間,對(duì)接時(shí)間越短,對(duì)接精度越高,表示分子對(duì)接軟件性能越好。

在相同的實(shí)驗(yàn)環(huán)境中比較QVina2與QVina2-GPU分子對(duì)接軟件的性能,兩種算法的介紹如下:

1)QVina2(Quick Vina2):是一種快速準(zhǔn)確的分子對(duì)接軟件,并提出了啟發(fā)式算法改進(jìn)了局部搜索算法,與傳統(tǒng)的分子對(duì)接軟件相比提高了運(yùn)行速度。

2)QVina2-GPU:是我們?cè)赒Vina2對(duì)接軟件基礎(chǔ)上,提出的一種異構(gòu)并行架構(gòu)并實(shí)現(xiàn)了基于Wolfe-Powell準(zhǔn)則的蒙特卡羅搜索,利用GPU硬件高度并行體系加速分子對(duì)接過程。

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

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

QVina2-GPU利用GPU并行化處理能力大量增加初始化分子構(gòu)象數(shù)量N,使得大量線程分別對(duì)應(yīng)不同的分子初始構(gòu)象,使得每個(gè)線程獨(dú)立承擔(dān)蒙特卡羅搜索算法子任務(wù),獨(dú)立搜索最佳的分子構(gòu)象。以減少每個(gè)子任務(wù)的每次搜索深度D,在保證全局優(yōu)化搜索結(jié)果一致的情況下,顯著減少了分子對(duì)接軟件運(yùn)行時(shí)間。下面將重點(diǎn)討論初始化分子構(gòu)象數(shù)量N和搜索深度D對(duì)性能的影響。

3.3 兩種準(zhǔn)則對(duì)接精度比較

比較對(duì)象QVina2使用的是基于Armijo準(zhǔn)則的蒙特卡羅搜索算法,我們提出的QVina2-GPU使用的是基于Wolfe-Powell準(zhǔn)則的蒙特卡羅搜索算法。因而,本次實(shí)驗(yàn)討論基于Wolfe-Powell準(zhǔn)則的蒙特卡羅搜索算法和基于Armijo準(zhǔn)則的蒙特卡羅搜索算法的對(duì)接精度,即對(duì)接分?jǐn)?shù)Score和和RMSD指標(biāo)。

在140個(gè)配體小分子復(fù)合物上進(jìn)行對(duì)接實(shí)驗(yàn),對(duì)比兩種準(zhǔn)則下對(duì)對(duì)接分?jǐn)?shù)Score和RMSD的影響,實(shí)驗(yàn)結(jié)果如圖4所示,其中圖4a表示140個(gè)配體小分子在兩種準(zhǔn)則條件下的Score,其中Score的單位是kcal/mol,表示每摩爾分子千卡能量,圖4b表示兩種準(zhǔn)則條件下的RMSD指標(biāo)。兩幅圖的縱坐標(biāo)表示基于Armijo準(zhǔn)則的指標(biāo),橫坐標(biāo)表示基于Wolfe-Powell準(zhǔn)則的指標(biāo)。顏色條表示配體中的原子數(shù),顏色條從深到淺對(duì)應(yīng)原子數(shù)目由少至多,即配體分子復(fù)雜度越來越高。

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

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

3.4 QVina2-GPU中初始化分子構(gòu)象數(shù)和搜索深度對(duì)性能的影響

不失一般性,從實(shí)驗(yàn)數(shù)據(jù)集中隨機(jī)選取8個(gè)復(fù)合物,QVina2的初始化分子構(gòu)象數(shù)N和搜索深度D按照配置文件設(shè)置并固定不變,QVina2-GPU的搜索深度D設(shè)置為1,初始化分子構(gòu)象數(shù)N從100增大到8 000,討論N變化對(duì)QVina2-GPU的對(duì)接分?jǐn)?shù)Score、RMSD和對(duì)接時(shí)間的影響。

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

表1 初始化分子構(gòu)象數(shù)N 對(duì)Score的影響(Score越小越好,標(biāo)粗表示性能達(dá)到或優(yōu)于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 初始化分子構(gòu)象數(shù)N 對(duì)于RMSD的影響(RMSD越小越好,標(biāo)粗表示性能達(dá)到或優(yōu)于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)

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

表4 搜索深度D 對(duì)Score的影響(Score越小越好,標(biāo)粗表示性能達(dá)到或優(yōu)于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 對(duì)RMSD的影響(RMSD越小越好,標(biāo)粗表示性能達(dá)到或優(yōu)于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 對(duì)于對(duì)接時(shí)間的影響(對(duì)接時(shí)間越小越好,標(biāo)粗表示性能達(dá)到或優(yōu)于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)

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

3.5 QVina2-GPU與QVina2的性能對(duì)比

首先對(duì)比了QVina2-GPU和QVina2在140個(gè)復(fù)合物上的對(duì)接分?jǐn)?shù)Score和RMSD對(duì)接精度指標(biāo),兩幅圖的縱坐標(biāo)表示QVina2的性能指標(biāo),橫坐標(biāo)表示QVina2-GPU的性能指標(biāo)。顏色條表示配體中的原子數(shù),顏色條從深到淺對(duì)應(yīng)原子數(shù)目由少至多,即配體分子復(fù)雜度越來越高。

對(duì)接分?jǐn)?shù)Score結(jié)果如圖5a所示,大多數(shù)復(fù)合物分布在對(duì)角線周圍,其對(duì)接分?jǐn)?shù)的皮爾森相關(guān)系數(shù)為0.914,表示它們之間存在極強(qiáng)相關(guān)性,即QVina2-GPU和QVina2對(duì)接分子Score的指標(biāo)基本相同;圖5b顯示大多數(shù)復(fù)合物落入左下角區(qū)域,QVina2(水平紅色虛線下方區(qū)域內(nèi))和QVina2-GPU(垂直紅色虛線下方區(qū)域內(nèi))分別有58.28%和55%的復(fù)合物的預(yù)測(cè)結(jié)果可被接受,結(jié)果表明本文提出的QVina2-GPU在初始化分子構(gòu)象數(shù)N為800和搜索深度D為1的條件下,達(dá)到了QVina2的對(duì)接精度。

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

下面重點(diǎn)討論一下QVina2-GPU和QVina2在140個(gè)復(fù)合物上的對(duì)接時(shí)間指標(biāo),這里我們使用兩個(gè)對(duì)接時(shí)間指標(biāo),對(duì)接時(shí)間加速比Acc和蒙特卡羅時(shí)間加速比Accd,定義如下:

(10)

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

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

圖6 140個(gè)復(fù)合物的對(duì)接時(shí)間加速比Fig.6 Docking time acceleration ratio of 140 composites

圖7 140個(gè)復(fù)合的物蒙特卡羅時(shí)間加速比Fig.7 Monte Carlo time acceleration ratio of 140 compounds

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

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

3.6 QVina2-GPU與Vina-GPU性能對(duì)比

為了區(qū)分Vina-GPU[26]和QVina2-GPU兩種方法在性能上的不同,我們從Vina-GPU論文里獲取可執(zhí)行程序,在相同的硬件平臺(tái)實(shí)驗(yàn)環(huán)境下,對(duì)上文提到的隨機(jī)8個(gè)復(fù)合物對(duì)比QVina2-GPU與Vina-GPU的對(duì)接分?jǐn)?shù)Score和RMSD對(duì)接精度指標(biāo)和對(duì)接時(shí)間,如表8所示,實(shí)驗(yàn)表明,對(duì)于簡(jiǎn)單分子和中等復(fù)雜分子QVina2-GPU可以在對(duì)接精度基本相同的情況下加速Vina-GPU,對(duì)于復(fù)雜分子,我們提出的QVina2-GPU在對(duì)接精度略微下降的條件下,運(yùn)行時(shí)間上有明顯的減少,相對(duì)于Vina-GPU的加速比最大可達(dá)32.92,這對(duì)于實(shí)際藥物篩選應(yīng)用具有很大的意義。

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

4 結(jié) 論

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

由于QVina2-GPU計(jì)算構(gòu)象的能量使用的依舊是QVina2的能量計(jì)算函數(shù),是一種經(jīng)驗(yàn)公式,篩選出能量較低的復(fù)合物用于生物實(shí)驗(yàn)中,存在較高的偏差[27],近年來已有學(xué)者研究基于深度學(xué)習(xí)的能量計(jì)算函數(shù)。未來,我們的工作將通過研究新的能量計(jì)算函數(shù)來提高QVina2-GPU的對(duì)接性能。