CN103714265A - 一种基于蒙特卡洛局部抖动和片段组装的蛋白质三维结构预测方法 - Google Patents

Abstract

一种基于蒙特卡洛局部抖动和片段组装的蛋白质三维结构预测方法，包括以下步骤：首先针对蛋白质高维构象空间搜索空间复杂难题，在Rosetta力场模型下，根据蛋白质数据库构建片段库，利用蒙特卡洛统计方法，判断片段替换的有效性；在差分进化群体算法框架下，片段组装使得搜索空间的复杂度降低，同时，利用蒙特卡洛统计方法剔除错误的片段组装，通过进化算法的多样性，逐步减小构象搜索空间以提高搜索效率；同时，利用粗粒度的模型，忽略侧链，有效减小了搜索的代价。本发明有效得到局最优稳定构想、预测效率较高、收敛正确性较好。

Description

一种基于蒙特卡洛局部抖动和片段组装的蛋白质三维结构预测方法

技术领域

本发明涉及生物信息学中蛋白质三维结构预测技术领域，特别涉及一种蛋白质三维结构预测方法，属于将现代智能优化方法应用到蛋白质三维结构预测。

背景技术

生物信息学通过综合利用生物学，计算机科学和信息技术而揭示大量而复杂的生物数据所赋有的生物学奥秘。是当前研究的一个热点。生物信息学研究成果已经被广泛应用于序列比对，蛋白质比对，基因识别分析，分子进化，序列重叠群装配，遗传密码，药物设计，生物系统，蛋白质结构预测等。其中蛋白质结构预测是生物信息学领域的一个重要的分支。著名的诺贝尔化学奖获得者Anfinsen通过对核糖核酸酶A的经典研究表明去折叠的蛋白质在体外可以自发的进行再折叠，仅仅是序列本身已经包括了蛋白质正确折叠的所有信息，并提出蛋白质折叠的热力学假说，这一理论认为结构同源的蛋白质可以通过不同的折叠途径形成相似的天然构象，蛋白质的一级结构决定了其三维结构，即蛋白质的氨基酸序列决定了蛋白质的三维空间结构，从热力学分析，蛋白质的天然构型对应着其物理能量最小的状态。

随着2000年基因组工作草图的完成，基因组的工作重心已经从测序转向了基因功能的识别，这一问题的解决显得尤为重要。人们通过研究蛋白质的内在结构来了解其内在机理，并找出蛋白质折叠的真正驱动力，这样可以按照人们的设想设计出符合特定需求的非天然蛋白质，对疾病进行有效的预测和控制。蛋白质结构从头预测也因此具有重大的意义。然后蛋白质的结构测定远远赶不上基因组测序的速度，单纯的依靠实验的方法测定所有的蛋白质的三维结构已经远远不能满足后基因时代对于蛋白质数据分析的需求，有必要发展一种可靠的理论预测方法，从而迅速、简便的获取蛋白质的结构信息。

蛋白质的从头预测方法就是选择一种力场模型，将蛋白质结构抽象成一个数学模型，使其具有可计算性，然后进行算法的设计，从算法的时间和空间复杂度出发，利用全局以及局部的优化算法在能量曲面上搜索全局最优解，在数学上属于一种极其复杂的非凸多极值优化问题。利用蛋白质的一级序列，利用热力学方法，计算得出蛋白质的三维结构。

针对力场模型的优化问题，2005年Bradley在结合蒙特卡洛优化方法开发了分子建模Rosetta程序，成功预测CASP6中T0281目标蛋白(包含70个残基)的C_α-RMSD为对某些小目标蛋白(<85个残基)，甚至得到误差精度小于

的高分辨率三维解析结构；2006年，Zhan采用Basin-Paving(BP)方法得到了脑啡肽在ECEPP/2、ECEPP/3两种力场模型下最小能量构象；2009年Roosi提出几种BP改进策略，并给出了比较结果。结合随机全局优化算法，如遗传算法(GeneticAlgorithms,GA)、差分进化(Differential Evolution,DE)算法也应用于力场模型的优化，但是由于蛋白质是一种强柔性的大分子体系，其势函数表达式极其复杂，存在极多局部极小点，这些算法极容易陷入某个局优解；同时由于随机算法缺乏全局收敛理论依据及解的不确定性，使用了全局选择因子，算法只能收敛到全局最优解，而忽略了众多局部极值解，进一步限制了它们在实际问题中的应用。

因此，现有的技术在确定蛋白质稳定构象方面存在着缺陷，需要进行一系列的改进。

发明内容

