CN105205347A - 基于bsa-ts算法的蛋白质三维结构预测方法 - Google Patents

Abstract

本发明涉及生物计算蛋白质结构预测领域，设计了一种基于BSA算法与TS算法的混合搜素优化方法。该方法将BSA算法应用到蛋白质三维结构预测中，并对BSA算法与TS算法分别进行了一定的改进。在目前广泛使用的斐波纳契序列和真实蛋白质序列上进行实验，从实验得出的数据和与其他方法的比较结果来看，该方法较好的克服了传统算法对于控制参数的初值比较敏感，并且易陷入局部最优的缺点，具有良好的性能和精度。

Description

基于BSA-TS算法的蛋白质三维结构预测方法

技术领域

本发明涉及BSA算法，TS算法和三维AB非晶格模型，具体讲的是在蛋白质三维AB非晶格模型中，通过对BSA算法及TS算法的改进与混合，以提高单一算法的算法精度，预测到更稳定的蛋白质三维结构，其属于生物信息学蛋白质结构预测领域。

背景技术

研究蛋白质结构，有助于了解蛋白质作用，了解其如何行使其生物功能，认识蛋白质与蛋白质(或其它分子)之间的相互作用，这对于生物学，医学和药学都非常重要。而近年来蛋白质结构测定的试验方法虽然得到很好发展，但仍比较耗时和昂贵，且对于某些不易结晶的蛋白质来说不适用，因此需要发展理论分析方法。随着计算机技术的发展，其逐步成为处理蛋白质分子超大数据的重要工具。

蛋白质三维结构的理论预测方法主要分为以下三步：一是提出能反映氨基酸残基间相互作用和环境等的数学模型；二是根据热力学假设建立一个计算上简单又能正确区分蛋白质天然结构与其它结构的能量函数，三是找到相应模型对应的能量函数最小值的全局优化方法。由于研究人员的努力，关于上面三个步骤分别得到了很多的成果。对于蛋白质结构模型，目前两个应用非常广的模型是HP晶格模型与AB非晶格模型。HP晶格模型使用两种残基——疏水性残基(H)与亲水性残基(P)表示氨基酸链，且各残基分布在堆叠的立方晶格顶点。这个模型的缺点是忽略了在蛋白质折叠过程中很重要的局部反应。AB非晶格模型比HP晶格模型更精确，因为其中氨基酸间键角是任意的。由于蛋白质折叠过程的复杂性，准确的蛋白质能量函数的建立是困难的，目前有很多简化的能量函数被提出。对于全局优化方法，基于HP晶格模型与AB非晶格模型的许多优化蛋白质结构预测算法已经被提出。近年来，基于HP晶格模型的优化算法比如MOSE,ACO,MCMPSO-TS等算法。MOSE算法用一个启发式偏差函数帮助理论的蛋白质结构形成疏水内核，但是这仅仅在简单折叠的蛋白质中高效。ACO算法提出的灵感来源于对蚂蚁寻找食物的行为的观察，其算法框架可以被分成三部分：蚂蚁解的构造，更新信息素，进程守护。MCMPSO-TS算法是一种联接了粒子群算法(PSO)和禁忌算法(TS)以提高全局优化能力的混合算法。基于AB非晶格模型的蛋白质结构预测算法比如PSO,TPSO,LPSO,PGATS,ABC等算法。PSO算法的灵感来源于鸟群迁徙过程中的信息交流，每只鸟不但能记住自己当前找到的距离食物最佳的位置，还可以知道种群中所有鸟儿当前所能找到的最佳位置，由这两个最优值判断食物的最佳位置。TPSO算法联接了PSO算法和TS算法，利用两个算法的优势提高搜索精度，跳出局部最优解。LPSO算法在PSO算法中加入一种称为利维飞行的随机过程以提高算法全局搜索能力。PGATS算法是在GA-PSO基础上提出的，同时加入了改进的禁忌搜索算法的混合算法。ABC算法的提出源于对蜜蜂采蜜行为的观察，此算法的主要特征是不用知道问题的特殊信息，只需要比较解的优劣性。既然设计算法的目的是预测蛋白质结构，那么最重要的判断标准就是每个算法所能预测到的对应能量函数的最低值。因此本文主要关注得到更低的能量值及对应的空间结构，以蛋白质三维结构预测为研究目的，并基于改进BSA算法与禁忌算法，针对AB非晶格模型提出了一种混合优化算法。通过在AB非晶格模型上的实验，结果表明此混合算法能得到比几乎所有以上算法得到的最优结果更低的能量值。

