一种计算机预测蛋白功能的方法

摘要

本发明涉及一种计算机预测蛋白功能的方法，包括以下步骤：对生物大分子的原子之间的距离进行半粗粒化统计，得到生物大分子之间相互作用的半粗粒化统计势函数；利用半粗粒化统计势函数在目标蛋白周围空间找到相互作用的热点区域；结合分子碎片生长方法，在热点区域寻找可有与目标蛋白紧密结合的蛋白质、DNA或RNA序列片段；计算蛋白质、DNA或RNA序列片段与目标蛋白结合的强度并排序；选取规定数量的序列片段及其结构，在相关的生物功能数据库中进行搜索，查找相应的生物功能，预测到目标蛋白可能具有的功能。本发明具有已有粗粒化统计势函数计算速度快的优点，而且也将与分子碎片生长方法结合的更加紧密、具有目前已有粗粒化统计势所不具备的准确性。

说明书

技术领域

本发明涉及的是一种生物信息技术，具体的说是一种计算机预测蛋白功能的方法。

背景技术

预测蛋白功能的计算方法主要有四种：

一类是基于结构相似性的方法。这类方法利用三维结构的相似程度，分析目标蛋白最可能有功能的结构单元，在已知蛋白质数据库中搜索与这些结构单元相似的蛋白质的功能数据，从而达到预测目标蛋白功能的目标。有些算法可以不用考虑拓扑结构，而只考虑空间结构。

一类是基于三维基序的方法。这些方法更多的是试图找到蛋白质中具有功能的结构元素，比如在很多DNA结合蛋白中都有发现的Helix-Turn-Helix(HTH)基序。这样，人们就可以通过类似的数据库近似得到具有这些基序的蛋白质的功能。比如用TESS算法在酶的活性位点建立起来的PROCAT数据库；直接利用蛋白质结构数据库(Protein Database Bank，简称PDB)文件中的site信息建立的PDBSite库等。

一类是基于表面的方法。这类方法通常会根据各原子的三维坐标和互相之间的距离来产生蛋白质的表面，而认为这些表面的形状及相互之间吻合的程度是决定蛋白质功能的主要因素。当然，也有研究表明某些蛋白质的功能是由表面的疏水性质或者静电性质决定的。这类方法通常被叫做分子对接方法，在计算机辅助药物设计中有着广泛的应用。

还有一类是基于学习的方法。这类方法包括数据挖掘、人工智能，各种方法层出不穷，比如使用支持向量机，决策树，人工神经网络等等等等。他们用各种方式试图找到蛋白质的功能信息和蛋白质的序列结构等信息之间的关联(显示的或隐式的)，以达到预测的目的。

基于知识的统计势函数是以已知稳定结构的一些特性统计平均为基础的势能函数，基本假设是，测到的稳定结构的势能最低。这样在稳定结构中统计到的各种特性数量，并通过一些类似波尔兹曼分布函数的方法得到各种特性的势能，然后对其他结构进行打分。与传统的以物理为基础势能函数相比，它们不太关心计算势能时候的具体物理规律，而是更关心它们需要用来计算势能的特性是否能被统计平均。比如有的统计势以二面角做统计，有的以溶液的表面做统计，有以接触对做统计的，有以距离对做统计的。以距离对做统计的基于知识的统计势函数其基本计算方法如下：对假定的目标蛋白对目标蛋白把所有原子在以r_cut为半径的球体内部分为很多的球壳(r，r+Δr]，这样统计目标蛋白中原子类型i的球体中原子类型为j原子在不同球壳中出现的数量为N_obs(i，j，r)。通过它与期望的分布值比较，可以得到原子对i，j在(r，r+Δr)距离的能量值：

$>u(i,j,r)=-\mathrm{RT}\ln \frac{N_{\mathrm{obs}}(i,j,r)}{N_{\exp }(i,j,r)}$>

基于知识的统计势函数按照参照态N_exp(i，j，r)的不同而不同。

