CN103077226B - 一种多模态蛋白质构象空间搜索方法 - Google Patents

Abstract

本发明提出了一种多模态蛋白质构象空间搜索方法，即在排挤差分进化算法的基础上，综合了空间局部性原理和集结过程思想，并采用能量极小化过程对实验所得蛋白质构象进行处理。空间局部性原理提升了算法的收敛速度，并有效地平衡了多模态优化问题的局部收敛性和模态多样性；集结过程通过随机选取不同的交叉策略，避免了构象中较好片段被算法破坏，改善了蛋白质构象群体的多样性；能量极小化过程降低了蛋白质构象空间求解的复杂度，有效缩小了其可行域的搜索空间。本发明以脑啡肽为例，不仅得到了其公认的全局最稳定结构，还获得了一系列高质量的局部稳定结构，解决了传统针对单一靶点的单模态研究方法所不能解决的疾病多基因和多靶点通路的问题，满足了现阶段计算机辅助药物设计中对于蛋白质结构多模态的需求。

Description

一种多模态蛋白质构象空间搜索方法

技术领域

本发明设计涉及蛋白质构象空间搜索技术领域，特别是涉及一种基于排挤差分进化算法的多模态蛋白质构象空间搜索方法，属于生物信息技术、现代智能优化方法和计算机虚拟现实技术综合交叉的技术。

背景技术

2003年 4 月 14日，美国人类基因组研究项目首席科学家 Collins F博士在华盛顿隆重宣布人类基因组序列图绘制成功，人类基因组计划(Human Genome Project, HGP)的所有目标全部实现，这标志着人类基因组计划胜利完成和后基因组时代 (Post Genome Era，PGE)已来临。在当今分子生物学领域中，蛋白质分子空间结构与功能的研究无疑是最具有挑战性的问题，虽然三联密码的破译让人们了解了分子生物学的中心法则，但这些有一定的氨基酸顺序排列的多肽链如何形成具有一定空间结构的蛋白质分子的，这是生物学中心法则中目前还没解决的问题。基因组测序计划产生了大量的氨基酸序列，然而要了解蛋白质在生物体中的具体作用，仅有这些序列是不够的，还必须找到蛋白质的结构和功能。蛋白质的结构和功能是统一的，一种特定的蛋白质结构是真正理解其机理所不可缺少的。因此，根据蛋白质分子的氨基酸序列预测其三维空间结构，它将会使人们能够更系统和完整的理解生物信息从DNA到具有生物活性蛋白质的传递全过程，使中心法则得到更完整的阐明，从而对生命过程中的各种现象有进一步的深刻认识，最终推动生命科学的快速发展。

蛋白质结构预测除了其自身的理论意义外，还具有很重要的实际应用意义，知道蛋白质功能并找到其致病的分子机理，只靠氨基酸的测序是不够的，必须知道它们的空间结构，即药物设计是基于蛋白质空间结构的。现代药物工程为了抑制某些酶或蛋白质的变异，在已知其三级结构的基础上，利用分子对接算法，在计算机上设计抑制剂分子作为候选药物，一些常见病如疯牛病、白内障等都是由于蛋白质结构发生变异而引起的。人类基因组为药物开发提供了新源泉，迄今已应用的人类药物靶标约500种，包括受体、酶、信号转导分子等，开发成功的药物约2000种，估计人类基因组中3-4万个基因中，约5000个基因产物可成为潜在的药物靶标。基因组药物的种类大致可以分为：基因工程重组蛋白质药物；以人类基因编码蛋白为靶标的化学药物；以人类基因编码蛋白为靶标的人源化抗体；反义核酸类和RNA类药物等。因此，人们希望通过研究蛋白质的内在结构来了解其内在机理，并找出蛋白质折叠的真正驱动力，这样不仅可以对疾病进行有效的预测和控制，还可以按照人们的设想设计出符合特定需求的非天然蛋白质。

