一种基于树搜索和片段组装的蛋白质结构预测方法

摘要

一种基于树搜索和片段组装的蛋白质结构预测方法，包括以下步骤：A1、获取蛋白质的pdb格式的文件并清洗出所需数据；A2、生成片段库；A3、选取力场模型；A4、采用Rosetta的score3能量函数；A5、将整个能量图景离散化，并将每层进一步离散化，分成单独的分块区，每次搜索根据能量的权重随机选取一个能量层，并在该能量层依概率选取一个分块区，如果分块区中包含有构象，采用片段组装方法，随机选取序列上的一个片段，然后在片段库中随机选取一个片段，替换序列上的目标片段，用蒙特卡洛准则来判断是否接受这个构象，如果接受则将它放入一个集合当中。本发明可以大大减少计算量和计算时间，保证搜索到能量较低的构象。

说明书

技术领域

本发明涉及计算机应用领域，生物信息学，最优化理论，分子生物学等，特别涉及一种蛋白质三维结构预测方法，属于将现代智能优化方法应用到蛋白质三维结构预测。

背景技术

生物信息学通过综合利用生物学，计算机科学和信息技术而揭示大量而复杂的生物数据所赋有的生物学奥秘。是当前研究的一个热点。生物信息学研究成果已经被广泛应用于序列比对，蛋白质比对，基因识别分析，分子进化，序列重叠群装配，遗传密码，药物设计，生物系统，蛋白质结构预测等。其中蛋白质结构预测是生物信息学领域的一个重要的分支。著名的诺贝尔化学奖获得者Anfinsen通过对核糖核酸酶A的经典研究表明去折叠的蛋白质在体外可以自发的进行再折叠，仅仅是序列本身已经包括了蛋白质正确折叠的所有信息，并提出蛋白质折叠的热力学假说，这一理论认为结构同源的蛋白质可以通过不同的折叠途径形成相似的天然构象，蛋白质的一级结构决定了其三维结构，即蛋白质的氨基酸序列决定了蛋白质的三维空间结构，从热力学分析，蛋白质的天然构型对应着其物理能量最小的状态。

随着2000年基因组工作草图的完成，基因组的工作重心已经从测序转向了基因功能的识别，这一问题的解决显得尤为重要。人们通过研究蛋白质的内在结构来了解其内在机理，并找出蛋白质折叠的真正驱动力，这样可以按照人们的设想设计出符合特定需求的非天然蛋白质，对疾病进行有效的预测和控制。蛋白质结构从头预测也因此具有重大的意义。然后蛋白质的结构测定远远赶不上基因组测序的速度，单纯的依靠实验的方法测定所有的蛋白质的三维结构已经远远不能满足后基因时代对于蛋白质数据分析的需求，有必要发展一种可靠的理论预测方法，从而迅速、简便的获取蛋白质的结构信息。

蛋白质的从头预测方法就是选择一种力场模型，将蛋白质结构抽象成一个数学模型，使其具有可计算性，然后进行算法的设计，从算法的时间和空间复杂度出发，利用全局以及局部的优化算法在能量曲面上搜索全局最优解，在数学上属于一种极其复杂的非凸多极值优化问题。利用蛋白质的一级序列，利用热力学方法，计算得出蛋白质的三维结构。

针对力场模型的优化问题，2005年Bradley在结合蒙特卡洛优化方法开发了分子建模Rosetta程序，成功预测CASP6中T0281目标蛋白(包含70个残基)的C_α-RMSD为1.6对某些小目标蛋白(<85个残基)，甚至得到误差精度小于1.5的高分辨率三维解析结构；2006年，Zhan采用Basin-Paving(BP)方法得到了脑啡肽在ECEPP/2、ECEPP/3两种力场模型下最小能量构象；2009年Roosi提出几种BP改进策略，并给出了比较结果。结合随机全局优化算法，如遗传算法(GeneticAlgorithms,GA)、差分进化(Differential Evolution,DE)算法也应用于力场模型的优化，但是由于蛋白质是一种强柔性的大分子体系，其势函数表达式极其复杂，存在极多局部极小点，这些算法极容易陷入某个局优解；同时由于随机算法缺乏全局收敛理论依据及解的不确定性，使用了全局选择因子，算法只能收敛到全局最优解，而忽略了众多局部极值解，进一步限制了它们在实际问题中的应用。

