CN105912885A - 一种基于群体Lipschitz下界估计的蛋白质结构预测方法 - Google Patents

Abstract

一种基于群体Lipschitz下界估计的蛋白质结构预测方法，首先，对整个初始构象种群构建Lipschitz下界估计支撑面，从而建立原能量函数模型的下界估计模型；然后，基于片段组装技术产生测试构象，进而根据下界估计模型获取测试构象的能量下界估计值，根据能量下界估计值判断是否需要对测试构象进行实际能量函数评价，并指导种群更新，从而有效减少目标函数评价次数，降低计算代价；最后，根据进化信息对较优的测试构象构建Lipschitz下界估计支撑面，使得下界估计模型向原能量函数模型不断收紧，从而获得更加精确的下界估计信息。

Description

一种基于群体Lipschitz下界估计的蛋白质结构预测方法

技术领域

本发明涉及一种生物学信息学、智能优化、计算机应用领域，尤其涉及的是，一种基于群体Lipschitz下界估计的蛋白质预测方法。

背景技术

生物细胞中包含许多蛋白质(由20多种氨基酸所形成的长链)，这些大分子对于完成生物功能至关重要。蛋白质分子从分子水平上揭示了其结构和功能之间的重要关系。不同的蛋白质在生物体中有着不同的功能，而蛋白的功能通常由其空间结构直接决定。尤其是蛋白质的三维结构是理解和转换生物功能的关键。因此，蛋白质三维结构预测对新蛋白的设计、药物制造、蛋白质之间的相互作用建模以及蛋白质稳定性预测至关重要。

蛋白质结构预测问题是计算生物学领域的一项挑战性工作。蛋白质三维结构可以通过一些实验手段测得，例如，核磁共振和X光晶体衍射。然而由于实验测定比较耗时和昂贵，对于某些不易结晶的蛋白质来说这两种实验方法不适用。蛋白质结构从头预测方法仅通过所给氨基酸序列就可以预测得到三维结构，因此被广泛应用于各种蛋白质的结构预测。在从头预测中，能量函数被用来评价构象的质量，并且通过有效的算法来搜索构象。由于从头预测方法遵循Anfinsen的热力学假说，即能量函数的全局最小解所对应的构象被认为就是天然态结构，因此能量函数在搜索过程中引导算法向查询序列的低能量构象搜索。可以看出，蛋白质结构从头预测方法涉及能量函数的优化问题。如何设计有效的方法来搜索能量函数的全局最小解是一项挑战性工作。

为了在构象空间中采样得到低能量构象，国内外专家相继提出了各种优化技术。这些技术包括进化算法(EAs)、蒙特卡洛(MC)、分子动力学(MD)和构象空间退火(CSA)等。这些方法在蛋白质结构预测中得到了广泛应用，同时也达到了一定的效果。然而，蛋白质结构预测的能量函数计算代价通常很昂贵，该能量函数自由度很高，能量函数曲面极其复杂，通常包含了大量的局部最优解和不可行构象区域。有时对其评价还需要调用到第三发能量包，导致评价一次达到数秒，甚至数分钟。从而导致算法计算代价极高，搜索效率低。

因此，现有的构象空间优化方法在计算代价和搜索效率方面存在着缺陷，需要改进。

发明内容

为了克服现有的蛋白质结构预测方法在计算代价和搜索效率方面的不足，本发明提出一种计算代价较低、搜索效率较高的基于群体Lipschitz下界估计的蛋白质结构预测方法。

本发明解决其技术问题所采用的技术方案是：

一种基于群体Lipschitz下界估计的蛋白质结构预测方法，所述方法包括以下步骤：

1)给定输入序列信息；

2)设置参数：种群大小NP，交叉因子CR，常数M，最大迭代次数；

3)种群初始化：由输入序列产生初始构象种群 其中，N表示维数，表示第i个构象Cⁱ的第N维元素，并初始化迭代次数G＝0；

4)对初始种群中的每一个构象Cⁱ构建Lipschitz下界估计支撑面lⁱ：

$l^i=(\frac{f\left(C^i\right)}{M}-x_1^i,\frac{f\left(C^i\right)}{M}-x_2^i,...,\frac{f\left(C^i\right)}{M}-x_{N+1}^i)---\left(1\right)$

