阅读说明：本技术 一种基于混合烟花算法的分子构象搜索方法 (Molecular conformation search method based on mixed firework algorithm ) 是由 李瑾 刘伟超 杨佳艳 于 2020-03-19 设计创作，主要内容包括：本发明提供了一种基于混合烟花算法的分子构象搜索方法,其特征在于：包含以下步骤：S1、设置受体分子的对接区域,用对接盒子表示,对接盒子用于存放配体构象；S2、初始化若干个初始烟花,每个烟花代表一个配体构象；将配体构象表示为解向量,设置受体-配体的结合亲和力打分函数为适应度函数；S3、由解向量构建解空间,解空间包含若干个层次；S4、构建烟花算法的算子；S5、结合烟花算法与局部搜索算法构建混合烟花算法,利用混合烟花算法搜索对接盒子内的近似最优配体构象。本发明减少了分子对接中测试复合物上对接的平均用时,提高了分子对接的速度；同时能够找到适应度函数的近似最优值,提高了分子对接的精度。(The invention provides a molecular conformation search method based on a mixed firework algorithm, which is characterized by comprising the following steps of: comprises the following steps: s1, setting a docking area of the receptor molecule, and representing the docking area by a docking box, wherein the docking box is used for storing ligand conformation; s2, initializing a plurality of initial fireworks, wherein each fireworks represents a ligand conformation; expressing the ligand conformation as a solution vector, and setting a receptor-ligand binding affinity scoring function as a fitness function; s3, constructing a solution space by the solution vector, wherein the solution space comprises a plurality of layers; s4, constructing operators of the firework algorithm; s5, a mixed firework algorithm is constructed by combining the firework algorithm and the local search algorithm, and the mixed firework algorithm is used for searching for the approximately optimal ligand conformation in the butt joint box. The invention reduces the average time consumption of the butt joint on the test compound in the molecular butt joint, and improves the speed of the molecular butt joint; meanwhile, an approximate optimal value of the fitness function can be found, and the precision of molecular docking is improved.)

一种基于混合烟花算法的分子构象搜索方法

技术领域

本发明涉及计算机辅助药物设计领域，主要涉及分子对接，具体涉及一种构象搜索方法。

背景技术

分子对接是计算机辅助药物研究领域的一项重要技术，被广泛应用于新药研发的多个环节，如药物发现阶段的早期虚拟筛选、药物作用靶点的发现、药物潜在作用机制研究、药物代谢位点的预测等方面。分子对接的操作过程是逐个将已知三维结构的分子放在靶标分子的活性位点上，预测二者的结合模式和亲和力，通过打分函数挑选出接近天然构象的受体亲和力最佳的配体。

分子对接中的构象搜索是极其复杂的问题，分子构象是指一个分子中各原子的连接方式不发生变化，而相对中心原子的各原子的位置发生变化而形成的分子的不同结构。从计算上看，分子对接是一个优化任务，采用何种分子构象搜索算法对系统的搜索速度和命中率有极大的影响，目前采用的算法存在预测精度低且分子对接用时较长的问题。

发明内容

为实现上述目的，本发明采用了如下的技术方案：一种基于混合烟花算法的分子构象搜索方法，包含以下步骤：

S1、设置受体分子的对接区域，用对接盒子表示，对接盒子用于存放配体构象；

S2、初始化若干个初始烟花，每个初始烟花代表一个配体构象；将配体构象表示为解向量，设置受体-配体的结合亲和力打分函数为适应度函数；

S3、由解向量构建解空间，解空间包含若干个层次，若干个层次分别表示烟花的不同的爆炸范围；

S4、构建烟花算法的算子；

S5、结合烟花算法与局部搜索算法构建混合烟花算法，利用混合烟花算法搜索对接盒子内的近似最优配体构象。

相比于现有技术，本发明具有如下有益效果：本发明将烟花算法作为全局优化器在解空间中快速定位有希望的区域，并且烟花算法内部融入局部搜索算法，在局部进行精细搜索，增强了局部搜索能力，减少了分子对接中测试复合物上对接的平均用时，提高了分子对接的速度；同时由于局部搜索算法能够将近似解定位到局部极小值，从而找到适应度函数的近似最优值，提高了分子对接的精度。

附图说明

图1为本发明的流程图。

图2为实施例1的程序执行逻辑框图。

图3为实施例1的结合亲和力结果图。

图4为实施例1的RMSD结果图。

具体实施方式

下面结合附图及实施例对本发明中的技术方案作进一步说明。

如图1所示，本发明提出了一种基于混合烟花算法的分子构象搜索方法，包含以下步骤：

S1、设置受体分子的对接区域，用对接盒子表示，对接盒子用于存放配体构象；