为了克服现有技术中采用单一的现代进化算法进行蛋白质三维结构预测时，会面对容易陷入极小值而不易得到全局最优稳定构想、预测效率低、错误收敛的不足。本发明结合蒙特卡洛统计模拟方法与现有的已知知识库，提出了一种基于蒙特卡洛局部抖动和片段组装的蛋白质三维结构预测方法，有效得到局最优稳定构想，预测效率较高，收敛正确性较好，弥补了现有搜索空间过于复杂以及以上的缺陷。

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

一种基于蒙特卡洛局部抖动和片段组装的蛋白质三维结构预测方法，包括以下步骤：

1）选取力场模型

采用Rosetta力场模型能量函数的表示形式如下：

E_protein＝W_inter repE_inter rep+W_inter atrE_inter atr+W_solvationE_solvation+W_bb/sc hbE_bb/sc hb+W_bb/bb hbE_bb/bb hb+W_sc/sc hbE_sc/sc hb+W_pairE_pair+W_dunbrackE_dunbrack+W_ramaE_rama+W_referenceE_reference

式中，E_protein表示蛋白质的总能量，E_inter rep表示范德华排斥力作用，W_inter rep为E_inter rep在整体的权重，E_inter atr表示范德华吸引力作用，W_inter atr为E_inter atr在整体的权重，E_solvation为Lazarids和Karplus描述的隐含的溶解作用，W_solvation为E_solvation在整体的权重，E_bb/sc hb、E_bb/bb hb、E_sc/sc hb为依赖方向的氢键能量，W_bb/sc hb、W_bb/bb hb、W_sc/sc hb分别为其能量在整体的权重，E_pair为残基对静电作用，W_pair为E_pair在整体的权重，E_dunbrack为氨基酸基于旋转异构体库的内部的能量，W_dunbrack为E_dunbrack在整体的权重，E_rama为参考特定位置的Ramachandrin骨链扭力，W_rama为E_rama在整体的权重，E_reference为未折叠态的蛋白质的参考能量，W_reference为E_reference在整体的权重，Rosetta的能量函数就是将所有的能量项通过各自的权重线性相加；

2）设置权重比，采用Rosetta的score3权重比，采用一种粗粒度的能量函数；

3）片段库的构建，通过一条蛋白质序列，从已知的知识库，实验测得的蛋白质数据库获取各个片段的信息，采用n(n为大于2的整数)个氨基酸的片段，即将已有的一系列数据库中的蛋白质，分解成n(n为大于2的整数)个氨基酸的片段，通过判断各个原子之间的距离，分析四个原子构成的平面的二面角，以及均方根偏差等来分析各个片段的相似性，将每个位置相似性最大的m(m为大于1的整数)个片段信息记录下来，在使用过程中，通过索引使用此片段库；

4）初始化，设置群体规模popSize，n为3，m为200，进化代数G，变异因子F为0.5，随机产生popSize个蛋白质结构作为初始种群，每个个体表示为：X_i(i=1,2,…,popSize)，其中i为个体在种群中的序列，在进化过程中，种群规模保持不变；

5）对每一个目标个体X_i(i=1,2,…,popSize)作如下处理：

5.1)任意选三个个体{X_r1,X_r2,X_r3|r1,r2,r3∈{1,2,…,popSize},r1≠r2≠r3≠i}；

5.2)对选择出的种群的三个个体{X_r1,X_r2,X_r3}执行变异操作V_i＝X_r1+F*(X_r2-X_r3)，生成变异个体V_i；

5.3)设置蒙特卡洛方法的循环次数ncycle=50；

5.4)将变异个体V_i作为蒙特卡洛方法的初始构象，并且计算能量E_vi；

5.5)在变异个体V_i上分别随机选择三个不同的片段，分别从片段库中随机挑选相似的片段进行替换，即三种二面角的替换，交叉，生成测试个体

5.6)计算测试个体的能量

变异个体V_i与测试个体

的能量差

如果ΔU>0，则接受这个构象，将测试个体

作为变异个体V_i，i=i+1；如果i<ncycle，转至步骤5.4）

5.7)j=j+1；如果j<popSize，转至步骤5.1）

6）对目标个体X_i和测试个体

逐个更新操作：

6.1)计算目标个体X_i的能量和测试个体

的能量和均方根偏差，如果目标个体

的能量和均方根偏差比测试个体的大，则将目标个体X_i和测试个体

进行替换，更新目标个体；

6.2)n=n+1；如果n<popSize，转至6.1）

7）k=k+1；如果k<G，转至5）

8）k=G时的目标个体为结构接近实验测得的蛋白质结构。

本发明的有益效果为：在差分进化算法的框架下，降低了算法的搜索空间复杂度，同时利用蒙特卡洛统计方法的判别，提高了算法的收敛能力，弥补了容易陷入局部极小值，搜索空间过于复杂等缺陷。

附图说明

