基于CPU和GPU異構(gòu)的蛋白質(zhì)分子半柔性對(duì)接算法優(yōu)化

陸旭峰 陸振宇 梅向東 孔 韌 常 珊

(1.江蘇理工學(xué)院電氣信息工程學(xué)院生物信息與醫(yī)藥工程研究所, 常州，213001； 2.江蘇贊奇科技股份有限公司，常州，213022)

引 言

分子對(duì)接是用來(lái)預(yù)測(cè)分子之間相互作用和相互識(shí)別的模式，而蛋白質(zhì)分子對(duì)接是預(yù)測(cè)蛋白質(zhì)受體和蛋白質(zhì)配體的結(jié)合匹配模式[1]。通常，生物學(xué)實(shí)驗(yàn)是探究分子之間的相互作用的直接手段。但通過(guò)實(shí)驗(yàn)直接研究蛋白質(zhì)相互作用并不容易，獲得蛋白質(zhì)復(fù)合物晶體的生物實(shí)驗(yàn)條件比較嚴(yán)苛，并且實(shí)驗(yàn)的成本也非常昂貴，這使理論模擬方法得到了越來(lái)越多的重視和應(yīng)用。分子對(duì)接方法在研究分子間相互作用領(lǐng)域有著廣泛的應(yīng)用，可減少單純實(shí)驗(yàn)帶來(lái)的種種困難，為實(shí)驗(yàn)者們提供重要的理論依據(jù)。

分子對(duì)接主要分為3種，分別是剛性對(duì)接(Rigid docking)，半柔性對(duì)接(Semi flexible docking)和柔性對(duì)接(Flexible docking)[2-3]。柔性對(duì)接指的是在對(duì)接過(guò)程中，受體和配體都發(fā)生構(gòu)象的改變，這種對(duì)接在于精度上得到了很好的保證，但是隨之產(chǎn)生的負(fù)面效應(yīng)是計(jì)算量的大幅增加，而剛性對(duì)接則恰恰相反。本文主要采用的是半柔性對(duì)接，即在保證受體構(gòu)象不發(fā)生變化的情況下，配體的構(gòu)象發(fā)生一定的改變，在計(jì)算效率和精度上都能達(dá)到相應(yīng)的要求。

現(xiàn)在被廣泛使用的分子對(duì)接軟件主要有：AutoDock,Molegro Vitual Docker,3D_Dock,DOCK以及FlexX等[4]。其中，3D_Dock程序具有很大的優(yōu)越性[5]，能夠快速進(jìn)行全局對(duì)接和復(fù)合物結(jié)構(gòu)數(shù)據(jù)的篩選。但是，隨著對(duì)接數(shù)據(jù)量的不斷增大以及計(jì)算機(jī)硬件設(shè)備及處理能力的不斷更新，傳統(tǒng)的3D_Dock程序已經(jīng)不能勝任目前的情況。特別是針對(duì)半柔性對(duì)接計(jì)算量比較大，原有的計(jì)算方法和并行策略都明顯需要改進(jìn)。

NVIDIA公司發(fā)布的統(tǒng)一計(jì)算設(shè)備架構(gòu) (Compute unified device architecture，CUDA)[6-7]，讓GPU作為計(jì)算單元得到了強(qiáng)有力的開(kāi)發(fā)。GPU作為圖形處理器具有超強(qiáng)的浮點(diǎn)運(yùn)算能力，而將GPU從傳統(tǒng)的圖形處理器變?yōu)橛?jì)算單元來(lái)使用，最大化地利用其浮點(diǎn)運(yùn)算能力，有可能解決上述分子對(duì)接面臨的問(wèn)題。如何將GPU和3D_Dock程序結(jié)合起來(lái)是本文的研究重點(diǎn)。

本文在原有3D_Dock程序基礎(chǔ)上進(jìn)行改進(jìn)。首先，在CPU上，基于FFTW3對(duì)原有的程序進(jìn)行算法升級(jí)和程序改進(jìn)。然后，在GPU上，基于CUFFTW和CUFFT對(duì)改進(jìn)后的程序進(jìn)行再升級(jí)。通過(guò)3種方法對(duì)程序進(jìn)行改進(jìn)和算法升級(jí)，用改進(jìn)程序?qū)DB code 為1PPE的測(cè)試蛋白進(jìn)行結(jié)合態(tài)和自由態(tài)對(duì)接，結(jié)果驗(yàn)證了程序的正確性并給出了改進(jìn)后程序的加速比。

