CN109300505B - 一种基于有偏采样的蛋白质结构预测方法 - Google Patents

Abstract

一种基于有偏采样的蛋白质结构预测方法，在遗传算法的基本框架下，对每个目标个体进行基于二级结构的有偏交叉策略；然后，随机选择一个变异窗口，以每个片段的二级结构得分为适应值，利用轮盘赌的方式选择一个适应值较高的片段来替换构象中的片段信息；最后，以构象的结构属性为主，能量属性为辅共同来指导种群的更新，进而选出潜在的构象进入下一代种群。本发明提出一种预测精度高、计算代价低的基于有偏采样的蛋白质结构预测方法。

Description

一种基于有偏采样的蛋白质结构预测方法

技术领域

本发明涉及一种生物学信息学、智能优化、计算机应用领域，尤其涉及的是一种基于有偏采样的蛋白质结构预测方法。

背景技术

生物领域高通量的基因数据以及蛋白质测序技术的飞速发展，使蛋白质结构测定技术面临前所未有的挑战。利用计算机的强大处理能力解决蛋白质结构预测问题是一项紧迫而现实的解决方案。而从头预测方法从氨基酸序列准确的预测蛋白质的三维结构是该领域中难度最高、意义最大、影响深远的问题之一。

在生物信息学中，蛋白质结构预测是一个组合优化问题。而由于能量函数的不准确性，导致优化目标难以准确量化，增加了该问题的求解复杂性。同时，蛋白质结构预测中解的构成不完整，解的长度不唯一，更加增大了求解的难度。目前常见的预测方法往往侧重能量函数的作用，而弱化搜索策略研究。而从头预测蛋白质结构最主要的问题就在于采样：天然结构和非天然结构相比具有最低自由能，但是我们却很难在巨大的构象空间中确定天然结构这一唯一的构象。

由于采样算法在从头预测中的重要性，所以许多研究者针对这一问题开发了多种针对蛋白质三维结构的采样算法，比如遗传算法、结构空间模拟退火、和分子动力学模拟采样。近几年的CASP比赛结果表明，Rosetta可以说是目前蛋白质结构从头预测领域最成功的方法，其核心的算法在于使用已知结构的蛋白质片段，通过蒙特卡罗的方法进行组装，产生接近天然结构的蛋白质构象。对于序列长度较短的，Rosetta已经能够很好的采样出非常接近天然结构的构象。但是对于更长的序列，更复杂的结构，却预测精度不高。

因此，现有的构象空间优化方法在预测精度和采样效率方面存在着缺陷，需要改进。

发明内容

为了克服现有的蛋白质构象优化方法中采样效率低、预测精度较低的不足，本发明提出一种采样效率较高、预测精度高的基于有偏采样的蛋白质结构预测方法。

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

一种基于有偏采样的蛋白质结构预测方法，所述方法包括以下步骤：

1)给定输入序列信息；

2)利用PSIPRED平台预测目标蛋白的二级结构信息，利用ROSETTA平台构建9片段的片段库；

3)选取Rosetta score3能量函数E(x)；

4)参数初始化：设置种群规模Ps，迭代计数器g，最大遗传代数G_max，初始种群搜索轨迹长度N，交叉片段长度c；

5)初始化种群：启动Ps条Monte Carlo轨迹，每条轨迹搜索N次，即生成Ps个初始个体；

6)对每个目标个体x_i,i∈{i＝1,.2..,Ps}进行如下操作：

6.1)随机选择一个个体x_j,j∈{1,2,...,Ps}且j≠i，进行交叉操作：

6.1.1)产生一个随机整数r∈[1,l-c]，确定出交叉区域[r,c+r]，其中l为序列长度；

6.1.2)根据预测的二级结构知识，确定出目标蛋白在区域[r,c+r]中对应的二级结构S_rc＝{q_k|r≤k≤r+c}，其中，q_k∈{H,E,L}是在目标蛋白中预测的第k氨基酸对应的二级结构类型；

6.1.3)利用DSSP测定个体x_i的二级结构，并确定出交叉区域对应的二级结构序列

其中，

是构象x_i中第k个残基对应的二级结构类型；

6.1.4)根据公式

