一种基于差分进化局部扰动的蛋白质构象空间优化方法与流程

本发明涉及生物信息学、计算机应用领域，尤其涉及的是一种基于差分进化局部扰动的蛋白质构象空间优化方法。

背景技术：

蛋白质是以氨基酸为基本单位构成的生物大分子。蛋白质是生命的物质基础，是构成细胞的基本有机物。生物体内的氨基酸种类有20种，20种氨基酸排列组合形成了各种各样具有特定功能的蛋白质。蛋白质的功能由其空间结构决定，许多疾病就是由于蛋白质的空间结构发生变化导致的。因此，确定蛋白质的空间结构将有助于相关疾病的治疗。

由于通过实验测定蛋白质三维结构的方法存在局限，因此利用计算机技术根据蛋白质的氨基酸序列预测蛋白质的三维结构成为测定蛋白质空间结构的另一种选择。根据氨基酸序列预测蛋白质三维结构的方法主要分为同源建模法和从头预测法。其中从头预测法不依赖目标蛋白的同源信息。目前比较成功的从头蛋白质结构预测方法有baker团队开发的rosetta和张阳团队开发的quark等。

根据氨基酸序列预测蛋白质三维结构实质上是在能量模型引导下的构象空间优化问题。目前的构象空间优化方法对α螺旋和β折叠的探索能力较强，往往能形成较高精度的α螺旋和β折叠；但是这些方法对loop区域的探索能力不够，不能形成很好的loop结构，从而影响整体的预测精度。

因此，目前的构象空间优化方法对loop区域的探索存在搜索效率低和预测精度不足等问题，需要改进。

技术实现要素：

为了克服现有的构象空间优化方法对loop区域的探索存在搜索效率低和预测精度不足等问题，本发明提供一种预测精度较高的一种基于差分进化局部扰动的蛋白质构象空间优化方法，本方法在大规模搜索构象空间的同时增强对loop区域的探索，利用差分进化算法实现对loop区域的扰动，从而提高loop区域结构的多样性，在已有结构的基础上进一步优化loop区域的结构，从而提高整体的预测精度。

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

一种基于差分进化局部扰动的蛋白质构象空间优化方法，所述方法包括以下步骤：

1)输入预测蛋白质的序列信息；

2)设置参数，过程如下：

2.1)构象搜索差分进化参数：种群规模np，迭代次数g；

2.2)二面角扰动差分进化参数：种群规模np′，迭代次数g′，交叉概率cr，变异算子f；

3)种群初始化：迭代rosetta协议第一、二、三阶段，产生具有np个个体的种群p＝{p1,p2,...,pnp}，其中pn表示种群p中的第n个个体，n∈{1,2,...,np}；

4)设g＝1，其中g∈{1,2,...,g}；

5)设n＝1，其中n∈{1,2,...,np}；

6)变异操作，过程如下：

6.1)从种群p中随机选择三个互不相同的个体pselect1、pselect2、pselect3；

6.2)在[0,l-3]内生成两个不同的均匀随机整数rand1和rand2，其中l表示氨基酸序列的长度；

6.3)将pselect1的第rand1至rand1+2号残基的二面角值替换成pselect2对应残基的二面角值，将pselect1的第rand2至rand2+2号残基的二面角值替换成pselect3对应残基的二面角值，形成变异后的个体p^mutate；

7)交叉操作，过程如下：

7.1)在[0,l-3]内生成与rand1和rand2不同的均匀随机整数rand3；

7.2)将p^mutate的第rand3至rand3+2号残基的二面角值替换成pn对应残基的二面角值，形成交叉后的个体p^cross；

8)利用psipred算法获取p^cross的二级结构信息，随机选择一段二级结构为loop且长度大于4的区域，记作[b,e]，其中b和e分别表示该区域第一个残基和最后一个残基的残基号；

9)利用rosetta协议第四阶段使用的片段组装技术对p^cross的[b,e]区域执行一次片段组装，得到片段组装后的个体p^assem；

10)对p^assem的[b,e]区域进行局部扰动，过程如下：

