本發(fā)明涉及一種生物學(xué)信息學(xué)、智能優(yōu)化、計算機應(yīng)用領(lǐng)域，尤其涉及的是，一種基于Rosetta局部增強的群體蛋白質(zhì)結(jié)構(gòu)預(yù)測方法。

背景技術(shù)：

蛋白質(zhì)是細(xì)胞功能的核心，與大部分核心生命過程息息相關(guān)。事實上，蛋白質(zhì)只有折疊成特定的三維結(jié)構(gòu)(即蛋白質(zhì)三級結(jié)構(gòu))之后才能產(chǎn)生其特定的生物學(xué)功能。因此，為了解蛋白質(zhì)的功能，就必須獲得其三維空間結(jié)構(gòu)，從而通過了解蛋白質(zhì)的三維結(jié)構(gòu)推動功能材料設(shè)計和新型藥物研制的發(fā)展，幫助人們理解生命的基本過程，包括對阿爾茲海默癥、帕金森病以及II型糖尿病等蛋白質(zhì)折疊病的認(rèn)識。

目前常用的蛋白質(zhì)結(jié)構(gòu)測定方法有X射線衍射和核磁共振(NMR)，這兩種方法雖然預(yù)測精度高，但是對于X射線衍射來說，難以培養(yǎng)晶體且晶體結(jié)構(gòu)測定的周期較長，核磁共振對樣品的需要量大、純度要求高，目前只能用于小分子蛋白質(zhì)結(jié)構(gòu)的測定。因此，以計算機為工具，利用適當(dāng)?shù)膬?yōu)化算法，直接通過氨基酸序列預(yù)測蛋白質(zhì)三維結(jié)構(gòu)，進而設(shè)計具有潛在藥物價值的新功能蛋白質(zhì)與多肽分子是生命科學(xué)領(lǐng)域需要解決的一個根本問題。該問題的最終解決關(guān)鍵在于：如何利用現(xiàn)有技術(shù)，設(shè)計一種高效的蛋白質(zhì)構(gòu)象空間優(yōu)化算法。

經(jīng)過40多年的發(fā)展，尤其是進入21世紀(jì)以來，分子動力學(xué)模擬(MD)、蒙特卡羅(MC)、構(gòu)象空間退火(CSA)、進化類優(yōu)化算法(EA)等隨機優(yōu)化算法在從頭預(yù)測領(lǐng)域得到了成功應(yīng)用；格點系統(tǒng)搜索(SGS)、分枝定界(BB)等確定性全局優(yōu)化算法，理論研究超前于其數(shù)值應(yīng)用，其極高的計算復(fù)雜度，限制了它們在中等規(guī)模以上蛋白構(gòu)象優(yōu)化方面的應(yīng)用。基于MC及CSA系列改進算法，Baker團隊開發(fā)的Rosetta從頭預(yù)測服務(wù)器、Zhang團隊開發(fā)的I-TASSER及QUARK從頭預(yù)測服務(wù)器目前已經(jīng)成為國際領(lǐng)先的預(yù)測軟件。上述方法在預(yù)測一些序列長度較短的小蛋白來說，能夠有效的得到三維結(jié)構(gòu)。然而，由于蛋白質(zhì)能量模型考慮分子體系成鍵作用以及范德華力、靜電、氫鍵、疏水等非成鍵作用，致使其形成的能量曲面極其粗糙，構(gòu)象對應(yīng)局部極小解數(shù)目隨序列長度的增加呈指數(shù)增長，對于這些傳統(tǒng)方法進行預(yù)測顯得力不從心，其原因在于極大的構(gòu)象搜索空間會導(dǎo)致算法在預(yù)測過程中搜索能力漸漸下降，同時群體的多樣性也變得越來越小，從而導(dǎo)致算法失去搜索的動力，影響最終的預(yù)測精度。

因此，現(xiàn)有的群體蛋白質(zhì)結(jié)構(gòu)預(yù)測方法在搜索能力和種群多樣性保持方面存在著缺陷，需要改進。

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