1 3D_Dock主要流程

分子對(duì)接首先要對(duì)受體和配體分子結(jié)合模式進(jìn)行全空間搜索[8-9]，但是，蛋白質(zhì)的自身特點(diǎn)決定了空間搜索的困難，即蛋白質(zhì)對(duì)接所產(chǎn)生的復(fù)合物數(shù)量十分龐大。傳統(tǒng)的用于分子對(duì)接的空間搜索算法主要有遺傳算法、模擬退火算法、快速傅里葉變換算法(Fast Fourier transform，F(xiàn)FT)、片段生長(zhǎng)法以及蒙特卡羅算法等。3D_Dock程序主要采用FFT算法來(lái)進(jìn)行復(fù)合物構(gòu)象的計(jì)算，只要所需評(píng)價(jià)函數(shù)滿足∑ijpiqj的形式，就可以通過(guò)FFT算法進(jìn)行程序的加速，提高采樣速度。其主要理論依據(jù)在于通過(guò)離散函數(shù)表現(xiàn)分子特性，具體的流程如下：

假設(shè)需要進(jìn)行對(duì)接的受體和配體分子的特性，以及通過(guò)對(duì)接得到的復(fù)合物的特性都通過(guò)格點(diǎn)表現(xiàn)出來(lái)，則可以通過(guò)格點(diǎn)特性來(lái)求解特性匹配值。匹配值計(jì)算公式如下

(1)

步驟3求頻域的點(diǎn)乘。

步驟4根據(jù)步驟3得到了所需的復(fù)合物結(jié)構(gòu)數(shù)據(jù)Co,p,q，但是得到的值為頻域下的值，因此，下一步需要進(jìn)行的是傅里葉逆變換。

主要計(jì)算公式如下

(2)

(3)

(4)

(5)

經(jīng)過(guò)式(5)計(jì)算后得到復(fù)合物的相關(guān)特性匹配值，根據(jù)排序篩選，取得該輪計(jì)算的極小值并進(jìn)行保存，至此便完成了一輪計(jì)算。接下來(lái)通過(guò)改變配體的不同角度繼續(xù)進(jìn)行對(duì)接，重復(fù)以上步驟，得到一系列對(duì)接復(fù)合物結(jié)構(gòu)。需要特別注意的是，因?yàn)槭前肴嵝詫?duì)接，即受體在構(gòu)象上不發(fā)生變化，因此受體只需要進(jìn)行一次傅里葉變換即可。完成所有對(duì)接情況后，對(duì)保留下來(lái)的數(shù)據(jù)繼續(xù)進(jìn)行排序，得到打分項(xiàng)最高的復(fù)合物結(jié)構(gòu)。

2 優(yōu)化方案介紹

2.1 基于FFTW的優(yōu)化方案一

3D_Dock的核心程序是FFTW。FFTW主要有兩個(gè)經(jīng)典版本FFTW2和FFTW3，但是FFTW3并不兼容FFTW2，而原程序的版本為FFTW2[10-11]。因此，需要有在CPU上優(yōu)化程序的方案，即將原來(lái)FFTW2的函數(shù)全部升級(jí)為FFTW3的函數(shù)。優(yōu)化方案一基于CPU，將原程序的FFTW2進(jìn)行升級(jí)為FFTW3，其主要流程如下：

步驟1準(zhǔn)備工作：主要涉及內(nèi)存分配、默認(rèn)參數(shù)的設(shè)置、角度計(jì)算、原子帶電量設(shè)定、離散后網(wǎng)格劃分、動(dòng)靜態(tài)分子結(jié)構(gòu)的讀入，以及建立傅里葉變換和逆變換策略。

步驟2將靜態(tài)分子進(jìn)行平面化和離散化，動(dòng)態(tài)分子先采用分子動(dòng)力學(xué)模擬生成3種構(gòu)象，然后進(jìn)行離散化。