我们了解到蛋白质迅速折叠是因为随机热运动使得构象发生改变，而根据漏斗状的能量图景中所体现的原理，构象变化会使蛋白质沿着能量下降的方向搜索到天然结构。

然而，如果采用穷举搜索的方式来遍历整个能量图景，那是非常不科学的方式，首先要想遍历整个能量图景，其计算量非常庞大，这个值与离散化能量图景后的最小单位值相关，离散得越精细，计算量相应越大，相应的时间成本也越高，需要的计算机资源也越多。

发明内容

为了克服现有的蛋白质结构预测方式的计算量较大、计算时间很长的不足，本发明提供一种基于树搜索和片段组装的蛋白质结构预测方法，采用将能量图景适度离散化并依概率随机选择能量层和分块区的方式，可以大大减少计算量，计算时间，同时保证搜索到能量较低的构象，达到预测蛋白质结构的目的。

本发明解决其技术问题所采用的技术方案是：

一种基于树搜索和片段组装的蛋白质结构预测方法，所述预测方法包括以下步骤：

A1、获取蛋白质的pdb格式的文件并清洗出所需数据；

A2、生成片段库；

A3、选取力场模型，力场模型表示形式如下：

E_protein＝W_inter repE_inter rep+W_inter atrE_inter atr+W_solvationE_solvation

+W_bb/sc hbE_bb/sc hb+W_bb/bb hbE_bb/bb hb+W_sc/sc hbE_sc/sc hb

+W_pairE_pair+W_dunbrackE_dunbrack+W_ramaE_rama+W_referenceE_reference    （1）

式中，E_protein表示蛋白质的总能量，E_inter rep表示范德华排斥力作用，W_inter rep为E_inter rep在整体的权重，E_inter atr表示范德华吸引力作用，W_inter atr为E_inter atr在整体的权重，E_solvation为Lazarids和Kar_plus描述的隐含的溶解作用，W_solvation为E_solvation在整体的权重，E_bb/sc hb、E_bb/bb hb、E_sc/sc hb为依赖方向的氢键能量，W_bb/sc hb、W_bb/bb hb、W_sc/sc hb分别为其能量在整体的权重，E_pair为残基对静电作用，W_pair为E_pair在整体的权重，E_dunbrack为氨基酸基于旋转异构体库的内部的能量，W_dunbrack为E_dunbrack在整体的权重，E_rama为参考特定位置的Ramachandrin骨链扭力，W_rama为E_rama在整体的权重，E_reference为未折叠态的蛋白质的参考能量，W_reference为E_reference在整体的权重，Rosetta的能量函数就是将所有的能量项通过各自的权重线性相加；

用于计算能量层的权重W(l)的公式：

$>W\left(l\right)=E_{\mathrm{protein}}^2+\frac{1.0}{2^{22}}---\left(2\right)$>

W(l)为第L层的权重值，E_protein是蛋白质构象的能量值，该能量层被选中的概率P(l)为：

$>P\left(l\right)=\frac{W\left(l\right)}{\underset{l=1}{\overset{100}{\Sigma }}W\left(l\right)}---\left(3\right)$>

某一个分块区被选中的的概率P(cell)为：

$>P\left(\mathrm{cell}\right)=\frac{1.0}{[(1.0+\mathrm{nsel})*\mathrm{nconfs}]}---\left(4\right)$>

其中nsel记录的是该分块区被选中的次数，nconfs记录的是绑定到该分块区上的构象的数目；

A4、选取能量函数，采用Rosetta的score3能量函数；

A5、蛋白质结构预测过程如下：

Step1：初始化工作，导入各种程序运行所需的包和模块，氨基酸序列的初始化，创建相应的文件夹存放低能量的构象和相关数据和图表，创建MC对象和Mover对象等；

Step2：进入迭代循环，给定迭代次数，由迭代次数控制循环次数；

Step3：根据权重函数选择相应的能量层，权重函数列表由公式（2）动态生成并扩展，E_protein由score3函数计算得出并放入一个能量列表中；

Step4：依概率随机选中特定的分块区，概率由公式（4）计算得出；