其中，f(Cⁱ)为构象Cⁱ的能量函数值，为辅助变量，M为常数；

5)针对种群中的每个构象个体Cⁱ，i∈{1,2,3,…,NP}，令C_target＝Cⁱ，C_target表示目标构象个体，执行以下操作得到变异构象C_mutant：

5.1)随机生成正整数rand1,rand2,rand3∈{1,2,3,......NP}，且rand1≠rand2≠rand3≠i；再生成4个随机整数randrange1,randrange2,randrange3,randrange4；其中randrange1≠randrange2，randrange3≠randrange4∈{1,2,…,L},L为序列长度；

5.2)令a＝min(randrange1,randrange2)，b＝max(randrange1,randrange2)，k∈[a,b]；令c＝min(randrange3,randrange4)，d＝max(randrange3,randrange4)，p∈[c,d]；其中min表示取两个数的最小值，max表示取两个数的最大值；

5.3)用C_rand2上位置a到位置b的片段的氨基酸k所对应的二面角phi、psi、omega替换C_rand1的相同位置所对应的二面角phi、psi、omega；再使用C_rand3上位置c到位置d的片段的氨基酸p所对应的二面角phi、psi、omega替换C_rand1上相同位置所对应的二面角phi、psi、omega，然后将所得C_rand1进行片段组装得到变异构象个体C_mutant；

6)对变异构象C_mutant执行交叉操作生成测试构象C_trial：

6.1)生成随机数rand4，rand5，其中rand4∈(0,1)，rand5∈(1,L)；

6.2)根据式(2)执行交叉过程：

$C_{trail}=\left\{\begin{array}{cc}{C}_{mu\tan t,rand5}\←C_{t\arg et,rand5},& if(rand4\≤CR)\\ C_{mu\tan t,rand5},& otherwise\end{array}\right.---\left(2\right)$

其中C_mutant,rand5表示变异构象中的rand5片段，若随机数rand4≤CR，变异构象C_mutant的片段rand5替换为目标构象C_target中对应的片段，否则C_trial直接等于变异构象C_mutant；

7)计算测试构象C_trial的下界估计值

$\stackrel{\‾}{f}\left(C_{trial}\right)=\underset{k\≤NP}{\max }\underset{i=1,\mathrm{...}N+1}{\min }M(l_j^k+x_{trial}^j)---\left(3\right)$

其中max表示求最大值，min表示求最小值，为测试构象向量的第j维元素，为支撑向量l^k的第j维元素；

8)根据如下操作决定测试构象C_trial是否替换目标构象C_target：

8.1)计算目标构象的能量函数值f(C_target)；

8.2)如果则目标构象C_target保持不变，并直接进入下一次迭代；

8.3)如果则计算测试构象C_trial的实际能量函数值f(C_trial)，如果f(C_trial)＜f(C_target)，则测试构象C_trial替换目标构象C_target并继续步骤8.4)；

8.4)根据式(1)对测试构象C_trial构建Lipschitz下界估计支撑面；

9)判断是否满足终止条件，若满足则输出结果并退出，否则返回步骤5)。

进一步，所述步骤9)中，对种群中的每一个个体都执行完步骤5)—8)以后，迭代次数G＝G+1，终止条件为迭代次数G达到预设最大迭代次数。

本发明的技术构思为：首先，对整个初始构象种群构建Lipschitz下界估计支撑面，从而建立原能量函数模型的下界估计模型；然后，基于片段组装技术产生测试构象，进而根据下界估计模型获取测试构象的能量下界估计值，根据能量下界估计值判断是否需要对测试构象进行实际能量函数评价，并指导种群更新，从而有效减少目标函数评价次数，降低计算代价；最后，根据进化信息对较优的测试构象构建Lipschitz下界估计支撑面，使得下界估计模型向原能量函数模型不断收紧，从而获得更加精确的下界估计信息。

本发明的有益效果表现在：一方面，根据下界估计值指导种群更新，减少不必要的能量函数评价，降低算法计算代价，提高搜索效率；另一方面，根据进化信息更新下界支撑面，使得下界估计模型向原能量函数模型不断收紧。

附图说明

图1是基于群体Lipschitz下界估计的蛋白质结构预测方法对蛋白质1FD4进行结构预测时的构象更新示意图。

图2是基于群体Lipschitz下界估计的蛋白质结构预测方法对蛋白质1FD4进行结构预测时得到的构象分布图；

