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

摘要

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

说明书

技术领域

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

背景技术

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

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

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

针对力场模型的优化问题，2005年Bradley在结合蒙特卡洛优化方法开发了 分子建模Rosetta程序，成功预测CASP6中T0281目标蛋白(包含70个残基)的 C_α-RMSD为对某些小目标蛋白(<85个残基)，甚至得到误差精度小于的高分辨率三维解析结构；2006年，Zhan采用Basin-Paving(BP)方法得到了脑啡 肽在ECEPP/2、ECEPP/3两种力场模型下最小能量构象；2009年Roosi提出几 种BP改进策略，并给出了比较结果。结合随机全局优化算法，如遗传算法(Genetic  Algorithms,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为实例所得出的优良效果，并非限定本发明的实 施范围，在不偏离本发明基本内容所涉及范围的的前提下对其做各种变形和改 进，不应排除在本发明的保护范围之外。