一种基于Rosetta局部增强的群体蛋白质结构预测方法

摘要

一种基于Rosetta局部增强的群体蛋白质结构预测方法，首先，将结构预测中的整个算法搜索过程分为四个阶段，对每个阶段设置片段长度进行片段组装，并选用不同的能量函数来衡量构象个体的质量；然后，基于二级结构信息，采用不同的变异策略利用loop区域信息来生成测试构象，并通过随机交换loop区域信息实现交叉过程，保持种群多样性，同时对每个阶段的测试构象和目标构象执行Rosetta局部增强；最后，提取构象的特征向量来衡量各构象个体的多样性，从而以能量函数为主要衡量指标，并以多样性为辅助衡量指标来指导构象种群更新。本发明搜索能力强、能量保持种群多样性、预测精度高。

说明书

技术领域

本发明涉及一种生物学信息学、智能优化、计算机应用领域，尤其涉及的是，一种基于Rosetta局部增强的群体蛋白质结构预测方法。

背景技术

蛋白质是细胞功能的核心，与大部分核心生命过程息息相关。事实上，蛋白质只有折叠成特定的三维结构(即蛋白质三级结构)之后才能产生其特定的生物学功能。因此，为了解蛋白质的功能，就必须获得其三维空间结构，从而通过了解蛋白质的三维结构推动功能材料设计和新型药物研制的发展，帮助人们理解生命的基本过程，包括对阿尔兹海默症、帕金森病以及II型糖尿病等蛋白质折叠病的认识。

目前常用的蛋白质结构测定方法有X射线衍射和核磁共振(NMR)，这两种方法虽然预测精度高，但是对于X射线衍射来说，难以培养晶体且晶体结构测定的周期较长，核磁共振对样品的需要量大、纯度要求高，目前只能用于小分子蛋白质结构的测定。因此，以计算机为工具，利用适当的优化算法，直接通过氨基酸序列预测蛋白质三维结构，进而设计具有潜在药物价值的新功能蛋白质与多肽分子是生命科学领域需要解决的一个根本问题。该问题的最终解决关键在于：如何利用现有技术，设计一种高效的蛋白质构象空间优化算法。

经过40多年的发展，尤其是进入21世纪以来，分子动力学模拟(MD)、蒙特卡罗(MC)、构象空间退火(CSA)、进化类优化算法(EA)等随机优化算法在从头预测领域得到了成功应用；格点系统搜索(SGS)、分枝定界(BB)等确定性全局优化算法，理论研究超前于其数值应用，其极高的计算复杂度，限制了它们在中等规模以上蛋白构象优化方面的应用。基于MC及CSA系列改进算法，Baker团队开发的Rosetta从头预测服务器、Zhang团队开发的I-TASSER及QUARK从头预测服务器目前已经成为国际领先的预测软件。上述方法在预测一些序列长度较短的小蛋白来说，能够有效的得到三维结构。然而，由于蛋白质能量模型考虑分子体系成键作用以及范德华力、静电、氢键、疏水等非成键作用，致使其形成的能量曲面极其粗糙，构象对应局部极小解数目随序列长度的增加呈指数增长，对于这些传统方法进行预测显得力不从心，其原因在于极大的构象搜索空间会导致算法在预测过程中搜索能力渐渐下降，同时群体的多样性也变得越来越小，从而导致算法失去搜索的动力，影响最终的预测精度。

因此，现有的群体蛋白质结构预测方法在搜索能力和种群多样性保持方面存在着缺陷，需要改进。

发明内容

为了克服现有的群体蛋白质结构预测方法在搜索能力和种群多样性方面的不足，本发明提出一种搜索能力强，且能够保持种群多样性的基于Rosetta局部增强的群体蛋白质结构预测方法。

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

一种基于Rosetta局部增强的群体蛋白质结构预测方法，所述方法包括以下步骤：

1)输入待测蛋白质的氨基酸序列信息；

2)初始化：设置种群规模NP，交叉概率CR，策略选择因子CS，多样性接受概率RS，Rosetta轨迹长度T，片段长度L₁，L₂；

3)根据序列信息以片段长度L₁进行随机片段组装生成初始构象种群P＝{C¹,C²,...,C^NP}，其中，Cⁱ的表示当前种群中的第i个构象个体，并根据能量函数RosettaScore0计算各构象个体的能量，同时初始化迭代次数G＝0；

