基于蚁群算法的蛋白质折叠优化方法

专利名称：基于蚁群算法的蛋白质折叠优化方法
技术领域：
本发明涉及生物信息技术和智能计算两大领域，主要涉及一种运用蚁群算法优化蛋白质折叠的方法。
技术背景 蛋白质的功能和它们的结构是紧密相连的，在蛋白质分析中最重要的一步就是找出蛋白质的结构。在自然界中，水溶液里的蛋白质会在一个非常短的时间内自然地折叠到其天然的结构。虽然现有的实验技术，例如X射线晶体学、核磁共振、蛋白质工程、质谱分析、荧光共振能量转移和原子力显微镜，可以用于生成蛋白质的影像，但是它们普遍的缺点是十分昂贵且耗时。更重要的是，每一种实验技术都有其不同的应用限制。由于很多简单的蛋白质序列只有唯一的天然结构，这使得运用计算方法预测蛋白质结构成为可能。计算方法可以在较大的范围内应用，而且预测的精确性比实验方法要高。这种基于序列预测蛋白质结构的问题，称为蛋白质折叠问题。
蛋白质折叠的一个公认的假设是蛋白质序列会折叠到自由能最小的平衡结构。于是，解决蛋白质折叠问题就是要寻找满足自由能最小状态的蛋白质折叠结构，其中亲疏水(HP)格点模型是一种非常重要的计算模型。它将蛋白质序列放置在一个方格阵中，其中每一个氨基酸无重叠地占据一格。在这种模型中，蛋白质内的氨基酸被分成两大类，一类是疏水性的氨基酸(H)，另一类是亲水性的氨基酸(P)。两个在蛋白质序列中非相连而在格点中邻接的疏水性氨基酸将带有能量值-1，这种邻接方式称为H-H键。寻找这种HP格点模型的最小自由能状态，就是要找出带有最多H-H键的蛋白质结构。
除了一些特殊的蛋白质序列外，蛋白质折叠问题是一个NP完全问题，也就是说无法在多项式的计算时间内得到问题的解。迄今为止，有很多基于HP格点模型的计算方法被提出来，例如遗传算法、蒙特卡罗算法、免疫算法和蚁群算法等。蚁群算法最早由Dorigo提出，其后发展出精华蚂蚁系统、最大最小蚂蚁系统和蚁群系统等一系列算法，它们统称为蚁群算法。从发明以来，蚁群算法已被成功运用于求解各种应用问题，例如车辆调度问题、车间调度问题和网络路由问题等。蚂蚁释放信息素和构造启发式信息的特点也在蛋白质折叠问题上得到了成功运用。


