基于遺傳算法預(yù)測2D三向的蛋白質(zhì)結(jié)構(gòu)

夏慧芳，郭雨珍，江宏昊

(南京航空航天大學(xué) 理學(xué)院數(shù)學(xué)系，南京 211106)

蛋白質(zhì)是生命活動(dòng)的重要承擔(dān)者，其空間結(jié)構(gòu)在很大程度上決定了它所具有的生物學(xué)功能，因此蛋白質(zhì)結(jié)構(gòu)的預(yù)測對于理解蛋白質(zhì)的結(jié)構(gòu)與功能之間的關(guān)系，并在此基礎(chǔ)上進(jìn)行蛋白質(zhì)復(fù)性、突變體設(shè)計(jì)以及基于結(jié)構(gòu)的藥物設(shè)計(jì)有著極其重要的意義[1]。蛋白質(zhì)分子是由二十多種氨基酸通過共價(jià)鍵連接而成的肽鏈形成，這些肽鏈?zhǔn)且罁?jù)什么原則形成具有一定空間結(jié)構(gòu)的蛋白質(zhì)分子，仍然是目前沒有解決的生物學(xué)問題[2]。隨著基因組測序工作的完成，生物學(xué)研究領(lǐng)域迫切需要找到一種從氨基酸序列出發(fā)，以此來預(yù)測蛋白質(zhì)結(jié)構(gòu)和功能的方法。在進(jìn)行蛋白質(zhì)結(jié)構(gòu)預(yù)測過程中，研究者提出了許多模型，最簡單的是Dill等人提出的HP格點(diǎn)模型[3-4],該模型將所有的氨基酸分為親水性(H)氨基酸和疏水性(P)氨基酸兩類，不考慮側(cè)鏈的影響，于是氨基酸序列被定義為一個(gè)由H和P組成的序列，這個(gè)序列遵循自回避原則，可以顯示在網(wǎng)格上。蛋白質(zhì)的天然構(gòu)象是吉布斯自由能最低的構(gòu)象是解決蛋白質(zhì)結(jié)構(gòu)預(yù)測問題的基礎(chǔ)。截止到現(xiàn)在，已經(jīng)有許多近似算法應(yīng)用在HP模型中，如粒子群算法[5-6]、神經(jīng)網(wǎng)絡(luò)算法[7]、遺傳算法[8-9]等,這些算法各有各的優(yōu)缺點(diǎn)，但至今還未發(fā)現(xiàn)一種算法完全好于其它算法。HP模型是一個(gè)偏理想化的模型，它需要將氨基酸鏈限制在正方形或矩形區(qū)域中，并且最大限度的將所有氨基酸只分為親水氨基酸和疏水氨基酸，但是有十幾種氨基酸并不能夠明確區(qū)分其疏水性及親水性，因此憑借HP模型來預(yù)測蛋白質(zhì)結(jié)構(gòu)并不符合實(shí)際。

疏水氨基酸相互作用，共價(jià)鍵和范德華力等會影響蛋白質(zhì)結(jié)構(gòu)的穩(wěn)定性，自然狀態(tài)下的蛋白質(zhì)有一個(gè)很緊湊的內(nèi)部結(jié)構(gòu)，范德華力在短程效應(yīng)中扮演著一個(gè)不可替代的角色，由范德華力方程式所產(chǎn)生的能量越大，蛋白質(zhì)結(jié)構(gòu)將會越緊湊。因此可以考慮基于范德華力勢能解決蛋白質(zhì)結(jié)構(gòu)的預(yù)測問題。

遺傳算法(Genetic Algorithm ,GA)是由美國密西根大學(xué)的Holland教授和他的學(xué)生在20 世紀(jì)60年代創(chuàng)立的[10]，該算法以遺傳機(jī)理和自然進(jìn)化為基礎(chǔ)，模擬了自然界中發(fā)生的自適應(yīng)現(xiàn)象，該算法被創(chuàng)立之后就被廣泛引用到工程問題中，現(xiàn)在已經(jīng)發(fā)展成為一種“自適應(yīng)啟發(fā)式概率性迭代式全局搜索算法”。目前，已被廣泛應(yīng)用到功能優(yōu)化、神經(jīng)網(wǎng)絡(luò)、機(jī)器學(xué)習(xí)、模式識別以及圖像處理[11]等領(lǐng)域。