目前现有的能够同时在更大的序列和结构范围内寻找与目标蛋白可能结合的蛋白质、DNA和RNA片段方法主要有两类，一类是基于随机性的方法，比如蒙特卡洛、遗传算法等，另外一类是基于确定性的方法，例如dead-end-elimination，branch-and-bound等，但是它们都只能在给定的序列长度以及固定的蛋白质骨架下进行序列和构象搜索，而且不能够对DNA和RNA进行序列和构象搜索。

分子碎片生长法是比较经典的计算机辅助药物设计的方法之一，被广泛用来进行小分子的全新设计和改造。分子碎片生长方法具备所有已有蛋白质序列和构象搜索方案所不具备的独特的优点。一个就是能够同时在序列的长度和组成以及结构空间进行搜索，另外一个优点就是完全不依赖于目标蛋白质的构象数目，克服了基于结构相似性分析的蛋白质功能预测方法的主要缺陷，而且蛋白质的构象可以是实验得到的，也可以是通过分子动力学模拟或者其他理论计算方法产生的。

聚类分析指将物理或抽象对象的集合分组成为由类似的对象组成的多个类的分析过程。它是一种重要的人类行为。聚类分析的目标就是在相似的基础上收集数据来分类。基于网格的聚类方法首先将对象空间划分为有限个单元以构成网格结构，然后利用网格结构完成聚类。

同源建模是利用序列相似的蛋白质(已知结构)的结构作为“模板”，以得到未知结构蛋白存在时比较可能的结构。如果能够获得比较好的所要预测结构的蛋白和已知蛋白序列的比对结果，则能够获得精确的结构模型。

使用全原子的统计势能函数计算量大，不利于构象空间的搜索，而粗粒化的统计势能函数丢失的信息比较多。

发明内容

针对现有技术中预测蛋白功能的计算方法存在计算量大，不利于构象空间的搜索以及丢失的信息比较多等不足之处，本发明要解决的技术问题是提供一种减少计算量、将与分子碎片生长方法结合的更加紧密的计算机预测蛋白功能方法。

为解决上述技术问题，本发明采用的技术方案是：

本发明计算机预测蛋白功能的方法包括以下步骤：

对生物大分子的原子之间的距离进行半粗粒化统计，得到生物大分子之间相互作用的半粗粒化统计势函数；

利用半粗粒化统计势函数在目标蛋白周围空间找到相互作用的热点区域；

结合分子碎片生长方法，在热点区域寻找可有与目标蛋白紧密结合的蛋白质、DNA或RNA序列片段；

计算蛋白质、DNA或RNA序列片段与目标蛋白结合的强度并排序；

选取规定数量的序列片段及其结构，在相关的生物功能数据库中进行搜索，查找相应的生物功能，预测到目标蛋白可能具有的功能。

所述半粗粒化统计势函数通过以下过程得到：

读入蛋白质三维结构，按照相应原子在其所在氨基酸残基或DNA碱基或RNA碱基类型上的拓扑位置给定其原子类型

利用给定类型统计原子对之间的距离分布，通过与期望的分布值比较，得到半粗粒化统计势函数，通过以下公式将分布值转化成统计势函数：

$>\bar{u}(i,j,r)=\left(\begin{array}{cc}-\mathrm{RT}\ln \frac{P^c(i,j,r)}{\frac{\left(r^{\alpha }\mathrm{\Delta r}\right)f_V\left(r\right)}{{\Sigma }_r\left(r^{\alpha }\mathrm{\Delta r}\right)f_V\left(r\right)}},& r<r_{\mathrm{cut}},\\ 0,& r\ge r_{\mathrm{cut}},\end{array}\right)$>

其中，R为理想气体常数，T为温度取值300K，r为原子对之间的距离，Δr为统计区间宽度，r_cut为截断距离。α和N₀为经验参数分别取值1.61和75。把距离在区间内的原子对按照其相应目标蛋白中的原子类型i和相应预测蛋白、DNA或RNA中的粗粒化原子类型j统计为N_obs(i，j，r)。为N_obs(i，j，r)中原子类型对相应的两个原子一个是全原子模型并且另一个是粗粒化模型的部分。则中间变量

