羅　升，呂　強

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

?

距離約束的HMC采樣算法在蛋白質結構預測中的運用

羅升1，呂強2*

(1.蘇州大學計算機科學與技術學院，蘇州 2150062.蘇州大學江蘇省計算機信息處理技術重點實驗室，蘇州 215006)

摘要：蛋白質結構預測中，采樣是指在構象空間中生成具有最小自由能的狀態(tài)。傳統(tǒng)的采樣方法是對自由度直接賦值。這種方法在處理較少的殘基時能取得好的效果。但是對于包含100個殘基以上的蛋白質結構，由于構象空間的急劇增長，難以得到理想的結構。本文引入深度學習中的HMC(Hybrid Monte Carlo)采樣方法，以概率分布為依據(jù)對蛋白質的自由度進行采樣，能夠對包含100、200甚至更多個殘基的蛋白質結構進行采樣。并且，在采樣的過程中加入殘基間的距離約束，使得一個結構中，相對于Rosetta的ab initio最多有75%(平均40% )的殘基對得到優(yōu)化，滿足距離約束。

關鍵詞：距離約束；HMC；采樣；結構預測；蛋白質

近年來，從頭預測的方法在蛋白質結構預測領域取得了不錯的成績。但是，就目前的技術手段而言，從頭預測的方法仍然面臨兩大基本難點[1]：第一，由于蛋白質內部原子相互作用的復雜性，人們難以找到足夠準確的能量函數(shù)來描述一個蛋白質構象；第二，蛋白質結構的構象空間相當大，尤其在殘基序列較長的時候，如果沒有合適的、高效的采樣算法，蛋白質構象采樣將是一個“災難性”的計算問題。由此可見，采樣算法在從頭預測的方法中有著舉足輕重的地位[2]，因此，許多研究者針對這一問題開發(fā)了多種適用于蛋白質三維結構采樣的算法，比如改進的遺傳算法[3]、構象空間退火[4-5]和分子動力學模擬采樣[6]。

在蛋白質結構預測的語境下，采樣是指產生一個滿足要求的構象，而根據(jù)這個要求的不同，我們可以將蛋白質構象采樣分為狹義的采樣和廣義的采樣：

(1) 狹義的采樣：一個蛋白質構象可以通過不同的自由度來表示，常見的有二面角、原子坐標等。如果使用二面角表示蛋白質構象，那么這里的采樣就是指對二面角的一次賦值；如果使用原子坐標來表示蛋白質構象，那么這里的采樣就是指對原子坐標的一次賦值。Rosetta[7]中的每一個Mover都是這種意義下的采樣。

(2) 廣義的采樣：對于表示蛋白質構象的自由度，如果我們認為它們滿足某一特定的概率分布，比如正態(tài)分布N(d,u)，那么就可以對這個特定的概率分布進行一次采樣，同樣可以得到該自由度表示下的蛋白質的構象，而這樣的采樣就是廣義上的采樣。

對于狹義的采樣，由于是在自由度上進行不斷的賦值嘗試，再通過能量函數(shù)進行打分評估。不可避免的問題就是隨著殘基數(shù)量的增加，自由度的數(shù)量會成指數(shù)級增長，即使在現(xiàn)有的計算能力下，這樣的計算量增長對采樣算法來說，仍然是不可行的。而且，能量函數(shù)本身就是一個值得研究的問題。本文提出的基于距離約束模型的HMC(Hybrid Monte Carlo)采樣算法則是對概率模型的采樣，即使殘基數(shù)量增加，蛋白質序列變長，但是在參數(shù)數(shù)量上面臨的問題會比直接對自由度賦值的方法好許多。

HMC采樣[8-9]是一種對概率模型的采樣算法，其采樣機制通過引入Hamiltonian動力學系統(tǒng)[10]，避免在巨大的采樣空間中隨機采樣，并且可以在采樣的過程中保留概率分布所具有的統(tǒng)計規(guī)律。例如，在本文的案例中，如果我們使用原子坐標作為自由度來表示一個蛋白質結構，每個殘基選取4個原子，即N、Cα、C和O；每個原子有3個笛卡爾坐標，即x、y和z坐標。因此，表示一個殘基的就是12個自由度。對于一個長度為10殘基的蛋白質構象，用來表示它的自由度就有120個。進一步認為這120個自由度滿足一個120個隨機變量的概率分布，那么就可以用HMC算法對這個概率分布進行一次采樣，可以認為這是對該結構的一次“理性”的采樣。因為采樣的過程是模擬的Hamiltonian動力學系統(tǒng)，即粒子在空間中運行的狀態(tài)，那么如何保證蛋白質結構在采樣的過程中不會因為隨機的走動而散架呢？本文引入了距離約束模型。