10.1)以p^assem第b至e号残基二面角和ψ对应的键作为旋转轴，取键方向上的单位向量得到旋转轴坐标，记作ω＝{ω1,ω2,...,ωs}，其中s表示旋转轴数量，ωs＝(xs,ys,zs)表示第s个旋转轴的坐标，s∈{1,2,...,s}；

10.2)记旋转角度φ＝{φ1,φ2,...,φs}，其中φs表示第s个旋转轴对应的旋转角度，即p^assem的[b,e]区域第s个二面角的扰动角度；

10.3)以p^assem第e至l号残基的cα原子作为旋转点，记作x＝{χ1,χ2,...,χk}，其中k表示旋转点数量，χk＝(xk,yk,zk)表示第k个旋转点的坐标，k∈{1,2,...,k}；

10.4)模型评估，过程如下：

10.4.1)将x中的旋转点χk依次按照旋转角φs绕旋转轴ws旋转，其中k∈{1,2,...,k}，s∈{s,s-1,...,1}；得到旋转后的x′＝{χ′1,χ′2,...,χ′k}，χ′k＝(x′k,y′k,z′k)；

10.4.2)打分函数，按如下公式计算φ对应的得分：

10.5)求解最优的φ，使其得分最低，过程如下：

10.5.1)初始化：生成初始种群γ＝{φ1,φ2,...,φnp′}，其中

10.5.2)设g′＝1，其中g′∈{1,2,...,g′}；

10.5.3)设n′＝1，其中n′∈{1,2,...,np′}；

10.5.4)变异操作，过程如下：

10.5.4.1)从种群γ中选择得分最低的个体，记作φbest；再随机选择两个互不相同且与φbest不同的个体，记作φselect1、φselect2；

10.5.4.2)按如下方式生成变异后的个体φ^mutate：

10.5.5)交叉操作，过程如下：

10.5.5.1)在[1,s]内生成均匀随机整数srand；

10.5.5.2)生成s个均匀随机小数rs，rs∈[0,1]，s∈{1,2,...,s}；

10.5.5.3)按如下方式生成交叉后的个体φ^cross：

10.5.6)选择操作，过程如下：

10.5.6.1)根据步骤10.4)计算φn′和φ^cross的得分d^n′和d^cross；

10.5.6.2)若e^cross＜e^n′，用φ^cross替换种群γ中的个体φn′；

10.5.7)n′＝n′+1；若n′≤np′，转至步骤10.5.4)；

10.5.8)g′＝g′+1；若g′≤g′，转至步骤10.5.3)；

10.5.9)从种群中选择五个能量最低的个体作为候选扰动角度，记作h＝{φ1,φ2,φ3,φ4,φ5}；

10.6)选择扰动，过程如下：

10.6.1)设i＝1，其中i∈{1,2,3,4,5}；

10.6.2)将h中的第i个个体φi的扰动角度添加到p^assem对应的二面角上，形成扰动后的个体p^perturb；

10.6.3)用rosettascore3能量函数计算p^assem和p^perturb的能量e^assem和e^perturb，并按照如下公式计算接收扰动的概率：

10.6.4)生成均匀随机小数rand4，rand4∈[0,1]；若pacc＜rand4，令p^trial＝p^perturb，转至步骤11)；否则，转至步骤10.6.5)；

10.6.5)i＝i+1；若i≤5，转至步骤10.6.2)；否则，令p^trial＝p^assem；

11)用rosettascore3能量函数计算p^trial和pn的能量，根据metropolis准则决定是否用p^trial替换种群p中的个体pn；

12)n＝n+1；若n≤np，转至步骤6)；

13)g＝g+1；若g≤g，转至步骤5)；

14)利用spicker聚类算法输出最终的预测结果。

本发明的有益效果为：在差分进化算法的框架下，利用种群中个体间的信息交换增强算法的探索能力；同时，利用差分进化算法实现对loop区域的微调，增加loop区域结构的多样性，从而在已有的结构的基础上进一步增强对loop区域的探索，从而提高整体的探索效率和预测精度。

附图说明