然而试验方法测定蛋白质结构的速度远远落后于测序速度，X射线晶体学方法和多维核磁共振技术是目前测定蛋白质结构的主要物理方法。X射线衍射分析方法能得到高分辨率的蛋白质晶体结构，但不能测定溶液中蛋白质分子的三维结构，同时该方法对蛋白质晶体之辈要求苛刻，所以带来成本高、周期长的缺点；多维核磁共振技术能够测定溶液中蛋白质分子结构，但是对样本的需求量大，样本蛋白质分子量也受到一定的限制。相比之下，蛋白质氨基酸序列被测定的数目由于DNA测序的实现而高速增长，截至目前，UniProtKB/TrEMBL数据库存放蛋白质序列已达到10,891,661条，而其中只有63,559条序列的结构被实验测定，仅占序列总数0.58%，这意味着已知序列的蛋白质数量和已测定结构的蛋白质数量的差距将会越来越大。因此，单纯依靠试验方法测定所有蛋白质的空间结构远远不能满足后基因时代对于蛋白质数据分析的需求，有必要发展一种可靠的理论预测方法，从而迅速、简便的获取蛋白质的结构信息。

现有预测蛋白质构象的方法主要有两种：同源建模法和从头预测法。同源建模法是通过比较目标蛋白质与另外一种已知结构，且其氨基酸序列与目标蛋白质联系密切的蛋白质进行预测，如果不存在与目标蛋白质氨基酸序列相似的蛋白质，则无法进行预测。而从头预测是假设折叠后的蛋白质取能量最低的构象，通过计算得到蛋白质每一部分不同的卷曲状态对应的构象，直到发现最低能量状态。由于它仅利用一级序列信息进行预测，而不需要任何其他已知蛋白质构象信息，所以该方法是比较理想的预测方法。但是在蛋白质分子设计过程中，很多情况下，预测算法所得到的蛋白质全局稳定构象可能并不满足实际的需求，例如使用计算机等技术进行蛋白质药物设计时，由于蛋白质构象预测所选取的能量模型的复杂性和实验的误差性，同时现在许多疾病涉及多基因、多靶点通路的问题，所以传统的、针对单一靶点的单模态研究方法预测所得的蛋白质全局稳定构象和实测靶点的结构并不能很好的吻合，这就需要设计一种多模态蛋白质构象空间搜索方法能够给出蛋白质其他的高质量的局部稳定构象。

近年来，许多随机全局优化算法陆续提出来解决多模态优化问题，如遗传算法(Genetic Algorithms, GA、差分进化(Differential Evolution, DE)算法，但是面对多模态函数优化时，智能优化算法要么只能寻找到问题的全局最优解，要么极易陷入局部极值解，必须与拥挤模型、适应度共享机制等小生境模型结合使用，才能找到尽可能多的极值解。2012年，K-C Wong提出的基于空间局部性原理的排挤差分进化算法(CrowdingDE-SL)，较好的解决了多模态优化问题。然而蛋白质是一个大分子的复杂体系，它的能量势能面存在着大量的局部极小，如何避免陷入局部极小的陷阱，找到全局能量极小构象成为解决多模态蛋白质构象空间搜索问题的关键。理想的蛋白质构象空间搜索方法是基于能量极小化的理论计算方法，因为它建立在热力学基础上，通常势能面可得到许多极小值的位置，其中对应于最低能量的点称为全局最小值，相当于蛋白质分子最稳定的构象。通过能量极小化过程进行结构优化时，避免了陷入局部相对稳定的构象，而力求得到全局最小值，从而实现全局优化。差分进化算法虽然具有很强的全局搜索能力，但是其局部搜索的能力较弱。本发明采用能量极小化过程，能提高差分进化算法的局部搜索能力，大大降低算法的搜索空间。特别是当群体进入一个全局最小的区域时，采用能量极小化过程可快速找到最优解。

尽管现阶段多模态蛋白质构象空间搜索取得了一些成果，但是要想精确、稳定、有效的搜索到蛋白质的稳定构象，同时还保证得到良好的模态分布性，蛋白质构象空间搜索仍然是一项艰巨的任务。其原因在于合理的多模态蛋白质构象空间搜索方法必须满足以下三个条件：一是要有一个合理的势函数，能够将蛋白质结构抽象转化成一个数学模型；二是蛋白质能量模型是一个高维的非凸函数，要保证算法在有效的计算时间找到势能函数的全局最优；三是在蛋白质分子设计过程中，可能算法预测所得的全局稳定构象并不满足实际的需求，那就要求新的算法不仅更快的得到蛋白质的全局稳定构象，还要尽可能的找到一系列高质量的局部最优构象。

发明内容

本发明克服现有技术出现的问题，提供了一种基于排挤差分进化算法的多模态蛋白质构象空间搜索方法，其综合了空间局部性原理和集结过程、能量极小化过程，不但能够快速、有效的定位到蛋白质的全局稳定构象，还能得到一系列高质量的局部稳定构象。

本发明所述的多模态蛋白质构象空间搜索方法，包括如下步骤：

1)、确定优化目标，选择蛋白质构象能量函数E(x)；

