預測蛋白質折疊結構的剪枝算法

陳昊，王代萍

(1.湖北大學 數學與計算機科學學院，湖北 武漢 430062；2.湖北大學 知行學院，湖北 武漢 430011)

0 引言

蛋白質折疊結構問題是生物信息學領域的核心問題之一.目前確定蛋白質折疊結構的一種途徑是通過物理化學的方法直接觀測，這種途徑能較容易測出蛋白質鏈構成，但要觀測蛋白質的空間結構卻十分困難.隨著計算機科學技術的進步，人們開始尋求理論計算的方法直接預測蛋白質的空間折疊結構.

針對蛋白質折疊結構，近年來科學家提出了一些簡化模型，目前國際上研究最為廣泛的是Dill和他的合作者提出的二維HP(hydrophobic-hydrophilic)格點模型[1].盡管蛋白質由20類氨基酸組成，但根據每類氨基酸的親水疏水特性的不同，分為兩大類，一類是疏水氨基酸(用H表示)，一類是親水氨基酸(用P表示)，因此任意的蛋白質鏈可表示為{H,P}上一個有限長的字符串，為表述方便，將P稱為白球，H稱為黑球，球的直徑均為1.一個合法的蛋白質空間構形必須滿足以下三個約束條件：

①序列中每個小球的球心必須擺放在二維歐氏空間整數坐標上.

②任意兩個在鏈上相鄰的小球擺放后必須在二維空間上也相鄰，即編號相鄰小球球心距離為1.

③二維空間上的每個整數格點至多只能擺放一個小球，即任意兩個小球都不能重疊.

滿足以上三個條件的蛋白質折疊結構稱為一個合法構形.任意合法構形都有其能量.能量定義只考慮疏水作用力，每一對在序列中不相鄰而在二維空間中相鄰的黑球產生能量為-1，其它情況能量計為0，計算合法構形中所有序列不相鄰而物理位置相鄰的H對的能量之和即是整個構形能量.一個鏈長為N的蛋白質合法構形能量E的形式化描述如下：

預測蛋白質折疊結構問題就是給定任意的氨基酸序列，找出滿足3個約束條件的最低能量構形.例如圖1是一個鏈長為20的蛋白質最優能量構形(圖中黑色方塊表示親水氨基酸H，白色方塊表示疏水氨基酸P)，能量為-9.二維格點模型是一種簡化模型，但這種簡化模型在一定程度上反映了蛋白質系統的一些基本特征，主要考慮了氨基酸的親水和疏水作用力，以最小能量作為優化指標，且該模型已被證實對預測蛋白質α螺旋結構有極高的可信度[2].

圖1 蛋白質構形(N=20,E=-9)

文獻[3]已證明基于HP格點模型的蛋白質折疊問題求解是NP難度的，因此依HP格點模型的求解蛋白質折疊問題大都尋求啟發式搜索算法.目前已提出的著名的近似算法有：基于重要性抽樣的SISPER算法[4]；基于Monte Carlo的MSOE算法[5]；基于蟻群優化思想的New ACO算法[6]；基于裁減復制策略的PERM算法[7-8]；基于模擬退火思想的SA算法[9].其中PERM算法是在目前公開發表的文獻中依據格點模型求解蛋白質折疊問題效率最高的近似求解算法，本文在PERM算法基礎上，結合擬物擬人思想，重新定義了權重、權重預測值和上下門限，得到了一種新的基于剪枝策略的啟發式搜索算法.

1 基于剪枝策略的啟發式搜索算法

圖2 剪枝策略示意圖

1.1剪枝算法(pruning algorithm)總體思想框架依據二維格點模型，可以將長度為N的蛋白質折疊構形生成過程理解為蛋白質鏈逐漸生長的過程.具體操作如下：考慮到空間對稱性，可先將第一，二號球固定，然后在第二號球周圍合法位置上放置第三號球，再在第三號球周圍合法位置上放置第四號球，依次類推，逐一放置黑白球直至放完N個球.可以用一棵搜索樹描述蛋白質構形的生成過程，由于從第三號小球開始，每個小球至多有3個合法的生長方向，因此一棵蛋白質折疊構形完全搜索樹的葉子結點個數大約為3N-2，顯然搜索空間隨鏈長的增加成指數級增長.