本文剩余部分按如下安排：第二章中我們介紹了范德華力勢能預(yù)測蛋白質(zhì)結(jié)構(gòu)問題的數(shù)學(xué)模型，第三章中介紹了基本遺傳算法以及定義了調(diào)整策略，第四章中執(zhí)行數(shù)值實(shí)驗(yàn)并對結(jié)果進(jìn)行分析，最后在第五章中對整篇論文做了總結(jié)并對未來的研究做了展望。

1 數(shù)學(xué)模型

范德華力是分子間作用力，是由分子(原子)間相互接近造成的極化耦合引起的。范德華力勢能可由Lennard-Jones勢能函數(shù)如下表示：

對于蛋白質(zhì)結(jié)構(gòu)折疊問題，給出一個(gè)氨基酸序列，它被抽象為一個(gè)C原子鏈,氨基酸之間通過范德華力相互作用。本文中，只考慮基于范德華力的蛋白質(zhì)結(jié)構(gòu)折疊。我們知道范德華力產(chǎn)生的勢能越大，蛋白質(zhì)結(jié)構(gòu)就越緊湊，也即最大范德華勢能對應(yīng)的結(jié)構(gòu)是最優(yōu)蛋白質(zhì)構(gòu)象。

為了能找到穩(wěn)定的蛋白質(zhì)結(jié)構(gòu)，基于范德華勢能的蛋白質(zhì)結(jié)構(gòu)預(yù)測問題將會轉(zhuǎn)化為數(shù)學(xué)問題。因?yàn)橄噜彽膬蓚€(gè)氨基酸之間的距離不是零，所以任意兩個(gè)原子之間的范德華力勢能可由L-J勢能方程計(jì)算得出。于是，本問題的數(shù)學(xué)模型按如下表示：

化為標(biāo)準(zhǔn)形式：

氨基酸序列被抽象為C原子鏈，查閱文獻(xiàn)[12]得知εi=0.12 kcal/mol，rmin= 0.21nm，同時(shí)，我們規(guī)定在鏈上相鄰的兩個(gè)氨基酸之間的距離rij=0.52 nm。

模型中約束條件意味著任意兩個(gè)氨基酸之間的距離不是零，也即任何兩個(gè)氨基酸不會處于同一個(gè)位置并且只有一個(gè)氨基酸能占有該位置?，F(xiàn)在，一個(gè)生物問題轉(zhuǎn)化成了數(shù)學(xué)優(yōu)化問題，并且范德華力勢能拋棄了HP模型的局限性，能夠更真實(shí)的反映出蛋白質(zhì)的空間結(jié)構(gòu)。

2 遺傳算法

遺傳算法(Genetic Algorithm , GA)最初是由美國的Holland教授提出的模擬自然界中生物進(jìn)化機(jī)制的一種算法，它把達(dá)爾文進(jìn)化論和孟德爾遺傳學(xué)說作為基礎(chǔ)，仿照生物的進(jìn)化與遺傳過程，遵循適者生存和優(yōu)勝劣汰的規(guī)則，通過復(fù)制、交叉和變異等一系列操作，將需要解決的問題從初始解逐代逼近最優(yōu)解。

2.1 調(diào)整算子

氨基酸序列經(jīng)過交叉、變異操作后，后代可能會出現(xiàn)循環(huán)狀態(tài)，即相同的位置同時(shí)被兩個(gè)氨基酸占據(jù)。為了克服這個(gè)缺點(diǎn)，我們構(gòu)造了調(diào)整算子。

由于對氨基酸序列的編碼代表了方向，所以先根據(jù)初始點(diǎn)與編碼將每一個(gè)氨基酸的坐標(biāo)確定下來，接著從序列中第一個(gè)氨基酸開始檢驗(yàn)，若遇到序列中重復(fù)的氨基酸，則從當(dāng)前重復(fù)的氨基酸開始，向后調(diào)整直到最后一個(gè)點(diǎn)的無重復(fù)坐標(biāo)確定。

