本發(fā)明涉及一種生物學信息學、智能優(yōu)化、計算機應用領域，尤其涉及的是，一種基于粒子群優(yōu)化算法的蛋白質結構從頭預測方法。

背景技術：

蛋白質是所有生物中最重要的大分子之一，一個細胞的干重的一半以上是蛋白質。蛋白質是由單一的肽鍵連接在一起的氨基酸的順序鏈，這些連接的鏈折疊成能體現(xiàn)它們功能的三維結構，并調(diào)節(jié)生物體中的細胞活動，以保持生命。蛋白質分子在分子水平上具體體現(xiàn)了蛋白質結構和功能之間的顯著關系。蛋白質在生物學中有許多不同的功能。蛋白質的三維天然結構決定它的生物學機制，而反過來，蛋白質氨基酸單體的一維鏈又對它的天然結構進行編碼。蛋白質的三維結構在藥物設計、蛋白質工程和生物技術中具有決定性的重要，因此，蛋白質結構預測是一個重要的研究問題。

蛋白質結構實驗測定方法包括x射線晶體學、核磁共振光譜、和電子顯微鏡等，這些方法被廣泛用于蛋白質結構測定。x射線晶體學被認為是這些方法中相對可行且準確的測定方法之一。然而，x射線晶體學需要進行一個復雜的結晶過程，而對于一些不易結晶的蛋白(如膜蛋白)，此方法無法用于結構測定。此外，這些實驗測定方法極其耗時，代價昂貴，而且容易發(fā)生錯誤。

anfinsen熱力學假說表明：蛋白質的空間結構由氨基酸序列唯一確定；蛋白質的空間結構是穩(wěn)定的；蛋白質的天然構象處于自由能最低點。根據(jù)anfinsen原則，以計算機為工具，運用適當?shù)乃惴?，從氨基酸序列出發(fā)直接預測蛋白質的三維結構，成為當前生物信息學中一種主要的研究課題。在上世紀文獻中，通常將預測方法劃分為以下三類：針對高相似序列的同源建模方法；針對較低相似性序列的折疊識別方法；以及不依賴模板的從頭預測方法。從頭預測方法直接基于anfinsen假說建立蛋白質物理或知識能量模型，然后設計適當優(yōu)化算法求解最小能量構象。一方面，從生物學意義上來講，有助于揭示蛋白質折疊機理，進而能夠最終闡明生物學中心法則中的第二遺傳密碼理論部分；另一方面，從現(xiàn)實意義來講，該方法具有普遍性，對于序列相似度<20％或寡肽(<10個殘基的小蛋白)來說，從頭預測方法是唯一的選擇。

目前主要的技術瓶頸在于兩個方面，第一方面在于采樣方法，現(xiàn)有技術對構象空間采樣能力不強，另一方面在于構象更新方法，現(xiàn)有技術對構象的更新精度仍然不足。因此，現(xiàn)有的構象空間搜索方法存在不足，需要改進。

技術實現(xiàn)要素：

為了克服現(xiàn)有的蛋白質結構預測方法在預測精度和計算代價方面的不足，本發(fā)明提供一種基于粒子群優(yōu)化算法的蛋白質結構從頭預測方法。根據(jù)當前個體的二面角與歷史最優(yōu)個體的二面角的差值，以及當前個體的二面角與種群最優(yōu)個體的二面角的差值引導下一代構象搜索，并利用片段組裝進行隨機擾動，提高了全局探測能力，加快了收斂速度，降低了計算代價。

本發(fā)明解決其技術問題所采用的技術方案是：

一種基于粒子群優(yōu)化算法的蛋白質結構從頭預測方法，所述方法包括以下步驟：

1)選取rosettascore3作為蛋白質力場模型，即能量函數(shù)e(x)；

2)給定輸入序列信息；

3)參數(shù)初始化：設置粒子群規(guī)模np、慣性權重ω、加速度加速度和最大迭代次數(shù)gmax；

4)用二面角對表示蛋白質構象cⁱ(g)，其中，l為蛋白序列長度，表示第g代種群中第i個個體的第j個二面角對，為便于描述，將其記為其中表示構象cⁱ(g)的第k維分量，k∈{1,...,2l}；

