一种基于多模态构象空间采样的蛋白质结构预测方法与流程

本发明涉及生物信息学、计算机应用领域，尤其涉及的是一种基于多模态构象空间采样的蛋白质结构预测方法。

背景技术：

蛋白质结构预测问题也称为蛋白质折叠问题。蛋白质折叠结构的形状在很大程度上决定了其生物功能，精确地预测蛋白质的空间结构对生物制药和疾病研究都有重要的意义，利用生化手段测定蛋白质结构代价高，因此利用计算手段预测蛋白质三维结构逐渐成为计算生物学的重要课题。

从头蛋白质结构预测方法是一种常用的蛋白质结构预测方法，由于它仅利用一级序列信息进行预测，而不依赖于已知蛋白质结构模板，所以该方法也是一种较理想的预测方法。从头蛋白质结构预测的理论依据是一定环境中天然蛋白质的三维结构是整个系统自由能最小的结构。因此，从头蛋白质结构预测有两个关键：一是要有一个合理的势能函数，势能函数的全局最小点对应于蛋白质的天然结构；二是要有一个高效的构象空间搜索算法，保证在有效的计算时间内找到势能函数的全局最小值。在蛋白质结构从头预测的过程中，能量函数的不精确以及采样能力的不足导致预测结果并不理想。

实际优化问题很多都属于多模态函数优化问题,通常需要求出多个解,包括全局最优解和局部最优解。在蛋白质结构与预测问题中，局部最优解极有可能对应于更高精度的结构，使用多模态优化方法可以缓解能量函数的不准确性。在使用多模态优化方法时涉及到一个关键问题：如何将一个种群划分为多个模态？

通常采用聚类的方法将种群划分为多个模态，但每种聚类方法多少都存在一定的缺陷。k-均值聚类算法是最常用的聚类算法之一，其最大的缺点是需要人为设定k值；然而很多时候很难确定数据集应该分成多少个类别才最合适。在蛋白质结构预测问题中，需要根据两个结构的相似度去划分模态；使用rmsd计算两个结构的相似度过多的关注于局部结构，不能很好的从全局结构出发去评判相似度，所以如何设计精确的相似度指标是非常关键的。

技术实现要素：

为了克服现有的蛋白质结构预测方法能量函数的不准确性，本发明提出一种基于多模态构象空间采样的蛋白质结构预测方法。本发明设计了一种新的相似度评判指标，使用接触图对两个结构的整体拓扑相似度进行打分，忽略局部结构的一些微小差异，旨在从全局结构评判相似性，从而更精确的对种群进性划分。得到种群个体相似度得分矩阵后采用n-1阶近邻方法确定k值，n-1阶近邻方法通过分析群体个体间的n-1阶最短近邻距离来计算种群的全局分布,并利用阶跃信息自动确定k值，弥补k-均值聚类算法的缺陷，从而有效的把种群划分为多个模态，求解出多个可能的最优解。

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

一种基于多模态构象空间采样的蛋白质结构预测方法，所述方法包括以下步骤：

1)输入预测蛋白质的序列信息，读取序列长度l；设置参数：种群规模n，迭代次数g，阈值ε；

2)根据目标蛋白序列信息，利用robetta(http://robetta.bakerlab.org)构建片段库；

3)迭代rosetta第一、二阶段，生成具有n个个体的初始种群

4)按照如下方式计算两个个体之间的相似度：

4.1)对于每一个个体，计算每一对氨基酸之间的距离，当两个氨基酸之间的距离小于等于时，则定义为接触，记为1，否则记为0，最终构建出一个l×l的接触矩阵，记两个个体的接触矩阵分别为m¹、m²；

4.2)按照如下公式给形成的接触分配权重：

其中i和j均∈{1,2，…l}，表示氨基酸序号，|i-j|表示两个氨基酸在序列上的距离；

4.3)将两个个体的接触矩阵重叠，按如下公式计算两个接触矩阵形成交集的点的分数总和，记为soverlap：

按如下公式计算两个接触矩阵形成并集的点的分数总和，记为stotal：

将两个体的相似度分数为：

5)根据步骤4)计算种群中所有个体两两之间的相似度；

6)运用n-1阶近邻方法确定模态数目：

6.1)把两个个体之间的相似度分数作为两个个体之间的距离v，令v1ⁱ为个体pi^init的一阶最短近邻距离，为个体pi^init的二阶最短近邻距离，以此类推则有：

6.2)按如下公式分别计算j阶最短近邻距离的平方均值和平均值