在進(jìn)行調(diào)整操作的過程中，可能會碰到一個(gè)點(diǎn)的所有方向都不可以取的情況，在數(shù)值實(shí)驗(yàn)時(shí)，就要定義一個(gè)記憶函數(shù)，每一個(gè)氨基酸都會對應(yīng)一個(gè)集合，這個(gè)集合記錄了這個(gè)氨基酸除了當(dāng)前方向還可以改變的其它方向。如果有一個(gè)氨基酸所有方向都會造成重疊，就要返回上一個(gè)氨基酸，當(dāng)前方向不可行，改變上一個(gè)氨基酸的方向，并且改變對應(yīng)的集合。同時(shí)，其它方向也不是隨意選取的，選取時(shí)是存在優(yōu)先級的。由于和初始點(diǎn)距離越近，氨基酸的序列就會更緊致，所以首先取其它所有可行方向中，對應(yīng)坐標(biāo)和初始點(diǎn)的距離最近的方向，最后得到不會發(fā)生重疊的序列。

2.2 改進(jìn)遺傳算法的步驟

對于預(yù)測蛋白質(zhì)結(jié)構(gòu)的優(yōu)化問題，改進(jìn)的遺傳算法按照如下步驟進(jìn)行：

Step1隨機(jī)編碼產(chǎn)生初始種群。本文編碼方式為：將“沿x軸正方向”設(shè)置為1，“與x軸正方向成120°”設(shè)置為2，“與x軸正方向成240°”設(shè)置為3。種群中隨機(jī)設(shè)置五個(gè)個(gè)體(氨基酸序列)，檢驗(yàn)每個(gè)氨基酸序列的有效性，如果是不合理序列，就要通過調(diào)整算子把它變?yōu)楹侠淼男蛄?，?jì)算每個(gè)序列的適應(yīng)度，規(guī)定適應(yīng)度為每個(gè)序列的范德華勢能。

Step2選擇。采用輪盤賭選擇，進(jìn)行交叉的個(gè)體被選擇的概率與它的范德華勢能成正比，進(jìn)行變異的個(gè)體被選擇的概率與它的范德華勢能成反比:

在本文中，每次循環(huán)選擇三個(gè)準(zhǔn)備進(jìn)行交叉的個(gè)體，選擇兩個(gè)準(zhǔn)備進(jìn)行變異的序列，于是，交叉如果能夠進(jìn)行就會產(chǎn)生六個(gè)新個(gè)體，變異能夠進(jìn)行則會產(chǎn)生兩個(gè)新個(gè)體。

Step3交叉。采用單點(diǎn)交叉，確定交叉概率pc=0.8,之后產(chǎn)生一個(gè)隨機(jī)概率r，且0

Step4變異。采用均勻變異，設(shè)置變異概率pm=0.05，之后產(chǎn)生一個(gè)隨機(jī)概率，如果隨機(jī)概率小于變異概率，則執(zhí)行變異操作，即對被選擇進(jìn)行變異的兩個(gè)個(gè)體，隨機(jī)一個(gè)變異位點(diǎn)，只改變這一個(gè)位點(diǎn)的編碼，變異規(guī)則按照如下方式：1→2,2→3,3→1，生成兩個(gè)新個(gè)體。和交叉操作一樣，為了防止新生成的個(gè)體不符合規(guī)則，也要對新生成的個(gè)體進(jìn)行檢驗(yàn)，不合理則調(diào)用調(diào)整算子，合理則繼續(xù)。

Step5適應(yīng)度評價(jià)。一次循環(huán)下來，通常都會產(chǎn)生新個(gè)體，計(jì)算新產(chǎn)生個(gè)體的范德華力勢能。

Step6種群更新。將新生成的個(gè)體的適應(yīng)度與父代進(jìn)行比較，如果子代個(gè)體中有個(gè)體的適應(yīng)度大于父代的適應(yīng)度，保存子代的最優(yōu)個(gè)體，淘汰父代中差的個(gè)體，總之要始終保持種群中有五個(gè)個(gè)體，在迭代過程中不斷更新種群。