具体实施方式

下面结合附图对本发明作进一步描述。

参照图1和图2，一种基于群体Lipschitz下界估计的蛋白质结构预测方法，包括以下步骤：

1)给定输入序列信息；

2)设置参数：种群大小NP，交叉因子CR，常数M，最大迭代次数；

3)种群初始化：由输入序列产生初始构象种群 其中，N表示维数，表示第i个构象Cⁱ的第N维元素，并初始化迭代次数G＝0；

4)对初始种群中的每一个构象Cⁱ构建Lipschitz下界估计支撑面lⁱ：

$l^i=(\frac{f\left(C^i\right)}{M}-x_1^i,\frac{f\left(C^i\right)}{M}-x_2^i,\mathrm{...},\frac{f\left(C^i\right)}{M}-x_{N+1}^i)---\left(4\right)$

其中，f(Cⁱ)维构象Cⁱ的能量函数值，为辅助变量，M为常数；

5)针对种群中的每个构象个体Cⁱ，i∈{1,2,3,…,NP}，令C_target＝Cⁱ，C_target表示目标构象个体，执行以下操作得到变异构象C_mutant：

5.1)随机生成正整数rand1,rand2,rand3∈{1,2,3,......NP}，且rand1≠rand2≠rand3≠i；再生成4个随机整数randrange1,randrange2,randrange3,randrange4；其中randrange1≠randrange2，randrange3≠randrange4∈{1,2,…,L},L为序列长度；

5.2)令a＝min(randrange1,randrange2)，b＝max(randrange1,randrange2)，k∈[a,b]；令c＝min(randrange3,randrange4)，d＝max(randrange3,randrange4)，p∈[c,d]；其中min表示取两个数的最小值，max表示取两个数的最大值；

5.3)用C_rand2上位置a到位置b的片段的氨基酸k所对应的二面角phi、psi、omega替换C_rand1的相同位置所对应的二面角phi、psi、omega；再使用C_rand3上位置c到位置d的片段的氨基酸p所对应的二面角phi、psi、omega替换C_rand1上相同位置所对应的二面角phi、psi、omega，然后将所得C_rand1进行片段组装得到变异构象个体C_mutant；

6)对变异构象C_mutant执行交叉操作生成测试构象C_trial：

6.1)生成随机数rand4，rand5，其中rand4∈(0,1)，rand5∈(1,L)；

6.2)根据式(2)执行交叉过程：

$C_{trail}=\left\{\begin{array}{cc}{C}_{mu\tan t,rand5}\←C_{t\arg et,rand5},& if(rand4\≤CR)\\ C_{mu\tan t,rand5},& otherwise\end{array}\right.---\left(5\right)$

其中C_mutant,rand5表示变异构象中的rand5片段，若随机数rand4≤CR，变异构象C_mutant的片段rand5替换为目标构象C_target中对应的片段，否则C_trial直接等于变异构象C_mutant；

7)计算测试构象C_trial的下界估计值

$\stackrel{\‾}{f}\left(C_{trial}\right)=\underset{k\≤NP}{\max }\underset{i=1,\mathrm{...}N+1}{\min }M(l_j^k+x_{trial}^j)---\left(6\right)$

其中max表示求最大值，min表示求最小值，为测试构象向量的第j维元素，为支撑向量l^k的第j维元素；

8)根据如下操作决定测试构象C_trial是否替换目标构象C_target：

8.1)计算目标构象的能量函数值f(C_target)；

8.2)如果则目标构象C_target保持不变，并直接进入下一次迭代；

8.3)如果则计算测试构象C_trial的实际能量函数值f(C_trial)，如果f(C_trial)＜f(C_target)，则测试构象C_trial替换目标构象C_target并继续步骤8.4)；

8.4)根据式(1)对测试构象C_trial构建Lipschitz下界估计支撑面；

9)判断是否满足终止条件，若满足则输出结果并退出，否则返回步骤5)。

进一步，所述步骤9)中，对种群中的每一个个体都执行完步骤5)—8)以后，迭代次数G＝G+1，终止条件为迭代次数G达到预设最大迭代次数。

本实施例序列长度为46的α/β折叠蛋白质1FD4为实施例，一种基于群体Lipschitz下界估计的蛋白质结构预测方法，其中包含以下步骤：