图1是一种基于差分进化局部扰动的蛋白质构象空间优化方法对蛋白质4uex进行结构预测时的构象更新示意图。

图2是一种基于差分进化局部扰动的蛋白质构象空间优化方法对蛋白质4uex进行结构预测得到的三维结构图。

具体实施方式

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

参照图1和图2，一种基于差分进化局部扰动的蛋白质构象空间优化方法，包括以下步骤：

1)输入预测蛋白质的序列信息；

2)设置参数，过程如下：

2.1)构象搜索差分进化参数：种群规模np，迭代次数g；

2.2)二面角扰动差分进化参数：种群规模np′，迭代次数g′，交叉概率cr，变异算子f；

3)种群初始化：迭代rosetta协议第一、二、三阶段，产生具有np个个体的种群p＝{p1,p2,...,pnp}，其中pn表示种群p中的第n个个体，n∈{1,2,...,np}；

4)设g＝1，其中g∈{1,2,...,g}；

5)设n＝1，其中n∈{1,2,...,np}；

6)变异操作，过程如下：

6.1)从种群p中随机选择三个互不相同的个体pselect1、pselect2、pselect3；

6.2)在[0,l-3]内生成两个不同的均匀随机整数rand1和rand2，其中l表示氨基酸序列的长度；

6.3)将pselect1的第rand1至rand1+2号残基的二面角值替换成pselect2对应残基的二面角值，将pselect1的第rand2至rand2+2号残基的二面角值替换成pselect3对应残基的二面角值，形成变异后的个体p^mutate；

7)交叉操作，过程如下：

7.1)在[0,l-3]内生成与rand1和rand2不同的均匀随机整数rand3；

7.2)将p^mutate的第rand3至rand3+2号残基的二面角值替换成pn对应残基的二面角值，形成交叉后的个体p^cross；

8)利用psipred算法获取p^cross的二级结构信息，随机选择一段二级结构为loop且长度大于4的区域，记作[b,e]，其中b和e分别表示该区域第一个残基和最后一个残基的残基号；

9)利用rosetta协议第四阶段使用的片段组装技术对p^cross的[b,e]区域执行一次片段组装，得到片段组装后的个体p^assem；

10)对p^assem的[b,e]区域进行局部扰动，过程如下：

10.1)以p^assem第b至e号残基二面角和ψ对应的键作为旋转轴，取键方向上的单位向量得到旋转轴坐标，记作ω＝{ω1,ω2,...,ωs}，其中s表示旋转轴数量，ωs＝(xs,ys,zs)表示第s个旋转轴的坐标，s∈{1,2,...,s}；

10.2)记旋转角度φ＝{φ1,φ2,...,φs}，其中φs表示第s个旋转轴对应的旋转角度，即p^assem的[b,e]区域第s个二面角的扰动角度；

10.3)以p^assem第e至l号残基的cα原子作为旋转点，记作x＝{χ1,χ2,...,χk}，其中k表示旋转点数量，χk＝(xk,yk,zk)表示第k个旋转点的坐标，k∈{1,2,...,k}；

10.4)模型评估，过程如下：

10.4.1)将x中的旋转点χk依次按照旋转角φs绕旋转轴ws旋转，其中k∈{1,2,...,k}，s∈{s,s-1,...,1}；得到旋转后的x′＝{χ′1,χ′2,...,χ′k}，χ′k＝(x′k,yk′,z′k)；

10.4.2)打分函数，按如下公式计算φ对应的得分：

10.5)求解最优的φ，使其得分最低，过程如下：

10.5.1)初始化：生成初始种群γ＝{φ1,φ2,...,φnp′}，其中

10.5.2)设g′＝1，其中g′∈{1,2,...,g′}；

10.5.3)设n′＝1，其中n′∈{1,2,...,np′}；

10.5.4)变异操作，过程如下：

10.5.4.1)从种群γ中选择得分最低的个体，记作φbest；再随机选择两个互不相同且与φbest不同的个体，记作φselect1、φselect2；

10.5.4.2)按如下方式生成变异后的个体φ^mutate：

