基于差分进化和构象空间退火的蛋白质三维结构预测方法

摘要

针对蛋白质三维结构预测这个高维复杂的优化问题，本发明提出了一种综合差分进化算法和构象空间退火算法的一种蛋白质三维结构预测方法，将差分进化算法融入到构象空间退火算法中，形成混合优化算法。在构象空间退火算法的框架下，差分进化算法能够很好的提高种群的多样性，同时差分进化算法中的全局优化因子能够增加种群竞争，提高了算法的收敛能力。仿真结果表明，该方法具有较好的收敛性能和预测精度，能够有效的找到脑啡肽的全局能量最小对应的结构。

说明书

技术领域

本发明涉及蛋白质三维结构预测技术领域，特别涉及一种基于差分进化和构象空间退火相结合的一种蛋白质三维结构预测方法,属于将现代智能优化方法应用到蛋白质三维结构的工程预测方法。 

背景技术

随着人类基因组计划的完成和蛋白质的研究日益深入，在当今分子生物学领域中，蛋白质分子空间结构与功能的研究无疑是最具有挑战性的问题。三联密码的破译使人们掌握了遗传信息从DNA到氨基酸序列的复制规律，然而仅有氨基酸序列是不够的，氨基酸序列还必须形成一定的空间结构，才真正完成了蛋白质的合成，也才能行使其特定的生物功能，因此蛋白质结构的确定将是研究其生物功能的基础。同时蛋白质结构预测对于全新蛋白质分子设计、药物分子设计、生物和化学催化剂、光能-化学能-动能的转换、生物传感器、以及环境科学等许多领域都有深远的意义。特别是随着2000年基因组工作草图的完成，基因组的工作重心已经从测序转向了基因功能的识别，这一问题的解决显得尤为重要。人们通过研究蛋白质的内在结构来了解其内在机理，并找出蛋白质折叠的真正驱动力，这样可以按照人们的设想设计出符合特定需求的非天然蛋白质，对疾病进行有效的预测和控制。 

然而蛋白质结构的测定远远赶不上基因组测序速度，X射线晶体学方法和多维核磁共振技术是目前测定蛋白质结构的主要物理方法。用X射线晶体学方法测定蛋白质结构不仅需要花相当长的时间，在技术上也受到相当大的限制；多维核磁共振技术能够有效测定蛋白质的结构，但仅适用于小蛋白。因此，单纯依靠试验方法测定所有蛋白质的空间结构远远不能满足后基因时代对于蛋白质数据分析的需求，有必要发展一种可靠的理论预测方法，从而迅速、简便的获取蛋白质的结构信息。 

早在1961年，Anfinsen提出蛋白质一级结构完全决定其空间结构的著名论断，并因此获得诺贝尔奖，从此揭开了理论预测蛋白质结构的序幕。解决蛋白质空间结构预测问题可以分为三部曲：即首先进行形式化，将蛋白质结构抽象转化成一个数学模型；然后分析模型的可计算性；最后进行算法设计，从算法的时间和空间复杂度出发，寻找最优算法。理论计算方法（也称热力学方法）是一种常用的蛋白质结构预测方法，由于它仅利用一级序列信息进行预测，而不需要任何其他已知蛋白质结构信息，所以该方法是比较理想的预测方法。其基本假设是：一定环境中的天然蛋白质的三维结构是整个系统自由能最小的结构。要实现这一方法有两个关键：一是要有一个合理的势函数，势函数的全局极小值对应蛋白质的天然结构； 二是要有个好的算法，保证在有效的计算时间找到势能函数的全局最小。 

目前，根据Anfinsen假设，直接从氨基酸序列出发，基于分子力场势能模型，采用全局优化方法，在其势能面上搜索蛋白质分子系统的最小能量状态，从而能够高通量、廉价地预测肽链的天然构像，已经成为生物信息学最重要的研究课题之一。由于蛋白质是一种强柔性的大分子体系，其势函数表达式极其复杂，存在极多局部极小点。据估计，N个残基的氨基酸序列，其局部极小点在10^N以上，而迄今发现的蛋白质中肽链平均含300-500个氨基酸残基，所以有效解决全局优化问题是蛋白质结构预测中的一个关键。 

近年来，许多随机全局优化算法陆续提出来解决蛋白质三维结构预测问题，如遗传算法、差分进化算法、构象空间退火算法。差分进化算法虽然具有较强的全局搜索能力和求解问题的能力，还具有简单、通用和并行处理等特点。但是用差分进化这种群体优化算法处理多模态优化问题时，由于使用了全局选择因子，算法只能收敛到全局最优解，而忽略了众多局部极值解；其次，模型的复杂性造成这些算法极易陷入某个局优解；同时差分这种随机算法缺乏全局收敛理论依据及解的不确定性，进一步限制了它们在实际问题中的应用。 