2)、初始化种群：采用随机法生成初始种群，其中N为优化问题的维数，P为种群大小，并设定算法的各控制参数；

3)、采用能量极小化过程对初始种群预处理，并按能量从小到大排序，选取种群前M个蛋白质构象作为种子个体；

4)、变异操作：对每个目标个体（其中i、m表示第m个目标个体的第i维，i∈{1,2,……,N}）做以下处理：

4.1)初始化：令G=1(G为迭代次数)，i=1；

4.2)计算种群中非亲代和亲代个体的距离，其中m≠n∈{1,2,……,P}，i∈N，并由大到小排序，再通过转换函数f(x)将距离转化为轮盘赌形式，其中转换函数f(x)形式为：，其中d为和亲代目标个体的距离，SD等于，d_MAX为距亲代最远的距离；

4.3)首先选取距亲代最近的个体作为变异的基准矢量，然后在轮盘中随机选取个体和，最后对执行变异操作得到变异测试个体，变异过程可表示为：，其中i表示目标个体的第i维，且a≠b≠c∈{1,2,……,P}，i∈N，F为变异因子；

5)、为了保证种群的多样性，算法以不同的概率（X、Y、Z）(其中X、Y、Z为 (0,1)之间不同的随机数，且X+Y+Z=1)从下列三种不同的交叉策略选取一种执行交叉操作：

5.1)以概率X执行基本的DE的交叉策略，得到交叉测试个体，其过程可表示为： $v_m^i=\left\{\begin{array}{cc}{u}_m^i& \mathrm{ifrand}\left(i\right)\≤\mathrm{CRori}=\mathrm{rnbr}\left(i\right)\\ x_m^i& \mathrm{otherwise}\end{array}\right.$ ，其中i表示目标个体的第i维，i∈N，rand(i)为(0,1)之间的随机数，CR为交叉因子，rnbr(i)为(0,N)之间的随机整数，N为优化目标的维数；

5.2)  以概率Y随机选取一个小组，并随机从M个种子中选取一个种子个体，将种子中和小组相对应的局部片段直接复制给变异测试个体相对应的片段中，得到交叉测试个体；

5.3)  以概率Z随机选取一个集合组，并随机从M个种子中选取一个种子个体，将种子中和集合组相对应的局部片段直接复制给变异测试个体相对应的片段中，得到交叉测试个体；

6)、对所得的测试个体进行能量极小化处理得到测试个体；

7)、令i=i+1，若i<=N，继续循环执行步骤4）~7），否则循环结束；

8)、对所得的测试个体执行基本DE的选择操作，如果，则替换亲代个体，否则保持种群不变；

9)、令G=G+1，若G<=CG，继续循环执行3）~7），否则循环结束；

10)、根据阀值标准D_IJ（两个构象解I和J的二面角距离）对算法所得的蛋白质构象解进行筛选，最终得到一系列高质量的蛋白质稳定构象，其中D_IJ的表达式为：

$D_{\mathrm{IJ}}=\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}\min \left[\mathrm{mod}\right\{({\mathrm{\&theta;}}_k^I-{\mathrm{\&theta;}}_k^J),\mathrm{sym}\left(k\right)\},\{\mathrm{sym}\left(k\right)-\mathrm{mod}\{({\mathrm{\&theta;}}_k^I-{\mathrm{\&theta;}}_k^J),\mathrm{sym}\left(k\right)\}\left\}\right]$

式中N为优化目标维数，θ为蛋白质构象的二面角解集，sym(k)取360o、180o、120o，表示蛋白质构象关于二面角k的对称重复结构；

11)、算法结束，输出实验所得的蛋白质稳定构象。