6.3)计算n-1阶方差：

6.4)根据计算出的n-1个方差值按如下公式统计方差的阶跃点数目，当公式成立时，阶跃点数目加一：

6.5)统计得出阶跃点数目为u，则模态数目为k＝u+1；

7)运用k-均值聚类算法对种群进行聚类，形成k个模态；根据这k个模态将种群划分成k个子种群，记作

8)设g＝1，g∈{1,2,...,g}：

9)对每个模态中的每个个体依次进行交叉、变异和排挤操作，过程如下：

9.1)交叉操作：从目标个体pi^k所在的模态中随机选择两个互不相同且与pi^k不同的个体，记作在[1,l-2]内生成两个不同的随机整数rand1和rand2；将pi^k的第rand1至rand1+2号残基的二面角值替换成对应残基的二面角值，将pi^k的第rand2至rand2+2号残基的二面角值替换成对应残基的二面角值，形成交叉后的个体pi^k′；

9.2)变异操作：在[1,l-2]内生成与rand1和rand2不同的随机整数rand3；在rand3对应的滑动窗口的片段库中随机选一个片段，替换pi^k′的第rand3至rand3+2号残基，形成变异后的个体pi^k″；

9.3)排挤操作：用rosettascore3能量函数计算个体pi^k″的能量值，按照步骤4)的方法计算个体pi^k″与该模态中所有个体的相似度，找到与个体pi^k″最相似的个体，若该个体能量值比pi^k″高，则替换掉该个体；

10)g＝g+1；若g≤g，转至步骤9)；

11)选出每个模态中能量最低的个体作为最终输出结果。

本发明的有益效果为：设计了一种新的相似度评判指标，根据接触图匹配度评判全局结构的相似性，从而更精确的对种群进行簇划分。采用n-1阶近邻方法通过分析群体个体间的n-1阶最短近邻确定聚类簇数目，有效的把种群划分为多个模态，从而求解出多个可能的近天然态结构。

附图说明

图1是一种基于多模态构象空间采样的蛋白质结构预测方法对蛋白质1dhn进行结构预测得到的三维结构图。

图2是一种基于多模态构象空间采样的蛋白质结构预测方法对蛋白质1dhn进行结构预测的流程图。

具体实施方式

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

参照图1和图2，一种基于多模态构象空间采样的蛋白质结构预测方法，包括以下步骤：

1)输入预测蛋白质的序列信息，读取序列长度l；设置参数：种群规模n，迭代次数g，阈值ε；

2)根据目标蛋白序列信息，利用robetta(http://robetta.bakerlab.org)构建片段库；

3)迭代rosetta第一、二阶段，生成具有n个个体的初始种群

4)按照如下方式计算两个个体之间的相似度：

4.1)对于每一个个体，计算每一对氨基酸之间的距离，当两个氨基酸之间的距离小于等于时，则定义为接触，记为1，否则记为0，最终构建出一个l×l的接触矩阵，记两个个体的接触矩阵分别为m¹、m²；

4.2)按照如下公式给形成的接触分配权重：

其中i和j均∈{1,2，…l}，表示氨基酸序号，|i-j|表示两个氨基酸在序列上的距离；

4.3)将两个个体的接触矩阵重叠，按如下公式计算两个接触矩阵形成交集的点的分数总和，记为soverlap：

按如下公式计算两个接触矩阵形成并集的点的分数总和，记为stotal：

将两个体的相似度分数为：

5)根据步骤4)计算种群中所有个体两两之间的相似度；

6)运用n-1阶近邻方法确定模态数目：

6.1)把两个个体之间的相似度分数作为两个个体之间的距离v，令v1ⁱ为个体pi^init的一阶最短近邻距离，为个体pi^init的二阶最短近邻距离，以此类推则有：

6.2)按如下公式分别计算j阶最短近邻距离的平方均值和平均值

6.3)计算n-1阶方差：

6.4)根据计算出的n-1个方差值按如下公式统计方差的阶跃点数目，当公式成立时，阶跃点数目加一：

6.5)统计得出阶跃点数目为u，则模态数目为k＝u+1；

7)运用k-均值聚类算法对种群进行聚类，形成k个模态；根据这k个模态将种群划分成k个子种群，记作

8)设g＝1，g∈{1,2,...,g}：

9)对每个模态中的每个个体依次进行交叉、变异和排挤操作，过程如下：