剪枝算法總體思想是控制分支繁殖的規模，只對有潛力的分支進行繁殖.具體是通過制定一定的評判準則，讓那些有前途的個體得以繁衍，而不具備發展前景的個體終止繁殖，從而有效地模擬了生物繁殖過程中優勝劣汰的自然規律，如圖2所示.對個體的評判準則可以理解為剪枝策略，如何制定合適的剪枝策略直接決定了算法的執行效率.

1.2剪枝算法PA中的主要術語為便于描述剪枝算法，先給出算法中幾個主要術語的定義.

圖3 6個球的格局

定義1 格局：已放入前n(2≤n≤N)個球時蛋白質鏈所處的狀態稱為一個格局，并且將含有n個球的格局簡稱為格局n.如圖3，稱為格局6.

定義2 自由度：對于格局n，當第n+1個球放入時，其所有的合法位置數目稱為格局n的自由度，記作kfree.例如圖3中的格局6的kfree=3.

定義3 權重：對于格局n，其優劣程度稱為該格局的權重，記作Wn.剪枝算法中對權重的采用如下遞歸定義：

(1)

其中，ΔEn為第n個球放入后新增的能量值；T為溫度常數，一般取值范圍為T∈[0.2,0.4].上述Wn實際上是表達前n個球所組成的部分鏈構形的能量值.

(2)

1.3剪枝算法PA描述以上分析可以看出鏈長為N的蛋白質構形生長過程對應于剪枝搜索樹的生成過程，要描述剪枝算法只需給出如何從任一格局n-1(2≤n-1≤N-1)，通過選擇分支將第n個球放下，從而發展出若干新格局n.

(3)～(5)

其中C>、C<為(0,1)中的常數，Zn為所有已產生格局n的權重Wn算術平均值.用Cn表示算法運行到某一時刻所有長度為n的格局的個數，則在計算過程中每得出一個新的格局，都需依據其權重Wn更新Zn和Cn，Zn?(CnZn+Wn)/(Cn+1),Cn?Cn+1.

剪枝算法運行時，需首先給Zn和Cn賦初始值，一般可隨機生成一個鏈長為N的合法構形，計算相應的Wn(1≤n≤N)，則Zn和Cn可初始化為Z1=W1,…,Zn=WN；C1=C2=…=CN=1，然后采用深度優先方式遞歸生成剪枝搜索樹，最后在所有生成的長度為N的構形中找出能量最低構形作為算法的最終輸出.

由于權重Wn直接決定了繁殖門限的大小，因此權重Wn是引導算法走向最優解的實質因素，考慮到構形能量大小是格點模型的唯一優化指標，剪枝算法中按公式(1)定義的權重Wn只與前n個球構成的部分鏈構形的能量值有關，這一權重的引導使得算法更加趨于問題的本質.

2 計算結果與分析

2.1實驗結果用VC語言實現了剪枝算法PA，并對表1中列舉的10個典型算例進行了實算，該組算例是當今國際文獻中預測蛋白質折疊結構問題的公認算例.剪枝算法PA運行環境為Celeron 1.3 GHzCPU，512 M內存微型機.其他算法實驗環境分別是：MSOE使用的是500 MHz Dec21164A Chip機器，SISPER使用的是 600 MHz Sony VAIO laptop機器，New ACO使用的機器是1 GHz CPU，PERM算例8使用的機器是Sun Ultra 1，算例9、10使用的機器是500 MHz Dec21164A，可以看出運行5個算法的硬件性能相當.

圖4 N=64,E=-42

國際上具有代表性的4個算法以及剪枝算法均在較短時間得到了鏈長N<64的算例的最優構形，因此本文只對當今公認的3個典型難例進行比較，計算結果如表2所示.從計算結果的質量相比較，只有PA和PERM找到了所有算例的最低能量構形，而MOSE、SISPER和NEW ACO未能計算出部分算例的最優構形.從計算時間效率比較，PA算法總體優于PERM，尤其對于鏈8，PA計算速度比PERM算法快90余倍.