发明内容

鉴于已有生产方法存在的缺陷，本发明提供一种基于三维AB非晶格模型的蛋白质三维结构预测算法BSA-TS混合算法，该方法是将BSA算法用于蛋白质结构三维预测且加以改进，同时在算法后期使用加入变异算子的TS算法，利用其全局搜索能力较强的优点对BSA算法运行得到的优化结果进行进一步全局搜索，有效提高了算法的搜索精度，得到更稳定的蛋白质三维结构。

该方法首先通过改进的BSA算法运行3000次得到较稳定的优化结果，然后利用改进的TS算法进行2000次全局局部邻域搜索得到最终预测结果。在蛋白质三维结构预测过程中,对使用的BSA算法所作的改进为:将BSA算法选择-Ⅱ策略替换为惩罚机制，判断种群中个体两两间欧氏距离，如果距离小于预设临界值D(经验值，与序列长度相关)时则对其中适应值较高的个体实施惩罚操作，以此增加种群中个体的多样性；对TS算法所作的改进为：引进了变异算子，使邻域产生前期对当前个体产生较大扰乱度，保证变异个体多样性；邻域产生后期对当前解产生较小扰乱度，避免错过全局最优，以提高TS算法的全局搜索能力和算法精度。

为实现上述目的，本发明所采用的技术方案是：

步骤1：设置参数，初始化种群P与历史种群P^old；

步骤2：进入改进的BSA算法迭代循环部分，判断终止条件是否满足，满足则执行步骤4,不满足则执行步骤3；

步骤3：对种群依次执行选择-Ⅰ策略、执行变异操作、执行交叉策略、执行边界控制、执行惩罚机制；

步骤4：对BSA部分得到的优化结果用改进的TS算法，通过引入变异算子进行全局邻域搜索；当满足终止条件时，退出循环，输出优化结果。

所述的方法包括以下具体步骤:

1)随机初始化种群P和历史种群P^old；

2)所述步骤2中终止条件为判断循环代数epk值是否小于3000，若大于3000，则满足终止条件，执行9)，否则执行3)；

3)执行选择-I策略：在每次迭代初始时重新确定将被用来计算搜索方向的历史种群P^old：

Ifa＜bthenP^old＝P

其中，a，b为随机数，a,b～U(0,1)，从之前迭代得到的种群中选择一个作为历史种群P^old，并记住它直到再次发生变化；

4)在P^old确定之后，对其中个体进行随机排序，并重新赋值给P^old：

P^old＝P^old(randperm(popsize),:)

5)执行变异操作：产生实验种群Mutant的初始值：

Mutant＝P+F*(P^old-P)

其中，F是变异尺度，控制搜索方向矩阵(P^old-P)的变异幅度，F取值为2.5*randn，其中randn是一个符合均值为0，方差为1的标准正态分布的随机数。

6)执行交叉策略：通过判断两个0到1之间的随机数大小而随机调用其两种交叉策略；分别进行两种交叉时，首先确定交叉长度u以及对应交叉位置，交叉位置是先判断交叉长度u，交叉位置在程序中是通过随机函数确定的u个交叉位置；然后将种群P与Mutant中相同位置个体的同维元素互换(此处相同位置指的是两个200个个体的种群中都选择第i个个体，同维指的是每个个体中的dim个维度中都选择第j个元素，以此形成相同大小矩阵的对应元素互换)，生成新种群P’与Mutant’，并将P’赋给Mutant。在每次交叉时，u的两种选择方式如下：

a)交叉策略I：u(i)＝1；

b)交叉策略II：其中mixrate为交叉概率，dim为每个染色体长度。

7)执行边界控制：

在执行交叉策略后，对种群Mutant个体中超出搜索范围的基因执行边界控制，通过判断两个0到1之间的随机数大小而等概率调用以下两种边界控制方法产生新基因，新生成的实验种群记为Q：

a)当种群Mutant中某基因P_ij超过基因搜索范围的上界，则重新赋值为上界up_j,此处上界为角度大小的上界π，也是问题定义域的上界。

低于基因搜索范围的下界时，则重新赋值为下界low_j，此处为-π；

b)当种群Mutant中某基因P_ij超过基因搜索范围的上界或下界时，P_ij＝rand*(up(j)-low(j))+low(j)；