发明内容
本发明提出一种有效的蚁群算法用于求解蛋白质折叠问题。发明的蚁群算的执行步骤如下 (1)读入蛋白质序列并初始化各个参数。其中蛋白质序列以亲水性氨基酸(P)和疏水性氨基酸(H)的序列表示。
(2)每一只蚂蚁都是在网格中央从蛋白质序列的中间开始构造路径的。对于带有n个氨基酸的蛋白质序列{s0，s1，...，sn-1}(sj∈{P，H}，j＝0，...，n-1)，每只蚂蚁都是从(n+2)×(n+2)大小的网格模板的中央的两个水平方格处开始路径的构造的。位于网格模板中央靠左的水平方格称为左起始方格，靠右的水平方格称为右起始方格。蚂蚁开始构造路径时，首先将氨基酸sn/2放置在左起始方格中，氨基酸sn/2+1放置在右起始方格中。蛋白质的子序列{s0，...，sn/2}从左起始方格开始构建，称为“左路径”，而另一个子序列{sn/2+1，...，sn-1}则从右起始方格开始构建，称为“右路径”。蚂蚁随机选择从蛋白质序列的左端或者右端开始走下一步。
(3)所有蚂蚁开始根据蛋白质序列建立可行的折叠结构。如果蚂蚁当前处于方格i，下一个要放置的氨基酸为sj，那么选择下一个方向d的概率如下
其中，τid和ηjd分别为信息素和启发式信息的取值，β为启发式信息的加强参数。蚂蚁每移动一步都要执行局部信息素的更新，降低刚刚经过的两个相邻方格之间的信息素 τid←δ·τid 其中i是蚂蚁当前所在方格的序号，d是蚂蚁选择放置下一个氨基酸的方格所在的方向。δ＝(m-1)/m＜1称为局部蒸发率，m是蚂蚁的总数。设置信息素的最低取值为τmin，当信息素低于这个值时将会重置。
(4)统计每只蚂蚁构建出来的蛋白质在网格模板中形成的处于相邻方格但不在蛋白质序列相邻的H-H键数量，并选出本次迭代中的最优蚂蚁。
(5)执行全局信息素的更新。首先蒸发所有相邻方格之间的信息素 τid←ρ·τid 其中，ρ是全局信息素蒸发率，i＝0，1，2，...，(n+2)2-1。然后，加强本次循环中找到的最好的折叠路径上的信息素 其中，i′∈{最优蚂蚁所经过的方格}，d为蚂蚁从方格i′移向下一个方格的方向，ε是在本次循环中最大的H-H键数量，是在HP格点模型中近似的最小自由能。
(6)如果达到结束条件就终止整个算法，否则返回步骤(2)。
提出的方法带有三种特别的机制，包括路径的修补、信息素的吸引和折叠的启发式信息。在运用HP格点模型解决蛋白质折叠问题上，这些新的机制使得提出的方法能够更快速有效地求解蛋白质折叠问题。
路径的修补HP格点模型中的蛋白质序列是一个由疏水氨基酸(H)和亲水氨基酸(P)组成的串。蚂蚁在网格中一个接一个地放置氨基酸。如果氨基酸周围的方格都已经被占用，那么就无法防止下一个氨基酸了，这种情况被称为停滞。这时候运用路径修补策略，对发生停滞的路径进行修补。
信息素在发明的方法中，信息素是放置在连接相邻两个方格的虚边上的。在一个方格中，这样的虚边实质上就是指向与其相邻的方格的四个方向，分为上(U)、下(D)、左(L)和右(R)。这种信息素放置方法与之前提出的用于求解蛋白质折叠问题的蚁群算法是不同的，其好处在于信息素被固定在网格中，因此蚂蚁可以从网格中直接感知蛋白质折叠的全局状况。
启发式信息在发明的方法中，蚂蚁在信息素和启发式信息的共同作用下构造蛋白质的结构。疏水性的氨基酸的启发式信息定义为 ηjd＝hjd+1 其中，hjd是把氨基酸sj放置相邻的方格中将会新产生的H-H键的数量，d代表蚂蚁移动的方向。
亲水性的氨基酸的启发式信息为，下一个可能被选择的方格周围的空格和亲水性氨基酸的总和加一，也就是ηjd＝vjd+hjd′+1 其中，vjd和hjd′分别为sj下一个可能位置周围的空格和亲水性氨基酸的数量。两种氨基酸的启发式信息是不同的，这也是提出的方法与之前的算法的一个区别。