以鏈8(N=64,E=-42)為例，該算例雖然鏈長不算太長，但是公認的難例，尤其對生長型算法.主要表現在該鏈的最優構形在生長的初期能量減少較慢，直到后期才減少較多能量.剪枝算法中諸參數取值：T=0.2,C>=0.01,C<=0.2.算法運行了93 min結束共找到了5個最低能量構形，平均18.6 min找到一個最優構形，如圖4所示，運行效率明顯優于PERM、MSOE和New ACO，SISPER則未能找到最低能量構形.

表1 國際公認算例

E*指目前已發現的最低能量值

表2 計算結果

①②③：計算時間指算法運行過程中首次出現表中所得到能量構形的運算時間,

④⑤：計算時間指平均得到一個最低能量構形所用的時間,--表示未報道運算時間.

表3 剪枝算法得到的能量分布

2.2剪枝算法PA特征分析PA是一種近似求解算法，通過對一棵完全搜索樹的部分分支進行剪枝從而使得算法時間復雜度降為多項式時間.表3以算例1(N=20,E=-9)為例，比較了完全搜索與剪枝算法的搜索空間大小，可以看出對于算例1剪枝算法搜索空間大小不到完全搜索空間大小的萬分之六.隨著鏈長的增長，剪枝算法搜索空間所占完全搜索空間的比重將越來越小.

同時我們發現剪枝算法中參數C>、C<的設定對算法運行效率有較大的影響，參數的設定需綜合考慮解的質量和運算時間.一般而言，參數設定較大，則算法運算時間較短，但卻很難找到最優解;參數設定較小，則解的質量一般很好，但算法運算時間較長.

3 結論

預測蛋白質折疊結構問題的HP格點模型是一個典型的NP難問題.目前求解該類問題大都采用一些通用的啟發式搜索算法，如局部搜索、遺傳算法、蟻群算法、退火算法等，但隨著問題規模增大，計算優度和效率都會明顯下降.本文學習了PERM算法中的剪枝思想，通過定義全新的權重、上下門限制定了一套符合自然界優勝劣汰規律的剪枝算法，實驗結果表明剪枝算法是高效的，對目前公認算例組中最難的3個算例，剪枝算法均在較短時間內找到了當前已知的最低能量構形.同時剪枝算法還為求解NP難問題提供了一種新的思路，尤其對于離散空間中NP難度的搜索問題更是有很好的借鑒作用.

參考文獻：

[1] Dill K A,Bromberg S,Yue K. Principles of protein folding:a perspective from simple exact models [J]. Protein Science,1995(4):561-602.

[2] Shih C T,Su Z Y,Gwan J F. The HP model,designability,and alpha:helices in protein structures[J]. Physical Review Letter,2000,84:384-389.

[3] Berger B,Leighton T. Protein folding in the hydrophobic-hydrophilic(HP) model is NP complete[J]. Journal of Computational Biology,1998,5(1):27-40.

[4] Junni L Z,Jun S L. A new sequential importance sampling method and its application to the two-dimensional hydrophobic:hydrophilic Model[J]. Journal of Chemical Physics,2002,117(7):3492-3498.

[5] Chikenji G,Kikuchi M,Iba Y. Mutlti-self-overlap ensemble for protein folding:ground state search and thermodynamics[J]. Physical Review Letter,1999,83(9):1886-1889.

[6] Shmygelska A,Hoos H H. An improved ant colony optimization algorithm for the 2D HP protein folding problem[M]//Canadian Conference on AI,Lecture Notes in Computer Science. Berlin:Springer Verlag,2003:2671-2683.

[7] Hsu H P,Mehra V,Nadler W,et al. Growth algorithms for lattice heteropolymers at low temperatures[J]. Journal of Chemical Physics,2003,118(1):444-452.

[8] Hsu H P,Mehra V,Nadler W,et al. A growth-based optimization algorithm for lattice heteropolymers [J]. Physical Review E,2003,68(2):1007-1011.

[9] 陳矛,黃文奇,呂志鵬.求解蛋白質折疊問題的模擬退火算法[J].小型微型計算機系統,2007,28(1):75-78.