步驟3分別將靜態(tài)分子和動(dòng)態(tài)分子得到的網(wǎng)格數(shù)據(jù)進(jìn)行傅里葉變換，并進(jìn)行點(diǎn)乘計(jì)算，得到復(fù)合物投影離散函數(shù)，并對(duì)其進(jìn)行傅里葉逆變換，得到復(fù)合物相關(guān)函數(shù)。

步驟4對(duì)得到的復(fù)合物相關(guān)函數(shù)進(jìn)行匹配因子計(jì)算，并進(jìn)行打分排序，篩選出最好的3個(gè)結(jié)果進(jìn)行保留。

步驟5將動(dòng)態(tài)分子進(jìn)行旋轉(zhuǎn)，重復(fù)步驟3，4，直到所有旋轉(zhuǎn)結(jié)束。

步驟6將得到的所有結(jié)果進(jìn)行排序，并寫(xiě)入文件中，上限為10 000個(gè)數(shù)據(jù)。

步驟7按上述步驟完成受體與另外兩個(gè)動(dòng)態(tài)分子構(gòu)象的對(duì)接測(cè)試。

2.2 基于CUFFTW的優(yōu)化方案二

隨著GPU技術(shù)的不斷革新，NVIDIA公司提出了CUDA架構(gòu)，同時(shí)也針對(duì)傅里葉變換發(fā)布了兩種加速方案，一種是CUFFT加速，另一種是CUFFTW加速。

CUFFTW和CUFFT的區(qū)別主要在于，CUFFTW是針對(duì)FFTW庫(kù)開(kāi)發(fā)的，可以最快速、最便捷地為使用FFTW庫(kù)(方案一采用的加速方法)的用戶帶來(lái)GPU的加速效果，即只要原程序采用的是FFTW庫(kù)，那么通過(guò)CUFFTW可以直接在GPU端進(jìn)行運(yùn)行，不需要太多的改動(dòng)。而CUFFT則不同，它是針對(duì)所有傅里葉變換開(kāi)發(fā)的，是一種完全的CUDA接口。

優(yōu)化方案二在優(yōu)化方案一的基礎(chǔ)上，采用CUFFTW通過(guò)部分CUDA GPU接口進(jìn)行加速。換言之，方案二是將原來(lái)在CPU端運(yùn)行的程序移植到GPU端運(yùn)行，通過(guò)CUDA GPU架構(gòu)來(lái)對(duì)方案一進(jìn)行再加速。具體流程與方案一相同，只是運(yùn)行環(huán)境發(fā)生了改變。

圖1 基于GPU的分子對(duì)接流程圖 Fig.1 Flow chart of molecular docking based on GPU

2.3 基于CUFFT的優(yōu)化方案三

優(yōu)化方案一雖然程序確實(shí)得到了加速，但是加速效果并不理想；優(yōu)化方案二專門針對(duì)FFTW庫(kù)開(kāi)發(fā)，不是采用純CUDA GPU接口，若換成其他函數(shù)庫(kù)則無(wú)法使用。因此，考慮到CUFFT版本的加速效果，本文又提出優(yōu)化方案三，即在GPU端進(jìn)行純CUDA接口傅里葉變換，將耗時(shí)的快速傅里葉變換運(yùn)算和復(fù)合物點(diǎn)乘計(jì)算移植到GPU端，將所得結(jié)果返回CPU進(jìn)行篩選排序處理，并利用GPU的高速浮點(diǎn)運(yùn)算能力加速程序的運(yùn)行[12-13]。具體運(yùn)輸步驟如下：

步驟1根據(jù)所做的工作制定一個(gè)計(jì)劃。因?yàn)楸疚母道锶~變換是在三維空間進(jìn)行的，因此需要調(diào)用的函數(shù)是cufftPlan3d。其他低維情況可以參考CUFFT的官方文檔做相應(yīng)變動(dòng)。

步驟2調(diào)用相關(guān)函數(shù)開(kāi)始執(zhí)行計(jì)劃。因?yàn)楸疚乃龅氖请p精度的傅里葉變換和傅里葉逆變換，所以需要調(diào)用的函數(shù)分別是cufftExecR2C( )和cufftExecuC2R( )。

