【摘要】 根据氨基酸的序列预测蛋白质的空间结构在基因治疗药物分子设计等方面有巨大的潜在应用价值。本研究基于hp格子模型利用改进的遗传算法预测了蛋白质的三维空间结构。改进的遗传算法引入了克隆体数量限制策略、巢穴竞争选择策略及局部优化策略等。实验结果表明，改进的遗传算法显著地提高了蛋白质结构的预测效率，模拟的蛋白质结构紧凑，更接近真实蛋白质的构型。
【关键词】 遗传算法 蛋白质折叠 三维hp模型
1 引言
蛋白质是生命活动的重要承担者，蛋白质所具有的功能在很大程度上取决于其空间结构，掌握蛋白质的空间结构在基因治疗和药物分子设计方面有极大的潜在应用价值[1]。目前,测定蛋白质空间结构的方法主要是核磁共振和x射线衍射技术，这些技术消耗巨大，测定效率低下，远远满足不了日益增加的待测定海量蛋白质的需要[2]。根据氨基酸的序列,从理论上预测蛋白质的空间结构有助于提高测定蛋白质结构的效率，对生物医学的发展有重要的意义[3]。
蛋白质结构预测是个典型的“np问题”（算法的复杂性随着规模的增长成指数增长），也就是蛋白质的结构不能用一个多项式来明确表示，其能量的最小值必须通过启发式算法来搜索[4]。目前发展的启发式算法主要有蒙特卡罗模拟算法、禁忌算法、蚁群算法、链增长算法、模拟退火算法和遗传算法等[5～8]。Www.133229.COM其中遗传算法由于高效地搜索效率得到了广泛的应用。
本研究针对hp模型（疏水性和亲水性格子模型）采用改进遗传算法模拟了蛋白质的三维空间折叠行为，改进的遗传算法主要引入了克隆体数量限制策略、巢穴竞争选择策略和局部优化策略等。
2 原理和方法
2.1 3dhp模型
最简单的蛋白质分子模型是hp格子模型，该模型把所有的氨基酸残基按疏水性和亲水性分成两类：疏水性残基（h）和亲水性残基（p）。因此，蛋白质序列被抽象为一个由h和p组成的序列[9]。hp格子模型在三维空间中的折叠简称3dhp模型，每个残基的折叠方向可以向左、向右、向上、向下90°或者向前,折叠的残基不能重叠在其它残基上，整个蛋白质序列在一个三维方格上折叠。3dhp模型的理论基础是氨基酸的疏水性是球蛋白形成的主要驱动力[2]。该模型忽略了侧链的影响，符合真实蛋白的基本特征。疏水性的氨基酸为了减小与水分子的接触面积而彼此靠近并进入分子的内部，形成了疏水相互作用；亲水性氨基酸则形成了分子的表面，形成紧密的团状构象。3dhp模型虽然过于粗糙，与真实蛋白分子相差甚远，但是它能模拟真实蛋白的折叠行为，且计算简单，有利于对比不同折叠搜索算法。
hp格子模型中，一个构象的能量计算规则如下：当两个在序列上不相邻的节点在空间上相邻时，便提供给构象一个相互作用能量.对于一个特定的序列结构，它的总能量e为：ｅ＝∑i＜n,j＜ni=1,j=i+1δreij，式中n为蛋白质序列的长度。如果i与j在空间中拓扑相邻但序列不相邻，那么δr等于1，否则等于0。eij表示在序列中第i个氨基酸与第j个氨基酸之间的能量。三维空间中拓扑相邻的残基有３种情形：hh、hp、pp，３种拓扑关系的能量规定如下[5]：ｅｈｈ＝－１．０， ｅｈｐ＝０．０， ｅｐｐ＝０．０（1）由此，蛋白质三维折叠模拟的命题表述为：搜索蛋白质序列在空间中的结构，使该结构中拓扑相邻的hh数量最多。
上述的模型得到了广泛的应用，然而这种模型只考虑hh间的相互作用，而未考虑hp间的相互作用。实际的蛋白质结构是亲水性残基包裹疏水性残基形成球状结构，忽略hp间的相互作用将导致虽然找到了最多的hh接触数量，但是末端的疏水性分子p没有任何约束而随意折叠，蛋白质空间结构的自由度太大，甚至形成与真实蛋白质结构相差太远的结构。
实际３种拓扑关系的能量大小关系为：ｅｈｈ＜ｅｈｐ＜ｅｐｐ，本研究对３种拓扑关系做如下修正：ｅｈｈ＝１．０， ｅｈｐ＝－０．４， ｅｐｐ＝０．０（2）这种修正考虑了氨基酸残基应满足的物理制约条件，不同类型的氨基酸残基趋向于分离，满足关系式[1１]：２ｅηρ＞ｅηη＋ｅｐｐ（３）本研究中，个体适应度规定为：ｆi＝－ei＋０．０１（4） 分 析 化 学第37卷第1期李绍新等：基于改进遗传算法的蛋白质三维折叠模拟 该规定保证了适应度总为正数，个体能量越低，适应度越大。增加的常量（0.01）保证了个别个体能量为零时适应度不为零，也有机会参与遗传操作。修正后的蛋白质三维折叠模拟命题表述为：寻求给定蛋白质序列具有最大适应度的三维空间结构。
2.2 遗传算法 遗传算法首先是由美国的holland教授提出来的启发式优化组合方法[1２]。它基于达尔文进化论和孟德尔遗传学说，仿效生物的进化与遗传，根据“生存竞争”和“优胜劣汰”的原则，借助复制、交换、突变等操作，使所要解决的问题从初始解一步步逼近最优解。与其他搜索方法相比，ga具有随机性、鲁棒性、并行性、全局搜索等优越性[1３]。
遗传算法运行时首先编码建立解的初始群体，编码一般采用二进制或浮点，每个解用特定的基因串表示，突变算子独立作用在串上，在最初的方案中，突变算子就是改变串上的一个位。在执行完一定数量的突变后，由交叉操作产生新的串：选择集团中的两个串，并确定串中的断点，两个新的集团成员由一个串的左边部分连接到另一个串的右边而形成。这样的操作进行到一个由可接受串组成的新的群体形成为止。接着进行下一阶段的循环。这个步骤重复进行直到集团收敛于一个串，适应值函数则用来评估突变和交叉所产生新串的质量。
2.2.1 遗传算法的编码 蛋白质折叠模拟的遗传算法个体编码与一般问题的遗传算法编码不同，其编码表示氨基酸残基的折叠方向，一般有绝对编码和相对编码，本研究采用绝对编码。一个蛋白质个体由数字0～5组成，表示氨基酸残基的折叠方向分别为东南西北上下。这种编码容易产生无效个体，也就是个体中的多个氨基酸残基会落在同一位置。但是这种编码产生的个体容易做局部的调整，个体中的氨基酸残基变异的时候对其他残基的折叠影响较小。
2.2.2 标准遗传算法 蛋白质折叠模拟的标准遗传算法如下：
（１）随机产生初始群体，计算每个个体的适应度；（２）生存选择：根据个体适应度大小选择生存个体，一般采用轮盘赌选择，适应度越大的个体被选中的概率越大；（３）交叉：采用单点或两点交叉。根据交叉概率随机选择一对交叉个体，在选中的个体上随机选择交叉位点，形成两个新个体；（４）变异：根据变异概率随机选择变异位点实施基因突变，一般采用均匀变异；（５）适应度评价：根据能量法则计算每个个体的适应度大小；（６）群体更新：如果子代个体中最优个体的适应度大于父代最优个体，则保存子代的最优个体，通过遗传操作后的所有个体代替父代个体，重复步骤2～6直到产生满足要求的最优个体。
2.2.3 改进策略 本研究对上述的标准遗传算法进行了改进，引进了新的遗传算子和计分方法，改进的策略如下：
（ⅰ）按公式（４）计算适应度，该计算方法可以保证所有个体都有机会参与遗传操作；当群体中出现无效个体（几个氨基酸残基重叠在同一位置），对该个体给予惩罚扣分，不是简单的丢弃该个体； （ⅱ）生存选择阶段引进克隆体数量限制策略。在用轮盘赌选择个体时候，个别个体的竞争力很强，会被大量的繁殖，群体逐渐同质化。该策略限制了在进化中个别个体被克隆的数量，保持了群体的多样性，防止群体的早熟收敛；（ⅲ）交叉阶段引进多点交叉，巢穴竞争选择策略。一般的进化算法是两个亲代个体交叉后产生两个子代个体。巢穴选择策略是两个亲代个体杂交后产生多个子代个体，子代个体与亲代个体竞争选择最好的两个个体遗传进化。根据氨基酸的长度，采取3点交叉，每对随机选择的亲代个体随机交叉2次产生4个后代个体；（ⅳ）局部优化策略。当算法搜索到一定阶段后，染色体进化速度骤然降低，最优个体往往停止进化。因此，对最优个体进行局部优化有利于算法跳出‘局部陷阱’。局部优化策略操作如下：首先选择群体中最优个体；再从第二个位开始，对最优个体进行随机变异操作；再计算变异个体的适应度。如果变异后个体的适应度f2大于等于变异前的适应度f1，接受变异后的新个体，最后对新个体的下一位继续进行变异操作，重复步骤（ⅲ）和（ⅳ）直到个体的所有位变异操作结束。
3 结果与讨论
３．１ 改进的遗传算法的性能比较
利用改进的遗传算法对含27个残基的标准hp序列进行了三维折叠模拟，序列如表1所示。该序列在许多文献中多次应用[1４～1６]。程序采用matlab语言编写，优化后的参数为群体规模200，交叉概率0.75，变异概率0.05。对每个序列折叠模拟20次。实验结果发现，改进后的算法性能得到了显著的提高，不仅能以较小的代价搜索到最低能量构型，而且搜索到的构型紧凑，更接近真实蛋白的结构。
表1 测试的蛋白质序列（略）
table 1 peptide length test cases
为了便于比较，对最后搜索到的最优个体采用公式（1）重新计算能量。表2是改进算法后的结果与其它标准算法结果比较。由２表可以看出，搜索到最低能量时，unger需要的能量评价函数较多，patton需要的能量评价函数大为减少，本研究需要的函数评价数目比patton算法有所减少，但是个别序列有所增多。
表2 测试能量评价数结果比较（略）
table 2 result comparison of energy evaluation
unger采用的遗传算法在初始阶段所有个体从一条直线开始变异[15]，变异后的个体用蒙特卡罗方法过滤。在交叉阶段算法实行单点交叉，交叉后的个体也用蒙特卡罗方法过滤。当产生的后代个体出现无效个体时抛弃该个体，重新产生新的个体。这种算法类似于模拟退火算法，抑制了遗传算法的搜索性能，所以该算法的能量评价数目非常多。patton的遗传算法采用相对编码，两点交叉，当出现无效个体时候，对每个重叠位置采取惩罚性扣分[1４]。patton 的算法性能得到了较大的提高，但是patton的格子模型没有考虑hp的相互作用。
3.2 改进策略的影响
本研究中的克隆体数量限制策略对维持种群的多样性起了很重要的作用。实验发现，当群体遗传一定代数后，群体进化陷入停滞，群体中的个别个体大量繁殖， 甚至占了近群体20%～50%，这样的群体很难有新的进化。采用克隆体数量限制策略有效地解决了过度繁殖的问题，该策略规定群体中相同个体不能超过一定数量，超过的部分用随机产生新的个体来代替。该策略不必频繁使用，每遗传10代使用一次比较节约资源。实验发现克隆体限制数量设定为3～6比较合适，本研究将克隆体数量限定为4。
本研究中的多点交叉策略也有利于保持个体的有效性。对于一个染色体，改变其中一个氨基酸的折叠方向将对整个个体产生巨大的影响，而多点交换策略只改变其中一段染色体的结构，降低了单点交叉带来的压力。巢穴竞争选择策略使得新的个体不仅面临与同辈个体间的竞争，也面临与父辈个体的竞争，提高了繁殖优秀个体的能力。
本研究采用的局部优化策略是系统变异，类似于monte carlo搜索方法[17]。本策略对最优染色体进行二次寻优，在算法的初期阶段爬山能力较强，但是在后期基本上失去了对染色体的改造能力，产生有效个体数不多。
3.3 改进能量关系的影响
本研究的适应度的规定与其它文献有所差别[5,7,9]，一般的适应度都是直接用hh间的接触数量表示适应度的高低，没有hh接触的个体适应度为0，没有机会参与遗传，这种个体中也存在优秀基因。本实验增加了一个常量0.01，个体适应度都不为零，所有的个体都有被选中的机会，这种策略不仅保持了群体的多样性，也使更多的优秀基因有机会参与遗传。
图1是序列27.０9的折叠模拟结果，图中黑球表示非极性分子h，白球表示极性分子p，图1a为采用公式（1）能量关系得到的结果，图1b为采用公式（2）的结果，两种结果能量值相同，但是形状不同，图1b呈球形，更接近天然蛋白的形状。原来的能量关系公式（1）只考虑hh的能量相互作用，而修正的能量关系公式（2）不仅考虑了hh的相互作用，也考虑hp的相互作用，这种模型更接近真实蛋白质内部分子相互作用机制。采用（1）的能量关系能够搜索到最低能量，但不一定是最优的构型，极性分子p不能保证一定要包围非极性分子h。为了更好的说明这个问题，考察序列p8h8p8，该序列中的h单体都集中在一起，容易折叠，分别采用两种能量关系对该序列折叠模拟，图2a、2b分别是采用公式(1)和（2）得到的结构图，两种构型都达到了它的最低hh能量5。图2a结构松散，两端的极性分子p没有包裹非极性分子h的选择压力而展开，该结构与真实蛋白分子结构相差甚远；图2b的极性分子p包围h分子，构型紧凑，接近真实蛋白质分子球形结构。说明对能量关系的修正是有效的。
图1 序列27.０9的两种不同结构（a）为未改进算法得到的结构，（b） 为改进算法后的结构，两种结构的hh键数量都是7，图中黑球表示非极性分子h，白球表示极性分子p （略）
fig.1 structures from case 27.０9 under the different models （a）: original model, (b): modified model. the number of hydrophobic contacts of both structure is 7. dark spheres are hydrophobic monomers h and white spheres are polar monomers p)。
图2 序列p8h8p8的两种不同结构a为未改进算法得到的结构，b为改进算法后的结构，两种结构的hh键数量都是5，图中黑球表示非极性分子h，白球表示极性分子p（略）
fig.2 structures from sequence p8h8p8 under the different models. (a): original model, (b): modified model. the number of hydrophobic contacts of both structure is 5. dark spheres are hydrophobic monomers h and white spheres are polar monomers p
结果表明，改进的遗传算法维持了种群的多样性，增强了算法寻优能力，提高了搜索效率，模拟的蛋白质结构紧凑，更接近真实蛋白质的构型。
【参考文献】
1 baker d & sali a. science, 2001, 294(5540): 93～96
2 anfinsen c． science, 1973, 181(96): 223～230
3 chris t, alena s, holger h h. bmc bioinformatics, 2007, 8: 342～362
4 hart w e, istrail s. journal of computational biology, 1997, 4(1): 1～22
5 unger r, moult j. j. mol. biol., 1993, 231(1): 75～81
6 andrzej s， piotr r. acta biochimica polonica, 2001, 48(1): 77～81
7 lu b z, wang b h, chen w z. protein engineering, 2003, 16(9) :659～663
８ wang yinghua(王英华), ouyang liqun(欧阳立群), li nan(李 楠), wang hongyan(王洪艳), yang bing(杨 兵). chinese j. anal. chem.(分析化学), 2006, 34(12): 1787～1790
9 jiang t, cui q, shi g, ma s. journal of chemical physics, 2003, 119 (8): 4592～4595
１０ lau k f, dill k a. macromolecules, 1989, 22(10): 3986～3997
1１ li h, helling r, tang c. science, 1996, 273(2): 666～669
1２ holland j. adaptation in natural and artificial systems. ann arbor: uni. mich. press, 1975
1３ li shaoｘin (李绍新), xing da (邢 达), qin huaｍing (秦华明), yang xiangｂo (杨湘波)，tan shiｃi (谭石慈). chinese j. anal. chem.(分析化学), 2004, 32(4): 481～484
1４ patton a l, punch ⅲ w, goodman e ｄ. ｐroceedings of the ｓixth ｉnternational conference on genetic algorithms, 1995: 574～581
1５ yu z g, en m. the journal of chemical physics, 2006, 125（２３）: 234703
1６ unger r, moult j. proceedings of the fifth international conference on genetic algotithms, 1993: 581～588
1７ andrzej s， piotr r. acta biochimica polonica, 2001, 48(1): 77～81