為了克服現(xiàn)有的群體蛋白質(zhì)結(jié)構(gòu)預(yù)測方法在搜索能力和種群多樣性方面的不足，本發(fā)明提出一種搜索能力強，且能夠保持種群多樣性的基于Rosetta局部增強的群體蛋白質(zhì)結(jié)構(gòu)預(yù)測方法。

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

一種基于Rosetta局部增強的群體蛋白質(zhì)結(jié)構(gòu)預(yù)測方法，所述方法包括以下步驟：

1)輸入待測蛋白質(zhì)的氨基酸序列信息；

2)初始化：設(shè)置種群規(guī)模NP，交叉概率CR，策略選擇因子CS，多樣性接受概率RS，Rosetta軌跡長度T，片段長度L₁，L₂；

3)根據(jù)序列信息以片段長度L₁進行隨機片段組裝生成初始構(gòu)象種群P＝{C¹,C²,...,C^NP}，其中，Cⁱ的表示當(dāng)前種群中的第i個構(gòu)象個體，并根據(jù)能量函數(shù)Rosetta Score0計算各構(gòu)象個體的能量，同時初始化迭代次數(shù)G＝0；

4)采用能量函數(shù)Rosetta Score0評價構(gòu)象的質(zhì)量，以片段長度L₁對初始種群中的每個構(gòu)象個體執(zhí)行軌跡長度為T的Rosetta局部增強，并計算每個構(gòu)象的特征向量；

5)對步驟4)中增強后的每個構(gòu)象個體Cⁱ,i∈{1,2,…,NP}作如下處理：

5.1)設(shè)置能量函數(shù)和片段長度：

5.1.1)如果當(dāng)前迭代次數(shù)0＜G≤G_max/3，則片段長度l＝L₁，且選用能量函數(shù)Rosetta Score1；

5.1.2)如果當(dāng)前迭代次數(shù)G_max/3＜G≤2G_max/3，則片段長度l＝L₁，且選用能量函數(shù)Rosetta Score2；

5.1.3)如果當(dāng)前迭代次數(shù)G＞2G_max/3，則片段長度l＝L₂，且選用能量函數(shù)Rosetta Score3

5.2)如果當(dāng)前迭代次數(shù)G為G_max/3的整數(shù)倍，則對以片段長度l對構(gòu)象個體Cⁱ執(zhí)行軌跡長度為T的Rosetta局部增強，并根據(jù)步驟5.1)中設(shè)置的能量函數(shù)進行評價；

5.3)計算目標(biāo)構(gòu)象Cⁱ的特征向量，以及Cⁱ與當(dāng)前種群中其他構(gòu)象之間的特征向量歐氏距離，并以最小距離為Cⁱ的多樣性值Dⁱ；

5.4)根據(jù)序列信息，利用DSSP得到待測蛋白的loop區(qū)域，并隨機生成一個0到1之間的隨機數(shù)p；

5.5)如果p<CS，則從當(dāng)前種群中選取三個互不相同的構(gòu)象個體C^a、C^b和C^c，其中，a≠b≠c≠i，a,b,c∈[1,NP]，從構(gòu)象個體C^b和C^c中各隨機選取一個片段替換C^a中對應(yīng)位置的片段，并從C^a中隨機選取一個不包含loop區(qū)

域的窗口進行片段組裝生成變異構(gòu)象C^mutant；

5.6)如果p≥1-CS，則選出當(dāng)前能量值最低的構(gòu)象個體C^best，并從當(dāng)前種群中選取兩個互不相同的構(gòu)象個體C^a和C^b，其中，a≠b≠i，a,b∈[1,NP]，從構(gòu)象個體C^a和C^b中各隨機選取一個片段替換C^best中對應(yīng)位置的片段，并從C^best中隨機選取一個不包含loop區(qū)域的窗口進行片段組裝生成變異構(gòu)象C^mutant；