Step5：首先判断该分块区中是否存在构象，如果存在则取出该构象；

Step6：对取出的蛋白质构象进行扰动，也就是片段组装，随机选取该构象中的一段长度为3的片段，然后随机选取片段库中与其序列相同的片段，其本质就是替换该片段的3个二面角的值；

Step7：用蒙特卡洛评判标准判断构象是否接受；如果接受则放入集合，不接受则跳转到Step2，进入下一次迭代，直到循环次数达到设定的循环值则停止。

进一步，所述步骤A1中，从RCSB官网上下载所需的蛋白质pdb文件，使用python脚本语言，编写一个脚本程序，选出pdb文件中包含蛋白质原子详细信息的那段信息，截取出来并保存成一个新的pdb文件。

本发明的有益效果主要表现在：本发明降低了算法的搜索空间复杂度，同时利用蒙特卡洛统计方法的判别，提高了算法的收敛能力，弥补了容易陷入局部极小值、搜索空间过于复杂等缺陷。

附图说明

图1是一种基于树搜索和片段组装的蛋白质结构预测方法的流程图。

具体实施方式

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

参照图1，一种基于树搜索和片段组装的蛋白质结构预测方法，其中包含以下步骤：

A1、获取蛋白质的pdb格式的文件并清洗出所需数据。从RCSB官网上下载所需的蛋白质pdb文件。由于pdb格式的文件中包含了许多在预测中不需要的信息，所以需要进一步“清洗”才能用于下一步的结构预测。实现“清洗”pdb文件的方法：使用python脚本语言，编写一个脚本程序，选出pdb文件中包含蛋白质原子详细信息的信息，截取出来并保存成一个新的pdb文件。

A2、从相关网站下载生成片段库。

A3、选取力场模型。采用Rosetta力场模型能量函数。选取力场模型。力场模型的选取直接关系到计算蛋白质能量的准确性，目前已有的能量模型都不能绝对精确的计算其能量值，因此我们选用了一种非常接近真实值的力场模型，模型表示形式如下：

E_protein＝W_inter repE_inter rep+W_inter atrE_inter atr+W_solvationE_solvation

+W_bb/sc hbE_bb/sc hb+W_bb/bb hbE_bb/bb hb+W_sc/sc hbE_sc/sc hb

+W_pairE_pair+W_dunbrackE_dunbrack+W_ramaE_rama+W_referenceE_reference

A4、选取能量函数，采用Rosetta的score3能量函数，score3能量函数属于粗粒度能量函数中计算最精确的一种；

A5、这个预测算法采用python语言编写，算法思想是将整个能量图景离散化，分成100层，并将每层进一步离散化，分成一个个单独的分块区，每次搜索根据能量的权重随机选取一个能量层，并在该能量层依概率选取一个分块区，如果这个分块区中包含有构象，则给它一个扰动，用蒙特卡洛准则来判断是否接受这个构象，如果接受则将他放入一个集合当中。第二阶段采用片段组装方法，随机选取序列上的一个片段，然后在片段库中随机选取一个片段，替换序列上的目标片段，从而得到一个新的蛋白质构象；

蛋白质结构预测过程如下：Step1：初始化工作，导入各种程序运行所需的包和模块，氨基酸序列的初始化，创建相应的文件夹存放低能量的构象和相关数据和图表，创建MC对象和Mover对象等，并将片段组装的次数设置为10；

Step2：进入迭代循环，迭代次数设置成12000次，由迭代次数控制循环次数；

Step3：根据权重函数选择相应的能量层，权重函数列表由公式（2）动态生成并扩展，E_protein由score3函数计算得出并放入一个能量列表中；

Step4：依概率随机选中特定的分块区，这里的概率由公式（4）计算得出；

Step5：首先判断该分块区中是否存在构象，如果存在则取出该构象；

Step6：对选定的蛋白质构象进行扰动，也就是片段组装，随机选取该构象中的一段长度为3的片段，然后随机选取片段库中与其序列相同的片段，替换该片段的3个二面角；

Step7：用蒙特卡洛评判标准判断构象是否接受；如果接受则放入集合，不接受则跳转到Step2，进入下一次迭代，直到循环次数达到设定的循环值则停止。