构象空间退火算法结合了遗传算法、模拟退火等算法的优点，具有很好的搜索全局最优的能力。在构象空间退火法中，首先以一组随机选择的蛋白质构象开始，然后对这些构想进行能量极小化，所得到的局部最小蛋白质构象作为种群库中的新成员，这些构象要求是随机分布的，并且在相互独立的条件下进一步被优化，以覆盖种群库中尽可能大的区域。单种智能算法总存在一些算法固有的缺点，而将两种或多种智能算法按照某种规则融合在一起或在单种智能算法中引入其他优化思想，形成混合优化算法，则可以有效地扬长避短，发挥智能算法的优点，大大提高算法的全局和局部收敛能力。构象空间退火法虽然能够很好的解决蛋白质三维结构预测问题，但是由于蛋白质结构极其复杂，其势能函数曲面分布大量的局部极值点，而构象空间退火算法优化种群较小，优化过程中种群多样性较小，同时种群缺乏竞争机制，所以算法容易陷入局部极值解。 

发明内容

本发明针对现有的蛋白质三维结构预测在采用单一现代进化算法时，容易陷入极小值、预测效率差等缺点，提出了一种综合构象空间退火算法和差分进化算法的一种蛋白质三维结构预测方法。 

实现本发明而采取的技术方案概括如下： 

1、选择优化目标，即蛋白质力场模型函数E(X)。 

2、初始化进化种群：随机产生N_bank个蛋白质结构作为初始种群POP，同时设定算法各参数，并使用拟牛顿法对初始种群进行能量极小化处理。 

3、将处理过的N_bank个蛋白质个体作为初始种群库firstbank，并设定初始D_cut=1/2D_ave, 其中$D_{\mathrm{ave}}=\frac{1}{N_{\mathrm{bank}}*(N_{\mathrm{bank}}-1)}\underset{i=1}{\overset{N_{\mathrm{bank}}}{\Sigma }}\underset{j=1}{\overset{N_{\mathrm{bank}}}{\Sigma }}\underset{k=1}{\overset{N_{\mathrm{dire}}}{\Sigma }}|{\theta }_k^i-{\theta }_k^j|,$其中i≠j,N_dre表示的是优化目标的维数，表示的是第i个蛋白质结构的第k个二面角。 

4、将firstbank中的个体直接复制给Bank中，并更新Bank中的个体。 

5、随机从Bank中选取M个个体作为种子，每个种子按照以下四种扰动策略产生10个新的个体，总共产生M*10个新个体。同时检查Bank中的个体是否全部都作为过种子，若库中的个体全部做过种子，而循环仍然没有结束，则重新随机产生N_bank个蛋白质结构加入Bank中。 

5.1以基本差分进化算法中变异操作产生3个新的蛋白质结构，其中具体操作公式为：u＝x_a+F·(x_b-x_c)a≠b≠c，其中公式中基准矢量x_a为种子个体；个体x_b和x_c是随机从Bank中选取的蛋白质结构，F为变异因子。 

5.2以改进差分进化算法中变异操作产生3个新的蛋白质结构，其中具体操作公式为：u＝x_best+F·(x_b-x_c)best≠b≠c，其中x_best是Bank库中能量最低个体；x_b和x_c是随机从Bank中选取的蛋白质结构，F为变异因子。 

5.3随机选取一种蛋白质结构片段（小组），同时随机从M个种子中和Bank中各选取一个种子和待测试个体，将待测试个体中和小组相对应的局部片段直接复制给种子相应的片段。以这种扰动策略产生2个新的蛋白质结构。 

5.4随机选取一种蛋白质结构片段（组合组），同时随机从M个种子中和Bank中各选取一个种子和待测试个体，将待测试个体中和组合组相对应的局部片段直接复制给种子相应的片段。以这种扰动策略产生2个新的蛋白质结构。 

6、将新产生的个体经拟牛顿法进行能量极小化处理，然后将这些新产生的蛋白质结构更新给Bank各种子生境中。其中更新的具体操作： 

6.1、通过计算新产生个体B和Bank中的种子A之间的距离d（A,B），再将d（A,B）和D_cut进行比较。其中表示的是蛋白质结构A的第s个二面角，N_dre表示的是优化目标的维数。 

6.2、若d（A,B）<D_cut，则将新产生个体更新至此种子的生境中，同时若个体B所代表的能量值小于若小于A，则B成为这个种子生境的代表。 

6.3、若新产生个体B不属于Bank中任何的种子生境，则重新成为一个种子生境，同时为了保证Bank中生境总数不变，则移除Bank中能量值最大的种子生境。 

7、检测此时种群中是否出现蛋白质全局能量最小结构。如若找到循环结束；若没找到循环继续。每循环一轮，通过减小D_cut大小更新Bank，使D_cut减小至D_ave。 

8、循环执行步骤4-7，直至达到终止条件或找到蛋白质的全局能量最小结构。 

针对蛋白质三维结构预测这个高维复杂的优化问题，本发明的优点是：本发明提出了一种基于构象空间退火和差分进化相结合的一种蛋白质三维结构预测方法,在构象空间退火算法的框架下，差分进化算法能够很好的提高种群的多样性，同时差分进化算法中的全局优化因子能够增加种群竞争，提高了算法的收敛能力。 

附图说明

图1为本发明中算法的算法流程图； 

图2为算法对蛋白质构象相对应小组的划分图； 

图3为本发明中算法所得到全局能量最低的脑啡肽PDB图； 

