基于pose共享的蛋白質(zhì)-配體構(gòu)象并行搜索算法

楊 偉，呂 強(qiáng)

(1.蘇州大學(xué)計算機(jī)科學(xué)與技術(shù)學(xué)院，江蘇 蘇州215006;2.蘇州大學(xué)江蘇省計算機(jī)信息處理技術(shù)重點實驗室，江蘇蘇州215006)

隨著X射線衍射以及NMR等技術(shù)的發(fā)展，越來越多的蛋白質(zhì)的三維結(jié)構(gòu)被測定出來，使基于受體結(jié)構(gòu)的計算機(jī)輔助藥物設(shè)計更具現(xiàn)實意義。分子對接(Molecular docking)方法是基于受體結(jié)構(gòu)的藥物設(shè)計(Structure-based drug design，SBDD)中的重要方法之一。分子對接是指兩個或多個分子通過幾何匹配和能量匹配相互識別的過程，通過分子對接可以從已有藥物分子庫篩選出有希望的先導(dǎo)化合物，避免了繁瑣的化學(xué)合成實驗過程，縮短了藥物研發(fā)周期，降低了研發(fā)成本［1］。分子對接還被用于重新評估已知藥物，發(fā)現(xiàn)已知藥物新的適應(yīng)癥［2］。近年來，科研人員把注意力轉(zhuǎn)向天然藥物，希望從中藥中開發(fā)出更為安全有效的新藥，中藥成分復(fù)雜，分子對接成為科學(xué)闡述中藥藥效物質(zhì)基礎(chǔ)和作用機(jī)理成為中藥研究的重要工具［3］。分子對接操作就是尋找配體與受體結(jié)合在受體活性位點處的低能構(gòu)象的過程。對接算法的性能依賴于兩個因素:能量函數(shù)的精度和復(fù)合物構(gòu)象搜索算法的性能。分子對接所涉及的搜索空間非常巨大，即便對柔性小分子，粗略估計其搜索空間至少含1030個解，要從中找出低能構(gòu)象必須借助于各種優(yōu)化算法。目前，已經(jīng)有許多優(yōu)化算法用來解決分子對接問題，典型的模擬退火算法、遺傳算法、禁忌搜索算法、蒙特卡羅方法以及這些算法的各種修正變種［4-6］。

常用的分子對接模擬軟件有 AUTODOCK［10］，RosettaLigand［7］，GLIDE［12］， DOCK［13］，F(xiàn)IexX［14］，GOLD［15］等。Rosetta是華盛頓大學(xué)開發(fā)的開源生物大分子建模軟件包，其中包含用于對蛋白質(zhì)和核酸進(jìn)行結(jié)構(gòu)預(yù)測、設(shè)計和重建模的工具，在Rosetta中，分子的結(jié)構(gòu)用pose數(shù)據(jù)結(jié)構(gòu)表示［4］。RosettaLigand是其中利用模擬退火算法進(jìn)行蛋白質(zhì)-配體柔性對接的程序工具。本文的研究基于Rosetta v3.4版，下文中提及RosettaLigand處均指Rosetta v3.4中的RosettaLigand對接程序。Rosetta源代碼可通過http://www.rosettacommons.org/網(wǎng)站獲取。

1 材料與方法

RosettaLigand的對接算法流程如圖1所示［6］。

圖1 RosettaLigand對接算法Fig.1 Original dock algorithm in RosettaLigand