图1蛋白质序列的二维HP折叠结构示意图 图2蚂蚁在网格中构造路径的方法示意图 图3蚂蚁选择下一步移动的方法示意图 图4路径停滞状态示意图 图5蚁群算法优化蛋白质折叠问题的整体流程图
具体实施例方式 以下结合附图对发明的方法作进一步的描述。
在自然界中一共存在二十种氨基酸，用不同的分类方法，它们可以分为酸性、碱性或中性；带正电、带负电或不带电；疏水的或亲水的等等。对于在水溶液中的球状蛋白，它的亲水性氨基酸由于受到水分子的吸引会倾向于分布在球的表面。而疏水性氨基酸由于受到水分子的排斥，除了一些特殊的分布在表面的疏水氨基酸区域外，其余大部分的疏水性氨基酸将会聚集在球状蛋白的内部形成一个核心。
含有48个氨基酸的蛋白质序列的一些二维HP折叠结构如图1所示，其中实心正方形表示疏水性氨基酸，空心正方形表示亲水性氨基酸。对于格点HP模型，蛋白质结构的优劣是根据H-H键(图1左图中的虚线)的数量判断的。在图1中，每一种结构的H-H键的数量都是23，因此它们的最小自由能状态均为E*＝-23。从图中还可以看到疏水性氨基酸在蛋白质结构中形成内核，而亲水性氨基酸则分布在表面。
给定二维的格点模板，求解蛋白质折叠问题就是将蛋白质序列放置在网格中形成一条自回避的路径。蚁群算法的目标就是要寻找一条含有最多H-H键的路径。
为了不超出网格的边界，每一只蚂蚁都是在网格中央从蛋白质序列的中间开始构造路径的。对于带有n个氨基酸的蛋白质序列{s0，s1，...，sn-1}(sj∈{P，H}，j＝0，...，n-1)，每只蚂蚁都是从(n+2)×(n+2)大小的网格模板的中央的两个水平方格处开始路径的构造的，如图2所示。模板中的方格从左上到右下依次编号为0到(n+2)2-1。于是，左右两个起始方格的编号分别为(n/2+1)·(n+2)+n/2和(n/2+1)·(n+2)+n/2+1。这种对网格中的方格进行编号的方法有两个好处。其一是方格的坐标可以用一维的数字来表示。其二是可以很容易地得到四个相邻的方格。如图3所示，蚂蚁在方格i中，它可以向左、右、上或者下移动，相对应的方格的编号为i-1、i+1、i-(n+2)、或者i+(n+2)。由于蚂蚁已经经过了右边的方格，所以它只能选择向剩下的三个方格移动。
蚂蚁开始构造路径时，首先将氨基酸sn/2放置在左起始方格中，氨基酸sn/2+1放置在右起始方格中。蛋白质的子序列{s0，...，sn/2}从左起始方格开始构建，称为“左路径”，而另一个子序列{sn/2+1，...，sn-1}则从右起始方格开始构建，称为“右路径”。蚂蚁随机选择从蛋白质序列的左端或者右端开始走下一步。如果蚂蚁当前处于方格i，下一个要放置的氨基酸为sj，那么选择下一个方向d的概率如下
其中，τid和ηjd分别为信息素和启发式信息的取值，β为启发式信息的加强参数。
在发明的方法中，信息素被放置在连接两个相邻方格的虚边上，并以变量τid表示，其中i＝0，1，2，...，(n+2)2-1，d＝L，R，U，D。
在启发式信息的设置方面，不同于之前提出的算法只考虑疏水性氨基酸的启发式信息，本文提出的方法不但考虑疏水性氨基酸，而且还考虑亲水性氨基酸的启发式信息。由于蛋白质折叠问题的目标是寻找最小自由能的结构，因此疏水性氨基酸(H)的启发式信息由H-H键的数量决定的。如果一种结构能够生成更多的H-H键，那么它就应该有更高的概率被选择。如果蚂蚁要放置的下一个氨基酸sj是疏水性的，那么启发式信息定义为 ηjd＝hjd+1 其中，hjd是把氨基酸sj放置在相邻的方格中将会新产生的H-H键的数量，d代表蚂蚁移动的方向。
如果下一个要放置的氨基酸是亲水性的，则启发式信息是下一个可能被选择的方格周围的空格和亲水性氨基酸的总和加一 ηjd＝vjd+hjd′+1 其中，vjd和hjd′分别为sj下一个可能位置周围的空格和亲水性氨基酸的数量。值得注意的是，该方法不考虑亲水性氨基酸在蛋白质序列中的连续性。
为了避免早熟收敛，蚂蚁每移动一步都要执行局部信息素的更新，降低刚刚经过的两个相邻方格之间的信息素τid←δ·τid 其中i是蚂蚁当前所在方格的序号，d是蚂蚁选择放置下一个氨基酸的方格所在的方向。δ＝(m-1)/m＜1称为局部蒸发率，m是蚂蚁的总数。设置信息素的最低取值为τmin，当信息素低于这个值时将会重置。
当蚂蚁在放置氨基酸sj的时候，如果它已经访问过所有相邻的方格，那么蛋白质就不能进行折叠了。在这种情况下，就需要对折叠路径进行修补。假设已经构造了子序列{sleft，...，sstartL，sstartR，...，sright}，其中left代表最左端的氨基酸的序号，right代表最右边的氨基酸的序号，startL＝n/2，startR＝n/2+1。对右路径进行修补的方法为，根据下式随机生成一个数jj＝rand％(right-startR)+startR 其中rand是一个随机的非负整数。释放序号从sj+1到sright的氨基酸，并且清空相对应的方格。同样的，对左路径进行修补方法也是首先生成一个随机数jj＝rand％(startL-left)+left+1 释放从sleft到sj-1的氨基酸并清空相对应的方格。
需要注意的是，即使停滞只发生在蛋白质的右路径，也并不意味着只有右路径需要修补，因为一些停滞状况并不能简单地通过相应边的路径修补得到清除。图4举例说明了在右路径上的两种停滞状态。其中，两个空心圆分别代表左右两个起始氨基酸，实心正方形代表在左路径上的氨基酸，空心三角形代表在右路径上的氨基酸。在图4的左图中，当j＝21的时候停滞状态可以通过对右路径修补得到解决，之后蚂蚁可以选择往上走。但是在图4的右图中，停滞却不能通过对右路径修补得到解决，而必须对左路径进行修补。因此，运用两个布尔变量RightRetrievalBool和LeftRetrievalBool来判断，以保证同一个方向上不会进行两次连续的修补。
把对右路径上的停滞进行修补的程序称为“右边修补(RightSideRetrieval)”，对左路径上的停滞进行修补的程序称为“左边修补(LeftSideRetrieval)”。分别定义两个函数“RightRetrievalSequence()”和“LeftRetrievalSequence()”执行左右的修补步骤，程序的伪代码如下 /*startL＝n/2，startR＝n/2+1 left已构造的最左端氨基酸的序号，right已构造的最右端氨基酸的序号*/ Procedure RightSideRetrieval({sleft，...，sstartL，sstartR，...，sright}) If startR＜right && RightRetrievalBool＝＝falsej＝rand％(right-startR)+startR；RightRetrievalSequence({sleft，...，sstartL，sstartR，...，sright}，j)；LeftRetrievalBool＝false；RightRetrievalBool＝true； Else If startL＞leftj＝rand％(startL-left)+left+1；LeftRetrievalSequence({sleft，...，sstartL，sstartR，...，sright}，j)；RightRetrievalBool＝false； End Procedure LeftSideRetrieval({sleft，...，sstartL，sstartR，...，sright}) If startL＞left && LeftRetrievalBool＝＝falsej＝rand％(startL-left)+left+1；LeftRetrievalSequence({sleft，...，sstartL，sstartR，...，sright}，j)；LeftRetrievalBool＝true；RightRetrievalBool＝false； Else If startR＜rightj＝rand％(right-startR)+startR；RightRetrievalSequence({sleft，...，sstartL，sstartR，...，sright}，j)；LeftRetrievalBool＝false； End 以图4中的右图为例，由于停滞发生在右路径，因此进行右边修补。根据startR＝21，right＝24和RightRetrievalBool＝false，可以生成一个随机数j。假设j＝22，然后执行函数“RightRetrievalSequence()”释放编号为23至24的氨基酸，并将RightRetrievalBool设置为true，LeftRetrievalBool设置为false。根据之前得到的结论，这并不能解决停滞问题。随着路径构建的继续进行，将再次发生停滞。这时，再次执行右边修补程序。由于RightRetrievalBool的取值为true，因此只能执行“LeftRetrievalSequence()”函数去释放左路径上的氨基酸，此时如果生成的随机数j在5和20之间，停滞问题就可以得到解决。
当所有蚂蚁都完成蛋白质折叠路径的构造后，执行全局信息素的更新，所有相邻方格之间的信息素都会被蒸发 τid←ρ·τid 其中，ρ是全局信息素蒸发率，i＝0，1，2，...，(n+2)2-1，d＝L，R，U，D。
然后，加强本次循环中找到的最好的折叠路径上的信息素 其中，i′∈{最优蚂蚁所经过的方格}，d为蚂蚁从方格i′移向下一个方格的方向，ε是在本次循环中最大的H-H键数量，是在HP格点模型中近似的最小自由能。
蚁群算法优化蛋白质折叠问题的完整的流程图如图5所示。
一些在二维HP格点模型中的基准蛋白质序列如下表所示。其中，l表示氨基酸的数量而E*表示最小自由能水平。字母“H”代表疏水性氨基酸，字母“P”代表亲水性氨基酸。
以这八种蛋白质的折叠问题为例，对本文提出的优化算法进行测试。其中，参数的设置为，τ0＝1/3，τmin＝0.05，ρ＝0.9，对于蛋白质序列1至7，m＝10，β＝2，对于蛋白质序列8，m＝100，β＝3。将发明的方法求得的结果与遗传算法算法相比，提出的算法对所有八个蛋白质折叠问题求解速度都更快。与免疫算法相比，发明的蚁群算法在所有问题中都能以100％的概率得到最优的结构，而免疫算法求解蛋白质序列8的成功率仅为56.67％。综上可得，发明的方法在蛋白质折叠问题的求解中是高效并且快速的。
权利要求
1、一种基于蚁群算法的蛋白质折叠优化方法，其特征在于，该方法包括以下步骤
(1)读入蛋白质序列并初始化各个参数。其中蛋白质序列以亲水性氨基酸(P)和疏水性氨基酸(H)的序列表示。
(2)每一只蚂蚁都是在网格中央从蛋白质序列的中间开始构造路径的。对于带有n个氨基酸的蛋白质序列{s0，s1，...，sn-1}(sj∈{P，H}，j＝0，...，n-1)，每只蚂蚁都是从(n+2)×(n+2)大小的网格模板的中央的两个水平方格处开始路径的构造的。位于网格模板中央靠左的水平方格称为左起始方格，靠右的水平方格称为右起始方格。蚂蚁开始构造路径时，首先将氨基酸sn/2放置在左起始方格中，氨基酸sn/2+1放置在右起始方格中。蛋白质的子序列{s0，...，sn/2}从左起始方格开始构建，称为“左路径”，而另一个子序列{sn/2+1，...，sn-1}则从右起始方格开始构建，称为“右路径”。蚂蚁随机选择从蛋白质序列的左端或者右端开始走下一步。
(3)所有蚂蚁开始根据蛋白质序列建立可行的折叠结构。如果蚂蚁当前处于方格i，下一个要放置的氨基酸为sj，那么选择下一个方向d的概率如下
其中，τid和ηjd分别为信息素和启发式信息的取值，β为启发式信息的加强参数。蚂蚁每移动一步都要执行局部信息素的更新，降低刚刚经过的两个相邻方格之间的信息素
τid←δ·τid
其中i是蚂蚁当前所在方格的序号，d是蚂蚁选择放置下一个氨基酸的方格所在的方向。δ＝(m-1)/m＜1称为局部蒸发率，m是蚂蚁的总数。设置信息素的最低取值为τmin，当信息素低于这个值时将会重置。
(4)统计每只蚂蚁构建出来的蛋白质在网格模板中形成的处于相邻方格但不在蛋白质序列相邻的H-H键数量，并选出本次迭代中的最优蚂蚁。
(5)执行全局信息素的更新。首先蒸发所有相邻方格之间的信息素
τid←ρ·τid
其中，ρ是全局信息素蒸发率，i＝0，1，2，...，(n+2)2-1。然后，加强本次循环中找到的最好的折叠路径上的信息素
其中，i′∈{最优蚂蚁所经过的方格}，d为蚂蚁从方格i′移向下一个方格的方向，ε是在本次循环中最大的H-H键数量，是在HP格点模型中近似的最小自由能。
(6)如果达到结束条件就终止整个算法，否则返回步骤(2)。
2、基于权利要求1所述的一种基于蚁群算法的蛋白质折叠优化方法，其特征在于，蚂蚁在网格中逐个放置氨基酸时，如果氨基酸周围的方格都已经被占用，则对相应路径执行修补程序。
3、基于权利要求1所述的一种基于蚁群算法的蛋白质折叠优化方法，其特征在于，信息素是放置在连接相邻两个方格的虚边上的。在一个方格中，这样的虚边实质上就是指向与其相邻的方格的四个方向，分为上(U)、下(D)、左(L)和右(R)。
4、基于权利要求1所述的一种基于蚁群算法的蛋白质折叠优化方法，其特征在于，疏水性和亲水性两种氨基酸的启发式信息是不同的。疏水性的氨基酸的启发式信息定义为
ηjd＝hjd+1
其中，hjd是把氨基酸sj放置在相邻的方格中将会新产生的H-H键的数量，d代表蚂蚁移动的方向。亲水性的氨基酸的启发式信息为，下一个可能被选择的方格周围的空格和亲水性氨基酸的总和加一，也就是
ηjd＝vjd+hjd′+1
其中，vjd和hjd′分别为sj下一个可能位置周围的空格和亲水性氨基酸的数量。
全文摘要
本发明公开了一种基于蚁群算法的蛋白质折叠优化方法。运用亲疏水格点模型，首先将蚂蚁放置于左右两个起始方格中。然后让蚂蚁根据信息素和启发式信息将蛋白质中的氨基酸放置在网格中，蚂蚁每走一步执行一次局部信息素的更新，如果发生了停滞状况，则对相应的路径进行修补。最后，当所有蚂蚁都完成蛋白质折叠路径的构造后，执行全局信息素的更新。不断重复这些步骤以得到自由能最小的蛋白质结构。提出的方法与用于求解蛋白质折叠问题的已有算法相比，具有高效、快速等优点。
文档编号G06N3/00GK101236616SQ20081002647
公开日2008年8月6日 申请日期2008年2月27日 优先权日2008年2月27日
发明者军 张, 胡晓敏, 韬 黄 申请人:中山大学