$>f_V\left(r\right)={\Sigma }_{i,j}{N}_{\mathrm{obs}}^{\mathrm{Protein}-\mathrm{CG}}(i,j,r)/{\Sigma }_{i,j}{N}_{\mathrm{obs}}(i,j,r),$>中间变量

$>P(i,j,r)=N_{\mathrm{obs}}^c(i,j,r)/{\Sigma }_r{N}_{\mathrm{obs}}^c(i,j,r),$>中间变量

$>N_{\mathrm{obs}}^c(i,j,r)=N_{\mathrm{obs}}(i,j,r)+N_0\frac{{\Sigma }_{i,j}{N}_{\mathrm{obs}}^{\mathrm{Protein}-\mathrm{CG}}(i,j,r)}{{\Sigma }_{i,j,r}{N}_{\mathrm{obs}}^{\mathrm{Protein}-\mathrm{CG}}(i,j,r)}$>

所述分子碎片生长方法为：

以构成天然蛋白质的20种氨基酸，构成DNA的4种脱氧核苷酸以及构成RNA的4种核苷酸做为片段前体；

对片段前体进行修改后形成生长使用的片段；

父代的一个分子选择一个片段，连接片段，调整新连接片段的构象，使每一个构象都对应一个子代分子；

通过计算每个子代分子和目标蛋白质之间的结合能量，按照带权随机法选择一批子代进行下一代的繁殖。

在目标蛋白周围把空间划分为以立方体为结构单元的格点，以格点位置为假想各种粗粒化粒子的坐标，在每个格点上使用各种粗粒化类型并使用导出的半粗粒化统计势能函数评价格点和目标蛋白之间的能量，其中能量最低的值留做此格点的能量；

使用能量阀值删掉能量比较低的点，再删掉相对孤立的点，对剩下的点进行聚类分析——将剩下的点按照空间位置关系分成彼此不相连的子集。能量最低的点集作为最优预测位点。以得到的最优预测位点为第一个种子片段的位置，对种子的空间取向进行搜索来找到能量最低的方向，然后指定生长分子的长度以及想要得到的最终生长出来的分子数量等信息，在目标蛋白周围长出所有可能的分子；

分析长出来的分子的特征，通过已有生物功能数据库预测目标蛋白可能具有的功能。

本发明具有以下有益效果及优点：

1.本发明使用半粗粒化的统计势能函数，在半粗粒化的统计势能函数里，目标蛋白将使用全原子模型，而与之结合的蛋白质、DNA或RNA将使用粗粒化模型，这样统计出原子-残基(碱基)对在距离上的分布以计算相应的能量；采用半粗粒化统计势函数不仅减少计算量、具有已有粗粒化统计势函数计算速度快的优点，而且也将与分子碎片生长方法结合的更加紧密、具有目前已有粗粒化统计势所不具备的准确性。

附图说明

图1为本发明方法流程图；

图2为新分子形成过程图示意图。

具体实施方式

如图1所示，本发明计算机预测蛋白功能的方法包括以下步骤：

对生物大分子的原子之间的距离进行半粗粒化统计，得到生物大分子之间相互作用的半粗粒化统计势函数；

利用半粗粒化统计势函数在目标蛋白周围空间找到相互作用的热点区域；

结合分子碎片生长方法，在热点区域寻找可有与目标蛋白紧密结合的蛋白质、DNA或RNA序列片段；

计算蛋白质、DNA或RNA序列片段与目标蛋白结合的强度并排序；

选取规定数量的序列片段及其结构，在相关的生物功能数据库中进行搜索，查找相应的生物功能，预测到目标蛋白可能具有的功能。

自然界一共存在20种氨基酸，DNA链包括4种碱基，RNA也包括4种碱基。DNA与RNA不同的地方在于每个DNA的主链5碳糖结构上比RNA少一个氧原子，而且在DNA中有胸腺嘧啶而RNA中没有胸腺嘧啶但是有尿嘧啶。