Step7重復(fù)步驟 Step 1～Step 6，一直循環(huán)到5 000代，最后得到最優(yōu)解。

由于在遺傳算法的過程中，可能會出現(xiàn)局部最優(yōu)解的情況出現(xiàn)，所以為了克服這個(gè)缺陷，在進(jìn)行數(shù)值實(shí)驗(yàn)的過程中要重復(fù)進(jìn)行五次以上的實(shí)驗(yàn)取最優(yōu)解。

3 數(shù)值模擬

為了驗(yàn)證模型和改進(jìn)算法的有效性，進(jìn)行數(shù)值實(shí)驗(yàn)，分別預(yù)測氨基酸序列長度為15，17，20，25，30，35的蛋白質(zhì)結(jié)構(gòu)。

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

在進(jìn)行數(shù)值實(shí)驗(yàn)時(shí)，對于不同長度的氨基酸序列，我們都重復(fù)預(yù)測了五次，比較得出一個(gè)范德華勢能最大的構(gòu)象，結(jié)果見圖1所示。

進(jìn)行數(shù)值實(shí)驗(yàn)時(shí)，累加五次實(shí)驗(yàn)所得構(gòu)象的范德華勢能，求出平均勢能，同時(shí)記錄不同長度序列計(jì)算每代的運(yùn)行時(shí)間以及得到最優(yōu)解時(shí)平均運(yùn)行時(shí)間，結(jié)果分別如表1和表2所示。

從表1中可以看出，平均勢能與最大范德華勢能的誤差比較小，完全在可接受的范圍內(nèi)，這也反映出改進(jìn)后的遺傳算法的有效性。觀察表1和表2中的數(shù)據(jù)，我們推測：(1)序列越長，運(yùn)行時(shí)間會越長。(2)范德華勢能隨著序列長度的增加而增大。(3)應(yīng)用本文的方法，可以在可接受的時(shí)間里得到較長的序列的構(gòu)象。(4)蛋白質(zhì)的構(gòu)象越緊致，結(jié)構(gòu)會更穩(wěn)定。

3.2 最大范德華勢能擬合函數(shù)及誤差分析

通過觀察表1中范德華勢能與序列長度的數(shù)據(jù)，擬合得出能量與序列長度的關(guān)系函數(shù)及其函數(shù)圖像(見圖2)：

y=0.0052x2+0.9285x-0.0678

其中，x表示氨基酸序列長度，y表示對應(yīng)的范德華勢能。

分別用擬合函數(shù)和改進(jìn)遺傳算法計(jì)算了表1中序列的最大范德華勢能，比較結(jié)果如表3所示：

圖1 不同長度的氨基酸序列的二維拆疊構(gòu)象
Fig.1 Amino acid sequences with different lengths

表1 不同長度序列對應(yīng)的范德華勢能Table 1 Van der Waals potential energy corresponding to sequences of different lengths

表2 不同長度序列對應(yīng)的運(yùn)行時(shí)間Table 2 Running time corresponding tosequences of different lengths

圖2 能量與長度的關(guān)系Fig.2 Relationship between energy and length

氨基酸序列長度/個(gè)實(shí)際能量/(kJ·mol-1)擬合能量/(kJ·mol-1)誤差/(kJ·mol-1)1514.895 015.029 70.134 71716.817 017.219 50.402 52020.660 920.582 20.078 72527.387 726.394 70.993 03031.712 032.467 20.755 23538.919 338.799 70.119 6

從表3的誤差來看，擬合的效果非常接近程序的結(jié)果，這說明擬合函數(shù)是可以接受的。于是我們采用擬合函數(shù)分別預(yù)測序列長度為500,1 000和2 000的氨基酸序列的最大范德華勢能，結(jié)果見表4。

表4 能量擬合函數(shù)的預(yù)測結(jié)果Table 4 Prediction results of energy fitting function