蛋白質的空間結構并不僅僅是由一個個的殘基的獨立位置決定的。事實上，殘基間的相互作用對蛋白質結構有著巨大的影響[11]。在蛋白質結構中，按照殘基的間隔，一般可以分為三類距離約束[12]： 1) 短距離約束(Short-range)，此類約束的殘基間隔相對較短，范圍通常在6-12個殘基； 2) 中等距離約束(Medium-range)，此類約束的殘基間隔一般在12-24個殘基； 3) 長距離約束(Long-range)，此類約束的殘基間隔一般在24個殘基以上。同時，在蛋白質結構中，一般認為殘基空間歐氏距離在8 ?以內的殘基間是有相互作用的，因此本次實驗中使用的約束模型就是8 ?以內的長距離約束。

1材料與方法

1.1數(shù)據(jù)準備

蛋白質結構預測評估大會[13](CASP: Critical Assessment of techniques for protein Structure Prediction)是世界范圍內對蛋白質結構預測技術進行評估，評價當前蛋白質結構建模方法的能力及局限性的競賽，被稱為蛋白質結構預測領域的“奧運會”。本文從CASP8、CASP9和CASP10當中共選取了16個案例，它們的平均殘基長度為100～200，最短的包含101個殘基，最長的包含195個殘基，每個蛋白質都是提取的單鏈。

如上文所述，本文采用的自由度是原子三維坐標，每個殘基取其包括骨架原子Cα在內的4個原子的三維坐標，即每個殘基由12個值表示。每條蛋白質為一行，以300行為一個batch。因為本文的主要目的是驗證采樣算法的可行性，所以在數(shù)據(jù)準備的時候是從蛋白質的天然結構里獲取其距離約束的，而在實際應用的時候可以考慮使用同源蛋白質的距離約束。

從天然結構中獲取距離約束的過程如下：首先，遍歷天然結構的殘基，計算出每對殘基間隔大于24個殘基的殘基對的空間歐氏距離，然后找出小于8 ?的殘基對，記錄下這些殘基對以及距離值，作為約束。采樣的輸入數(shù)據(jù)是用Rosetta平臺ab initio得到的?？紤]到ab initio得到的結構具有較大的RMSD，可以更加顯著的看到采樣算法的效果，而且，采樣是為了得到更加符合距離約束的decoys，因此選用這樣的數(shù)據(jù)集作為采樣的輸入。

1.2HMC采樣

HMC采樣算法的核心在于通過引入Hamiltonian動力學系統(tǒng)，模擬粒子在空間的運動，從而對一個能量模型進行采樣。事實上，蛋白質的空間結構就是由空間中一個個粒子組成的，對蛋白質結構的采樣可以想象為讓這些粒子在空間中不斷的運動，以找到和天然結構最接近的構象。可見，如果不添加任何約束，那么采樣空間將接近于無窮大，因為粒子在空間里可以向任意方向以任意速度運動。面對如此巨大的采樣空間，任何采樣方法都是無能為力的。因此必須添加一定的約束從而一定程度上縮小采樣空間。

1.2.1Hamiltonian動力學系統(tǒng)

粒子在空間中的任何時刻都具有兩種屬性，即速度和位置，這兩種屬性對應于兩種能量：動能和勢能，我們把這樣的一個動力系統(tǒng)稱作Hamiltonian動力學系統(tǒng)。HMC采樣正是通過模擬這樣的一個系統(tǒng)，對一個高維空間進行采樣的。其能量函數(shù)的定義如公式(1)、(2)、(3)所示：

(1)

(2)

(3)

進行采樣，這是因為兩個概率分布p(s)和p(φ)是相互獨立的，邊緣化φ并不會產生很大的影響，并且可以重現(xiàn)原有的概率分布。

1.2.2蛙跳算法

因為模擬的時間是離散的，而粒子的運動是連續(xù)的，所以在實現(xiàn)過程中并不能很精確的模擬這樣一個系統(tǒng)。對于這個問題，有多種方案可以解決，其中比較好的是蛙跳算法[14-15]，可以通過三步操作解決這個問題，如公式(4)、(5)、(6)所示：

(4)

(5)

(6)

1.2.3Metropolis接受/拒絕

(7)

1.2.4HMC采樣流程

在本文的實驗中，勢能是根據(jù)蛋白質殘基的距離約束構建的模型，動能是一個隨機的均勻分布，整個采樣的流程如下：

(1) 從均勻分布中產生一個隨機的速度(動能)；