1)给定输入序列信息；

2)设置参数：种群大小NP＝50，算法的最大迭代次数＝10000，交叉因子CR＝0.5，常数M＝1000；

3)种群初始化：由输入序列产生初始构象种群 其中，N表示维数，表示第i个构象Cⁱ的第N维元素，并初始化迭代次数G＝0；

4)对初始种群中的每一个构象Cⁱ构建Lipschitz下界估计支撑面lⁱ：

$l^i=(\frac{f\left(C^i\right)}{M}-x_1^i,\frac{f\left(C^i\right)}{M}-x_2^i,\mathrm{...},\frac{f\left(C^i\right)}{M}-x_{N+1}^i)---\left(7\right)$

其中，f(Cⁱ)维构象Cⁱ的能量函数值，为辅助变量，M为常数；

5)针对种群中的每个构象个体Cⁱ，i∈{1,2,3,…,NP}，令C_target＝Cⁱ，C_target表示目标构象个体，执行以下操作得到变异构象C_mutant：

5.1)随机生成正整数rand1,rand2,rand3∈{1,2,3,......NP}，且rand1≠rand2≠rand3≠i；再生成4个随机整数randrange1,randrange2,randrange3,randrange4；其中randrange1≠randrange2，randrange3≠randrange4∈{1,2,…,L},L为序列长度；

5.2)令a＝min(randrange1,randrange2)，b＝max(randrange1,randrange2)，k∈[a,b]；令c＝min(randrange3,randrange4)，d＝max(randrange3,randrange4)，p∈[c,d]；其中min表示取两个数的最小值，max表示取两个数的最大值；

5.3)用C_rand2上位置a到位置b的片段的氨基酸k所对应的二面角phi、psi、omega替换C_rand1的相同位置所对应的二面角phi、psi、omega；再使用C_rand3上位置c到位置d的片段的氨基酸p所对应的二面角phi、psi、omega替换C_rand1上相同位置所对应的二面角phi、psi、omega，然后将所得C_rand1进行片段组装得到变异构象个体C_mutant；

6)对变异构象C_mutant执行交叉操作生成测试构象C_trial：

6.1)生成随机数rand4，rand5，其中rand4∈(0,1)，rand5∈(1,L)；

6.2)根据式(2)执行交叉过程：

$C_{trail}=\left\{\begin{array}{cc}{C}_{mu\tan t,rand5}\←C_{t\arg et,rand5},& if(rand4\≤CR)\\ C_{mu\tan t,rand5},& otherwise\end{array}\right.---\left(8\right)$

其中C_mutant,rand5表示变异构象中的rand5片段，若随机数rand4≤CR，变异构象C_mutant的片段rand5替换为目标构象C_target中对应的片段，否则C_trial直接等于变异构象C_mutant；

7)计算测试构象C_trial的下界估计值

$\stackrel{\‾}{f}\left(C_{trial}\right)=\underset{k\≤NP}{\max }\underset{i=1,\mathrm{...}N+1}{\min }M(l_j^k+x_{trial}^j)---\left(9\right)$

其中max表示求最大值，min表示求最小值，为测试构象向量的第j维元素，为支撑向量l^k的第j维元素；

8)根据如下操作决定测试构象C_trial是否替换目标构象C_target：

8.1)计算目标构象的能量函数值f(C_target)；

8.2)如果则目标构象C_target保持不变，并直接进入下一次迭代；

8.3)如果，则计算测试构象C_trial的实际能量函数值f(C_trial)，如果f(C_trial)＜f(C_target)，则测试构象C_trial替换目标构象C_target并继续步骤8.4)；

8.4)根据式(1)对测试构象C_trial构建Lipschitz下界估计支撑面；

9)对种群中的每一个个体都执行完步骤5)—8)以后，迭代次数G＝G+1，判断是否满足终止条件G＝10000，若满足则输出结果并退出，否则返回步骤5)。

以序列长度为46的α/β折叠蛋白质1FD4为实施例，运用以上方法得到了该蛋白质的近天然态构象，最小均方根偏差为平均均方根偏差为构象更新图如图1所示，构象分布图如图2所示。

以上阐述的是本发明给出的一个实施例表现出来的优良优化效果，显然本发明不仅适合上述实施例，而且可以应用到实际工程中的各个领域(如电力系统，路径规划等优化问题)，同时在不偏离本发明基本精神及不超出本发明实质内容所涉及内容的前提下可对其做种种变化加以实施。

