本发明涉及生物信息学、计算机应用领域，尤其涉及的是一种基于局部扰动的群体蛋白质构象空间优化方法。
背景技术：
蛋白质分子是组成人体一切细胞、组织的重要成分。机体所有重要的组成部分都需要有蛋白质的参与。蛋白质的功能十分丰富，对机体的正常运行起着至关重要的作用。而蛋白质的三维结构决定着蛋白质的功能，蛋白质只有正确的折叠成特定的三维结构才能产生特有的生物功能。由于蛋白质错误折叠引起的疾病并不少见。因此，要了解蛋白质的功能、治愈与蛋白质有关的多种疾病，就必须获得蛋白质的三维结构。
二十世纪末以来，生命科学领域飞速发展，作为生物体内分布最广、功能最复杂的一类大分子，蛋白质受到尤为广泛的关注和研究。蛋白质的三维结构预测是生物信息学的一个重要任务。不同的蛋白质拥有不同的氨基酸序列，所有蛋白质都会在其一维序列的基础上折叠形成特定的三维结构，了解蛋白质的三维结构是研究其生物功能的基础。目前蛋白质的三维结构主要通过x射线晶体衍射和核磁共振成像技术得到，但这两种通过实验的方法获取蛋白质结构的成本均不菲，并且有各自的应用局限。因此，通过计算机技术结合生物信息学的方法模拟蛋白质从氨基酸序列折叠成特定的空间结构的过程，从而预测蛋白质的三维结构就显得尤为重要。根据序列信息预测蛋白质三级结构依然是一个重大的难题。基于片段组装原理的技术是解决这项难题的最先进的方法，但是这种方法对较大、较难的蛋白质的预测精度不够。
目前比较成功的从头蛋白质结构预测方法有davidbaker及其团队设计的rosetta方法、张阳及其团队开发的quark方法等。但是至今还没有一种非常完善的蛋白质三维结构预测方法。目前主要的技术难点来源于两个方面，一方面是现有的构象空间采样方法的采样能力不够强，另一方面在于能量函数的不精确使得仅以能量值作为构象更新的选择标准并不能产生十分令人满意的结果。
因此，目前的蛋白质结构预测方法存在不足，需要改进。
技术实现要素：
为了克服目前的蛋白质结构预测方法的不足，本发明提供一种预测精度较高的基于局部扰动的群蛋白质体构象空间优化方法。首先，利用rosetta协议的第一、二、三阶段对种群进行初始化；片段局部扰动阶段，迭代rosetta协议的第四阶段进行构象搜索，并对构象的片段进行旋转操作；二级结构扰动阶段，对种群中的每一个个体的一个二级结构进行旋转操作，迭代此过程，直到达到给定的迭代次数。
本发明解决其技术问题所采用的技术方案是：
一种基于局部扰动的群体蛋白质构象空间优化方法，所述方法包括以下步骤：
1)输入预测蛋白质的序列信息和二级结构信息，并记录每一个二级结构的第一个和最后一个残基的残基号；
2)设置参数：种群规模np，二级结构扰动迭代次数g；
3)种群初始化：迭代rosetta协议第一、二、三阶段，产生具有np个个体的种群p＝{p1,p2,...,pnp}；
4)片段局部扰动，过程如下：
4.1)设n＝1，其中n∈{1,2,...,np}；
4.2)对种群中第n个个体pn执行rosetta协议第四阶段的3片段组装；
4.3)随机选择构象pn的一个3片段窗口，记该窗口第一个残基的cα原子坐标为(x1,y1,z1)、最后一个残基的cα原子坐标为(x2,y2,z2)；
4.4)令
4.5)以矢量
4.5.1)设r＝1，其中需要进行旋转的原子的编号r∈{1,2,...,r}，r表示需要旋转的原子的总数；
4.5.2)生成[-5,5]的均匀分布随机数μ，作为旋转角；
4.5.3)根据旋转轴和旋转角确定旋转矩阵：
其中，c＝cosμ，s＝sinμ；
4.5.4)记需要旋转的原子的坐标为(xold,yold,zold)，则旋转后的坐标为(xnew,ynew,znew)＝t(xold,yold,zold)，即：
其中ti,j表示矩阵t的第i行第j列的元素，i和j均∈{1,2,3}；
4.5.5)计算构象ptrial和pn的能量，并根据metropolis准则决定是否用ptrial替换pn；
4.5.6)r＝r+1；
4.5.7)若r≤r，转至步骤4.5.2)执行下一个原子的旋转；否则，转至步骤4.6)；
4.6)n＝n+1；
4.7)若n≤np，转至步骤4.2)对下一个个体进行扰动；否则，转至步骤5)；
5)二级结构扰动，过程如下：
5.1)开始迭代，设g＝1，其中迭代次数g∈{1,2,...,g}；
5.2)对当前种群p中的每一个个体进行局部扰动，过程如下：
5.2.1)设m＝1，其中m∈{1,2,...,np}；
5.2.2)随机选择种群p中第m个个体pm的一个二级结构，记该二级结构第一个残基的cα原子坐标为(x1,y1,z1)、最后一个残基的cα原子坐标为(x2,y2,z2)；
5.2.3)令
5.2.4)以矢量
5.2.4.1)设h＝1，其中需要进行旋转的原子的编号h∈{1,2,...,h}，h表示需要旋转的原子的总数；
5.2.4.2)生成[-5,5]的均匀分布随机数γ，作为旋转角；
5.2.4.3)根据旋转轴和旋转角确定旋转矩阵：
其中，c＝cosγ，s＝sinγ；
5.2.4.4)记需要旋转的原子的坐标为(xold,yold,zold)，则旋转后的坐标为(xnew,ynew,znew)＝t(xold,yold,zold)，即：
其中ti,j表示矩阵t的第i行第j列的元素，i和j均∈{1,2,3}；
5.2.4.5)计算构象ptrial和pm的能量，并根据metropolis准则决定是否用ptrial替换pm；
5.2.4.6)h＝h+1；
5.2.4.7)若h≤h，转至步骤5.2.4.2)执行下一个原子的旋转；否则，转至步骤4.2.6)；
5.2.5)m＝m+1；
5.2.6)若m≤np，转至步骤5.2.2)对下一个个体进行扰动；否则，转至步骤5.3)
5.3)g＝g+1；
5.4)若g≤g，转至步骤5.2)执行下一次迭代；否则，转至步骤6)；
6)根据rosetta协议输出预测结果。
本发明的有益效果为：首先利用rosetta协议进行大范围的构象搜索，并在rosetta协议的第四阶段结合对片段的局部扰动，增强对构象空间的探索，此时已形成较为完整的二级结构；然后对二级结构进行扰动，使二级结构的空间排列更加接近蛋白质的真实结构，形成更接近天然蛋白质的三维结构，从而提高蛋白质结构预测的精度。
附图说明
图1是基于局部扰动的群体蛋白质构象空间优化方法对蛋白质1acf进行结构预测时的构象更新示意图。
图2是基于局部扰动的群体蛋白质构象空间优化方法对蛋白质1acf进行结构预测得到的三维结构图。
具体实施方式
下面结合附图对本发明作进一步描述。
参照图1和图2，一种基于局部扰动的群体蛋白质构象空间优化方法，包括以下步骤：
1)输入预测蛋白质的序列信息和二级结构信息，并记录每一个二级结构的第一个和最后一个残基的残基号；
2)设置参数：种群规模np，二级结构扰动迭代次数g；
3)种群初始化：迭代rosetta协议第一、二、三阶段，产生具有np个个体的种群p＝{p1,p2,...,pnp}；
4)片段局部扰动，过程如下：
4.1)设n＝1，其中n∈{1,2,...,np}；
4.2)对种群中第n个个体pn执行rosetta协议第四阶段的3片段组装；
4.3)随机选择构象pn的一个3片段窗口，记该窗口第一个残基的cα原子坐标为(x1,y1,z1)、最后一个残基的cα原子坐标为(x2,y2,z2)；
4.4)令
4.5)以矢量
4.5.1)设r＝1，其中需要进行旋转的原子的编号r∈{1,2,...,r}，r表示需要旋转的原子的总数；
4.5.2)生成[-5,5]的均匀分布随机数μ，作为旋转角；
4.5.3)根据旋转轴和旋转角确定旋转矩阵：
其中，c＝cosμ，s＝sinμ；
4.5.4)记需要旋转的原子的坐标为(xold,yold,zold)，则旋转后的坐标为(xnew,ynew,znew)＝t(xold,yold,zold)，即：
其中ti,j表示矩阵t的第i行第j列的元素，i和j均∈{1,2,3}；
4.5.5)计算构象ptrial和pn的能量，并根据metropolis准则决定是否用ptrial替换pn；
4.5.6)r＝r+1；
4.5.7)若r≤r，转至步骤4.5.2)执行下一个原子的旋转；否则，转至步骤4.6)；
4.6)n＝n+1；
4.7)若n≤np，转至步骤4.2)对下一个个体进行扰动；否则，转至步骤5)；
5)二级结构扰动，过程如下：
5.1)开始迭代，设g＝1，其中迭代次数g∈{1,2,...,g}；
5.2)对当前种群p中的每一个个体进行局部扰动，过程如下：
5.2.1)设m＝1，其中m∈{1,2,...,np}；
5.2.2)随机选择种群p中第m个个体pm的一个二级结构，记该二级结构第一个残基的cα原子坐标为(x1,y1,z1)、最后一个残基的cα原子坐标为(x2,y2,z2)；
5.2.3)令
5.2.4)以矢量
5.2.4.1)设h＝1，其中需要进行旋转的原子的编号h∈{1,2,...,h}，h表示需要旋转的原子的总数；
5.2.4.2)生成[-5,5]的均匀分布随机数γ，作为旋转角；
5.2.4.3)根据旋转轴和旋转角确定旋转矩阵：
其中，c＝cosγ，s＝sinγ；
5.2.4.4)记需要旋转的原子的坐标为(xold,yold,zold)，则旋转后的坐标为(xnew,ynew,znew)＝t(xold,yold,zold)，即：
其中ti,j表示矩阵t的第i行第j列的元素，i和j均∈{1,2,3}；
5.2.4.5)计算构象ptrial和pm的能量，并根据metropolis准则决定是否用ptrial替换pm；
5.2.4.6)h＝h+1；
5.2.4.7)若h≤h，转至步骤5.2.4.2)执行下一个原子的旋转；否则，转至步骤4.2.6)；
5.2.5)m＝m+1；
5.2.6)若m≤np，转至步骤5.2.2)对下一个个体进行扰动；否则，转至步骤5.3)
5.3)g＝g+1；
5.4)若g≤g，转至步骤5.2)执行下一次迭代；否则，转至步骤6)；
6)根据rosetta协议输出预测结果。
本实施例以序列长度为125的蛋白质1acf为实施例，一种基于局部扰动的群体蛋白质构象空间优化方法，其中包含以下步骤：
1)输入预测蛋白质1acf的序列信息和二级结构信息，并记录每一个二级结构的第一个和最后一个残基的残基号；
2)设置参数：种群规模np＝100，二级结构扰动迭代次数g＝500；
3)种群初始化：迭代rosetta协议第一、二、三阶段，产生具有np个个体的种群p＝{p1,p2,...,pnp}；
4)片段局部扰动，过程如下：
4.1)设n＝1，其中n∈{1,2,...,np}；
4.2)对种群中第n个个体pn执行rosetta协议第四阶段的3片段组装；
4.3)随机选择构象pn的一个3片段窗口，记该窗口第一个残基的cα原子坐标为(x1,y1,z1)、最后一个残基的cα原子坐标为(x2,y2,z2)；
4.4)令
4.5)以矢量
4.5.1)设r＝1，其中需要进行旋转的原子的编号r∈{1,2,...,r}，r表示需要旋转的原子的总数；
4.5.2)生成[-5,5]的均匀分布随机数μ，作为旋转角；
4.5.3)根据旋转轴和旋转角确定旋转矩阵：
其中，c＝cosμ，s＝sinμ；
4.5.4)记需要旋转的原子的坐标为(xold,yold,zold)，则旋转后的坐标为(xnew,ynew,znew)＝t(xold,yold,zold)，即：
其中ti,j表示矩阵t的第i行第j列的元素，i和j均∈{1,2,3}；
4.5.5)计算构象ptrial和pn的能量，并根据metropolis准则决定是否用ptrial替换pn；
4.5.6)r＝r+1；
4.5.7)若r≤r，转至步骤4.5.2)执行下一个原子的旋转；否则，转至步骤4.6)；
4.6)n＝n+1；
4.7)若n≤np，转至步骤4.2)对下一个个体进行扰动；否则，转至步骤5)；
5)二级结构扰动，过程如下：
5.1)开始迭代，设g＝1，其中迭代次数g∈{1,2,...,g}；
5.2)对当前种群p中的每一个个体进行局部扰动，过程如下：
5.2.1)设m＝1，其中m∈{1,2,...,np}；
5.2.2)随机选择种群p中第m个个体pm的一个二级结构，记该二级结构第一个残基的cα原子坐标为(x1,y1,z1)、最后一个残基的cα原子坐标为(x2,y2,z2)；
5.2.3)令
5.2.4)以矢量
5.2.4.1)设h＝1，其中需要进行旋转的原子的编号h∈{1,2,...,h}，h表示需要旋转的原子的总数；
5.2.4.2)生成[-5,5]的均匀分布随机数γ，作为旋转角；
5.2.4.3)根据旋转轴和旋转角确定旋转矩阵：
其中，c＝cosγ，s＝sinγ；
5.2.4.4)记需要旋转的原子的坐标为(xold,yold,zold)，则旋转后的坐标为(xnew,ynew,znew)＝t(xold,yold,zold)，即：
其中ti,j表示矩阵t的第i行第j列的元素，i和j均∈{1,2,3}；
5.2.4.5)计算构象ptrial和pm的能量，并根据metropolis准则决定是否用ptrial替换pm；
5.2.4.6)h＝h+1；
5.2.4.7)若h≤h，转至步骤5.2.4.2)执行下一个原子的旋转；否则，转至步骤4.2.6)；
5.2.5)m＝m+1；
5.2.6)若m≤np，转至步骤5.2.2)对下一个个体进行扰动；否则，转至步骤5.3)
5.3)g＝g+1；
5.4)若g≤g，转至步骤5.2)执行下一次迭代；否则，转至步骤6)；
6)根据rosetta协议输出预测结果。
以氨基酸序列长度为125的蛋白质1acf为实施例，运用以上方法得到了该蛋白质的近天然态构象，其构象更新示意图如图1所示,最小均方根偏差为
以上阐述是本发明给出的一个实施例表现出来的预测效果，显然本发明不仅适合上述实施例，在不偏离本发明基本思想及不超出本发明实质内容的前提下可对其做种种改进加以实施。