步驟3調(diào)用函數(shù)cufftDestroy( )釋放資源并銷毀計(jì)劃。

以一種動(dòng)態(tài)分子構(gòu)象計(jì)算為例，基于CUDA GPU的3D_Dock程序的具體流程如圖1所示。

3 分子對(duì)接實(shí)例及結(jié)果分析

3.1 對(duì)接實(shí)例

為了測(cè)試3D_Dock兩種優(yōu)化方案的有效性以及計(jì)算精度，在采用3D_Dock默認(rèn)對(duì)接參數(shù)的基礎(chǔ)上，將程序應(yīng)用于PDB code為1PEE，1B6C，4HX3和2SNI的結(jié)合態(tài)分子對(duì)接和自由態(tài)分子對(duì)接(圖2)，求取結(jié)合態(tài)和自由態(tài)模式下的前10名復(fù)合物的最小均方根偏差LRMSD值[14-15]，同時(shí)展示4個(gè)體系復(fù)合物預(yù)測(cè)的結(jié)構(gòu)圖，比較對(duì)接結(jié)果與原晶體結(jié)構(gòu)的差異，驗(yàn)證結(jié)果的正確性，如表1所示。

圖2 PDB code 分別為1PPE，1B6C，4HX3和2SNI的蛋白質(zhì)分子結(jié)構(gòu)(淡黃色為配體，藍(lán)色為受體)

Fig.2 Molecular structures of 1PPE, 1B6C, 4HX3 and 2SNI(PDB code)( Pale yellow shows the ligand and blue the receptor)
表1對(duì)接4個(gè)體系中具有最小LRMSD值的前10名結(jié)合態(tài)和自由態(tài)結(jié)果

Tab.1ThefirsttenboundandunboundresultsofdockinginstancewiththeminimumLRMSDvalue

PDB編號(hào)受體分子配體分子結(jié)合態(tài)LRMSD/?自由態(tài)LRMSD/?1PPEBovine trypsinCMTI-1 squash inhibitor2.093.771B6CFK506-BIND-INGTGF-B superfamily receptor type I2.182.244HX3Neutral pro-teinase inhibi-tor ScNPIExtracellular small neutral pro-tease0.933.732SNISubtilisinChymotrypsin in-hibitor 21.173.43

圖3 4個(gè)PDB 結(jié)合態(tài)結(jié)果在LRMSD最低值時(shí)與原晶體結(jié)構(gòu)對(duì)比圖 Fig.3 Comparison of four PDB bound results with the minimum values of LRMSD and the original crystal structure

評(píng)價(jià)方式采用蛋白質(zhì)復(fù)合物結(jié)構(gòu)預(yù)測(cè)競(jìng)賽(Critical assessment of predicted interactions, CAPRI)[16-17]組委會(huì)的指標(biāo)。在CAPRI實(shí)驗(yàn)中，正確的預(yù)測(cè)結(jié)果分為高精度結(jié)構(gòu)、中等精度結(jié)構(gòu)和可接受結(jié)構(gòu)[18-19]。根據(jù)LRMSD劃分情況為：當(dāng)LRMSD小于1 ?時(shí)為高精度結(jié)構(gòu)；當(dāng)LRMSD小于5 ?時(shí)為中等精度結(jié)構(gòu)；當(dāng)LRMSD小于10 ?時(shí)為可接受結(jié)構(gòu)。LRMSD表示預(yù)測(cè)的配體構(gòu)象與晶體結(jié)構(gòu)中配體結(jié)構(gòu)的均方根偏差。本文選取的標(biāo)準(zhǔn)為L(zhǎng)RMSD小于5 ?，即為中等精度下為近天然結(jié)構(gòu)。具體數(shù)據(jù)來(lái)源于CAPRI官方網(wǎng)站: http://www.ebi.ac.uk/msd-srv/capri/。