计算个体x_i中交叉片段的得分，其中，f(k)根据如下公式计算：

其中，p_k代表在预测的二级结构知识中第k个氨基酸的二级结构被预测为q_k的置信度；

6.1.5)与步骤6.1.3)同理，确定出个体x_j中区域[r,c+r]对应的二级结构序列

6.1.6)与步骤6.1.4)同理，计算个体x_j中交叉片段的二级结构得分

6.1.7)比较

与

的大小，如果

则x′_i＝x_i,并进入步骤6.2)，否则，执行步骤6.1.8)；

6.1.8)用个体x_j中[r,c+r]区域的扭转角度依次替换个体x_i中对应区域的扭转角，生成交叉个体x′_i；

6.2)对个体x′_i进行如下的变异操作：

6.2.1)随机产生一个整数u∈[1,l-9]，进而确定出个体x_i′的片段插入窗口[u,9+u]以及片段库中对应的片段

h∈[1,200]，其中h是片段库中对每一个片段的编号；

6.2.2)根据片段库中的信息确定出在第u个位置处，每个片段

的二级结构

其中，

是片段库中第h个片段的第t个氨基酸对应的二级结构类型；

6.2.3)与步骤6.1.4)同理计算得到在第u个位置处，片段库中每个片段

的二级结构得分

6.2.4)根据公式

计算在第u个位置处每个片段的适应值；

6.2.5)以

为每个片段的适应值，利用轮盘赌的方式随机选取一个片段，并用该片段的扭转角替换个体x′_i中的扭转角，生成变异分体x″_i；

6.3)根据能量和二级结构得分进行如下的选择操作：

6.3.1)利用DSSP分别测定出个体x_i和x″_i的二级结构序列

和

6.3.2)与步骤6.1.4)同理，分别计算个体x_i和x″_i的二级结构得分ss_i和ss″_i；

6.3.3)比较ss_i和ss″_i的大小，如果ss_i＜ss″_i，则变异个体x″_i换个体x_i，并作为下一代的目标个体，执行步骤7)，否则执行步骤6.3.4)和步骤6.3.5)；

6.3.4)根据Rosetta score3函数分别计算个体x_i和x″_i的能量E(x_i)和E(x″_i)；

6.3.5)比较E(x_i)和E(x″_i)的大小，如果E(x″_i)＜E(x_i)，则变异个体x″_i换个体x_i，并作为下一代的目标个体，执行步骤7)，否则个体x_i被保留，并作为下一代的目标个体，执行步骤7)；

7)对种群中的每一个个体都执行完步骤6)以后，迭代次数g＝g+1，判断g是否大于G_max，若g＞G_max，则停止迭代并退出，否则返回步骤6)。

本发明的技术构思为：在遗传算法的基本框架下，对每个目标个体进行基于二级结构的有偏交叉策略；然后，随机选择一个变异窗口，以每个片段的二级结构得分为适应值，利用轮盘赌的方式选择一个适应值较高的片段来替换构象中的片段信息；最后，以构象的结构属性为主，能量属性为辅共同来指导种群的更新，进而选出潜在的构象进入下一代种群。

本发明的有益效果表现在：一方面通过预测的二级结构设计有偏的构象空间采样策略，显著的降低了构象的采样空间，进而提高的搜索效率；另一方面，通过结构属性和能量属性共同指导种群更新，降低了由于能量函数不精确带来的误差，进而大大提高了预测精度。

附图说明

图1是基于有偏采样的蛋白质结构预测方法的基本示意图。

图2是基于有偏采样的蛋白质结构预测方法对蛋白质1fna进行结构预测时的构象更新示意图。

图3是基于有偏采样的蛋白质结构预测方法对蛋白质1fna进行结构预测得到的三维结构图。

具体实施方式

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

参照图1～图3，一种基于有偏采样的蛋白质结构预测方法，所述方法包括以下步骤：

1)给定输入序列信息；

2)利用PSIPRED平台预测目标蛋白的二级结构信息，利用ROSETTA平台构建9片段的片段库；

3)选取Rosetta score3能量函数E(x)；