Claims (2)

1.一种基于群体Lipschitz下界估计的蛋白质结构预测方法，其特征在于：所述蛋白质结构预测方法包括以下步骤：

1)给定输入序列信息；

2)设置参数：种群大小NP，交叉因子CR，常数M，最大迭代次数；

3)种群初始化：由输入序列产生初始构象种群 其中，N表示维数，表示第i个构象Cⁱ的第N维元素，并初始化迭代次数G＝0；

4)对初始种群中的每一个构象Cⁱ构建Lipschitz下界估计支撑面lⁱ：

$l^i=(\frac{f\left(C^i\right)}{M}-x_1^i,\frac{f\left(C^i\right)}{M}-x_2^i,...,\frac{f\left(C^i\right)}{M}-x_{N+1}^i)---\left(1\right)$

其中，f(Cⁱ)维构象Cⁱ的能量函数值，为辅助变量，M为常数；

5)针对种群中的每个构象个体Cⁱ，i∈{1,2,3,…,NP}，令C_target＝Cⁱ，C_target表示目标构象个体，执行以下操作得到变异构象C_mutant：

5.1)随机生成正整数rand1,rand2,rand3∈{1,2,3,......NP}，且rand1≠rand2≠rand3≠i；再生成4个随机整数randrange1,randrange2,randrange3,randrange4；其中randrange1≠randrange2，randrange3≠randrange4∈{1,2,…,L},L为序列长度；

5.2)令a＝min(randrange1,randrange2)，b＝max(randrange1,randrange2)，k∈[a,b]；令c＝min(randrange3,randrange4)，d＝max(randrange3,randrange4)，p∈[c,d]；其中min表示取两个数的最小值，max表示取两个数的最大值；

5.3)用C_rand2上位置a到位置b的片段的氨基酸k所对应的二面角phi、psi、omega替换C_rand1的相同位置所对应的二面角phi、psi、omega；再使用C_rand3上位置c到位置d的片段的氨基酸p所对应的二面角phi、psi、omega替换C_rand1上相同位置所对应的二面角phi、psi、omega，然后将所得C_rand1进行片段组装得到变异构象个体C_mutant；

6)对变异构象C_mutant执行交叉操作生成测试构象C_trial：

6.1)生成随机数rand4，rand5，其中rand4∈(0,1)，rand5∈(1,L)；

6.2)根据式(2)执行交叉过程：

$C_{trail}=\left\{\begin{array}{cc}{C}_{mu\tan t,rand5}\←C_{t\arg et,rand5},& if(rand4\≤CR)\\ C_{mu\tan t,rand5},& otherwise\end{array}\right.---\left(2\right)$

其中C_mutant,rand5表示变异构象中的rand5片段，若随机数rand4≤CR，变异构象C_mutant的片段rand5替换为目标构象C_target中对应的片段，否则C_trial直接等于变异构象C_mutant；

7)计算测试构象C_trial的下界估计值

$\stackrel{\‾}{f}\left(C_{trial}\right)=\underset{k\≤NP}{max}\underset{i=1,...N+1}{min}M(l_j^k+x_{trial}^j)---\left(3\right)$

其中max表示求最大值，min表示求最小值，为测试构象向量的第j维元素， 为支撑向量l^k的第j维元素；

8)根据如下操作决定测试构象C_trial是否替换目标构象C_target：

8.1)计算目标构象的能量函数值f(C_target)；

8.2)如果则目标构象C_target保持不变，并直接进入下一次迭代；

8.3)如果则计算测试构象C_trial的实际能量函数值f(C_trial)，如果f(C_trial)＜f(C_target)，则测试构象C_trial替换目标构象C_target并继续步骤8.4)；

8.4)根据式(1)对测试构象C_trial构建Lipschitz下界估计支撑面；

9)判断是否满足终止条件，若满足则输出结果并退出，否则返回步骤5)。

2.如权利要求1所述的一种基于群体Lipschitz下界估计的蛋白质结构预测方法，其特征在于：所述步骤9)中，对种群中的每一个个体都执行完步骤5)—8)以后，迭代次数G＝G+1，终止条件为迭代次数G达到预设最大迭代次数。