9.1)交叉操作：从目标个体pi^k所在的模态中随机选择两个互不相同且与pi^k不同的个体，记作在[1,l-2]内生成两个不同的随机整数rand1和rand2；将pi^k的第rand1至rand1+2号残基的二面角值替换成对应残基的二面角值，将pi^k的第rand2至rand2+2号残基的二面角值替换成对应残基的二面角值，形成交叉后的个体pi^k′；

9.2)变异操作：在[1,l-2]内生成与rand1和rand2不同的随机整数rand3；在rand3对应的滑动窗口的片段库中随机选一个片段，替换pi^k′的第rand3至rand3+2号残基，形成变异后的个体pi^k″；

9.3)排挤操作：用rosettascore3能量函数计算个体pi^k″的能量值，按照步骤4)的方法计算个体pi^k″与该模态中所有个体的相似度，找到与个体pi^k″最相似的个体，若该个体能量值比pi^k″高，则替换掉该个体；

10)g＝g+1；若g≤g，转至步骤9)；

11)选出每个模态中能量最低的个体作为最终输出结果。

本实施例以序列长度为121的蛋白质1dhn为实施例，一种基于多模态构象空间采样的蛋白质结构预测方法，包括以下步骤：

1)输入预测蛋白质的序列信息，读取序列长度l＝121；设置参数：种群规模

n＝200，迭代次数g＝300，阈值ε＝0.1；

2)根据目标蛋白序列信息，利用robetta(http://robetta.bakerlab.org)构建片段库；

3)迭代rosetta第一、二阶段，生成具有n个个体的初始种群

4)按照如下方式计算两个个体之间的相似度：

4.1)对于每一个个体，计算每一对氨基酸之间的距离，当两个氨基酸之间的距离小于等于时，则定义为接触，记为1，否则记为0，最终构建出一个l×l的接触矩阵，记两个个体的接触矩阵分别为m¹、m²；

4.2)按照如下公式给形成的接触分配权重：

其中i和j均∈{1,2，…l}，表示氨基酸序号，|i-j|表示两个氨基酸在序列上的距离；

4.3)将两个个体的接触矩阵重叠，按如下公式计算两个接触矩阵形成交集的点的分数总和，记为soverlap：

按如下公式计算两个接触矩阵形成并集的点的分数总和，记为stotal：

将两个体的相似度分数为：

5)根据步骤4)计算种群中所有个体两两之间的相似度；

6)运用n-1阶近邻方法确定模态数目：

6.1)把两个个体之间的相似度分数作为两个个体之间的距离v，令v1ⁱ为个体pi^init的一阶最短近邻距离，为个体pi^init的二阶最短近邻距离，以此类推则有：

6.2)按如下公式分别计算j阶最短近邻距离的平方均值和平均值

6.3)计算n-1阶方差：

6.4)根据计算出的n-1个方差值按如下公式统计方差的阶跃点数目，当公式成立时，阶跃点数目加一：

6.5)统计得出阶跃点数目为u，则模态数目为k＝u+1；

7)运用k-均值聚类算法对种群进行聚类，形成k个模态；根据这k个模态将种群划分成k个子种群，记作

8)设g＝1，g∈{1,2,...,g}：

9)对每个模态中的每个个体依次进行交叉、变异和排挤操作，过程如下：

9.1)交叉操作：从目标个体pi^k所在的模态中随机选择两个互不相同且与pi^k不同的个体，记作在[1,l-2]内生成两个不同的随机整数rand1和rand2；将pi^k的第rand1至rand1+2号残基的二面角值替换成对应残基的二面角值，将pi^k的第rand2至rand2+2号残基的二面角值替换成对应残基的二面角值，形成交叉后的个体pi^k′；

9.2)变异操作：在[1,l-2]内生成与rand1和rand2不同的随机整数rand3；在rand3对应的滑动窗口的片段库中随机选一个片段，替换pi^k′的第rand3至rand3+2号残基，形成变异后的个体pi^k″；

9.3)排挤操作：用rosettascore3能量函数计算个体pi^k″的能量值，按照步骤4)的方法计算个体pi^k″与该模态中所有个体的相似度，找到与个体pi^k″最相似的个体，若该个体能量值比pi^k″高，则替换掉该个体；

10)g＝g+1；若g≤g，转至步骤9)；

11)选出每个模态中能量最低的个体作为最终输出结果。

以氨基酸序列长度为121的蛋白质1dhn为实施例，运用以上方法得到了该蛋白质在五个模态中的近天然态个体，预测的蛋白质的均方根偏差分别为预测结构如图1所示。

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