10.5.5)交叉操作，过程如下：

10.5.5.1)在[1,s]内生成均匀随机整数srand；

10.5.5.2)生成s个均匀随机小数rs，rs∈[0,1]，s∈{1,2,...,s}；

10.5.5.3)按如下方式生成交叉后的个体φ^cross：

10.5.6)选择操作，过程如下：

10.5.6.1)根据步骤10.4)计算φn′和φ^cross的得分d^n′和d^cross；

10.5.6.2)若e^cross＜e^n′，用φ^cross替换种群γ中的个体φn′；

10.5.7)n′＝n′+1；若n′≤np′，转至步骤10.5.4)；

10.5.8)g′＝g′+1；若g′≤g′，转至步骤10.5.3)；

10.5.9)从种群中选择五个能量最低的个体作为候选扰动角度，记作h＝{φ1,φ2,φ3,φ4,φ5}；

10.6)选择扰动，过程如下：

10.6.1)设i＝1，其中i∈{1,2,3,4,5}；

10.6.2)将h中的第i个个体φi的扰动角度添加到p^assem对应的二面角上，形成扰动后的个体p^perturb；

10.6.3)用rosettascore3能量函数计算p^assem和p^perturb的能量e^assem和e^perturb，并按照如下公式计算接收扰动的概率：

10.6.4)生成均匀随机小数rand4，rand4∈[0,1]；若pacc＜rand4，令p^trial＝p^perturb，转至步骤11)；否则，转至步骤10.6.5)；

10.6.5)i＝i+1；若i≤5，转至步骤10.6.2)；否则，令p^trial＝p^assem；

11)用rosettascore3能量函数计算p^trial和pn的能量，根据metropolis准则决定是否用p^trial替换种群p中的个体pn；

12)n＝n+1；若n≤np，转至步骤6)；

13)g＝g+1；若g≤g，转至步骤5)；

14)利用spicker聚类算法输出最终的预测结果。

本实施例以序列长度为81的蛋白质4uex为实施例，一种基于差分进化局部扰动的蛋白质构象空间优化方法，包括以下步骤：

1)输入预测蛋白质4uex的序列信息；

2)设置参数，过程如下：

2.1)构象搜索差分进化参数：种群规模np＝100，迭代次数g＝100；

2.2)二面角扰动差分进化参数：种群规模np′＝100，迭代次数g′＝1000，交叉概率cr＝0.1，变异算子f＝0.5；

3)种群初始化：迭代rosetta协议第一、二、三阶段，产生具有np个个体的种群p＝{p1,p2,...,pnp}，其中pn表示种群p中的第n个个体，n∈{1,2,...,np}；

4)设g＝1，其中g∈{1,2,...,g}；

5)设n＝1，其中n∈{1,2,...,np}；

6)变异操作，过程如下：

6.1)从种群p中随机选择三个互不相同的个体pselect1、pselect2、pselect3；

6.2)在[0,l-3]内生成两个不同的均匀随机整数rand1和rand2，其中l表示氨基酸序列的长度；

6.3)将pselect1的第rand1至rand1+2号残基的二面角值替换成pselect2对应残基的二面角值，将pselect1的第rand2至rand2+2号残基的二面角值替换成pselect3对应残基的二面角值，形成变异后的个体p^mutate；

7)交叉操作，过程如下：

7.1)在[0,l-3]内生成与rand1和rand2不同的均匀随机整数rand3；

7.2)将p^mutate的第rand3至rand3+2号残基的二面角值替换成pn对应残基的二面角值，形成交叉后的个体p^cross；

8)利用psipred算法获取p^cross的二级结构信息，随机选择一段二级结构为loop且长度大于4的区域，记作[b,e]，其中b和e分别表示该区域第一个残基和最后一个残基的残基号；

9)利用rosetta协议第四阶段使用的片段组装技术对p^cross的[b,e]区域执行一次片段组装，得到片段组装后的个体p^assem；

10)对p^assem的[b,e]区域进行局部扰动，过程如下：