4)参数初始化：设置种群规模Ps，迭代计数器g，最大遗传代数G_max，初始种群搜索轨迹长度N，交叉片段长度c；

5)初始化种群：启动Ps条Monte Carlo轨迹，每条轨迹搜索N次，即生成Ps个初始个体；

6)对每个目标个体x_i,i∈{i＝1,.2..,Ps}进行如下操作：

6.1)随机选择一个个体x_j,j∈{1,2,...,Ps}且j≠i，进行交叉操作：

6.1.1)产生一个随机整数r∈[1,l-c]，确定出交叉区域[r,c+r]，其中l为序列长度；

6.1.2)根据预测的二级结构知识，确定出目标蛋白在区域[r,c+r]中对应的二级结构S_rc＝{q_k|r≤k≤r+c}，其中，q_k∈{H,E,L}是在目标蛋白中预测的第k氨基酸对应的二级结构类型；

6.1.3)利用DSSP测定个体x_i的二级结构，并确定出交叉区域对应的二级结构序列

其中，

是构象x_i中第k个残基对应的二级结构类型；

6.1.4)根据公式

计算个体x_i中交叉片段的得分，其中，f(k)根据如下公式计算：

其中，p_k代表在预测的二级结构知识中第k个氨基酸的二级结构被预测为q_k的置信度；

6.1.5)与步骤6.1.3)同理，确定出个体x_j中区域[r,c+r]对应的二级结构序列

6.1.6)与步骤6.1.4)同理，计算个体x_j中交叉片段的二级结构得分

6.1.7)比较

与

的大小，如果

则x′_i＝x_i,并进入步骤6.2)，否则，执行步骤6.1.8)；

6.1.8)用个体x_j中[r,c+r]区域的扭转角度依次替换个体x_i中对应区域的扭转角，生成交叉个体x′_i；

6.2)对个体x′_i进行如下的变异操作：

6.2.1)随机产生一个整数u∈[1,l-9]，进而确定出个体x′_i的片段插入窗口[u,9+u]以及片段库中对应的片段

h∈[1,200]，其中h是片段库中对每一个片段的编号；

6.2.2)根据片段库中的信息确定出在第u个位置处，每个片段

的二级结构

其中，

是片段库中第h个片段的第t个氨基酸对应的二级结构类型；

6.2.3)与步骤6.1.4)同理计算得到在第u个位置处，片段库中每个片段

的二级结构得分

6.2.4)根据公式

计算在第u个位置处每个片段的适应值；

6.2.5)以

为每个片段的适应值，利用轮盘赌的方式随机选取一个片段，并用该片段的扭转角替换个体x′_i中的扭转角，生成变异分体x″_i；

6.3)根据能量和二级结构得分进行如下的选择操作：

6.3.1)利用DSSP分别测定出个体x_i和x″_i的二级结构序列

和

6.3.2)与步骤6.1.4)同理，分别计算个体x_i和x″_i的二级结构得分ss_i和ss″_i；

6.3.3)比较ss_i和ss″_i的大小，如果ss_i＜ss″_i，则变异个体x″_i换个体x_i，并作为下一代的目标个体，执行步骤7)，否则执行步骤6.3.4)和步骤6.3.5)；

6.3.4)根据Rosetta score3函数分别计算个体x_i和x″_i的能量E(x_i)和E(x″_i)；

6.3.5)比较E(x_i)和E(x″_i)的大小，如果E(x″_i)＜E(x_i)，则变异个体x″_i换个体x_i，并作为下一代的目标个体，执行步骤7)，否则个体x_i被保留，并作为下一代的目标个体，执行步骤7)；

7)对种群中的每一个个体都执行完步骤6)以后，迭代次数g＝g+1，判断g是否大于G_max，若g＞G_max，则停止迭代并退出，否则返回步骤6)。

本实施例序列长度为91的β折叠蛋白质1fna为实施例，一种基于有偏采样的蛋白质结构预测方法，包括以下步骤：

1)给定输入序列信息；

2)利用PSIPRED平台预测目标蛋白的二级结构信息，利用ROSETTA平台构建9片段的片段库；

3)选取Rosetta score3能量函数E(x)；