具体实施方式

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

脑啡肽（Try¹-Gly²-Gly³-Phe⁴-Met⁵）是由五个氨基酸组成的蛋白质分子，由75个原子组成，可用24个独立的主-侧链二面角描述，公认的脑啡肽稳定能量值为-11.707kcal/mol。参照图2，算法将脑啡肽为例，将其对应的24个二面角分为8个小组，小组中的ψ、ω代表脑啡肽主链中的二面角，χⁱ代表脑啡肽侧链中的二面角。在算法中，我们进一步的将8个小组分为7个集合组，如下表1所示。小组和集合组中的成员分别对应24个二面角中的某些片段，这些小组和组合组类似于原始集结过程中两种不同的片段。 

表1 

实现本发明的技术方案概括如下： 

1、以脑啡肽为优化目标，选取ECEPP/3能量模型E(x)为本发明的优化模型，搜索其能量全局最低所对应的脑啡肽构象。 

2、初始化进化种群：随机产生N_bank个蛋白质结构作为初始种群POP，同时设定算法各参数，并使用拟牛顿法对初始种群进行能量极小化处理。 

3、将处理过的N_bank个蛋白质个体作为初始种群库firstbank，并设定初始D_cut=1/2D_ave,其 中$D_{\mathrm{ave}}=\frac{1}{N_{\mathrm{bank}}*(N_{\mathrm{bank}}-1)}\underset{i=1}{\overset{N_{\mathrm{bank}}}{\Sigma }}\underset{j=1}{\overset{N_{\mathrm{bank}}}{\Sigma }}\underset{k=1}{\overset{N_{\mathrm{dire}}}{\Sigma }}|{\theta }_k^i-{\theta }_k^j|,$其中i≠j,N_dre表示的是优化目标的维数，表示的是第i个蛋白质结构的第k个二面角。 

4、将firstbank中的个体直接复制给Bank中，并更新Bank中的个体。 

5、随机从Bank中选取M个个体作为种子，每个种子按照以下四种扰动策略产生10个新的个体，总共产生M*10个新个体。同时检查Bank中的个体是否全部都作为过种子，若库中的个体全部做过种子，而循环仍没有结束，则重新随机产生N_bank个蛋白质结构加入Bank中。 

5.1以基本差分进化算法中变异操作产生3个新的蛋白质结构，其中具体操作公式为：u＝x_a+F·(x_b-x_c)a≠b≠c，其中公式中基准矢量x_a为种子个体；个体x_b和x_c是随机从Bank中选取的蛋白质结构，F为变异因子。 

5.2以改进差分进化算法中变异操作产生3个新的蛋白质结构，其中具体操作公式为：u＝x_best+F·(x_b-x_c)best≠b≠c，其中x_best是Bank库中能量最低个体；x_b和x_c是随机从Bank中选取的蛋白质结构，F为变异因子。 

5.3随机选取一种蛋白质结构片段（小组），同时随机从M个种子中和Bank中各选取一个种子和待测试个体，将待测试个体中和小组相对应的局部片段直接复制给种子相应的片段。以这种扰动策略产生2个新的蛋白质结构。 

5.4随机选取一种蛋白质结构片段（组合组），同时随机从M个种子中和Bank中各选取一个种子和待测试个体，将待测试个体中和组合组相对应的局部片段直接复制给种子相应的片段。以这种扰动方式策略2个新的蛋白质结构。 

6、将新产生的个体经拟牛顿法进行能量极小化处理，然后将这些新产生的蛋白质结构更新给Bank各种子生境中。其中更新的具体操作： 

6.1、通过计算新产生个体B和Bank中的种子A之间的距离d（A,B），再将d（A,B）和D_cut进行比较。其中表示的是蛋白质结构A的第s个二面角，N_dre表示的是优化目标的维数。 

6.2、若d（A,B）<D_cut，则将新产生个体更新至此种子的生境中，同时若个体B所代表的能量值小于若小于A，则B成为这个种子生境的代表。 

6.3、若新产生个体B不属于Bank中任何的种子生境，则重新成为一个种子生境，同时为了保证Bank中生境总数不变，则移除Bank中能量值最大的种子生境。 

7、检测此时种群中是否出现蛋白质全局能量最小结构。如若找到循环结束；若没找到循环继续。每循环一轮，使D_cut减小至D_ave，通过减小D_cut大小更新Bank。 

8、循环执行4-7步骤，直至达到终止条件或找到蛋白质全局能量最小结构。 

本发明中算法的参数设置如下：种群pop大小N_bank=50，选取的种子数目M=10，变异因子F=0.9，终止条件是算法循环10次。算法独立运行50，图3显示的是本发明50次运行所得到的所有能量值小于-10kcal/mol的脑啡肽结构分布图。如图3所示，算法在50次运行中，能够找到192个能量值小于-10kcal/mol的脑啡肽结构。其中能量值在-11.0至-11.8kcal/mol的脑啡肽结构有18个。表2显示的是本发明所得到的全局能量最小对应的脑啡肽的二面角解。图3显示的全局能量最小对应的脑啡肽的PDB图。 

表2 

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