针对蛋白质构象空间搜索这个高维复杂的优化问题，本发明的优点是：首先，改善的排挤差分进化算法采用全局搜索和内在并行的搜索方式，能够快速的定位优化问题的全局最优解和局部极值解；同时，局部空间原理对近亲个体进行操作既保证算法快速收敛至极值点，又能防止算法陷入局部极值点；其次，针对蛋白质预测问题中高维构象空间极其复杂的瓶颈，算法采用能量极小化过程对种群进行处理，大大降低蛋白质结构预测的搜索空间；最后，采用集结过程的思想直接继承种子个体中局部优秀片段，避免算法运行过程时蛋白质结构中优秀片段被算法破坏，而且不同的交叉策略又保证了种群具有好的模态分布。

附图说明

图1为本发明中算法的算法流程图；

图2为算法中变异操作的流程图；

图3为算法中交叉操作的流程图；

图4为算法对蛋白质构象相对应小组的划分图；

图5为四种算法50次运行结果的种群平均能量分布图；

图6为本发明中算法所得到6个高质量的蛋白质稳定构象的PDB图；

具体实施方式

下面结合附图和实施例对本发明进行详细的描述。

脑啡肽（Try¹-Gly²-Gly³-Phe⁴-Met⁵）是由五个氨基酸组成的蛋白质分子，由75个原子组成，可用24个独立的主-侧链二面角描述，公认的脑啡肽稳定能量值为-11.7073kcal/mol。参照图4，算法将脑啡肽为例，将其二面角搜索范围定在-180°到180°，并将其对应的24个二面角分为8个小组，小组中的代表脑啡肽主链中的二面角，χⁱ代表脑啡肽侧链中的二面角。在算法中，我们进一步的将8个小组分为7个集合组，如下表1所示。小组和集合组中的成员分别对应24个二面角中的某些片段，这些小组和集合组类似于原始集结过程中两种不同的片段。

表1

参照图1，为本发明中的算法流程图。

1)、以脑啡肽为优化目标，选取ECEPP/3能量模型E(x)为本发明的优化模型，搜索其能量全局最低和局部最低所对应的脑啡肽构象。

2)、初始化种群：采用随机法生成初始种群，其中N为优化问题的维数，P为种群大小，并设定算法的各控制参数。

3)、采用能量极小化过程对初始种群预处理，并按能量从小到大进行排序，选取种群前M个脑啡肽构象作为种子个体。

4)、变异操作：对每个目标个体（其中i、m表示第m个目标个体的第i维，i∈{1,2,……,N}）做以下处理（如图2所示）：

4.1)初始化：令G=1(G为迭代次数)，i=1。

4.2)计算种群中非亲代和亲代个体的距离，其中m≠n∈{1,2,……,P}，i∈N，并由大到小排序，再通过转换函数f(x)将距离转化为轮盘赌形式，其中转换函数f(x)形式为：

$f\left(x\right)=\exp \{-\left(\frac{d^2}{2\&times;{\mathrm{SD}}^2}\right)\}$

其中d为和亲代目标个体的距离，SD等于，d_MAX为距亲代最远的距离。

4.3)首先选取距亲代最近的个体作为变异的基准矢量，然后在轮盘中随机选取个体和，最后对执行变异操作得到变异测试个体，变异过程可表示为：

$u_m^i=x_a^i+F\&CenterDot;(x_b^i-x_c^i)$

其中i表示目标个体的第i维，且a≠b≠c∈{1,2,……,P}，i∈N，F为变异因子。

5)、为了保证种群的多样性，算法以不同的概率（X、Y、Z）(其中X、Y、Z为 (0,1)之间不同的随机数，且X+Y+Z=1)从下列三种不同的交叉策略选取一种执行交叉操作（如图3所示）：

5.1)以概率X执行基本的DE的交叉策略，得到交叉测试个体，其过程可表示为：

$v_m^i=\left\{\begin{array}{cc}{u}_m^i& \mathrm{ifrand}\left(i\right)\&le;\mathrm{CRori}=\mathrm{rnbr}\left(i\right)\\ x_m^i& \mathrm{otherwise}\end{array}\right.$

其中i表示目标个体的第i维，i∈N，rand(i)为(0,1)之间的随机数，CR为交叉因子，rnbr(i)为(0,N)之间的随机整数，N为优化目标的维数。

5.2)  以概率Y从图4所示的小组中随机选取一个小组，并随机从M个种子中选取一个种子个体，将种子中和小组相对应的局部片段直接复制给变异测试个体相对应的片段中，得到交叉测试个体。