S2、初始化若干个初始烟花，每个初始烟花代表一个配体构象；将配体构象表示为解向量，设置受体-配体的结合亲和力打分函数为适应度函数；

S3、由解向量构建解空间，解空间包含若干个层次，若干个层次分别表示烟花的不同的爆炸范围；

S4、构建烟花算法的算子；

S5、结合烟花算法与局部搜索算法构建混合烟花算法，利用混合烟花算法搜索对接盒子内的近似最优配体构象。

在步骤S1中，设置对接盒子的中心位置和边长，对接盒子的中心设为受体结合口袋中的所有原子的几何中心，对接盒子的边长是由受体结合口袋中任意两个原子的最大距离决定。

在步骤S2中，采用如下优化问题模型表示分子构象搜索过程：

min F(X)

s.t.g(X)≤0

约束函数g(X)将配体构象约束在蛋白质的对接盒子中，对于所有不在对接盒子中的解，均视为无效解；目标函数F(X)对应的是受体-配体结合亲和力打分函数即适应度函数。该优化问题的解表示为i＝1,2,……N，其中，表示配体构象的中心位置，是配体构象的旋转角，即配体构象的取向信息，的取值范围为[-π,π]；分别表示配体构象中每个可扭转键键角，取值范围为[-π,π]，其中1,2,…,b表示配体中的b个可扭转键的标号。

在步骤S3中，由解向量构建解空间，解空间包含三个层次，分别表示烟花的三个爆炸范围:第一个层次是指在配体构象的取向及可扭转键键角固定的条件下，配体构象的中心位置变化的范围；第二个层次是指在配体构象的中心位置及可扭转键键角固定的条件下，配体构象的取向旋转的范围；第三个层次是指在配体构象的中心位置及取向固定的条件下，配体构象的可扭转键旋转的范围。

步骤S4包括：

(1)构建爆炸算子

爆炸算子包括烟花爆炸范围和烟花爆炸产生的火花数量。

烟花爆炸范围计算方法：将烟花按适应度函数值由小到大进行排序后分为三个等份，优选地，第一个等份的烟花是适应度函数值较小的烟花，即适应度函数值属于第一预设范围的烟花，属于优质烟花，在第三个层次进行爆炸，随机改变烟花可扭转键键角，模拟烟花小范围爆炸现象；第二个等份的烟花包含适应度函数值属于第二预设范围的烟花，在第二个层次进行爆炸，可随机改变其取向模拟爆炸现象；第三个等份的烟花包含适应度函数值属于第三预设范围的烟花，在第一个层次进行爆炸，随机改变其中心位置，模拟烟花大范围爆炸现象。烟花爆炸后产生爆炸火花，每个爆炸火花代表一个配体构象。

烟花爆炸产生的火花数量计算公式：

S_i＝M-r_i

其中，S_i表示烟花X_i爆炸产生的火花数量，M为每个烟花最多爆炸的火花数目，r_i表示烟花X_i的适应度函数值的排序序号；

(2)构建变异算子

采用高斯变异的方式来产生变异烟花，以增加种群多样性，避免陷入局部最优；从当代烟花与爆炸火花中随机选择L个个体，从每个个体即解向量中随机选择若干个维度，将这些维度上的值乘以一个服从高斯分布的变量得到变异火花个体，每个变异火花代表一个配体构象。

第i个变异火花的第k个变量按如下公式进行变异：

其中，Gaussian(1,1)表示均值和方差均为1的高斯分布随机数，若变异火花为不可行解，则在解空间里随机生成一个变异火花。

(3)烟花选择策略

将当代烟花中的N个烟花、S个爆炸火花和L个变异火花作为候选集合，从中选择N个烟花个体作为下一代烟花。选择策略为：候选集合中适应度值最小的烟花被确定性地选择到下一代作为烟花，而对剩下的N-1个烟花用轮盘赌的方法进行选择,其中N为大于或等于2的整数。对于候选者a，其被选中概率的计算公式为：

上述公式中，RMSD_ab，用来衡量两个具有相似三维结构的小分子a和b之间的差异，计算的是结构a和b之间N个原子位置的均方根误差。其中，原子i在结构a中的位置为原子i在结构b中的位置为

上述公式中，RMSD_a为当前烟花个体a到候选者集合K中除了a以外所有烟花个体的RMSD之和，表示烟花a的差异性。若差异性较大，则表示与该个体相似的烟花个体较少，则该个体被选中的概率较大。

步骤S5包括：

S51、对当代烟花进行局部搜索，根据metropolis准则来决定是否接受这个局部最优值，这一局部搜索过程即是推举局部区域优秀个体的过程；

S52、根据烟花爆炸范围的计算公式和产生火花数量的计算方法，将当代烟花进行爆炸，产生爆炸火花；