在多肽链生长的过程中，碎片库由20种氨基酸残基构成。每种氨基酸残基的活性位点即为主链上羧基的C原子或者主链上的N原子。生成肽键的时候，通过重新设定某些原子的坐标，要将构成肽键的原子设置到一个平面上。每一个氨基酸中都有主链和侧链，大部分单键都是可以旋转的。构象的数量非常大。采取的办法是，主链上可以旋转的键采取均匀取值(比如10，20，30，40，...，360)，侧链中可以旋转的键直接到已经有的数据库(天然氨基酸构象出现几率的数据库)中搜索。选择的原则还是出现几率大的被选择的几率就大。种子一开始的方向决定开始了生长的方向。随机转动种子的方向，随机选择10000个方向，然后使用统计势能打分函数对所有的方向进行评价，保留分最高的一个。为了不漏掉结合能量最低位点相近的位点，采用多个位点放种子并列生长的办法，最后将所有种子的结构放在一起进行评价。然后按照分值进行排序，能量低的放在前面。

在RNA生长的过程，选择构成RNA的四个碱基作为碎片库基础。每种碱基的活性位点即为碱基上的P原子，新碎片的O3’会被接到旧分子的P原子上。主链的旋转按照转动的角度均匀取值，侧链转动的角度随机取值。生长失败的时候，可以适当调整生长点所在键的方向。

在DNA生长过程，选择构成DNA的四个碱基作为碎片库基础。每种碱基的活性位点即为碱基上的P原子，新碎片的O3’会被接到旧分子的P原子上。DNA是双螺旋的，所以不需要旋转。因为DNA的结构比较硬，不像蛋白质和RNA相对柔软。因此在生长的过程中，每个种子带上特定的方向都对应一个结果，将所有的结果综合在一起按照分值的排序来找到最可能的方向。

所述半粗粒化统计势函数通过以下过程得到：

读入蛋白质三维结构，按照相应原子在其所在氨基酸残基或DNA碱基或RNA碱基类型上的拓扑位置给定其原子类型

利用给定类型统计原子对之间的距离分布，通过与期望的分布值比较，得到半粗粒化统计势函数，通过以下公式将分布值转化成统计势函数：

$>\bar{u}(i,j,r)=\left(\begin{array}{cc}-\mathrm{RT}\ln \frac{P^c(i,j,r)}{\frac{\left(r^{\alpha }\mathrm{\Delta r}\right)f_V\left(r\right)}{{\Sigma }_r\left(r^{\alpha }\mathrm{\Delta r}\right)f_V\left(r\right)}},& r<r_{\mathrm{cut}},\\ 0,& r\ge r_{\mathrm{cut}},\end{array}\right)$>

本实施例中，读入蛋白质三维结构，以Protein Database Bank三维结构文件为例，按照相应原子在其所在氨基酸残基或DNA碱基或RNA碱基类型上的拓扑位置给定其原子类型。目标蛋白使用全原子模型，按照说明说表1给定其原子类型。预测蛋白、DNA、RNA使用粗粒化模型，其粗粒化原子类型为其侧链的几何中心。

对于需要预测功能的目标蛋白使用全原子结构模型，原子结构坐标数据主要从蛋白质结构数据库PDB获得，对于在PDB里不存在的目标蛋白结构，可以使用同源建模的方式获得。其相应原子类型如表1所示。

导出半粗粒化统计势函数：

对用于训练势函数的蛋白，遍历所有原子间距离在r_cut以内的原子对，把距离在(r，r+Δr]区间内的原子对按照其相应目标蛋白中的原子类型i和相应预测蛋白、DNA或RNA中的粗粒化原子类型j统计为N_obs(i，j，r)。通过它与我们期望的分布值比较，我们可以得到原子类型对i，j在(r，r+Δr]距离的能量值为：