4)参数初始化：设置种群规模Ps＝100，迭代计数器g＝0，最大遗传代数G_max＝200，初始种群搜索轨迹长度N＝2000，交叉片段长度c＝6；

5)初始化种群：启动Ps条Monte Carlo轨迹，每条轨迹搜索N次，即生成Ps个初始个体；

6)对每个目标个体x_i,i∈{i＝1,.2..,Ps}进行如下操作：

6.1)随机选择一个个体x_j,j∈{1,2,...,Ps}且j≠i，进行交叉操作：

6.1.1)产生一个随机整数r∈[1,l-c]，确定出交叉区域[r,c+r]，其中l为序列长度；

6.1.2)根据预测的二级结构知识，确定出目标蛋白在区域[r,c+r]中对应的二级结构S_rc＝{q_k|r≤k≤r+c}，其中，q_k∈{H,E,L}是在目标蛋白中预测的第k氨基酸对应的二级结构类型；

6.1.3)利用DSSP测定个体x_i的二级结构，并确定出交叉区域对应的二级结构序列

其中，

是构象x_i中第k个残基对应的二级结构类型；

6.1.4)根据公式

计算个体x_i中交叉片段的得分，其中，f(k)根据如下公式计算：

其中，p_k代表在预测的二级结构知识中第k个氨基酸的二级结构被预测为q_k的置信度；

6.1.5)与步骤6.1.3)同理，确定出个体x_j中区域[r,c+r]对应的二级结构序列

6.1.6)与步骤6.1.4)同理，计算个体x_j中交叉片段的二级结构得分

6.1.7)比较

与

的大小，如果

则x′_i＝x_i,并进入步骤6.2)，否则，执行步骤6.1.8)；

6.1.8)用个体x_j中[r,c+r]区域的扭转角度依次替换个体x_i中对应区域的扭转角，生成交叉个体x′_i；

6.2)对个体x′_i进行如下的变异操作：

6.2.1)随机产生一个整数u∈[1,l-9]，进而确定出个体x′_i的片段插入窗口[u,9+u]以及片段库中对应的片段

h∈[1,200]，其中h是片段库中对每一个片段的编号；

6.2.2)根据片段库中的信息确定出在第u个位置处，每个片段

的二级结构

其中，

是片段库中第h个片段的第t个氨基酸对应的二级结构类型；

6.2.3)与步骤6.1.4)同理计算得到在第u个位置处，片段库中每个片段

的二级结构得分

6.2.4)根据公式

计算在第u个位置处每个片段的适应值；

6.2.5)以

为每个片段的适应值，利用轮盘赌的方式随机选取一个片段，并用该片段的扭转角替换个体x′_i中的扭转角，生成变异分体x″_i；

6.3)根据能量和二级结构得分进行如下的选择操作：

6.3.1)利用DSSP分别测定出个体x_i和x″_i的二级结构序列

和

6.3.2)与步骤6.1.4)同理，分别计算个体x_i和x″_i的二级结构得分ss_i和ss″_i；

6.3.3)比较ss_i和ss″_i的大小，如果ss_i＜ss″_i，则变异个体x″_i换个体x_i，并作为下一代的目标个体，执行步骤7)，否则执行步骤6.3.4)和步骤6.3.5)；

6.3.4)根据Rosetta score3函数分别计算个体x_i和x″_i的能量E(x_i)和E(x″_i)；

6.3.5)比较E(x_i)和E(x″_i)的大小，如果E(x″_i)＜E(x_i)，则变异个体x″_i换个体x_i，并作为下一代的目标个体，执行步骤7)，否则个体x_i被保留，并作为下一代的目标个体，执行步骤7)；

7)对种群中的每一个个体都执行完步骤6)以后，迭代次数g＝g+1，判断g是否大于G_max，若g＞G_max，则停止迭代并退出，否则返回步骤6)。

以序列长度为91的β折叠蛋白质1fna为实施例，运用以上方法得到了该蛋白质的近天然态构象，最小均方根偏差为

平均均方根偏差为

预测结构如图3所示。

以上说明是本发明以1fna蛋白质为实例所得出的预测效果，并非限定本发明的实施范围，在不偏离本发明基本内容所涉及范围的前提下对其做各种变形和改进，不应排除在本发明的保护范围之外。