根據(jù)表1可以看出，結(jié)合態(tài)結(jié)構(gòu)的LRMSD值要比自由態(tài)結(jié)構(gòu)的LRMSD值小。其主要原因是自由態(tài)對(duì)接主鏈結(jié)構(gòu)有明顯變化，因此自由態(tài)的LRMSD值通常要比結(jié)合態(tài)對(duì)接的值大。而結(jié)合態(tài)結(jié)構(gòu)是直接由復(fù)合物晶體結(jié)構(gòu)拆分而來(lái)，因此在對(duì)接中更容易找到接近晶體結(jié)構(gòu)的復(fù)合物構(gòu)象。圖3是4個(gè)PDB結(jié)果在結(jié)合態(tài)結(jié)構(gòu)最低值時(shí)(根據(jù)表1結(jié)果)對(duì)應(yīng)的對(duì)接結(jié)果與原晶體結(jié)構(gòu)的對(duì)比圖，從圖中可以確認(rèn)對(duì)接結(jié)果的正確性。

3.2 性能分析

為了測(cè)試兩種方案對(duì)于程序的優(yōu)化效果，對(duì)上述4種蛋白均進(jìn)行了測(cè)試，并且加速比基本一致。為便于討論，選取受體Bovine trypsin和配體CMTI-1 squash inhibitor(PDB code:1PPE)的結(jié)合態(tài)和自由態(tài)對(duì)接作為優(yōu)化測(cè)試體系。在保持網(wǎng)格數(shù)一定的前提下，通過(guò)改變旋轉(zhuǎn)角度步長(zhǎng)來(lái)比較優(yōu)化的效果，分別取旋轉(zhuǎn)角度為12°，15°，18°進(jìn)行測(cè)試。表2，3為3種方案結(jié)合態(tài)和自由態(tài)結(jié)構(gòu)的測(cè)試結(jié)果，每個(gè)測(cè)試結(jié)果都是配體3種構(gòu)象的平均值，表中原方案在CPU上運(yùn)行，采用FFTW2函數(shù)庫(kù)。

表2基于結(jié)合態(tài)的4種方案運(yùn)行時(shí)間

Tab.2Runningtimeoffouroptimizationschemesbasedonboundstructure

旋轉(zhuǎn)角度步長(zhǎng)/(°)結(jié)合態(tài)分子對(duì)接方案運(yùn)行時(shí)間/s原方案優(yōu)化方案一優(yōu)化方案二優(yōu)化方案三129 805.993 553.281 207.79907.01154 909.901 807.44625.58476.89182 820.541 047.82365.27290.95

表3基于自由態(tài)的4種方案運(yùn)行時(shí)間

Tab.3Runningtimeoffouroptimizationschemesbasedonunboundstructure

旋轉(zhuǎn)角度步長(zhǎng)/(°)自由態(tài)分子對(duì)接方案運(yùn)行時(shí)間/s原方案優(yōu)化方案一優(yōu)化方案二優(yōu)化方案三128 842.643 150.291 058.62860.66154 717.491 691.28580.51474.74182 610.71949.71349.40288.23

表2中優(yōu)化方案一較原方案在運(yùn)行時(shí)間上有較大幅度的降低。這是由于原方案采用的快速傅里葉變換函數(shù)庫(kù)為FFTW2，而優(yōu)化方案一采用的函數(shù)庫(kù)為FFTW3。版本2的傅里葉變換主要有fftw.h和rfftw.h兩個(gè)頭文件，分別代表傅里葉變換和傅里葉逆變換；而版本3的傅里葉變換頭文件只有1個(gè)頭文件fftw3.h，即無(wú)論是正變換還是逆變換都集中于該頭文件。此外，在版本2的傅里葉變換中，原來(lái)的實(shí)數(shù)被雙精度浮點(diǎn)型替換，使得精度得到了較大的提高。而且原復(fù)數(shù)的實(shí)部和虛部是在結(jié)構(gòu)體中定義，而在版本3中直接進(jìn)行宏定義，這樣使用起來(lái)更加方便和快捷。同時(shí)，F(xiàn)FTW2和FFTW3的差別還在于策略的生成與執(zhí)行。FFTW3將傅里葉變換和反變換的重點(diǎn)從策略的制定轉(zhuǎn)移到策略的執(zhí)行。如果把計(jì)算都放在了策略的生成上面，那么每次生成策略都需要消耗一定的時(shí)間，而分子對(duì)接工作需要多次傅里葉變換計(jì)算，因此采用FFTW2函數(shù)庫(kù)運(yùn)行時(shí)間會(huì)大幅增加，而采用FFTW3則可以降低運(yùn)行時(shí)間。優(yōu)化方案二、三分別是在方案一的基礎(chǔ)上進(jìn)行的再升級(jí)，利用CUFFT的兩個(gè)子庫(kù)使得部分程序能夠在GPU端進(jìn)行運(yùn)算，從而降低了程序的運(yùn)行時(shí)間，較之方案一，兩者的提速比分別達(dá)到3倍和10倍。