图1为本发明中算法实例采用的蛋白质的实验测定结构图；

图2为本发明中算法所得的1GYZ最优PDB图；

图3为本发明算法流程图；

图4为未使用本发明算法得到的种群个体评价图；

图5为使用本发明得到的种群个体评价图。

具体实施方式

下面结合附图对本发明进行详细的描述。

参照图1～图5，一种基于蒙特卡洛局部抖动和片段组装的蛋白质三维结构预测方法，包括以下步骤：

1）选取力场模型

力场模型是依赖于原子三维坐标的经验势能函数，由于其忽略了电子的相互作用，使得分子力场模型结构相对简单，计算速度较快，辅以成熟的力场参数，可达到很高的精度，能够用于生物大分子的结构预测。力场模型总势能通常被经验性的划分成若干项，本发明采用Rosetta力场模型能量函数的表示形式如下：

E_protein＝W_inter repE_inter rep+W_inter atrE_inter atr+W_solvationE_solvation+W_bb/sc hbE_bb/sc hb+W_bb/bb hbE_bb/bb hb+W_sc/sc hbE_sc/sc hb+W_pairE_pair+W_dunbrackE_dunbrack+W_ramaE_rama+W_referenceE_reference

式中E_protein表示蛋白质的总能量，E_inter rep表示范德华排斥力作用，W_inter rep为E_inter rep在整体的权重，E_inter atr表示范德华吸引力作用，W_inter atr为E_inter atr在整体的权重，E_solvation为Lazarids和Karplus描述的隐含的溶解作用，W_solvation为E_solvation在整体的权重，E_bb/sc hb、E_bb/bb hb、E_sc/sc hb为依赖方向的氢键能量，W_bb/sc hb、W_bb/bb hb、W_sc/sc hb分别为其能量在整体的权重，E_pair为残基对静电作用，W_pair为E_pair在整体的权重，E_dunbrack为氨基酸基于旋转异构体库的内部的能量，W_dunbrack为E_dunbrack在整体的权重，E_rama为参考特定位置的Ramachandrin骨链扭力，W_rama为E_rama在整体的权重，E_reference为未折叠态的蛋白质的参考能量，W_reference为E_reference在整体的权重，Rosetta的能量函数就是将所有的能量项通过各自的权重线性相加。

2）设置权重比，本发明采用的是Rosetta的score3权重比，采用一种粗粒度的能量函数。Score3权重比的能量函数是一种粗粒度的能量函数，它忽略了侧链，将蛋白质质心化转化为骨干模型，使计算优化更加有效率。

3）片段库的构建，通过一条蛋白质序列，可以从已知的知识库，实验测得的蛋白质数据库等获取各个片段的信息，本发明采用n(n为大于2的整数)个氨基酸的片段，即将已有的一系列数据库中的蛋白质，分解成n(n为大于2的整数)个氨基酸的片段，通过判断各个原子之间的距离，分析四个原子构成的平面的二面角，以及均方根偏差等来分析各个片段的相似性，将每个位置相似性最大的m(m为大于1的整数)个片段信息记录下来，在使用过程中，可以通过索引使用此片段库。

4）初始化，设置群体规模popSize=200，n为3，m为200，进化代数G，变异因子F为0.5，随机产生popSize个蛋白质结构作为初始种群，终止条件是算法循环300次，算法独立运行50，每个个体表示为：X_i(i=1,2,…,popSize)，其中i为个体在种群中的序列，在进化过程中，种群规模保持不变。

5）对每一个目标个体X_i(i=1,2,…,popSize)作如下处理：

5.1)任意选三个个体{X_r1,X_r2,X_r3|r1,r2,r3∈{1,2,…,popSize},r1≠r2≠r3≠i}；

5.2)对选择出的种群的三个个体{X_r1,X_r2,X_r3}执行变异操作V_i＝X_r1+F*(X_r2-X_r3)，生成变异个体V_i；

5.3)设置蒙特卡洛方法的循环次数ncycle=50；

5.4)将变异个体V_i作为蒙特卡洛方法的初始构象，并且计算能量E_vi；

5.5)在变异个体V_i上分别随机选择三个不同的片段，分别从片段库中随机挑选相似的片段进行替换，即三种二面角的替换，交叉，生成测试个体

5.6)计算测试个体的能量

变异个体V_i与测试个体

的能量差

如果ΔU>0，则接受这个构象，将测试个体

作为变异个体V_i，i=i+1；如果i<ncycle，转至步骤5.4）

5.7)j=j+1；如果j<popSize，转至步骤5.1）

6）对目标个体X_i和测试个体逐个更新操作：

6.1)计算目标个体X_i的能量和测试个体

的能量和均方根偏差，如果目标个体

的能量和均方根偏差比测试个体