5.7)隨機生成一個0與1之間隨機數(shù)p′，如果p′＞CR，則隨機選取一個loop區(qū)域，替換目標(biāo)構(gòu)象個體Cⁱ與變異構(gòu)象個體C^mutant在該區(qū)域的二面角，從而生成測試構(gòu)象C^trial，否則C^trial直接等于變異構(gòu)象C^mutant；以片段長度l對測試構(gòu)象個體C^trial執(zhí)行軌跡長度為T的Rosetta局部增強；

5.8)計算增強后測試構(gòu)象的特征向量，并計算測試構(gòu)象的特征向量與當(dāng)前種群中各構(gòu)象個體的特征向量之間的距離，以最小距離為測試構(gòu)象的多樣性值D^trial；

5.9)計算測試構(gòu)象C^trial的能量函數(shù)值E^trial，并進行如下操作：

5.9.1)如果E^trial小于當(dāng)前目標(biāo)構(gòu)象個體Cⁱ的能量函數(shù)值Eⁱ，則測試構(gòu)象C^trial替換目標(biāo)構(gòu)象Cⁱ；

5.9.2)如果E^trial大于當(dāng)前目標(biāo)構(gòu)象個體Cⁱ的能量函數(shù)值Eⁱ，且測試構(gòu)象的多樣性值D^trial大于目標(biāo)構(gòu)象的多樣性值Dⁱ，則隨機生成一個0與1之間隨機數(shù)，如果p″>RS，則測試構(gòu)象C^trial替換目標(biāo)構(gòu)象Cⁱ；

6)判斷是否滿足終止條件，若滿足則輸出結(jié)果并退出，否則返回步驟5)。

進一步，所述步驟2)中，設(shè)置最大迭代次數(shù)G_max，所述步驟6)中，對種群中的每個構(gòu)象個體都執(zhí)行完步驟5)以后，迭代次數(shù)G＝G+1，終止條件為迭代次數(shù)G達到預(yù)設(shè)最大迭代次數(shù)G_max。

本發(fā)明的技術(shù)構(gòu)思為：首先，將結(jié)構(gòu)預(yù)測中的整個算法搜索過程分為四個階段，對每個階段設(shè)置片段長度進行片段組裝，并選用不同的能量函數(shù)來衡量構(gòu)象個體的質(zhì)量；然后，基于二級結(jié)構(gòu)信息，采用不同的變異策略利用loop區(qū)域信息來生成測試構(gòu)象，并通過隨機交換loop區(qū)域信息實現(xiàn)交叉過程，保持種群多樣性，同時對每個階段的測試構(gòu)象和目標(biāo)構(gòu)象執(zhí)行Rosetta局部增強；最后，提取構(gòu)象的特征向量來衡量各構(gòu)象個體的多樣性，從而以能量函數(shù)為主要衡量指標(biāo)，并以多樣性為輔助衡量指標(biāo)來指導(dǎo)構(gòu)象種群更新。

本發(fā)明的有益效果表現(xiàn)在：一方面，基于二級結(jié)構(gòu)信息，根據(jù)loop區(qū)域的殘基操作來實現(xiàn)不同策略的測試構(gòu)象生成，并對每個測試構(gòu)象進行Rosetta局部增強，從而提高算法的搜索能力；其次，針對不同階段的Rosetta局部增強，采用不同的片段長度進行片段組裝，并采用不同的能量函數(shù)衡量構(gòu)象的質(zhì)量，從而提高搜索效率；另一方面，在選擇過程中，基于各構(gòu)象個體之間的特征向量距離來衡量構(gòu)象的多樣性，并將其作為輔助指標(biāo)來衡量構(gòu)象的質(zhì)量，從而在搜索過程充分保持種群多樣性，進而提高預(yù)測精度。

附圖說明

圖1是基于Rosetta局部增強的群體蛋白質(zhì)結(jié)構(gòu)預(yù)測方法的流程圖。

圖2是基于Rosetta局部增強的群體蛋白質(zhì)結(jié)構(gòu)預(yù)測方法對蛋白質(zhì)1AIL進行結(jié)構(gòu)預(yù)測時的構(gòu)象更新示意圖。