由表3可知，自由態(tài)結(jié)構(gòu)和結(jié)合態(tài)結(jié)構(gòu)一樣，在改進(jìn)后的方案中速度都得到了提高，并且提高幅度相差不大。但是自由態(tài)結(jié)構(gòu)的運(yùn)行時(shí)間普遍比結(jié)合態(tài)結(jié)構(gòu)的運(yùn)行時(shí)間少，這主要是因?yàn)椴糠纸Y(jié)合態(tài)結(jié)構(gòu)計(jì)算所需的格點(diǎn)數(shù)要比自由態(tài)結(jié)構(gòu)多[20]，因此程序運(yùn)行消耗的時(shí)間也相應(yīng)會(huì)有所提高，但是增加的時(shí)間并不多，這也一定程度上證實(shí)了程序的正確性，以及優(yōu)化方法的正確性。

圖4，5是結(jié)合態(tài)結(jié)構(gòu)和自由態(tài)結(jié)構(gòu)中4種方案在旋轉(zhuǎn)角度為12°，15°，18°時(shí)的加速比結(jié)果，圖中所求得的加速比都是配體3種構(gòu)象測(cè)試結(jié)果的平均值。

通過(guò)表2，3及圖4,5可以得到兩方面結(jié)論。首先，當(dāng)旋轉(zhuǎn)步長(zhǎng)不斷增大時(shí)，程序的運(yùn)行速率得到了提高，這主要是因?yàn)楫?dāng)旋轉(zhuǎn)步長(zhǎng)增加時(shí)，需要旋轉(zhuǎn)采樣的次數(shù)隨之減少，這樣所需要的計(jì)算量便大幅降低，因此程序運(yùn)行時(shí)間自然減少。其次，如果旋轉(zhuǎn)采樣的次數(shù)減少，所需要使用的傅里葉變換次數(shù)便隨之減少，而本文所進(jìn)行的優(yōu)化大部分是針對(duì)傅里葉變換進(jìn)行的，一旦傅里葉變換計(jì)算次數(shù)減少，那么總體的加速比自然會(huì)下降。相反，如果旋轉(zhuǎn)步長(zhǎng)降低，程序的優(yōu)化效果則會(huì)得到較大的提高，原因在于旋轉(zhuǎn)次數(shù)增加，對(duì)接采樣的角度更加細(xì)化，使用快速傅里葉變換的計(jì)算次數(shù)也會(huì)增加。以上結(jié)論可以證明本文的優(yōu)化方案有效，通過(guò)對(duì)最主要的傅里葉變換算法進(jìn)行提速，使得程序運(yùn)行效率得到較大的提高。

Fig.5 Accelerated line chart of unbound structures using four schemes

4 結(jié)束語(yǔ)

本文主要針對(duì)3D_Dock中傅里葉變換提出了3種不同程度的優(yōu)化方案，分別是基于CPU的優(yōu)化和兩種基于GPU的優(yōu)化。從實(shí)驗(yàn)結(jié)果可以證明較之原程序而言，3種優(yōu)化方案都提高了程序的運(yùn)行速度，特別是基于GPU的方案，并行計(jì)算的優(yōu)化更大程度上提高了程序的運(yùn)行速度。這3種優(yōu)化方案對(duì)于蛋白質(zhì)結(jié)構(gòu)預(yù)測(cè)領(lǐng)域[21]有著重要的應(yīng)用，可以降低計(jì)算成本，提高預(yù)測(cè)效率。