(2) 執(zhí)行n次蛙跳算法，獲取新的狀態(tài)x′(勢能 + 動能)，即對殘基的三維坐標進行n次調整，模擬粒子的運動，每次賦值都是通過計算勢能函數(shù)的梯度得到的，并不是隨機的賦值；

(3) 對新的狀態(tài)x′和舊的狀態(tài)χ執(zhí)行一次Metropolis接受判斷，是否可以接受；

(4) 根據(jù)當前接受率和目標接受率對步長進行一次調整，調整采樣效率；

(5) 重復以上步驟，直到達到特定的步數(shù)；

HMC采樣的過程中，有幾個重要的參數(shù)，比如步長、蛙跳次數(shù)、步長增益幅度等等，這些參數(shù)直接涉及到采樣的效果和速度(見表1)。

表1　采樣參數(shù)

實驗中，步長的初始值和變化范圍設置的比較小，這是可以理解的，因為蛋白質原子空間的變化本身就不是很大，太大幅度的變化可能導致空間嚴重變形，或者直接越過最佳位置。而且本文僅僅加了距離約束，沒有添加其它約束，如果設置過大，那么將會出現(xiàn)不可預期的結構。

1.3模型構建和實驗流程

實驗中使用的是殘基間隔大于24并且殘基空間歐氏距離小于8 ?的殘基對作為約束。在實驗的過程中，對于每一次的采樣結果，都會與天然結構進行計算。主要分為三步：

(1) 計算標記的每個殘基對的空間歐氏距離；

(2) 將這些距離與天然結構中對應殘基對的空間歐氏距離作差；

(3) 計算這些差值的最小平方和。

用這個值作為HMC采樣的勢能函數(shù)。這里，使用最小平方和作為模型輸出，是因為該方法本身具有的優(yōu)點：它可以找出與數(shù)據(jù)最佳的匹配函數(shù)，并且對于曲線有很好的擬合能力。在本實驗中，該模型如公式(8)～(11)所示：

(8)

(9)

其中，ΓD和ΓN分別表示decoys和天然結構中，兩個標記的殘基對之間的空間歐氏距離。然后對它們求差值：

Δi=ΓDi-ΓNi

(10)

最后，再將這些差值求平方和：

(11)

φ的大小反映了采樣得到的構象是否滿足從天然結構中獲取的距離約束。其值越小，說明標記的每對殘基對的空間歐氏距離和天然結構中對應的殘基對的空間歐氏距離越接近，因此我們也就認為在采樣的過程中，距離約束起到了作用，得到了更加滿足距離約束的構象。

按照上述模型，在本地機器上(CPU：Intel(R) Xeon(R) E5-2620 v2 @ 2.10Hz 內存：16GB)跑完HMC算法用時3.5 h，每個案例在6～26 min不等。整個采樣實驗流程如下：

(1) 使用Rosetta平臺，對每個案例用ab initio的方法得到300個decoys作為采樣的輸入；

(2) 從天然結構中計算得到距離約束，即包括殘基對和對應的空間歐氏距離；

(3) 對于ab initio得到的decoys，提取每個殘基骨架原子的三維坐標，構成一個行數(shù)為300，列數(shù)為序列長度×12的輸入矩陣；

(4) 將做好的矩陣輸入到HMC當中，進行多次的迭代采樣，并輸出最終得到的構象。

之所以選用Rosetta平臺制作輸入集，也是為了作為對比結果。Rosetta的ab initio就是一種不帶距離約束的采樣，本文使用這樣的decoys作為輸入，經過帶距離約束的HMC采樣，然后比較輸出decoys的距離約束，分析HMC采樣的效果。

2結果與分析

實驗一共使用了16個案例，每個案例用300個decoys作為采樣的輸入，每個案例采樣1 001次，其中有1 000次作為HMC的burn-in階段，丟棄。保留最后一次的結果。對比采樣之前的input-decoys和采樣之后得到的output-decoys，分別對每個decoys中標記的距離小于8?的殘基對的個數(shù)變化做了統(tǒng)計并比較，變化比例的計算如公式(12)所示：

(12)

實驗統(tǒng)計結果如表2所示，其中第一列是所選案例的PDB ID，二、三兩列分別是殘基和約束對的個數(shù)，最后三列則是變化的百分比。

表2　采樣后距離小于8 ?的殘基對個數(shù)變化比例