8)执行惩罚机制：将种群P按适应度大小升序排列，取其前一半个体为种群R，将R与Q组成新种群Mu，求出Mu所有个体之间的欧氏距离d_ij，如果d_ij<D时，则将序列Mu_i与Mu_j之间适应度较大的值处以罚函数Penalty，Penalty为10¹³，然后按适应度值将Mu按升序排序，取前popsize个个体赋给种群P并参与下次迭代。

9)设定禁忌算法参数，将BSA部分得到的globalminimizer值作为初始解，禁忌表置空。

globalminimizer值指第8)步最后的种群P中的第一行(即BSA部分得到的适应度全局最小值对应的个体)。

10)所述步骤4中终止条件是判断循环代数T是否小于2000，若大于等于2000，则满足终止条件，执行15)，否则执行11)；

11)由当前解结合变异算子产生邻域解，并确定候选解；

12)判断是否满足藐视准则，判断候选解与目前得到的最优解对应的最低适应度函数值pbest的大小，如果某候选解适应度值比pbest低，则满足藐视准则)，若满足，则执行14)，否则执行13)；

13)判断候选解禁忌属性，将非禁忌对象对应的最优状态作为当前解，并用该对象替换最早进入禁忌表的对象。跳转10)。

14)将满足藐视准则的解作为当前解，其对应禁忌对象替换最早进入禁忌表的对象，更新最优状态。跳转10)。

15)输出优化结果。

将BSA算法引入到蛋白质结构预测中并加以改进，利用BSA算法对控制参数不过于敏感的优点可以降低优化结果的随机性，得到较稳定的适应度值。

所述变异算子为产生邻域的函数，对当前解随机选一个基因x用变异算子公式进行变异为x_new：

其中x^k为随机选择的粒子位置，c和r为0到1之间的随机数，rate是因子规模，本文取0.95，f(r)为相关系数。

所述适应度大小是能量函数值E，能量函数值E的计算公式：

其中r_ij为任意不相邻残基之间欧式距离，ξ_i,ξ_j为残基的二进制编码，如果第i个残基为疏水性残基A，则ξ_i＝1，如果为B，则ξ_i＝-1，C(ξ_i,ξ_j)公式在一定程度上反应真实蛋白质特性：AA残基对间有较强引力，AB残基对间有较弱斥力，BB残疾对间有较弱引力，求解适应度时需要用到的将每个个体的角度值转换为三维空间坐标的公式如下：

采用上述方法后，与现有技术方法相比具有以下有益效果：将BSA算法引入到蛋白质三维结构预测中并加以改进，利用其算法框架简单高效，控制参数少，且优化结果对控制参数的初值不过于敏感的优点可以得到较稳定的优化结果；在算法后期TS算法部分，使用变异算子的思想产生当前解的邻域，使邻域产生前期对当前个体产生较大扰乱度，保证变异个体多样性；邻域产生后期对当前解产生较小扰乱度，避免错过全局最优，以提高TS算法的全局搜索能力和算法精度。

附图表说明

图1三维AB非晶格模型结构示意图；

图2程序流程图；

图3实验用Fibonacci序列；

图4实验用真实蛋白质序列；

图5本发明与其他方法在Fibonacci序列上实验结果比较；

图6本发明与其他方法在真实蛋白质序列上实验结果比较。

具体实施方式

本发明的实施例是在以本发明技术方案为前提下进行实施的，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述实施例。

因为蛋白质序列种类很多，我们以图3中长度为13的Fibonacci序列为例。

步骤1：将蛋白质序列分为疏水性氨基酸(H)与亲水性氨基酸(P)，并且将蛋白质序列抽象为数字序列，以便计算机运算，其中疏水氨基酸对应数字1，亲水性氨基酸对应数字-1。因此长度为13的序列表示为[1-1-11-1-11-11-1-11-1]。

步骤2:根据序列长度随机化种群个体：序列长度为n时，dim(问题维度)为2n-5(n-2个相邻键之间的键角，n-3个平面扭转角)，因此n为13时，问题维度dim为21。因此可以初始化种群

P＝[P₁,P₂,...P_i...,P_popsize],popsize为种群大小，popsize取200，个体P_i维度dim为21，dim即每个染色体长度，是指种群P中每个个体的维度都为dim，每个个体代表基因算法中的一条染色体，其中每个个体的每个元素都表示随机的初始角度，范围为[-π，π]，并判断循环代数epk是否小于3000，若满足小于3000，则执行步骤3，否则执行步骤8。

步骤3：执行选择-Ⅰ，确定本次迭代中被用来计算搜索方向的历史种群P^old。