圖3是基于Rosetta局部增強的群體蛋白質(zhì)結(jié)構(gòu)預(yù)測方法對蛋白質(zhì)1AIL進行結(jié)構(gòu)預(yù)測時得到的構(gòu)象分布圖。

圖4是基于基于Rosetta局部增強的群體蛋白質(zhì)結(jié)構(gòu)預(yù)測方法對蛋白質(zhì)1AIL進行結(jié)構(gòu)預(yù)測得到的三維結(jié)構(gòu)圖。

具體實施方式

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

參照圖1～圖4，一種基于Rosetta局部增強的群體蛋白質(zhì)結(jié)構(gòu)預(yù)測方法，包括以下步驟：

1)輸入待測蛋白質(zhì)的氨基酸序列信息；

2)初始化：設(shè)置種群規(guī)模NP，交叉概率CR，策略選擇因子CS，多樣性接受概率RS，Rosetta軌跡長度T，片段長度L₁，L₂；

3)根據(jù)序列信息以片段長度L₁進行隨機片段組裝生成初始構(gòu)象種群P＝{C¹,C²,...,C^NP}，其中，Cⁱ的表示當(dāng)前種群中的第i個構(gòu)象個體，并根據(jù)能量函數(shù)Rosetta Score0計算各構(gòu)象個體的能量，同時初始化迭代次數(shù)G＝0；

4)采用能量函數(shù)Rosetta Score0評價構(gòu)象的質(zhì)量，以片段長度L₁對初始種群中的每個構(gòu)象個體執(zhí)行軌跡長度為T的Rosetta局部增強，并計算每個構(gòu)象的特征向量；

5)對步驟4)中增強后的每個構(gòu)象個體Cⁱ,i∈{1,2,…,NP}作如下處理：

5.1)設(shè)置能量函數(shù)和片段長度：

5.1.1)如果當(dāng)前迭代次數(shù)0＜G≤G_max/3，則片段長度l＝L₁，且選用能量函數(shù)Rosetta Score1；

5.1.2)如果當(dāng)前迭代次數(shù)G_max/3＜G≤2G_max/3，則片段長度l＝L₁，且選用能量函數(shù)Rosetta Score2；

5.1.3)如果當(dāng)前迭代次數(shù)G＞2G_max/3，則片段長度l＝L₂，且選用能量函數(shù)Rosetta Score3

5.2)如果當(dāng)前迭代次數(shù)G為G_max/3的整數(shù)倍，則對以片段長度l對構(gòu)象個體Cⁱ執(zhí)行軌跡長度為T的Rosetta局部增強，并根據(jù)步驟5.1)中設(shè)置的能量函數(shù)進行評價；

5.3)計算目標(biāo)構(gòu)象Cⁱ的特征向量，以及Cⁱ與當(dāng)前種群中其他構(gòu)象之間的特征向量歐氏距離，并以最小距離為Cⁱ的多樣性值Dⁱ；

5.4)根據(jù)序列信息，利用DSSP得到待測蛋白的loop區(qū)域，并隨機生成一個0到1之間的隨機數(shù)p；

5.5)如果p<CS，則從當(dāng)前種群中選取三個互不相同的構(gòu)象個體C^a、C^b和C^c，其中，a≠b≠c≠i，a,b,c∈[1,NP]，從構(gòu)象個體C^b和C^c中各隨機選取一個片段替換C^a中對應(yīng)位置的片段，并從C^a中隨機選取一個不包含loop區(qū)域的窗口進行片段組裝生成變異構(gòu)象C^mutant；

5.6)如果p≥1-CS，則選出當(dāng)前能量值最低的構(gòu)象個體C^best，并從當(dāng)前種群中選取兩個互不相同的構(gòu)象個體C^a和C^b，其中，a≠b≠i，a,b∈[1,NP]，從構(gòu)象個體C^a和C^b中各隨機選取一個片段替換C^best中對應(yīng)位置的片段，并從C^best中隨機選取一個不包含loop區(qū)域的窗口進行片段組裝生成變異構(gòu)象C^mutant；