4)采用能量函数Rosetta Score0评价构象的质量，以片段长度L₁对初始种群中的每个构象个体执行轨迹长度为T的Rosetta局部增强，并计算每个构象的特征向量；

5)对步骤4)中增强后的每个构象个体Cⁱ,i∈{1,2,…,NP}作如下处理：

5.1)设置能量函数和片段长度：

5.1.1)如果当前迭代次数0＜G≤G_max/3，则片段长度l＝L₁，且选用能量函数Rosetta>

5.1.2)如果当前迭代次数G_max/3＜G≤2G_max/3，则片段长度l＝L₁，且选用能量函数Rosetta>

5.1.3)如果当前迭代次数G＞2G_max/3，则片段长度l＝L₂，且选用能量函数RosettaScore3

5.2)如果当前迭代次数G为G_max/3的整数倍，则对以片段长度l对构象个体Cⁱ执行轨迹长度为T的Rosetta局部增强，并根据步骤5.1)中设置的能量函数进行评价；

5.3)计算目标构象Cⁱ的特征向量，以及Cⁱ与当前种群中其他构象之间的特征向量欧氏距离，并以最小距离为Cⁱ的多样性值Dⁱ；

5.4)根据序列信息，利用DSSP得到待测蛋白的loop区域，并随机生成一个0到1之间的随机数p；

5.5)如果p<CS，则从当前种群中选取三个互不相同的构象个体C^a、C^b和C^c，其中，a≠b≠c≠i，a,b,c∈[1,NP]，从构象个体C^b和C^c中各随机选取一个片段替换C^a中对应位置的片段，并从C^a中随机选取一个不包含loop区

域的窗口进行片段组装生成变异构象C^mutant；

5.6)如果p≥1-CS，则选出当前能量值最低的构象个体C^best，并从当前种群中选取两个互不相同的构象个体C^a和C^b，其中，a≠b≠i，a,b∈[1,NP]，从构象个体C^a和C^b中各随机选取一个片段替换C^best中对应位置的片段，并从C^best中随机选取一个不包含loop区域的窗口进行片段组装生成变异构象C^mutant；

5.7)随机生成一个0与1之间随机数p′，如果p′＞CR，则随机选取一个loop区域，替换目标构象个体Cⁱ与变异构象个体C^mutant在该区域的二面角，从而生成测试构象C^trial，否则C^trial直接等于变异构象C^mutant；以片段长度l对测试构象个体C^trial执行轨迹长度为T的Rosetta局部增强；

5.8)计算增强后测试构象的特征向量，并计算测试构象的特征向量与当前种群中各构象个体的特征向量之间的距离，以最小距离为测试构象的多样性值D^trial；

5.9)计算测试构象C^trial的能量函数值E^trial，并进行如下操作：

5.9.1)如果E^trial小于当前目标构象个体Cⁱ的能量函数值Eⁱ，则测试构象C^trial替换目标构象Cⁱ；

5.9.2)如果E^trial大于当前目标构象个体Cⁱ的能量函数值Eⁱ，且测试构象的多样性值D^trial大于目标构象的多样性值Dⁱ，则随机生成一个0与1之间随机数，如果p″>RS，则测试构象C^trial替换目标构象Cⁱ；

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

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

本发明的技术构思为：首先，将结构预测中的整个算法搜索过程分为四个阶段，对每个阶段设置片段长度进行片段组装，并选用不同的能量函数来衡量构象个体的质量；然后，基于二级结构信息，采用不同的变异策略利用loop区域信息来生成测试构象，并通过随机交换loop区域信息实现交叉过程，保持种群多样性，同时对每个阶段的测试构象和目标构象执行Rosetta局部增强；最后，提取构象的特征向量来衡量各构象个体的多样性，从而以能量函数为主要衡量指标，并以多样性为辅助衡量指标来指导构象种群更新。

本发明的有益效果表现在：一方面，基于二级结构信息，根据loop区域的残基操作来实现不同策略的测试构象生成，并对每个测试构象进行Rosetta局部增强，从而提高算法的搜索能力；其次，针对不同阶段的Rosetta局部增强，采用不同的片段长度进行片段组装，并采用不同的能量函数衡量构象的质量，从而提高搜索效率；另一方面，在选择过程中，基于各构象个体之间的特征向量距离来衡量构象的多样性，并将其作为辅助指标来衡量构象的质量，从而在搜索过程充分保持种群多样性，进而提高预测精度。