Claims (1)

1.一种基于有偏采样的蛋白质结构预测方法，其特征在于，所述方法包括以下步骤：

1)给定输入序列信息；

2)利用PSIPRED平台预测目标蛋白的二级结构信息，利用ROSETTA平台构建9片段的片段库；

3)选取Rosetta score3能量函数E(x)；

4)参数初始化：设置种群规模Ps，迭代计数器g，最大遗传代数G_max，初始种群搜索轨迹长度N，交叉片段长度c；

5)初始化种群：启动Ps条Monte Carlo轨迹，每条轨迹搜索N次，即生成Ps个初始个体；

6)对每个目标个体x_i,i∈{i＝1,2,. ..,Ps}进行如下操作：

6.1)随机选择一个个体x_j,j∈{1,2,...,Ps}且j≠i，进行交叉操作：

6.1.1)产生一个随机整数r∈[1,l-c]，确定出交叉区域[r,c+r]，其中l为序列长度；

6.1.2)根据预测的二级结构知识，确定出目标蛋白在区域[r,c+r]中对应的二级结构S_rc＝{q_k|r≤k≤r+c}，其中，q_k∈{H,E,L}是在目标蛋白中预测的第k氨基酸对应的二级结构类型；

6.1.3)利用DSSP测定个体x_i的二级结构，并确定出交叉区域对应的二级结构序列

其中，

是构象x_i中第k个残基对应的二级结构类型；

6.1.4)根据公式

计算个体x_i中交叉片段的得分，其中，f(k)根据如下公式计算：

其中，p_k代表在预测的二级结构知识中第k个氨基酸的二级结构被预测为q_k的置信度；

6.1.5)与步骤6.1.3)同理，确定出个体x_j中区域[r,c+r]对应的二级结构序列

6.1.6)与步骤6.1.4)同理，计算个体x_j中交叉片段的二级结构得分

6.1.7)比较

与

的大小，如果

则x_i′＝x_i,并进入步骤6.2)，否则，执行步骤6.1.8)；

6.1.8)用个体x_j中[r,c+r]区域的扭转角度依次替换个体x_i中对应区域的扭转角，生成交叉个体x_i′；

6.2)对个体x′_i进行如下的变异操作：

6.2.1)随机产生一个整数u∈[1,l-9]，进而确定出个体x′_i的片段插入窗口[u,9+u]以及片段库中对应的片段

其中h是片段库中对每一个片段的编号；

6.2.2)根据片段库中的信息确定出在第u个位置处，每个片段

的二级结构

其中，

是片段库中第h个片段的第t个氨基酸对应的二级结构类型；

6.2.3)与步骤6.1.4)同理计算得到在第u个位置处，片段库中每个片段

的二级结构得分

6.2.4)根据公式

计算在第u个位置处每个片段的适应值；

6.2.5)以

为每个片段的适应值，利用轮盘赌的方式随机选取一个片段，并用该片段的扭转角替换个体x′_i中的扭转角，生成变异分体x″_i；6.3)根据能量和二级结构得分进行如下的选择操作：

6.3.1)利用DSSP分别测定出个体x_i和x″_i的二级结构序列

和

6.3.2)与步骤6.1.4)同理，分别计算个体x_i和x″_i的二级结构得分ss_i和ss″_i；

6.3.3)比较ss_i和ss″_i的大小，如果ss_i＜ss″_i，则变异个体x″_i换个体x_i，并作为下一代的目标个体，执行步骤7)，否则执行步骤6.3.4)和步骤6.3.5)；

6.3.4)根据Rosetta score3函数分别计算个体x_i和x″_i的能量E(x_i)和E(x″_i)；

6.3.5)比较E(x_i)和E(x″_i)的大小，如果E(x″_i)＜E(x_i)，则变异个体x″_i换个体x_i，并作为下一代的目标个体，执行步骤7)，否则个体x_i被保留，并作为下一代的目标个体，执行步骤7)；

7)对种群中的每一个个体都执行完步骤6)以后，迭代次数g＝g+1，判断g是否大于G_max，若g＞G_max，则停止迭代并退出，否则返回步骤6)。

Images