在Rosetta框架內(nèi)，RosettaLigand對接協(xié)議被重復(fù)啟動N次(串行或并行)，生成指定數(shù)量的侯選結(jié)構(gòu)，稱為decoys。在這次N次軌跡中，RosettaLigand對接實例之間并不共享采樣信息，彼此之間完全獨立。我們認(rèn)為在多次對接實例之間共享pose信息可以幫助RosettaLigand更好的對接蛋白質(zhì)與配體。RosettaLigand搜索蛋白質(zhì)-配體復(fù)合物構(gòu)象的過程本質(zhì)是在一個在Rosetta全原子能量函數(shù)的指導(dǎo)下，在能量地形圖上搜索的過程。由于蛋白質(zhì)-配體復(fù)合物自由度非常大，而且已知能量函數(shù)并不足夠精確，所以不存在全局最優(yōu)結(jié)構(gòu)的精確表達(dá)式，通常對接實驗都會生成數(shù)目巨大的侯選結(jié)構(gòu)，以期能夠采樣到盡可能多的近天然結(jié)構(gòu)，然后通過基于能量若者機(jī)器學(xué)習(xí)的挑選方法將近天然結(jié)構(gòu)挑選出來［9］。所以對接程序的搜索構(gòu)象的能力依賴兩個因素，一是搜索的廣度，二是搜索的深度和充分性。RosettaLigand通過運行多次對接軌跡來解決搜索的廣度問題，我們現(xiàn)在希望通過在多個軌跡之間共享采樣信息來改善RosettaLigand的構(gòu)象采樣的深度和充分性問題，如圖2a，每一個點代表一次對接過程，在不共享采樣信息的情況下，每個對接實例只能在其附近空間進(jìn)行不充分的構(gòu)象搜索，而在圖2b中，大量的搜索實例某個進(jìn)行吸引到構(gòu)象空間最能較低處，并對此處進(jìn)行充分的采樣。

圖2 共享與非共享采樣信息搜索示意圖Fig.2 Searching pattern:sampling sharing vs non-sharing

RosettaLigand在對接的DockMCM階段使用MonteCarlo判斷接受或拒絕當(dāng)前采樣。在經(jīng)過修改的RosettaLigand算法中使用MonteCarlo判斷接受或拒絕采樣時，首先將其他對接實例以文件形式共享的最佳采樣結(jié)果讀入，與當(dāng)前采樣結(jié)果進(jìn)行比較，而不是與RosttaLigand內(nèi)置算法中實現(xiàn)的與本實例已采樣的最佳結(jié)果比較。經(jīng)過修改的DockCMC算法流程如圖3所示:

圖3 RosettaLigand DockCMC的改進(jìn)流程Fig.3 Modification to original RosettaLigand DockMCM

RosettaLigand原始對接算法在生成一個候選結(jié)構(gòu)的采樣過程中保有一個last_accepted_pose對象，用于存儲上一次模擬退火接受的采樣結(jié)果。在每一次新的采樣后，將當(dāng)前采樣結(jié)果與此last_accepted_pose比較:

如果當(dāng)前采樣結(jié)果優(yōu)于last_accepted_pose，則接受此次采樣，并用此采樣結(jié)果替換last_accepted_pose對象;如果如果當(dāng)前采樣結(jié)果不如last_accepted_pose，則按概率接受或拒絕當(dāng)前采樣結(jié)果，接受時同樣用當(dāng)前采樣結(jié)果替換last_accepted_pose，拒絕時用last_accepted_pose做為下次采樣的起點。

在共享pose的RosettaLigand對接算法中，同樣會保有一個 last_accepted_pose對象，與原始的RosettaLigand算法不同，在每一次新的采樣結(jié)束后，采樣結(jié)果不是與上一次接受的采樣結(jié)果比較，而是與共享的采樣結(jié)果進(jìn)行比較:如果當(dāng)前采樣被接受，則用當(dāng)前采樣結(jié)果替換共享的采樣結(jié)果;如果當(dāng)前采樣被拒絕，則用共享的采樣結(jié)果做為下一次采樣的起點。

修改過的采樣接受判定算法如下:

本文從meiler數(shù)據(jù)集［8］中選擇11個自對接的算例進(jìn)行實驗，這11個目標(biāo)是1AQ1，1DBJ，1DM2，1NJA，1NJE，1PB9，1PBQ，1Y1M，2AYR，2PRG 和4TIM。

配體的異構(gòu)體庫使用OMEGA2［11］生成，實驗組為16個并行進(jìn)程以共享采樣信息方式生成5 000個侯選結(jié)構(gòu)，對照組使用Rosetta平臺的內(nèi)置的MPI框架使用16個進(jìn)程生成5 000個侯選結(jié)構(gòu)。計算平臺為16核SMP Linux集群。

2 實驗結(jié)果