步骤4：对P^old中个体的顺序进行随机洗牌操作，并重新赋值给P^old：

步骤5：执行变异操作产生实验种群Mutant的初始值：

步骤6：通过判断两个0到1之间的随机数大小而随机调用上文介绍的两种交叉方式。交叉时，首先确定交叉长度L和交叉位置，然后将父代种群P与Mutant中相同顺序个体的相同维数元素进行互换，生成新种群P’与Mutant’，并将P’赋给Mutant。

交叉结束后通过判断两个0到1之间的随机数大小而等随机调用上文介绍的两种边界控制方法产生新元素，新生成的实验种群作为实验种群的最终形式，记为Q。

步骤7：实施惩罚操作，此处粒子边界距离D取1(经验值)，此处所用欧式距离公式：

步骤8：设定禁忌算法参数，将BSA部分得到的globalminimizer值作为TS算法初始解，禁忌表置空。

步骤9：判断循环代数T是否小于2000，若满足小于2000，则执行步骤10，否则跳转步骤14。

步骤:10：由当前解结合邻域产生函数产生800个邻域解，其中用于邻域产生的变异算子公式如下，并确定适应度最低的十个值为候选解：

步骤11：判断适应度最低的候选解是否满足藐视准则，若满足，则执行步骤13，否则执行步骤12。

步骤12：判断候选解禁忌属性，将非禁忌对象中适应度最低的个体作为当前解，并用该个体的适应度值替换最早进入禁忌表的对象。跳转步骤9。

步骤13：将满足藐视准则的解作为当前解并更新最优状态，其对应适应度值替换最早进入禁忌表的对象。跳转步骤9。

步骤14：输出优化结果。

根据以上所示步骤，我们在Fibonacci序列(见图3)和真实蛋白质序列(见图4)上进行相关操作，对所获得的数据与其他方法进行了比较(见图5和图6)。

综上所述，我们通过BSA-TS混合算法进行蛋白质三维结构预测的方法不但可以提高算法精度，得到更稳定的蛋白质空间结构，而且得到的结果随机性比较小，同时算法收敛速度也有一定提高，由此可以说明我们方法具有良好的性能。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明披露的技术范围内，根据本发明的技术方案及其发明构思加以等同替换或改变，都应涵盖在本发明的保护范围内。

Claims (5)

1.基于BSA-TS算法的蛋白质三维结构预测方法，其特征在于，包括以下步骤：

步骤1：设置参数，初始化种群P与历史种群P^old；

步骤2：进入改进的BSA算法迭代循环部分，判断终止条件是否满足，满足则执行步骤4,不满足则执行步骤3；

步骤3：对种群依次执行选择-Ⅰ策略、执行变异操作、执行交叉策略、执行边界控制、执行惩罚机制；

步骤4：对BSA部分得到的优化结果用改进的TS算法，通过引入变异算子进行全局邻域搜索；当满足终止条件时，退出循环，输出优化结果。

2.根据权利要求1所述的基于BSA-TS算法的蛋白质三维结构预测方法，其特征在于，所述的方法包括以下具体步骤：

1)随机初始化种群P和历史种群P^old；

2)所述步骤2中终止条件为判断循环代数epk值是否小于3000，若大于3000，则满足终止条件，执行9)，否则执行3)；

3)执行选择-I策略：在每次迭代初始时重新确定将被用来计算搜索方向的历史种群P^old：

Ifa＜bthenP^old＝P

其中，a，b为随机数，a,b～U(0,1)，从之前迭代得到的种群中选择一个作为历史种群P^old，并记住它直到再次发生变化；

4)在P^old确定之后，对其中个体进行随机排序，并重新赋值给P^old：

P^old＝P^old(randperm(popsize),:)

5)执行变异操作：产生实验种群Mutant的初始值：

Mutant＝P+F*(P^old-P)

其中，F是变异尺度，控制搜索方向矩阵(P^old-P)的变异幅度，F取值为2.5*randn，其中randn是一个符合均值为0，方差为1的标准正态分布的随机数；

6)执行交叉策略：通过判断两个0到1之间的随机数大小而随机调用其两种交叉策略；分别进行两种交叉时，首先确定交叉长度u以及对应交叉位置，然后将种群P与Mutant中相同位置个体的同维元素互换，生成新种群P’与Mutant’，并将P’赋给Mutant；在每次交叉时，u的两种选择方式如下：

a)交叉策略I：u(i)＝1；

b)交叉策略II：其中mixrate为交叉概率，dim为每个染色体长度；