S53、对爆炸火花进行局部搜索，根据metropolis准则来决定是否接受这个局部最优值；

S54、当代烟花与爆炸火花中随机选择L个个体根据变异算子进行变异，得到变异火花；

S55、根据下一代烟花选择策略，选择N个个体组成下一代烟花；

S56、烟花的适应度函数值收敛或者达到最大迭代次数，则搜索结束；否则，转入步骤S51。

局部搜索算法包含BFGS算法、DFP算法或L-BFGS算法中的至少一种。

实施例1中，将PDBbind标准测试集中核心集的195个受体-配体复合物对接作为测试对象，利用本发明搜索方法进行分子对接，每个复合物对接分别运行30次。图2是本发明实施例1的流程框图，以下结合图2进行详细阐述。

优选地，本实施例在分子对接软件Autodock Vina(简称为Vina)的框架中进行实现，采用Vina的打分函数作为适应度函数。具体步骤如下：

(1)将蛋白质的pdb文件以及配体的mol2文件分别转换为pdbqt文件；

(2)设置受体结合口袋范围，即设置对接盒子的中心位置(center_x,center_y,center_z)，长、宽、高(size_x,size_y,size_z)；

(3)产生初始烟花：初始化N个初始烟花(每个烟花代表一个配体的构象)，将配体构象表示为解向量，设置受体-配体的结合亲和力打分函数为适应度函数；

(4)将配体构象解空间的搜索范围分为三个层次来表示烟花不同的爆炸范围；

(5)构建烟花算法的爆炸算子、变异算子及下一代烟花选择策略；

(6)利用BFGS算法对每个当代烟花进行局部搜索，其中初始烟花为第一代烟花，按metropolis准则来决定是否接受这个局部最优值；

(7)根据烟花爆炸范围的计算公式和产生火花数量的计算方法，将当代烟花进行爆炸，产生爆炸火花；

(8)利用BFGS算法对每个爆炸火花进行局部搜索，按metropolis准则来决定是否接受这个局部最优值；

(9)烟花变异：当代烟花与爆炸火花中随机选择L个个体，每个个体随机选择若干个维度按变异算子进行变异，产生变异火花；

(10)根据下一代烟花选择策略，选择N个个体组成下一代烟花；

(11)适应度函数值是否收敛？是，转步骤(13)，否则，转步骤(12)；

(12)是否达到最大迭代次数？是，转步骤(13)，否则，转步骤(6)；

(13)输出近似最优配体构象以及结合亲和力。

利用测试集对预测模型在三个评价指标上对本发明的方法进行测试：

(1)分子对接运行时间t(秒)

(2)结合亲和力BA(binding affinity,kcal/mol)

(3)均方根误差RMSD

分子对接运行时间t是指在若干个测试复合物上对接的平均用时，其中单个受体-配体复合物对接一次的用时是指从开始对接到输出结果的时间，这个指标用于测试对接程序的运行效率。

结合亲和力BA是指搜索算法最终找到的适应度函数的近似最优值，这个数值越小说明搜索到的结果越好，这个指标用于测试搜索方法的准确度。

均方根误差RMSD是用来衡量对接软件预测的配体构象与天然构象结构的差异性，对于两个具有相似三维结构的小分子a和b，如下定义RMSD：

RMSD_ab计算的是a和b两个结构之间N个原子位置的均方根误差。其中，原子i在结构a中的位置为原子i在结构b中的位置为这个指标用于测试搜索方法预测配体构象的准确性。

在PDBbind数据集的195个复合物上用Vina对接软件中原配体构象搜索方法做分子对接，平均用时27.05秒，而将配体构象搜索方法替换为本发明搜索方法后，分子对接平均用时仅为12.86秒，说明本发明搜索方法能大大缩短分子对接时间。

图3所示为实施例1测试的结合亲和力结果示例图，对接运行30次，每次记录从第一次对接到当前对接的亲和力最小值。图4所示为实施例1测试的RMSD结果示例图，RMSD为图3中每次记录的亲和力值对应的配体与天然晶体构象结构的RMSD值。如图3所示，本发明搜索方法预测的最佳配体构象的结合亲和力在30次运行结果上均较小，说明本发明搜索方法搜索到的适应度函数的准确度更优；图4中，本发明搜索方法预测的最佳配体构象的RMSD结果在相同运行次数的结果上均小于Vina，说明本发明搜索方法搜索到的最佳结合配体构象更接近于天然构象，准确度更高。图3和图4的结果说明本发明搜索方法搜索到的最佳配体构象的平均结合亲和力值及RMSD更低，对接结果精确度更高。

最后说明的是，以上实施例仅用以说明本发明的技术方案而非限制，尽管参照较佳实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的宗旨和范围，其均应涵盖在本发明的权利要求范围当中。