圖4所示為11個目標(biāo)的候選結(jié)構(gòu)中Ligand-RMSD落在0～1和1～2兩個區(qū)間的累積直方圖。從圖中可以看出，1AQ1、1DBJ、1PBQ、1Y1M、2AYR這5個目標(biāo)，共享pose算法極大的提高了侯選結(jié)構(gòu)集合中 Ligand-RMSD小于 2.0的比例(Ligand-RMSD小于2.0的構(gòu)象被認(rèn)為是近天然構(gòu)象);1NJA和1NJE這兩個目標(biāo)，共享 pose算法和RosettaLigand內(nèi)置算法表現(xiàn)都很差;而 1DM2、1PB9、2PRG、4TIM這4個目標(biāo)，雖然共享pose算法不如RosettaLigand算法，但性能相近。表1統(tǒng)計了11個目標(biāo)的候選結(jié)構(gòu)中Ligand-RMSD小于2.0的構(gòu)象的數(shù)量、構(gòu)象最低能量以及實驗運行時間。從表中可以看出，除了目標(biāo)1Y1M，對于其他所有目標(biāo)，共享pose算法都采樣到了比對照組更低的能量，而計算時間平均僅增加了不到10%。

表1 對照組與實驗組候選結(jié)構(gòu)集合中Ligand-RMSD小于2.0的數(shù)量、最低能量值以及運行時間統(tǒng)計Table 1 Comparison of number of decoys with lrmsd ＜2.0，lowest energy and duration between experimental and control group

圖4 候選結(jié)構(gòu)集合中Ligand-RMSD落在0～1和1～2兩個區(qū)間的累積直方圖Fig.4 Stacked histogram of decoys with lrmsd fallen in interval 0～1 and 1～2

圖5 11個目標(biāo)候選結(jié)構(gòu)總體能量值分布直方圖Fig.5 Histogram of Rosetta energy of all decoys of all 11 targets

圖5所示為11個目標(biāo)的5 000個侯選結(jié)構(gòu)的能量分布的直方圖。從圖中可以看出，除1Y1M外，對大多數(shù)目標(biāo)來說，與RosettaLigand內(nèi)置算法相比，共享pose算法使降低了侯選結(jié)構(gòu)整體的平均能量。

3 結(jié)論

RosettaLigand以多次隨機(jī)啟動方式來解決構(gòu)象搜索的廣度問題，此方式類似于在崎嶇的能量地形圖上隨機(jī)分布起始點，并在每一個起始點附近使用模擬退火進(jìn)行構(gòu)象搜索，搜索的性能依賴隨機(jī)啟動的次數(shù)與隨機(jī)程度。我們的方法在N個并行的對接實例之間共享采樣信息，當(dāng)某個進(jìn)程采樣到一個近天然構(gòu)象時，其他的對接實例會被吸引到全局最優(yōu)點附近，并對該處進(jìn)行充分采樣，這就是在目標(biāo)1AQ1、1DBJ、1PBQ、1Y1M、2AYR 中看到的情況，近天然構(gòu)象的數(shù)量遠(yuǎn)遠(yuǎn)超過對照組。并且，由于每個對接實例在每一次采樣后都與會N個實例當(dāng)前所采樣到的能量最低的構(gòu)象進(jìn)行比較，所以絕大多數(shù)實驗組生成的候選結(jié)構(gòu)集合的平均能量低于對照組的平均能量。當(dāng)然，我們也看到，目標(biāo)1DM2、1PB9、2PRG、4TIM的實驗組采樣到的近天然構(gòu)象數(shù)量低于對照組，分析其原因是共享信息一方面強(qiáng)化了并行對接實例在能量地形圖上某處的采樣能力，但它同時另一方面也使得陷入局部最小的問題更加嚴(yán)重，原本在獨立采樣的情況下可能采樣到近天然結(jié)構(gòu)的實例可能被誤導(dǎo)到非近天然區(qū)域。共享采樣信息是有用的，但要避免陷入局部最小。如同樣啟動N個并行對接實例，在每個對接實例采樣都收斂之后，輸出采樣結(jié)果，再比較當(dāng)前N個實例采樣所得結(jié)果，以能量最低的結(jié)果當(dāng)作N個并行對接實例下一輪對接的起始結(jié)構(gòu)，以此循環(huán)。這樣在保證搜索廣度的同時，又可兼顧搜索的深度。更進(jìn)一步的，在共享整體pose的前提下，搜索的廣度不可能超出非共享條件下搜索的廣度，原因是無論如何選擇共享的時機(jī)與策略，所共享的采樣信息都是當(dāng)前N個并行對接實例可能采樣到的構(gòu)象。所以，為了拓寬搜索的廣度，以共享采樣信息為基礎(chǔ)，結(jié)合并行對接實例自身的采樣結(jié)果來構(gòu)造超出當(dāng)前所有對接實例本次對接軌跡可能覆蓋的構(gòu)象空間的構(gòu)象是一條可行的道路。