從表4中獲知，當(dāng)氨基酸序列長度是500時(shí)，最大范德華力勢能是1764.2 kJ·mol-1；當(dāng)序列長度是1 000時(shí)，最大范德華勢能是6128.4 kJ·mol-1；當(dāng)序列長度是2 000時(shí)，最大范德華勢能是22 657 kJ·mol-1。我們發(fā)現(xiàn)，隨著氨基酸序列變長，其最大范德華勢能也會增大，證實(shí)了之前的猜測。

3.3 時(shí)間與長度擬合函數(shù)及誤差分析

我們通過觀察表2中程序總運(yùn)行時(shí)間和序列長度的數(shù)據(jù)，可以擬合得出運(yùn)行時(shí)間和序列長度的函數(shù)及其函數(shù)圖像(見圖3)：

y=0.0282x2+0.8656x+0.2150

其中，x表示序列長度，y表示總運(yùn)行時(shí)間。

圖3 運(yùn)行時(shí)間與長度關(guān)系Fig.3 Relationship between running time and length

分別用擬合函數(shù)和改進(jìn)遺傳算法計(jì)算了長度為20和30的序列的運(yùn)行時(shí)間，比較結(jié)果如表5。

從表5的誤差和誤差率來看，擬合的效果非常接近程序的結(jié)果，這說明擬合函數(shù)是可以接受的，于是我們采用該擬合函數(shù)預(yù)測序列長度為500,1 000和2 000的氨基酸序列的平均運(yùn)行時(shí)間，結(jié)果見表6。

當(dāng)序列長度是500時(shí)，平均運(yùn)行時(shí)間大約是2.1 h；當(dāng)序列長度是1 000時(shí)，平均運(yùn)行時(shí)間是 8.1 h；當(dāng)序列長度是2 000時(shí)，平均運(yùn)行時(shí)間是31.8 h。此結(jié)果說明基于范德華勢能預(yù)測蛋白質(zhì)結(jié)構(gòu)是可行的。

表6 時(shí)間擬合函數(shù)的預(yù)測結(jié)果Table 6 Prediction results of time fitting function

通過對實(shí)驗(yàn)數(shù)據(jù)的分析可以看到，基于范德華勢能的數(shù)學(xué)模型，通過改進(jìn)的遺傳算法來預(yù)測蛋白質(zhì)的空間結(jié)構(gòu)具有很大的可行性，最后得到的氨基酸序列的構(gòu)象是很緊湊的，因此是比較符合真實(shí)結(jié)構(gòu)的。

4 總結(jié)與展望

本文討論了基于范德華力的蛋白質(zhì)結(jié)構(gòu)預(yù)測問題。選擇范德華勢能作為數(shù)學(xué)優(yōu)化模型，變量是任意兩個(gè) C 原子之間的距離。目標(biāo)函數(shù)要求范德華勢能最大，約束條件是兩個(gè)氨基酸不占據(jù)同一個(gè)位置。選擇遺傳算法來解決此數(shù)學(xué)模型，并且對遺傳算法做了改進(jìn)。為了防止氨基酸的位置重疊，引入了調(diào)整算子的概念，使氨基酸序列最大程度的符合其真實(shí)的生物學(xué)特性。在數(shù)值實(shí)驗(yàn)中，改進(jìn)的遺傳算法搜索能力和搜索效率都得到了提高，證明了模型和算法的可行性和有效性。

在未來有很多方向可以追求。首先，本文研究的是二維平面上蛋白質(zhì)結(jié)構(gòu)預(yù)測問題，而真實(shí)的蛋白質(zhì)結(jié)構(gòu)是三維的，在以后的研究中可以考慮將模型和改進(jìn)的算法擴(kuò)展到空間蛋白質(zhì)預(yù)測問題中去。其次，可以將模擬的結(jié)果與真實(shí)的蛋白質(zhì)結(jié)構(gòu)進(jìn)行比較，檢測模型和算法的有效性。第三，還可以比較蛋白質(zhì)結(jié)構(gòu)預(yù)測的疏水親水模型和范德華勢能模型的結(jié)果，分析出各自的優(yōu)缺點(diǎn)。

總而言之，本文的模型和方法為蛋白質(zhì)結(jié)構(gòu)預(yù)測問題提供了相當(dāng)大的潛力。