5.7)隨機生成一個0與1之間隨機數(shù)p′，如果p′＞CR，則隨機選取一個loop區(qū)域，替換目標(biāo)構(gòu)象個體Cⁱ與變異構(gòu)象個體C^mutant在該區(qū)域的二面角，從而生成測試構(gòu)象C^trial，否則C^trial直接等于變異構(gòu)象C^mutant；以片段長度l對測試構(gòu)象個體C^trial執(zhí)行軌跡長度為T的Rosetta局部增強；

5.8)計算增強后測試構(gòu)象的特征向量，并計算測試構(gòu)象的特征向量與當(dāng)前種群中各構(gòu)象個體的特征向量之間的距離，以最小距離為測試構(gòu)象的多樣性值D^trial；

5.9)計算測試構(gòu)象C^trial的能量函數(shù)值E^trial，并進行如下操作：

5.9.1)如果E^trial小于當(dāng)前目標(biāo)構(gòu)象個體Cⁱ的能量函數(shù)值Eⁱ，則測試構(gòu)象C^trial替換目標(biāo)構(gòu)象Cⁱ；

5.9.2)如果E^trial大于當(dāng)前目標(biāo)構(gòu)象個體Cⁱ的能量函數(shù)值Eⁱ，且測試構(gòu)象的多樣性值D^trial大于目標(biāo)構(gòu)象的多樣性值Dⁱ，則隨機生成一個0與1之間隨機數(shù)，如果p″>RS，則測試構(gòu)象C^trial替換目標(biāo)構(gòu)象Cⁱ；

6)判斷是否滿足終止條件，若滿足則輸出結(jié)果并退出，否則返回步驟5)。

進一步，所述步驟2)中，設(shè)置最大迭代次數(shù)G_max，所述步驟6)中，對種群中的每個構(gòu)象個體都執(zhí)行完步驟5)以后，迭代次數(shù)G＝G+1，終止條件為迭代次數(shù)G達到預(yù)設(shè)最大迭代次數(shù)G_max。

本實施例序列長度為56的α/β折疊蛋白質(zhì)1GB1為實施例，一種基于Rosetta局部增強的群體蛋白質(zhì)結(jié)構(gòu)預(yù)測方法，其中包含以下步驟：

1)輸入待測蛋白質(zhì)的氨基酸序列信息；

2)初始化：設(shè)置種群規(guī)模NP＝100，交叉概率CR＝0.5，策略選擇因子CS＝0.5，多樣性接受概率RS＝0.5，Rosetta軌跡長度T＝1000，最大迭代次數(shù)G_max＝1200，片段長度L₁＝3，L₂＝9；

3)根據(jù)序列信息以片段長度L₁進行隨機片段組裝生成初始構(gòu)象種群P＝{C¹,C²,...,C^NP}，其中，Cⁱ的表示當(dāng)前種群中的第i個構(gòu)象個體，并根據(jù)能量函數(shù)Rosetta Score0計算各構(gòu)象個體的能量，同時初始化迭代次數(shù)G＝0；

4)采用能量函數(shù)Rosetta Score0評價構(gòu)象的質(zhì)量，以片段長度L₁對初始種群中的每個構(gòu)象個體執(zhí)行軌跡長度為T的Rosetta局部增強，并計算每個構(gòu)象的特征向量；

5)對步驟4)中增強后的每個構(gòu)象個體Cⁱ,i∈{1,2,…,NP}作如下處理：

5.1)設(shè)置能量函數(shù)和片段長度：

5.1.1)如果當(dāng)前迭代次數(shù)0＜G≤G_max/3，則片段長度l＝L₁，且選用能量函數(shù)Rosetta Score1；

5.1.2)如果當(dāng)前迭代次數(shù)G_max/3＜G≤2G_max/3，則片段長度l＝L₁，且選用能量函數(shù)Rosetta Score2；

5.1.3)如果當(dāng)前迭代次數(shù)G＞2G_max/3，則片段長度l＝L₂，且選用能量函數(shù)Rosetta Score3