5)構象初始化：根據(jù)給定輸入序列，生成np個伸展鏈構象，作為初始粒子群p＝{c¹(g),c²(g),...,c^np(g)}，令g＝0；

6)對粒子群中的構象cⁱ(g)依次進行l(wèi)次片段組裝，用e(x)計算新生成構象的能量值，并將當前構象作為個體cⁱ(g)的歷史最優(yōu)構象表示pbestⁱ的第k維分量，k∈{1,2,...,2l}，選取當前粒子群中能量最低構象作為粒子群全局最優(yōu)構象gbest＝{y1,y2,...,yk,...,y2l}，yk表示gbest的第k維分量，k∈{1,2,...,2l}；

7)對粒子群中的每個個體迭代的執(zhí)行下述步驟：

7.1)根據(jù)計算構象的各個殘基二面角增量其中ω是慣性權重因子，是加速度常數(shù)，均為非負值。rand(0,a1)和rand(0,a2)為[0,a1]、[0,a2]范圍內(nèi)具有均勻分布的隨機數(shù)，a1與a2為控制參數(shù)，表示第g代種群中第i個個體的第k維分量，k∈{1,2,...,2l}；

7.2)更新構象：根據(jù)計算每個構象的各個殘基二面角；

7.3)采用片段組裝對生成的個體執(zhí)行隨機擾動；

7.4)根據(jù)e(x)計算新生成構象的能量，并與其歷史最優(yōu)個體的能量進行比較，若能量降低，則更新個體cⁱ(g)的歷史最優(yōu)個體pbestⁱ；

7.5)新生成的構象能量與全局最優(yōu)個體gbest的能量進行比較，若能量降低，則用當前個體替換gbest；

8)判斷是否滿足終止條件：

8.1)若g＝gmax，則輸出結果并退出；

8.2)若g<gmax，則令g＝g+1，返回步驟7)。

本發(fā)明的技術構思為：根據(jù)當前個體的二面角與歷史最優(yōu)個體的二面角的差值，以及當前個體的二面角與種群最優(yōu)個體的二面角的差值引導下一代構象搜索，并利用片段組裝進行隨機擾動，提高了全局探測能力，加快了收斂速度，降低了計算代價。

本發(fā)明的有益效果為：本發(fā)明在蛋白質結構預測中應用，可以得到預測精度較高、復雜度較低的構象。

附圖說明

圖1是基于粒子群優(yōu)化算法的蛋白質結構從頭預測方法對目標蛋白質1acf預測結構和實驗室測定結構的三維示意圖。

具體實施方式

下面結合附圖對本發(fā)明作進一步描述。

參照圖1，一種基于粒子群優(yōu)化算法的蛋白質結構從頭預測方法，包括以下步驟：

1)選取rosettascore3作為蛋白質力場模型，即能量函數(shù)e(x)；

2)給定輸入序列信息；

3)參數(shù)初始化：設置粒子群規(guī)模np、慣性權重ω、加速度加速度和最大迭代次數(shù)gmax；

4)用二面角對表示蛋白質構象cⁱ(g)，

其中，l為蛋白序列長度，表示第g代種群中第i個個體的第j個二面角對，為便于描述，將其記為其中表示構象cⁱ(g)的第k維分量，k∈{1,...,2l}；

5)構象初始化：根據(jù)給定輸入序列，生成np個伸展鏈構象，作為初始粒子群p＝{c¹(g),c²(g),...,c^np(g)}，令g＝0；

6)對粒子群中的構象cⁱ(g)依次進行l(wèi)次片段組裝，用e(x)計算新生成構象的能量值，并將當前構象作為個體cⁱ(g)的歷史最優(yōu)構象表示pbestⁱ的第k維分量，k∈{1,2,...,2l}，選取當前粒子群中能量最低構象作為粒子群全局最優(yōu)構象gbest＝{y1,y2,...,yk,...,y2l}，yk表示gbest的第k維分量，k∈{1,2,...,2l}；

