預(yù)測(cè)蛋白質(zhì)折疊結(jié)構(gòu)的剪枝算法

陳昊，王代萍

(1.湖北大學(xué) 數(shù)學(xué)與計(jì)算機(jī)科學(xué)學(xué)院，湖北 武漢 430062；2.湖北大學(xué) 知行學(xué)院，湖北 武漢 430011)

0 引言

蛋白質(zhì)折疊結(jié)構(gòu)問題是生物信息學(xué)領(lǐng)域的核心問題之一.目前確定蛋白質(zhì)折疊結(jié)構(gòu)的一種途徑是通過物理化學(xué)的方法直接觀測(cè)，這種途徑能較容易測(cè)出蛋白質(zhì)鏈構(gòu)成，但要觀測(cè)蛋白質(zhì)的空間結(jié)構(gòu)卻十分困難.隨著計(jì)算機(jī)科學(xué)技術(shù)的進(jìn)步，人們開始尋求理論計(jì)算的方法直接預(yù)測(cè)蛋白質(zhì)的空間折疊結(jié)構(gòu).

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

①序列中每個(gè)小球的球心必須擺放在二維歐氏空間整數(shù)坐標(biāo)上.

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

③二維空間上的每個(gè)整數(shù)格點(diǎn)至多只能擺放一個(gè)小球，即任意兩個(gè)小球都不能重疊.

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

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

圖1 蛋白質(zhì)構(gòu)形(N=20,E=-9)

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

1 基于剪枝策略的啟發(fā)式搜索算法

圖2 剪枝策略示意圖

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

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

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

圖3 6個(gè)球的格局

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

定義2 自由度：對(duì)于格局n，當(dāng)?shù)趎+1個(gè)球放入時(shí)，其所有的合法位置數(shù)目稱為格局n的自由度，記作kfree.例如圖3中的格局6的kfree=3.

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

(1)

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

(2)

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

(3)～(5)

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

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

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

2 計(jì)算結(jié)果與分析

2.1實(shí)驗(yàn)結(jié)果用VC語言實(shí)現(xiàn)了剪枝算法PA，并對(duì)表1中列舉的10個(gè)典型算例進(jìn)行了實(shí)算，該組算例是當(dāng)今國際文獻(xiàn)中預(yù)測(cè)蛋白質(zhì)折疊結(jié)構(gòu)問題的公認(rèn)算例.剪枝算法PA運(yùn)行環(huán)境為Celeron 1.3 GHzCPU，512 M內(nèi)存微型機(jī).其他算法實(shí)驗(yàn)環(huán)境分別是：MSOE使用的是500 MHz Dec21164A Chip機(jī)器，SISPER使用的是 600 MHz Sony VAIO laptop機(jī)器，New ACO使用的機(jī)器是1 GHz CPU，PERM算例8使用的機(jī)器是Sun Ultra 1，算例9、10使用的機(jī)器是500 MHz Dec21164A，可以看出運(yùn)行5個(gè)算法的硬件性能相當(dāng).

圖4 N=64,E=-42

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

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

表1 國際公認(rèn)算例

E*指目前已發(fā)現(xiàn)的最低能量值

表2 計(jì)算結(jié)果

①②③：計(jì)算時(shí)間指算法運(yùn)行過程中首次出現(xiàn)表中所得到能量構(gòu)形的運(yùn)算時(shí)間,

④⑤：計(jì)算時(shí)間指平均得到一個(gè)最低能量構(gòu)形所用的時(shí)間,--表示未報(bào)道運(yùn)算時(shí)間.

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

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

同時(shí)我們發(fā)現(xiàn)剪枝算法中參數(shù)C>、C<的設(shè)定對(duì)算法運(yùn)行效率有較大的影響，參數(shù)的設(shè)定需綜合考慮解的質(zhì)量和運(yùn)算時(shí)間.一般而言，參數(shù)設(shè)定較大，則算法運(yùn)算時(shí)間較短，但卻很難找到最優(yōu)解;參數(shù)設(shè)定較小，則解的質(zhì)量一般很好，但算法運(yùn)算時(shí)間較長.

3 結(jié)論

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

參考文獻(xiàn)：

[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] 陳矛,黃文奇,呂志鵬.求解蛋白質(zhì)折疊問題的模擬退火算法[J].小型微型計(jì)算機(jī)系統(tǒng),2007,28(1):75-78.