采用两阶段方法，第一阶段离散能量图景并离散每个能量层，依概率随机选中每个分块区，第二阶段片段组装，将搜索到的构象进行片段组装，搜索出能量更低的构象，可以大大减少计算量和计算时间，同时保证搜索到能量较低的构象，达到预测蛋白质结构的目的。

采用蒙特卡洛评判标准来评判是否接受构象，可以越过能量障碍，引导搜索方向指向能量更低的方向，从而有可能达到全局能量最小区域，找到更加接近天然构象的蛋白质结构。

以编号为1VII的蛋白质为例，它是由36个氨基酸组成的蛋白质分子，氨基酸序列为:MLSDEDFKAVFGMTRSAFANLPLWKQQNLKKEKGLF

一种基于树搜索和片段组装的蛋白质结构预测方法，其中包含以下步骤：

A1、获取蛋白质的pdb格式的文件并清洗出所需数据。从RCSB官网上下载所需的蛋白质pdb文件。由于pdb格式的文件中包含了许多在预测中不需要的信息，所以需要进一步“清洗”才能用于下一步的结构预测。实现“清洗”pdb文件的方法：使用python脚本语言，编写一个脚本程序，选出pdb文件中包含蛋白质原子详细信息的信息，截取出来并保存成一个新的pdb文件。

A2、从相关网站下载生成片段库。

A3、选取力场模型。采用Rosetta力场模型能量函数，选取力场模型。力场模型的选取直接关系到计算蛋白质能量的准确性，目前已有的能量模型都不能绝对精确的计算其能量值，因此我们选用了一种非常接近真实值的力场模型，模型表示形式如下：

E_protein＝W_inter repE_inter rep+W_inter atrE_inter atr+W_solvationE_solvation

+W_bb/sc hbE_bb/sc hb+W_bb/bb hbE_bb/bb hb+W_sc/sc hbE_sc/sc hb

+W_pairE_pair+W_dunbrackE_dunbrack+W_ramaE_rama+W_referenceE_reference

A4、选取能量函数，采用Rosetta的score3能量函数，score3能量函数属于粗粒度能量函数中计算最精确的一种；

A5、蛋白质结构预测过程如下：Step1：初始化工作，导入各种程序运行所需的包和模块，氨基酸序列的初始化，创建相应的文件夹存放低能量的构象和相关数据和图表，创建MC对象和Mover对象等，并将片段组装的次数设置为10；

Step2：进入迭代循环，迭代次数设置成12000次，由迭代次数控制循环次数；

Step3：根据权重函数选择相应的能量层，权重函数列表由公式（2）动态生成并扩展，E_protein由score3函数计算得出并放入一个能量列表中；

Step4：依概率随机选中特定的分块区，这里的概率由公式（4）计算得出；

Step5：首先判断该分块区中是否存在构象，如果存在则取出该构象；

Step6：对选定的蛋白质构象进行扰动，也就是片段组装，随机选取该构象中的一段长度为3的片段，然后随机选取片段库中与其序列相同的片段，替换该片段的3个二面角；

Step7：用蒙特卡洛评判标准判断构象是否接受；如果接受则放入集合，不接受则跳转到Step2，进入下一次迭代，直到循环次数达到设定的循环值则停止。

所述方法中，其中，在步骤A1中，对相关蛋白质实行数据清洗操作。

所述方法中，其中，在步骤A2中，从相关网站生成蛋白质片段库文件，用于后续的片段组装。

所述方法中，其中，在步骤A3中，力场模型的选取直接关系到计算蛋白质能量的准确性，目前已有的能量模型都不能绝对精确的计算其能量值，因此我们选用了一种非常接近真实值的力场模型，采用Rosetta软件提供的力场模型能量函数。

所述方法中，其中，在步骤A4中，选用粗粒度级别的能量函数，同时为了保证一定的精确度，因此选用score3能量函数来计算构型能量值。

所述方法中，其中，在步骤A5中，通过能量层权重列表和分块区权重列表，可以搜索变得更加高效，同时降低搜索空间的复杂度。

以上阐述的是本发明给出的一个实施例表现出来的优良效果，显然本发明不仅适合上述实施例，在不偏离本发明基本精神及不超出本发明实质内容所涉及内容的前提下可对其做种种变化加以实施。