附图说明

图1是基于Rosetta局部增强的群体蛋白质结构预测方法的流程图。

图2是基于Rosetta局部增强的群体蛋白质结构预测方法对蛋白质1AIL进行结构预测时的构象更新示意图。

图3是基于Rosetta局部增强的群体蛋白质结构预测方法对蛋白质1AIL进行结构预测时得到的构象分布图。

图4是基于基于Rosetta局部增强的群体蛋白质结构预测方法对蛋白质1AIL进行结构预测得到的三维结构图。

具体实施方式

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

参照图1～图4，一种基于Rosetta局部增强的群体蛋白质结构预测方法，包括以下步骤：

1)输入待测蛋白质的氨基酸序列信息；

2)初始化：设置种群规模NP，交叉概率CR，策略选择因子CS，多样性接受概率RS，Rosetta轨迹长度T，片段长度L₁，L₂；

3)根据序列信息以片段长度L₁进行随机片段组装生成初始构象种群P＝{C¹,C²,...,C^NP}，其中，Cⁱ的表示当前种群中的第i个构象个体，并根据能量函数RosettaScore0计算各构象个体的能量，同时初始化迭代次数G＝0；

4)采用能量函数Rosetta Score0评价构象的质量，以片段长度L₁对初始种群中的每个构象个体执行轨迹长度为T的Rosetta局部增强，并计算每个构象的特征向量；

5)对步骤4)中增强后的每个构象个体Cⁱ,i∈{1,2,…,NP}作如下处理：

5.1)设置能量函数和片段长度：

5.1.1)如果当前迭代次数0＜G≤G_max/3，则片段长度l＝L₁，且选用能量函数Rosetta>

5.1.2)如果当前迭代次数G_max/3＜G≤2G_max/3，则片段长度l＝L₁，且选用能量函数Rosetta>

5.1.3)如果当前迭代次数G＞2G_max/3，则片段长度l＝L₂，且选用能量函数RosettaScore3

5.2)如果当前迭代次数G为G_max/3的整数倍，则对以片段长度l对构象个体Cⁱ执行轨迹长度为T的Rosetta局部增强，并根据步骤5.1)中设置的能量函数进行评价；

5.3)计算目标构象Cⁱ的特征向量，以及Cⁱ与当前种群中其他构象之间的特征向量欧氏距离，并以最小距离为Cⁱ的多样性值Dⁱ；

5.4)根据序列信息，利用DSSP得到待测蛋白的loop区域，并随机生成一个0到1之间的随机数p；

5.5)如果p<CS，则从当前种群中选取三个互不相同的构象个体C^a、C^b和C^c，其中，a≠b≠c≠i，a,b,c∈[1,NP]，从构象个体C^b和C^c中各随机选取一个片段替换C^a中对应位置的片段，并从C^a中随机选取一个不包含loop区域的窗口进行片段组装生成变异构象C^mutant；

5.6)如果p≥1-CS，则选出当前能量值最低的构象个体C^best，并从当前种群中选取两个互不相同的构象个体C^a和C^b，其中，a≠b≠i，a,b∈[1,NP]，从构象个体C^a和C^b中各随机选取一个片段替换C^best中对应位置的片段，并从C^best中随机选取一个不包含loop区域的窗口进行片段组装生成变异构象C^mutant；

5.7)随机生成一个0与1之间随机数p′，如果p′＞CR，则随机选取一个loop区域，替换目标构象个体Cⁱ与变异构象个体C^mutant在该区域的二面角，从而生成测试构象C^trial，否则C^trial直接等于变异构象C^mutant；以片段长度l对测试构象个体C^trial执行轨迹长度为T的Rosetta局部增强；

5.8)计算增强后测试构象的特征向量，并计算测试构象的特征向量与当前种群中各构象个体的特征向量之间的距离，以最小距离为测试构象的多样性值D^trial；

5.9)计算测试构象C^trial的能量函数值E^trial，并进行如下操作：

5.9.1)如果E^trial小于当前目标构象个体Cⁱ的能量函数值Eⁱ，则测试构象C^trial替换目标构象Cⁱ；

5.9.2)如果E^trial大于当前目标构象个体Cⁱ的能量函数值Eⁱ，且测试构象的多样性值D^trial大于目标构象的多样性值Dⁱ，则随机生成一个0与1之间随机数，如果p″>RS，则测试构象C^trial替换目标构象Cⁱ；

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

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