的大，则将目标个体X_i和测试个体

进行替换，更新目标个体；

6.2)n=n+1；如果n<popSize，转至6.1）

7）k=k+1；如果k<G，转至5）

8）k=G时的目标个体为结构接近实验测得的蛋白质结构，将结构接近实验测得的蛋白质结构的目标个体输出。

实例：以编号为1GYZ的蛋白质为例，它是由60个氨基酸组成的蛋白质分子，氨基酸序列为:

WIARINAAVRAYGLNYSTFINGLKKAGIELDRKILADMAVRDPQAFEQV

VNKVKEALQVQ。

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

Claims (1)

1.一种基于蒙特卡洛局部抖动和片段组装的蛋白质三维结构预测方法，其特征在于：包括以下步骤：

1）选取力场模型

采用Rosetta力场模型能量函数的表示形式如下：

E_protein＝W_inter repE_inter rep+W_inter atrE_inter atr+W_solvationE_solvation+W_bb/sc hbE_bb/sc hb+W_bb/bb hbE_bb/bb hb+W_sc/sc hbE_sc/sc hb+W_pairE_pair+W_dunbrackE_dunbrack+W_ramaE_rama+W_referenceE_reference

式中，E_protein表示蛋白质的总能量，E_inter rep表示范德华排斥力作用，W_inter rep为E_inter rep在整体的权重，E_inter atr表示范德华吸引力作用，W_inter atr为E_inter atr在整体的权重，E_solvation为Lazarids和Karplus描述的隐含的溶解作用，W_solvation为E_solvation在整体的权重，E_bb/sc hb、E_bb/bb hb、E_sc/sc hb为依赖方向的氢键能量，W_bb/sc hb、W_bb/bb hb、W_sc/sc hb分别为其能量在整体的权重，E_pair为残基对静电作用，W_pair为E_pair在整体的权重，E_dunbrack为氨基酸基于旋转异构体库的内部的能量，W_dunbrack为E_dunbrack在整体的权重，E_rama为参考特定位置的Ramachandrin骨链扭力，W_rama为E_rama在整体的权重，E_reference为未折叠态的蛋白质的参考能量，W_reference为E_reference在整体的权重，Rosetta的能量函数就是将所有的能量项通过各自的权重线性相加；

2）设置权重比，采用Rosetta的score3权重比，采用一种粗粒度的能量函数；

3）片段库的构建，通过一条蛋白质序列，从已知的知识库，实验测得的蛋白质数据库获取各个片段的信息，采用n个氨基酸的片段，n为大于2的整数，即将已有的一系列数据库中的蛋白质，分解成n个氨基酸的片段，通过判断各个原子之间的距离，分析四个原子构成的平面的二面角，以及均方根偏差等来分析各个片段的相似性，将每个位置相似性最大的m个片段信息记录下来，m为大于1的整数，在使用过程中，通过索引使用此片段库；

4）初始化，设置群体规模popSize，进化代数G，变异因子F为0.5，随机产生popSize个蛋白质结构作为初始种群，每个个体表示为：X_i(i=1,2,…,popSize)，其中i为个体在种群中的序列，在进化过程中，种群规模保持不变；

5）对每一个目标个体X_i(i=1,2,…,popSize)作如下处理：

5.1)任意选三个个体{X_r1,X_r2,X_r3|r1,r2,r3∈{1,2,…,popSize},r1≠r2≠r3≠i}；

5.2)对选择出的种群的三个个体{X_r1,X_r2,X_r3}执行变异操作V_i＝X_r1+F*(X_r2-X_r3)，生成变异个体V_i；

5.3)设置蒙特卡洛方法的循环次数ncycle=50；

5.4)将变异个体V_i作为蒙特卡洛方法的初始构象，并且计算能量E_vi；

5.5)在变异个体V_i上分别随机选择三个不同的片段，分别从片段库中随机挑选相似的片段进行替换，即三种二面角的替换，交叉，生成测试个体

5.6)计算测试个体的能量

变异个体V_i与测试个体

的能量差

如果ΔU>0，则接受这个构象，将测试个体

作为变异个体V_i，i=i+1；如果i<ncycle，转至步骤5.4）

5.7)j=j+1；如果j<popSize，转至步骤5.1）

6）对目标个体X_i和测试个体

逐个更新操作：

6.1)计算目标个体X_i的能量和测试个体

的能量和均方根偏差，如果目标个体

的能量和均方根偏差比测试个体

的大，则将目标个体X_i和测试个体

进行替换，更新目标个体；

6.2)n=n+1；如果n<popSize，转至6.1）

7）k=k+1；如果k<G，转至5）

8）k=G时的目标个体为结构接近实验测得的蛋白质结构。

Images