5.3)  以概率Z从表1所示的集合组中随机选取一个集合组，并随机从M个种子中选取一个种子个体，将种子中和小组相对应的局部片段直接复制给变异测试个体相对应的片段中，得到交叉测试个体。

6)、对所得的测试个体进行能量极小化处理得到测试个体。

7)、令i=i+1，若i<=N，继续循环执行步骤4）~7），否则循环结束。

8)、对所得的测试个体执行基本DE的选择操作，如果，则替换亲代个体，否则保持种群不变。

9)、令G=G+1，若G<=CG，继续循环执行3）~7），否则循环结束。

10)、根据阀值标准D_IJ（两个构象解I和J的二面角距离）对算法所得的蛋白质构象解进行筛选，最终得到一系列高质量的蛋白质稳定构象，其中D_IJ的表达式为：

$D_{\mathrm{IJ}}=\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}\min \left[\mathrm{mod}\right\{({\mathrm{\&theta;}}_k^I-{\mathrm{\&theta;}}_k^J),\mathrm{sym}\left(k\right)\},\{\mathrm{sym}\left(k\right)-\mathrm{mod}\{({\mathrm{\&theta;}}_k^I-{\mathrm{\&theta;}}_k^J),\mathrm{sym}\left(k\right)\}\left\}\right]$

式中N为优化目标维数，θ为蛋白质构象的二面角解集，sym(k)取360o、180o、120o，表示蛋白质构象关于二面角k的对称结构。

11)、算法结束，输出实验所得的蛋白质稳定构象。

为了同本发明中算法的实验效果进行对比，文中给出了其他3种同类的算法，算法1为差分进化算法+能量极小化，算法2为差分进化算法+集结过程+能量极小化，算法3为基于空间局部性原理的排挤差分进化算法+能量极小化，算法4为本发明中的算法（基于空间局部性原理的排挤差分进化算法+集结过程+能量极小化）。算法的参数设置如下：种群大小P=100，变异因子F=0.9，交叉因子CR=0.1，最大迭代次数CG=400，选取的种子数目M=10，X=0.5，Y=0.2，Z=0.3，二面角距离D_IJ阀值为24?（若D_IJ?24?，则认为脑啡肽的两个构象解为对称重复结构）。四种算法独立运行50次，如表2是四种算法运行结果的各项性能的平均值，图5是四种算法50次运行结果的种群平均能量分布图。

表2

从表2中运行结果可以看出，算法1、2、3、4虽然基本都能搜索到脑啡肽的最稳态结构，但是加入了集结过程的算法2、3、4运行结果明显比较好，算法4在50次运行中有46次能够搜索到脑啡肽的全局最稳态结构，且能够搜索到最低能量值-11.7073 kcal/mol，可靠性相对最好。图5显示的是四种算法50次运行的种群平均能量分布图，如图所示，加入了空间局部性原理的算法4的平均能量曲线尽管开始阶段陡峭，平均迭代188.68次达到脑啡肽的稳态结构，但是相比其他几种算法，算法4能在保证收敛速度的前提下，其种群的平均能量曲线较为平缓，得到了更多的蛋白质结构。算法4平均每次运行就能得到96.64个不同的蛋白质结构，50次运行总共得到了679个能量值大于-10 kcal/mol的脑啡肽结构，经过阀值D_IJ的筛选，得到了84个独立的高质量脑啡肽稳定结构。图6是算法4所得到的6个高质量的蛋白质稳定构象的PDB图，从图中可以直观的看出，六个三维结构图有很大的差别。为了更清楚地区分不同结构之间的差异，文中固定选取各结构三维图中三个原子，并求出其角度，发现细微的能量差，三个原子所成角度是不同的。

以上说明是本发明以脑啡肽为实例所得出的优化效果，并非限定本发明的实施范围，在不偏离本发明基本内容所涉及范围的的前提下对其做各种变形和改进，不应排除在本发明的保护范围之外。

Claims (1)

1.一种多模态蛋白质构象空间搜索方法，其特征是：在排挤差分进化方法的基础上综合了空间局部性原理、集结过程、能量极小化过程，对蛋白质构象进行优化处理，具体实现如下：

