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

楊 偉，呂 強(qiáng)

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

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

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

1 材料與方法

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

圖1 RosettaLigand對(duì)接算法Fig.1 Original dock algorithm in RosettaLigand

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

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

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

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

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

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

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

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

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

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

2 實(shí)驗(yàn)結(jié)果

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

表1 對(duì)照組與實(shí)驗(yàn)組候選結(jié)構(gòu)集合中Ligand-RMSD小于2.0的數(shù)量、最低能量值以及運(yùn)行時(shí)間統(tǒng)計(jì)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兩個(gè)區(qū)間的累積直方圖Fig.4 Stacked histogram of decoys with lrmsd fallen in interval 0～1 and 1～2

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

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

3 結(jié)論

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