本实施例序列长度为56的α/β折叠蛋白质1GB1为实施例，一种基于Rosetta局部增强的群体蛋白质结构预测方法，其中包含以下步骤：

1)输入待测蛋白质的氨基酸序列信息；

2)初始化：设置种群规模NP＝100，交叉概率CR＝0.5，策略选择因子CS＝0.5，多样性接受概率RS＝0.5，Rosetta轨迹长度T＝1000，最大迭代次数G_max＝1200，片段长度L₁＝3，L₂＝9；

3)根据序列信息以片段长度L₁进行随机片段组装生成初始构象种群P＝{C¹,C²,...,C^NP}，其中，Cⁱ的表示当前种群中的第i个构象个体，并根据能量函数RosettaScore0计算各构象个体的能量，同时初始化迭代次数G＝0；

4)采用能量函数Rosetta Score0评价构象的质量，以片段长度L₁对初始种群中的每个构象个体执行轨迹长度为T的Rosetta局部增强，并计算每个构象的特征向量；

5)对步骤4)中增强后的每个构象个体Cⁱ,i∈{1,2,…,NP}作如下处理：

5.1)设置能量函数和片段长度：

5.1.1)如果当前迭代次数0＜G≤G_max/3，则片段长度l＝L₁，且选用能量函数Rosetta>

5.1.2)如果当前迭代次数G_max/3＜G≤2G_max/3，则片段长度l＝L₁，且选用能量函数Rosetta>

5.1.3)如果当前迭代次数G＞2G_max/3，则片段长度l＝L₂，且选用能量函数RosettaScore3

5.2)如果当前迭代次数G为G_max/3的整数倍，则对以片段长度l对构象个体Cⁱ执行轨迹长度为T的Rosetta局部增强，并根据步骤5.1)中设置的能量

函数进行评价；

5.3)计算目标构象Cⁱ的特征向量，以及Cⁱ与当前种群中其他构象之间的特征向量欧氏距离，并以最小距离为Cⁱ的多样性值Dⁱ；

5.4)根据序列信息，利用DSSP得到待测蛋白的loop区域，并随机生成一个0到1之间的随机数p；

5.5)如果p<CS，则从当前种群中选取三个互不相同的构象个体C^a、C^b和C^c，其中，a≠b≠c≠i，a,b,c∈[1,NP]，从构象个体C^b和C^c中各随机选取一个片段替换C^a中对应位置的片段，并从C^a中随机选取一个不包含loop区域的窗口进行片段组装生成变异构象C^mutant；

5.6)如果p≥1-CS，则选出当前能量值最低的构象个体C^best，并从当前种群中选取两个互不相同的构象个体C^a和C^b，其中，a≠b≠i，a,b∈[1,NP]，从构象个体C^a和C^b中各随机选取一个片段替换C^best中对应位置的片段，并从C^best中随机选取一个不包含loop区域的窗口进行片段组装生成变异

构象C^mutant；

5.7)随机生成一个0与1之间随机数p′，如果p′＞CR，则随机选取一个loop区域，替换目标构象个体Cⁱ与变异构象个体C^mutant在该区域的二面角，从而生成测试构象C^trial，否则C^trial直接等于变异构象C^mutant；以片段长度l对测试构象个体C^trial执行轨迹长度为T的Rosetta局部增强；

5.8)计算增强后测试构象的特征向量，并计算测试构象的特征向量与当前种群中各构象个体的特征向量之间的距离，以最小距离为测试构象的多样性值

D^trial；

5.9)计算测试构象C^trial的能量函数值E^trial，并进行如下操作：

5.9.1)如果E^trial小于当前目标构象个体Cⁱ的能量函数值Eⁱ，则测试构象C^trial替换目标构象Cⁱ；

5.9.2)如果E^trial大于当前目标构象个体Cⁱ的能量函数值Eⁱ，且测试构象的多样性值D^trial大于目标构象的多样性值Dⁱ，则随机生成一个0与1之间随机数，如果p″>RS，则测试构象C^trial替换目标构象Cⁱ；

6)当对种群中的每个构象都执行了步骤5)以后，G＝G+1，若G>G_max则输出结

果并退出，否则返回步骤5)。

以序列长度为56的α/β折叠蛋白质1GB1为实施例，运用以上方法得到了该蛋白质的近天然态构象，最小均方根偏差为平均均方根偏差为预测结构如图4所示。

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