1)、确定优化目标，选取蛋白质构象能量函数E(x)；

2)、设定差分进化方法各控制参数：种群大小P，变异因子F，交叉因子CR，种群迭代最大次数CG，并采用均匀分布函数随机法生成初始种群 其中N为优化问题的维数，P为种群大小；

3)、根据能量函数E(x)对初始种群进行能量极小化预处理，并按能量从小到大进行排序，选取种群前M个蛋白质构象作为种子个体；

4)、变异操作：对每个亲代目标个体做以下处理，其中i、m表示第m个目标个体的第i维，i∈{1,2,……,N}：

4.1)初始化：令G＝1，G为迭代次数，i＝1；

4.2)计算种群中非亲代个体和亲代目标个体的距离，其中m≠n,m、n∈{1,2,……,P}，i∈N，并由大到小排序，再通过转换函数f(x)将距离转化为轮盘赌形式，其中转换函数f(x)形式为：

$f\left(x\right)=\exp \{-\left(\frac{d^2}{2\&times;{\mathrm{SD}}^2}\right)\}$

其中d为非亲代个体和亲代目标个体的距离，SD等于d_MAX为距亲代目标个体最远的距离；

4.3)首先选取距亲代最近的个体作为变异的基准矢量，然后在轮盘中随机选取个体和最后对执行变异操作得到变异测试个体变异过程可表示为：

$u_m^i=x_a^i+F\&CenterDot;(x_b^i-x_c^i)$

其中i表示目标个体的第i维，且a≠b≠c，a、b、c∈{1,2,……,P}，i∈N，F为变异因子；

5)、以不同的概率X、Y、Z选取相应的交叉策略对变异测试个体执行交叉操作得到交叉测试个体

所述的交叉操作，在交叉操作中引入了集结过程，为了保证种群的多样性，以不同的概率X、Y、Z，其中X、Y、Z为(0,1)之间不同的随机数，且X+Y+Z＝1，从下列三种不同的交叉策略选取一种执行交叉操作；

5.1)以概率X执行基本的DE交叉策略，其过程可表示为：

$v_m^i=\left\{\begin{array}{cc}{u}_m^i& \begin{array}{cc}if& rand\left(i\right)\&le;CRori=rnbr\left(i\right)\end{array}\\ x_m^i& otherwise\end{array}\right.$

其中i表示亲代目标个体的第i维，i∈N，rand(i)为(0,1)之间的随机数，CR为交叉因子，rnbr(i)为(0,N)之间的随机整数，N为优化目标的维数；

5.2)以概率Y随机选取一个小组，并随机从M个种子中选取一个种子个体，将种子中和小组相对应的局部片段直接复制给测试个体相对应的小组片段中，得到交叉测试个体

5.3)以概率Z随机选取一个集合组，并随机从M个种子中选取一个种子个体，将种子中和集合组相对应的局部片段直接复制给测试个体相对应的集合组片段中，得到交叉测试个体

6)、对所得的交叉测试个体进行能量极小化处理得到测试个体

7)、令i＝i+1，若i<＝N，继续循环执行步骤4)～7)，否则循环结束；

8)、对所得的测试个体执行基本DE的选择操作，如果则替换亲代目标个体否则保持种群不变；

9)、令G＝G+1，若G<CG，继续循环执行3)～7)，否则循环结束；

10)、根据阀值标准D_IJ对所得的蛋白质构象解进行筛选，最终得到一系列高质量的蛋白质稳定构象，其中D_IJ为两个构象解I和J的二面角距离，其表达式为：

$D_{IJ}=\underset{k=1}{\overset{N}{\&Sigma;}}min\&lsqb;\mathrm{mod}\{({\mathrm{\&theta;}}_k^I-{\mathrm{\&theta;}}_k^J),sym\left(k\right)\},\{sym\left(k\right)-\mathrm{mod}\{({\mathrm{\&theta;}}_k^I-{\mathrm{\&theta;}}_k^J),sym\left(k\right)\left\}\right\}\&rsqb;$

式中N为优化目标维数，θ为蛋白质构象的二面角解集，sym(k)取360°、180°或120°，表示蛋白质构象关于二面角k的对称重复结构；

11)、输出实验所得的蛋白质稳定构象。