5.2)如果當(dāng)前迭代次數(shù)G為G_max/3的整數(shù)倍，則對以片段長度l對構(gòu)象個體Cⁱ執(zhí)行軌跡長度為T的Rosetta局部增強，并根據(jù)步驟5.1)中設(shè)置的能量

函數(shù)進行評價；

5.3)計算目標(biāo)構(gòu)象Cⁱ的特征向量，以及Cⁱ與當(dāng)前種群中其他構(gòu)象之間的特征向量歐氏距離，并以最小距離為Cⁱ的多樣性值Dⁱ；

5.4)根據(jù)序列信息，利用DSSP得到待測蛋白的loop區(qū)域，并隨機生成一個0到1之間的隨機數(shù)p；

5.5)如果p<CS，則從當(dāng)前種群中選取三個互不相同的構(gòu)象個體C^a、C^b和C^c，其中，a≠b≠c≠i，a,b,c∈[1,NP]，從構(gòu)象個體C^b和C^c中各隨機選取一個片段替換C^a中對應(yīng)位置的片段，并從C^a中隨機選取一個不包含loop區(qū)域的窗口進行片段組裝生成變異構(gòu)象C^mutant；

5.6)如果p≥1-CS，則選出當(dāng)前能量值最低的構(gòu)象個體C^best，并從當(dāng)前種群中選取兩個互不相同的構(gòu)象個體C^a和C^b，其中，a≠b≠i，a,b∈[1,NP]，從構(gòu)象個體C^a和C^b中各隨機選取一個片段替換C^best中對應(yīng)位置的片段，并從C^best中隨機選取一個不包含loop區(qū)域的窗口進行片段組裝生成變異

構(gòu)象C^mutant；

5.7)隨機生成一個0與1之間隨機數(shù)p′，如果p′＞CR，則隨機選取一個loop區(qū)域，替換目標(biāo)構(gòu)象個體Cⁱ與變異構(gòu)象個體C^mutant在該區(qū)域的二面角，從而生成測試構(gòu)象C^trial，否則C^trial直接等于變異構(gòu)象C^mutant；以片段長度l對測試構(gòu)象個體C^trial執(zhí)行軌跡長度為T的Rosetta局部增強；

5.8)計算增強后測試構(gòu)象的特征向量，并計算測試構(gòu)象的特征向量與當(dāng)前種群中各構(gòu)象個體的特征向量之間的距離，以最小距離為測試構(gòu)象的多樣性值

D^trial；

5.9)計算測試構(gòu)象C^trial的能量函數(shù)值E^trial，并進行如下操作：

5.9.1)如果E^trial小于當(dāng)前目標(biāo)構(gòu)象個體Cⁱ的能量函數(shù)值Eⁱ，則測試構(gòu)象C^trial替換目標(biāo)構(gòu)象Cⁱ；

5.9.2)如果E^trial大于當(dāng)前目標(biāo)構(gòu)象個體Cⁱ的能量函數(shù)值Eⁱ，且測試構(gòu)象的多樣性值D^trial大于目標(biāo)構(gòu)象的多樣性值Dⁱ，則隨機生成一個0與1之間隨機數(shù)，如果p″>RS，則測試構(gòu)象C^trial替換目標(biāo)構(gòu)象Cⁱ；

6)當(dāng)對種群中的每個構(gòu)象都執(zhí)行了步驟5)以后，G＝G+1，若G>G_max則輸出結(jié)

果并退出，否則返回步驟5)。

以序列長度為56的α/β折疊蛋白質(zhì)1GB1為實施例，運用以上方法得到了該蛋白質(zhì)的近天然態(tài)構(gòu)象，最小均方根偏差為平均均方根偏差為預(yù)測結(jié)構(gòu)如圖4所示。

以上說明是本發(fā)明以1GB1蛋白質(zhì)為實例所得出的優(yōu)化效果，并非限定本發(fā)明的實施范圍，在不偏離本發(fā)明基本內(nèi)容所涉及范圍的的前提下對其做各種變形和改進，不應(yīng)排除在本發(fā)明的保護范圍之外。