10.1)以p^assem第b至e号残基二面角和ψ对应的键作为旋转轴，取键方向上的单位向量得到旋转轴坐标，记作ω＝{ω1,ω2,...,ωs}，其中s表示旋转轴数量，ωs＝(xs,ys,zs)表示第s个旋转轴的坐标，s∈{1,2,...,s}；

10.2)记旋转角度φ＝{φ1,φ2,...,φs}，其中φs表示第s个旋转轴对应的旋转角度，即p^assem的[b,e]区域第s个二面角的扰动角度；

10.3)以p^assem第e至l号残基的cα原子作为旋转点，记作x＝{χ1,χ2,...,χk}，其中k表示旋转点数量，χk＝(xk,yk,zk)表示第k个旋转点的坐标，k∈{1,2,...,k}；

10.4)模型评估，过程如下：

10.4.1)将x中的旋转点χk依次按照旋转角φs绕旋转轴ws旋转，其中k∈{1,2,...,k}，s∈{s,s-1,...,1}；得到旋转后的x′＝{χ′1,χ′2,...,χ′k}，χ′k＝(x′k,y′k,z′k)；

10.4.2)打分函数，按如下公式计算φ对应的得分：

10.5)求解最优的φ，使其得分最低，过程如下：

10.5.1)初始化：生成初始种群γ＝{φ1,φ2,...,φnp′}，其中

10.5.2)设g′＝1，其中g′∈{1,2,...,g′}；

10.5.3)设n′＝1，其中n′∈{1,2,...,np′}；

10.5.4)变异操作，过程如下：

10.5.4.1)从种群γ中选择得分最低的个体，记作φbest；再随机选择两个互不相同且与φbest不同的个体，记作φselect1、φselect2；

10.5.4.2)按如下方式生成变异后的个体φ^mutate：

10.5.5)交叉操作，过程如下：

10.5.5.1)在[1,s]内生成均匀随机整数srand；

10.5.5.2)生成s个均匀随机小数rs，rs∈[0,1]，s∈{1,2,...,s}；

10.5.5.3)按如下方式生成交叉后的个体φ^cross：

10.5.6)选择操作，过程如下：

10.5.6.1)根据步骤10.4)计算φn′和φ^cross的得分d^n′和d^cross；

10.5.6.2)若e^cross＜e^n′，用φ^cross替换种群γ中的个体φn′；

10.5.7)n′＝n′+1；若n′≤np′，转至步骤10.5.4)；

10.5.8)g′＝g′+1；若g′≤g′，转至步骤10.5.3)；

10.5.9)从种群中选择五个能量最低的个体作为候选扰动角度，记作h＝{φ1,φ2,φ3,φ4,φ5}；

10.6)选择扰动，过程如下：

10.6.1)设i＝1，其中i∈{1,2,3,4,5}；

10.6.2)将h中的第i个个体φi的扰动角度添加到p^assem对应的二面角上，形成扰动后的个体p^perturb；

10.6.3)用rosettascore3能量函数计算p^assem和p^perturb的能量e^assem和e^perturb，并按照如下公式计算接收扰动的概率：

10.6.4)生成均匀随机小数rand4，rand4∈[0,1]；若pacc＜rand4，令p^trial＝p^perturb，转至步骤11)；否则，转至步骤10.6.5)；

10.6.5)i＝i+1；若i≤5，转至步骤10.6.2)；否则，令p^trial＝p^assem；

11)用rosettascore3能量函数计算p^trial和pn的能量，根据metropolis准则决定是否用p^trial替换种群p中的个体pn；

12)n＝n+1；若n≤np，转至步骤6)；

13)g＝g+1；若g≤g，转至步骤5)；

14)利用spicker聚类算法输出最终的预测结果。

以氨基酸序列长度为81的蛋白质4uex为实施例，运用以上方法得到了该蛋白质的近天然态构象，其构象更新示意图如图1所示，预测的蛋白质的均方根偏差为预测结构如图2所示。

以上阐述是本发明给出的一个实施的预测效果，本发明不仅适合上述实施例，在不偏离本发明基本思想及不超出本发明实质内容的前提下可对其做种种改进加以实施。