7)對粒子群中的每個個體迭代的執(zhí)行下述步驟：

7.1)根據(jù)計算構象的各個殘基二面角增量其中ω是慣性權重因子，是加速度常數(shù)，均為非負值。rand(0,a1)和rand(0,a2)為[0,a1]、[0,a2]范圍內(nèi)具有均勻分布的隨機數(shù)，a1與a2為控制參數(shù)，表示第g代種群中第i個個體的第k維分量，k∈{1,2,...,2l}；

7.2)更新構象：根據(jù)計算每個構象的各個殘基二面角；

7.3)采用片段組裝對生成的個體執(zhí)行隨機擾動；

7.4)根據(jù)e(x)計算新生成構象的能量，并與其歷史最優(yōu)個體的能量進行比較，若能量降低，則更新個體cⁱ(g)的歷史最優(yōu)個體pbestⁱ；

7.5)新生成的構象能量與全局最優(yōu)個體gbest的能量進行比較，若能量降低，則用當前個體替換gbest；

8)判斷是否滿足終止條件：

8.1)若g＝gmax，則輸出結果并退出；

8.2)若g<gmax，則令g＝g+1，返回步驟7)。

本實施例以目標蛋白質1acf為實施例，一種基于粒子群優(yōu)化算法的蛋白質結構從頭預測方法，所述方法包括以下步驟：

1)選取rosettascore3作為蛋白質力場模型，即能量函數(shù)e(x)；

2)給定輸入序列信息；

3)參數(shù)初始化：設置粒子群規(guī)模np＝100、慣性權重ω＝1.0、加速度加速度和最大迭代次數(shù)gmax＝100000；

4)用二面角對表示蛋白質構象cⁱ(g)，

其中，l＝125為蛋白序列長度，表示第g代種群中第i個個體的第j個二面角對，為便于描述，將其記為其中表示構象cⁱ(g)的第k維分量，k∈{1,...,2l}；

5)構象初始化：根據(jù)給定輸入序列，生成np個伸展鏈構象，作為初始粒子群p＝{c¹(g),c²(g),...,c^np(g)}，令g＝0；

6)對粒子群中的構象cⁱ(g)依次進行l(wèi)次片段組裝，用e(x)計算新生成構象的能量值，并將當前構象作為個體cⁱ(g)的歷史最優(yōu)構象表示pbestⁱ的第k維分量，k∈{1,2,...,2l}，選取當前粒子群中能量最低構象作為粒子群全局最優(yōu)構象gbest＝{y1,y2,...,yk,...,y2l}，yk表示gbest的第k維分量，k∈{1,2,...,2l}；

7)對粒子群中的每個個體迭代的執(zhí)行下述步驟：

7.6)根據(jù)計算構象的各個殘基二面角增量其中ω是慣性權重因子，是加速度常數(shù)，均為非負值。rand(0,a1)和rand(0,a2)為[0,a1]、[0,a2]范圍內(nèi)具有均勻分布的隨機數(shù)，a1＝1與a2＝1為控制參數(shù)，表示第g代種群中第i個個體的第k維分量，k∈{1,2,...,2l}；

7.7)更新構象：根據(jù)計算每個構象的各個殘基二面角；

7.8)采用片段組裝對生成的個體執(zhí)行隨機擾動；

7.9)根據(jù)e(x)計算新生成構象的能量，并與其歷史最優(yōu)個體的能量進行比較，若能量降低，則更新個體cⁱ(g)的歷史最優(yōu)個體pbestⁱ；

7.10)新生成的構象能量與全局最優(yōu)個體gbest的能量進行比較，若能量降低，則用當前個體替換gbest；

8)判斷是否滿足終止條件：

8.3)若g＝gmax，則輸出結果并退出；

8.4)若g<gmax，則令g＝g+1，返回步驟7)。

以上闡述的是本發(fā)明給出的一個實施例表現(xiàn)出來的優(yōu)良效果，顯然本發(fā)明不僅適合上述實施例，在不偏離本發(fā)明基本精神及不超出本發(fā)明實質內(nèi)容所涉及內(nèi)容的前提下可對其做種種變化加以實施。