從表2中可以看出，在所有的案例中，采樣后的300個decoys里，小于8?的殘基對個數(shù)有了明顯的增長。如案例3JV6里，最大的一個構象有75.56%的殘基對在采樣之后空間歐氏距離小于8?。即3JV6的天然結構中，一共有104個殘基，長距離殘基對個數(shù)為：1+2+3+…+79+80=3 240。其中，空間歐氏距離小于8?的有132個。變化最大的一個decoy，在采樣的過程中，這132個殘基對里有132×0.755 6≈100個殘基對，由原來大于8?的距離變成小于8?的距離，由不滿足距離約束變?yōu)榱藵M足距離約束。這就可以說明，本次實驗中，HMC采樣成功的加入了距離約束，使一些原本大于8?的殘基對的距離小于8?。這樣的意義在于使原來沒有滿足距離約束的殘基對變得滿足距離約束，如此一來，這樣的構象就是比較好的構象，對下一步的三級結構預測有很大的意義。

另外，表2是根據(jù)殘基數(shù)從少到多的順序排列的，通過觀察數(shù)據(jù)還可以看出，殘基數(shù)對采樣效果的影響很少。殘基數(shù)較少的案例其采樣的效果并不一定好，例如1WWV、1WWW等；而殘基數(shù)較多的案例其采樣的效果也不一定會差，例如1HH4、2G0N、2WKP。這也印證了之前所說的，HMC采樣避免了因序列增長，自由度數(shù)量急劇增長而導致的采樣空間巨大，采樣算法不可行的問題。

同時，本文也計算了采樣前后蛋白質結構RMSD的變化，發(fā)現(xiàn)RMSD的變化很小，分析可能的原因有兩點：

第一，采樣的步長取值很小，也就意味著每一次采樣對三維坐標的變化很小。模擬到Hamiltonian動力學系統(tǒng)中就是指粒子在空間中每次移動的速度很慢，因此原子移動距離很小，而且，實驗中僅僅取了Cα，并沒有全原子采樣。事實上，如果步長較大，那么反而會出現(xiàn)不可預期的結果。雖然RMSD變化不大，但是在經過多次移動之后，仍然可以將大于8 ?的殘基距離縮小到8 ?以內?？梢姴蓸拥倪^程還是符合預期的結果；

第二，在采樣的過程中，標記出的約束對里，雖然有部分殘基對空間歐氏距離從大于8 ?降到8?以內，但也有部分從8?以內拉長到8?以外。但是，從整體相互作用的殘基對數(shù)變化的情況來看，采樣還是可以獲得比較好的構象的。事實上，采樣得到的decoys并不是都要的，我們只需要保留比較好的構象。本文對這部分數(shù)據(jù)也做了統(tǒng)計，如表格3所示。

從表3中可以看出，對于好的采樣案例，空間歐氏距離小于8?的殘基對個數(shù)是增加的，而且增加的不少。在統(tǒng)計的16個案例中，有3個案例出現(xiàn)了超過50%的增幅(3JV5、3JV6、2WKP)。表中還統(tǒng)計了前1%平均和前10%平均，可以說，這些構象是比較好的構象，也正是我們需要的。結合表2和表3可見，基于距離約束的HMC采樣算法在蛋白質結構預測方面是可行的。并且，相對于目前常用的Rosetta方法，能夠得到更加滿足距離約束的采樣結果，取得更好效果，對于三級結構預測有很大的幫助。

表3　標記的殘基對中距離小于8 ?的殘基對增加個數(shù)及每個案例用時

3結束語

本文以蛋白質結構中的距離約束為模型，使用HMC采樣對蛋白質空間三維坐標進行采樣。實驗已可以證明HMC采樣可以用在蛋白質結構預測方面，并且相對于目前常用的Rosetta有不錯的效果。因為模擬了Hamiltonian動力學系統(tǒng)，避免了其它采樣的隨機性和盲目性，所以即使同時加入了上百對距離約束，仍然能在計算能力可及的范圍內取得結果。然而，實驗本身并不是完美的，仍然有需要改進的地方。比如，在加入距離約束的同時，如何能降低蛋白質結構的RMSD。因為結構預測的最終目的是獲得近天然結構，距離約束只是為了更好的預測結構。因此，如何優(yōu)化模型，加入其它信息，是本文的下一個工作重點。

參考文獻

[1]HARDIN C, POGORELOV T V, Luthey-Schulten Z. Ab initio protein structure prediction[J].Current Opinion in Structural Biology,2002,12(2):176-181.

[2]KIM D E, BLUM B, BRADLEY P, et al. Sampling bottlenecks in de novo protein structure prediction[J].Journal of Molecular Biology,2009,393(1):249-260.