7)执行边界控制：

在执行交叉策略后，对种群Mutant个体中超出搜索范围的基因执行边界控制，通过判断两个0到1之间的随机数大小而等概率调用以下两种边界控制方法产生新基因，新生成的实验种群记为Q：

a)当种群Mutant中某基因P_ij超过基因搜索范围的上界，则重新赋值为up_j,低于基因搜索范围的下界时，则重新赋值为low_j；

b)当种群Mutant中某基因P_ij超过基因搜索范围的上界或下界时，P_ij＝rand*(up(j)-low(j))+low(j)；

8)执行惩罚机制：将种群P按适应度大小升序排列，取其前一半个体为种群R，将R与Q组成新种群Mu，求出Mu所有个体之间的欧氏距离d_ij，如果d_ij<D时，则将序列Mu_i与Mu_j之间适应度较大的值处以罚函数Penalty，Penalty为10¹³，然后按适应度值将Mu按升序排序，取前popsize个个体赋给种群P并参与下次迭代；

9)设定禁忌算法参数，将BSA部分得到的globalminimizer值作为初始解，禁忌表置空；

10)所述步骤4中终止条件是判断循环代数T是否小于2000，若大于等于2000，则满足终止条件，执行15)，否则执行11)；

11)由当前解结合变异算子产生邻域解，并确定候选解；

12)判断是否满足藐视准则，判断候选解与目前得到的最优解对应的最低适应度函数值pbest的大小，如果某候选解适应度值比pbest低，则满足藐视准则)，若满足，则执行14)，否则执行13)；

13)判断候选解禁忌属性，将非禁忌对象对应的最优状态作为当前解，并用该对象替换最早进入禁忌表的对象，跳转10)；

14)将满足藐视准则的解作为当前解，其对应禁忌对象替换最早进入禁忌表的对象，更新最优状态，跳转10)；

15)输出优化结果。

3.根据权利要求1所述的基于BSA-TS算法的蛋白质三维结构预测方法，其特征在于：将BSA算法引入到蛋白质结构预测中并加以改进，利用BSA算法对控制参数不过于敏感的优点可以降低优化结果的随机性，得到较稳定的适应度值。

4.根据权利要求1所述的基于BSA-TS算法的蛋白质三维结构预测方法，其特征在于：所述变异算子为产生邻域的函数，对当前解随机选一个基因x用变异算子公式进行变异为x_new：

$x_{new}^k=x^k+2*\mathrm{\&pi;}*f\left(r\right)*c*{\mathrm{rate}}^i$

$f\left(r\right)=\left\{\begin{array}{c}1,r\&GreaterEqual;0.5\\ -1,r<0.5\end{array}\right.$

其中x^k为随机选择的粒子位置，c和r为0到1之间的随机数，rate是因子规模，f(r)为相关系数。

5.根据权利要求1所述的基于BSA-TS算法的蛋白质三维结构预测方法，其特征在于：所述适应度大小是能量函数值E，能量函数值E的计算公式：

$E=\underset{i=1}{\overset{n-2}{\&Sigma;}}\frac{1}{4}(1-{\mathrm{cos\&theta;}}_i)+\underset{i=1}{\overset{n-2}{\&Sigma;}}\underset{j=i+2}{\overset{n}{\&Sigma;}}4\&lsqb;r_{ij}^{-12}-C({\mathrm{\&xi;}}_i,{\mathrm{\&xi;}}_j)r_{ij}^{-6}\&rsqb;$

$C({\mathrm{\&xi;}}_i,{\mathrm{\&xi;}}_j)=\left\{\begin{array}{cc}+1& {\mathrm{\&xi;}}_i=1,{\mathrm{\&xi;}}_j=1\\ +0.5& {\mathrm{\&xi;}}_i=-1,{\mathrm{\&xi;}}_j=-1\\ -0.5& {\mathrm{\&xi;}}_i\&NotEqual;{\mathrm{\&xi;}}_j\end{array}\right.$

$r_{ij}=\left|\right|x_i-x_j\left|\right|=\sqrt{\underset{k=1}{\overset{\dim }{\&Sigma;}}{(x_{ik}-x_{jk})}^2}$

其中r_ij为任意不相邻残基之间欧式距离，ξ_i,ξ_j为残基的二进制编码，C(ξ_i,ξ_j)公式反应真实蛋白质特性：如果第i个残基为疏水性残基A，则ξ_i＝1，如果为B，则ξ_i＝-1。