$>\bar{u}(i,j,r)=\left(\begin{array}{cc}-\mathrm{RT}\ln \frac{P^c(i,j,r)}{\frac{\left(r^{\alpha }\mathrm{\Delta r}\right)f_V\left(r\right)}{{\Sigma }_r\left(r^{\alpha }\mathrm{\Delta r}\right)f_V\left(r\right)}},& r<r_{\mathrm{cut}},\\ 0,& r\ge r_{\mathrm{cut}},\end{array}\right)$>

其中，R为理想气体常数，T为温度取值300K，r为原子对之间的距离，Δr为统计区间宽度，r_cut为截断距离。α和N₀为经验参数分别取值1.61和75；把距离在区间内的原子对按照其相应目标蛋白中的原子类型i和相应预测蛋白、DNA或RNA中的粗粒化原子类型j统计为N_obs(i，j，r)。为N_obs(i，j，r)中原子类型对相应的两个原子一个是全原子模型并且另一个是粗粒化模型的部分。则中间变量

$>f_V\left(r\right)={\Sigma }_{i,j}{N}_{\mathrm{obs}}^{\mathrm{Protein}-\mathrm{CG}}(i,j,r)/{\Sigma }_{i,j}{N}_{\mathrm{obs}}(i,j,r),$>中间变量

$>P(i,j,r)=N_{\mathrm{obs}}^c(i,j,r)/{\Sigma }_r{N}_{\mathrm{obs}}^c(i,j,r),$>中间变量

$>N_{\mathrm{obs}}^c(i,j,r)=N_{\mathrm{obs}}(i,j,r)+N_0\frac{{\Sigma }_{i,j}{N}_{\mathrm{obs}}^{\mathrm{Protein}-\mathrm{CG}}(i,j,r)}{{\Sigma }_{i,j,r}{N}_{\mathrm{obs}}^{\mathrm{Protein}-\mathrm{CG}}(i,j,r)}$>

在得到了各种原子类型之间在某个距离区间的能量值以后，遍历所有原子间距离在r_cut以内的原子对，我们就可以得到整个结构的总能量：

(1)用于训练的蛋白由PDB获得，按照蛋白、DNA或RNA分为三类。每类都为PDB库里所有相应的结构，并去除35％同源性，其PDB ID列表如表2所示。

所述分子碎片生长方法为：

利用构成天然蛋白质的20种氨基酸，构成DNA的4种脱氧核苷酸，构成RNA的4种核苷酸为片段前体。对前体进行修改后形成生长使用的片段。整个过程按照进化选择方式进行控制。繁殖的方法是：父代的一个分子选择一个片段，然后连接片段，其次调整新连接片段的构象——每一个构象都对应了一个子代分子。然后通过计算每个子代分子和目标蛋白质之间的结合能量。按照带权随机法选择一批子代进行下一代的繁殖。

如图2所示，新的分子碎片(上右)半胱氨酸按照旧分子的取向(上左)接到旧分子的活性位点上，形成新的分子(下)。按照能量高低留下一定比例的新碎片，并在下一次生长中按照能量给予一定的几率选取。

所述选取规定数量的序列片段及其结构，在相关的生物功能数据库中进行搜索，查找相应的生物功能，预测到目标蛋白可能具有的功能包括以下步骤：

在目标蛋白周围以为边长，把空间划分为以立方体为结构单元的格点，以格点位置为假想各种粗粒化粒子的坐标，在每个格点上使用各种粗粒化格点类型并使用导出的半粗粒化统计势能函数评价格点和目标蛋白之间的能量，其中能量最低的值留做此格点的能量；

使用能量阀值删掉能量比较低的点，再删掉相对孤立的点，对剩下的点进行聚类分析——将剩下的点按照空间位置关系分成彼此不相连的子集。能量最低的点集作为最优预测位点。以得到的最优预测位点为第一个种子片段的位置，对种子的空间取向进行搜索来找到能量最低的方向，然后指定生长分子的长度以及想要得到的最终生长出来的分子数量等信息，在目标蛋白周围长出所有可能的分子；

分析长出来的分子的特征，通过已有生物功能数据库预测目标蛋白可能具有的功能。