[3]JUDSON R S, COLVIN M E, MEZA J C, et al.Do intelligent configuration search techniques outperform random search for large molecules?[J].International Journal of Quantum Chemistry,1992,44(2):277-290.

[4]LEE J, SCHERAGA H A, RACKOVSKY S. New optimization method for conformational energy calculations on polypeptides: conformational space annealing[J].Journal of Computational Chemistry,1997,18(9):1222-1232.

[5]LEE J, SCHERAGA H A.Conformational space annealing by parallel computations: extensive conformational search of Met-enkephalin and of the 20-residue membrane-bound portion of melittin[J].International Journal of Quantum Chemistry,1999,75(3):255-265.

[6]LEE M R, TSAI J, BAKER D, et al.Molecular dynamics in the endgame of protein structure prediction[J].Journal of Molecular Biology,2001,313(2):417-430.

[7]ROHL C A, STRAUSS C E M, MISURA K M S, et al. Protein structure prediction using Rosetta[J]. Methods in Enzymology,2004,383: 66-93.

[8]DUANE S, KENNEDY A D, PENDLETON B J, et al. Hybrid monte carlo[J]. Physics Letters B, 1987, 195(2): 216-222.

[9]SEXTON J C, WEINGARTEN D H. Hamiltonian evolution for the hybrid monte carlo algorithm[J]. Nuclear Physics B, 1992, 380(3): 665-677.

[10]DIRAC P A M. Generalized hamiltonian dynamics[C] //Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences. The Royal Society, 1958, 246(1246): 326-332.

[11]WANG Z, XU J. Predicting protein contact map using evolutionary and physical constraints by integer programming[J]. Bioinformatics, 2013, 29(13): i266-i273.

[12]MONASTYRSKYY B, D’ANDREA D, FIDELIS K, et al. Evaluation of residue-residue contact prediction in CASP10[J]. Proteins: Structure, Function, and Bioinformatics, 2014, 82(S2): 138-153.

[13]KRYSHTAFOVYCH A, MONASTYRSKYY B, FIDELIS K. CASP prediction center infrastructure and evaluation measures in CASP10 and CASP ROLL[J]. Proteins,2014,82(S2):7-13.

[14]VAN GUNSTEREN W F, BERENDSEN H J C. A leap-frog algorithm for stochastic dynamics[J].Molecular Simulation,1988,1(3): 173-185.

[15]SNYMAN J A. The LFOPC leap-frog algorithm for constrained optimization[J].Computers & Mathematics with Applications,2000,40(8):1085-1096.

[16]NEAL R M, EDU E R T. Probabilistic inference using markov chain monte carlo methods[J].Physics and Chemistry, Nato Science Series C: Mathematical and Physical Sciences,1993,92(440):497-537.

[17]ULYANOV N B,SCHMITZ U,JAMES T L.Metropolis monte carlo calculations of DNA structure using internal coordinates and NMR distance restraints:an alternative method for generating a high-resolution solution structure[J].Journal of Biomolecular NMR,1993,3(5):547-568.

Distance constrains model based hybrid monte carlo sampling algorithm in protein structure prediction

LUO Sheng1，Lü Qiang2*

(1.SchoolofComputerScience&Technology,SoochowUniversity,Suzhou215006,China;2.ProvincialKeyLaboratoryforComputerInformationProcessingTechnologyofJiangsu,Suzhou215006,China)

Abstract：Sampling is defined as searching the conformational space for the status with the minimum free energy in protein structure prediction. In this paper, the Hybrid Monte Carlo(HMC) method from deep learning algorithms is introduced to better sample the conformational space of protein structures with 100, 200, or even more residues according to the probability distributions, while traditional sampling methods succeed in cases that proteins usually have less residues by assigning each value of free degrees directly. But they often fail the situation in which proteins have more than 100 residues, because of the large conformational space. In addition, residue distance constrains are added to the sampling algorithm to optimize a maximum 75 percent (40 percent on average) of residue pairs in each structure compare with ab initio in Rosetta.

Keywords：Distance constrain; HMC; Sampling; Structure prediction; Protein

收稿日期：2016-03-04；修回日期：2016-04-25.

基金項目：國家自然科學基金項目(No.61170125)。

作者簡介：羅升，男，碩士生，蘇州大學。研究方向：生物信息計算；E-mail: great_luo@hotmail.com. *通信作者：呂強，男，教授，蘇州大學。研究方向：生物信息計算、元啟發(fā)搜索、并行計算；E-mail: qiang@suda.edu.cn.

doi:10.3969/j.issn.1672-5565.2016.02.09

中圖分類號：TP391

文獻標志碼：A

文章編號：1672-5565(2016)02-117-06