References)

［1］ KAR S，ROY K.How far can virtual screening take us in drug discovery［J］.Expert Opinion on Drug Discovery，2013，8(3):245-261.

［2］ MA D L，CHAN D S H，LEUNG C H.Drug repositioning by structure-based virtual screening［J］.Chemical Society Reviews，2013，42(5):2130-2141.

［3］ 任潔，魏靜.分子對接技術(shù)在中藥研究中的應(yīng)用［J］.中國中醫(yī)藥信息雜志 ，2014，(1):123-125.REN Jie，WEI Jing.Application of molecular docking techniques in Chinese herbal medicine［J］.Chinese Journal of Information on Traditional Chinese Medicine，2014，(1):123-125.

［4］ KAUFMANN K W，LEMMON G H，DELUCA S L，et al.Practically useful:What the Rosetta protein modeling suite can do for you［J］.Biochemistry，2006，49:2987-2998.

［5］ 劉敏，曾濤，徐開闊，等.一種基于免疫遺傳算法的分子對接構(gòu)象搜索策略［J］.計算機(jī)研究與發(fā)展，2009，46(z2):597-601.LIU Min，ZENG Tao，XU Kaikuo，et al.A molecular docking＇s conformational search strategy based on immune genetic algorithm［J］.Journal of Computer Research and Development，2009，46(z2):597-601.

［6］ 李純蓮，王希誠，趙金城，等.藥物分子對接中的構(gòu)象搜索策略［J］.計算機(jī)與應(yīng)用化學(xué)，2004，21(2):201-205.LI Chunlian，WANG Xicheng，ZHAO Jincheng，et al.Drug molecular docking design based on optimal conformation search ［J］.Computers and Applied Chemistry，2004，21(2):201-205.

［7］ DAVIS I W，BAKER D.RosettaLigand docking with full ligand and receptor flexibility［J］.Journal of molecular biology，2009，385:381-392.

［8］ MEILER J，Baker D.ROSETTALIGAND:protein-small molecule docking with full side-chain flexibility［J］.Proteins，2006，65，538-584.

［9］ 楊凌云，呂強(qiáng).一種基于SVR的分辨近天然GPCR―配體構(gòu)象的方法［J］.生物信息學(xué)，2011，9(2):167-170.YANG Lingyu，Lü Qiang.A SVR-Based method for indentifying near-native GPCR-Ligand confromation decoys［J］.China Journal of Bioinformatics，2011，09(2):167-170.

［10］ MORRIS G M，GOODSELL D S，HALLIDAY D S，et al.Automated docking using a lamarckian genetic algorithm and and empirical binding free energy function［J］.J Comp Chem，1998，19:1639-1662.

［11］ HAWKINS P C，SKILLMAN A G，WARREN G L，et al.Conformer generation with OMEGA:Algorithm and validation using high quality structures from the protein databank and cambridge structural database［J］.Journal of Chemical Information and Modeling，2010，50(4):572-584.

［12］ FRIESNER R A，MURPHY R B，REPASKY M P，et al.Extra precision glide:Docking and scoring incorporating a modelofhydrophobic enclosure for protein-Ligand complexes［J］.J.Med.Chem.，2006，49，6177-6196.

［13］ LANG P T，BROZELL S R，MUKHERJEE S，et al.DOCK 6:combining techniques to model RNA-small molecule complexes［J］.RNA ，2009，15(6):1219-1230.

［14］ RAREY M，KRAMER B，LENGAUER T，et al.A fast flexible docking method using an incremental construction algorithm［J］.J.Mol.Biol，1996，261:470-489.

［15］ JONES G，WILLETT P，GLEN R C，et al.Development and validation of a genetic algorithm for flexible docking［J］.J.Mol.Biol.，1997，267:727-748.