CN101794351A - 一种基于大间隔最近中心点的蛋白质二级结构的工程预测方法 - Google Patents

Abstract

一种基于大间隔最近中心点的蛋白质二级结构的工程预测方法，属于蛋白质二级结构的工程预测方法领域，解决了现有蛋白质二级结构的预测方法在采用机器学习算法时，存在的数据权值存在局部极小值、预测效率低的问题。本发明所述的蛋白质二级预测方法是：首先基于PDB数据库构造非冗余蛋白质二级结构训练数据集，然后基于NCBI nr数据库为目标蛋白质链构造多序列比对特征，最后调用大间隔最近中心点算法来构建蛋白质二级结构预测模型。所述的大间隔最近中心点算法利用欧氏距离的K-means聚类算法确定每类样本的中心点，通过最小化目标损失函数学习输入空间的线性变换。本发明实现了快速、高效率、高精度的蛋白质二级结构预测，适用于蛋白质的二级结构预测。

Description

一种基于大间隔最近中心点的蛋白质二级结构的工程预测方法

技术领域

本发明涉及的是一种基于机器学习的蛋白质二级结构的工程预测方法，特别涉及一种将大间隔最近中心点方法与多序列比对特征相结合的方法，属于蛋白质二级结构的工程预测方法领域。

背景技术

随着人类基因组图谱的公布，以及越来越多动物、植物、微生物的基因组全序列测定的完成，生物科学进入了“后基因组时代”。人类将在了解遗传物质DNA全部序列的基础上研究和认识生命的奥秘，阐明基因编码产物(蛋白质)的功能己成为主要研究目标。一系列的研究表明，一种蛋白质能够行使其特定的生物功能，是由它特定的结构所决定的。因此，了解蛋白质的结构是对其功能获得彻底理解的先决条件。尽管采用X-ray晶体衍射和核磁共振这两种实验方法可以测定蛋白质的三维结构，但由于其都具有周期长、成本高、技术难度大等缺点，使得采用实验方法测定蛋白质结构的速度远远低于由大规模测序工程生成的蛋白质序列的速度。为了缩小这一差距，采用理论预测蛋白质结构的方法势在必行。

蛋白质的氨基酸序列决定了它的三级结构，但是直接从氨基酸序列成功地预测蛋白质三级结构就目前的技术水平而言仍然是一个艰难的任务。于是，本领域技术人员采用一个折中却非常有效的步骤就是预测蛋白质二级结构，即将蛋白质序列中的每个氨基酸残基归约为螺旋构象、折叠构象、或是卷曲构象。蛋白质二级结构预测提供的信息，不仅可用于蛋白质三级结构的从头预测、蛋白质的互作与功能预测、蛋白质的分析和注释，还可以提高蛋白质折叠识别的敏感性。

自20世纪60年代以来，已经有越来越多的科研人员致力于蛋白质二级结构预测的研究并提出了各种各样的预测算法。按照算法的基本思想，它们大体上可以分为三类：1.基于规则的预测，包括生物规则、统计规则和组合规则。具有代表性的算法有Chou-Fasman、Lim和GGBSM等。由于这类方法的预测精度普遍低于70％，目前已极少采用。2.基于图模型的预测。由于图模型方法，包括隐马尔可夫模型、条件随机场和贝叶斯网等，能够成功地处理序列数据问题，因此此类方法在蛋白质的二级结构预测中得到了广泛应用。图模型方法是从生成模型的角度考虑蛋白质二级结构预测问题的。尽管图模型能融入标签的相关性和长距离信息，但由于蛋白质的二级结构主要是由氨基酸残基间的局部相互作用决定的，图模型并不适用于蛋白质二级结构的建模，长期的实验表明它们并不能给出精确的蛋白质二级结构预测。3.基于机器学习的预测，比如神经网络和支持向量机。早期的机器学习方法主要是基于单序列特征，它们的三态每残基预测精度只有65％左右。直到1993，Rost和Sander把神经网络方法和多序列比对特征相结合，使得蛋白质二级结构预测精度突破了70％。随后，Jones通过引入PSSM特征进一步改进了蛋白质二级结构预测精度。目前，预测精度较好的蛋白质二级结构预测方法都是基于机器学习算法，其中基于集成神经网络的蛋白质二级结构预测方法的预测精度已经达到了80％左右。此外，基于支持向量机的预测方法也获得了较好的预测结果。

尽管神经网络或是支持向量机的方法都能获得较高的蛋白质二级结构预测精度，但是它们都有不可克服的缺点。神经网络具有权值的训练存在局部极小值的风险，训练结果的好坏显著地依赖初始权值的选择。支持向量机的缺点是训练蛋白质二级结构预测模型会产生许多支持向量，因此不能快速地预测蛋白质二级结构。

发明内容

本发明的目的是针对现有蛋白质二级结构的预测方法在采用机器学习算法时，存在的数据权值存在局部极小值、预测效率低的问题，提供了一种基于大间隔最近中心点的蛋白质二级结构的工程预测方法。

本发明是通过下述方案予以实现的：一种基于大间隔最近中心点的蛋白质二级结构的工程预测方法，采用下列步骤实现：

步骤一、下载发布的NCBI nr数据库和PDB格式的蛋白质结构数据，基于PDB格式的蛋白质结构数据构造非冗余蛋白质二级结构训练数据集；

步骤二、给定目标蛋白质一级序列数据，根据步骤一提供的NCBI nr数据库为目标蛋白质一级序列中的每个残基构造多序列比对特征向量；

步骤三、基于步骤二中构造的目标蛋白质序列的多序列比对特征向量，调用大间隔最近中心点算法，获得目标蛋白质的二级结构预测数据，

在步骤三中，所述的大间隔最近中心点算法是通过以下步骤实现的：

步骤三·一、基于步骤二中为非冗余蛋白质二级结构训练数据集中的所有残基构造的多序列比对特征向量，以残基对应的二级结构作为特征向量的标签构造大间隔最近中心点算法的训练样本集；

步骤三·二、基于步骤三·一构造的训练样本集，利用欧氏距离的K-means聚类算法确定各类样本的中心点，其中，螺旋类样本、卷曲类样本和折叠类样本对应的K值分别为3，3和2；

步骤三·三、基于步骤三·二确定的各类样本的中心点和给定的初始超参数μ，利用子梯度投影算法，通过最小化目标损失函数求解大间隔最近中心点模型的线性变换矩阵，其中，目标损失函数形式化为凸半定规划问题。

本发明基于确定的超参数μ、从PDB数据库导出的非冗余蛋白质二级结构训练集以及多序列比对特征向量构造方法，构造一种基于大间隔最近中心点的蛋白质二级结构预测系统。

本发明采用大间隔最近中心点算法模型可以形式化为凸半定规划问题，因此参数求解不存在局部极小值问题。通过引入正则化技术，大间隔最近中心点算法克服了过学习的问题。此外，本发明的预测模型仅仅基于8个中心点和3个线性变换矩阵，因此可以快速高效地实现蛋白质二级结构预测。并且8个中心点分别对应二级结构的8种结构定义，因此，采用本发明所述的方法能够很好地拟合蛋白质氨基酸序列与其空间构象的关系，实现了快速高效的对蛋白质的二级结构进行预测。

附图说明

图1是一种基于大间隔最近中心点的蛋白质二级结构的工程预测方法的流程图；图2是采用大间隔最近中心点算法流程图；图3是构造非冗余蛋白质二级结构训练数据集的方法流程图；图4是为目标蛋白质一级序列中的每个残基构造多序列比对特征向量的方法的流程图。

具体实施方式

具体实施方式一：下面结合图1、图2具体说明本实施方式。一种基于大间隔最近中心点的蛋白质二级结构的工程预测方法，采用下列步骤实现：

步骤一、下载发布的NCBI nr数据库和PDB格式的蛋白质结构数据，基于PDB格式的蛋白质结构数据构造非冗余蛋白质二级结构训练数据集；

步骤二、给定目标蛋白质一级序列数据，根据步骤一提供的NCBI nr数据库为目标蛋白质一级序列中的每个残基构造多序列比对特征向量；

步骤三、基于步骤二中构造的目标蛋白质序列的多序列比对特征向量，调用大间隔最近中心点算法，获得目标蛋白质的二级结构预测数据，

在步骤三中，所述的大间隔最近中心点算法是通过以下步骤实现的：

步骤三·一、基于步骤二中为非冗余蛋白质二级结构训练数据集中的所有残基构造的多序列比对特征向量，以残基对应的二级结构作为特征向量的标签构造大间隔最近中心点算法的训练样本集；

步骤三·二、基于步骤三·一构造的训练样本集，利用欧氏距离的K-means聚类算法确定各类样本的中心点，其中，螺旋类样本、卷曲类样本和折叠类样本对应的K值分别为3，3和2；

步骤三·三、基于步骤三·二确定的各类样本的中心点和给定的初始超参数μ，利用子梯度投影算法，通过最小化目标损失函数求解大间隔最近中心点模型的线性变换矩阵，其中，目标损失函数形式化为凸半定规划问题。

本发明基于确定的超参数μ、从PDB数据库导出的非冗余蛋白质二级结构训练集以及多序列比对特征向量构造方法，构造一种基于大间隔最近中心点的蛋白质二级结构预测系统。

本发明采用大间隔最近中心点算法模型可以形式化为凸半定规划问题，因此参数求解不存在局部极小值问题。通过引入正则化技术，大间隔最近中心点算法克服了过学习的问题。此外，本发明的预测模型仅仅基于8个中心点和3个线性变换矩阵，因此可以快速高效地实现蛋白质二级结构预测。并且8个中心点分别对应二级结构的8种结构定义，因此，采用本发明所述的方法能够很好地拟合蛋白质氨基酸序列与其空间构象的关系，实现了快速高效的对蛋白质的二级结构进行预测。

本发明构建的基于大间隔最近中心点算法模型主要特征包括两方面：

一、鉴于蛋白质链氨基端和羧基端的二级结构形成方式与中间的形成方式不同，使用两个中心点模型：一个模拟蛋白质链两端的二级结构形成，一个模型蛋白质链中间的二级结构形成。

二、基于DSSP定义8种二级结构，构造8个中心点分别对应二级结构的8种结构定义，其中，螺旋构象中和卷曲构象中分别包括3个中心点，折叠构象中包括2个中心点。

本发明构建的基于大间隔最近中心点算法的功能：

基于构建的大间隔最近中心点算法模型和目标蛋白质一级序列，预测算法把距离目标蛋白质一级序列中残基的特征向量最近的中心点对应的蛋白质二级结构指派为目标蛋白质一级序列中残基的二级结构，作为预测的蛋白质二级结构输出。

步骤三·一中所述的大间隔最近中心点是采用K-means聚类算法进行确定的。

将步骤二中获得的目标蛋白质一级序列中的每个残基构造多序列比对特征向量与它所属的构象类别构造为训练数据集T，所述的训练数据集T表现形式记{(x₁，y₁)，(x₂，y₂)，...，(x_N，y_N)}，其中x_i是特征向量，y_i为类别标签，是从1到C的自然数，N是训练数据集中样本的数目，i是从1到N的自然数，(x_i，y_i)为训练样本点。

基于类别标签y_i把训练数据集划分为C个子集，采用欧氏距离的K-means聚类算法分别对每个子集进行聚类，聚类后的聚类中心点作为对应类别的中心点。对于每个聚类，把它的中心点称为其包含样本的目标中心点。一个样本的目标中心点，就是在度量学习过程中应该与其距离最近的中心点，目标中心点是在度量学习之前确定的并在学习过程中保持不变。

采用K-means聚类算法来确定中心点，计算效率高、易操作。

步骤三·二中所述的大间隔最近中心点模型是指基于步骤三·一中所述的大间隔最近中心点，学习一个线性变换矩阵L而实现的。

训练数据集T＝{(x₁，y₁)，(x₂，y₂)，...，(x_N，y_N)}中通过分别对每类的所有中心点从1进行依次编号，可以用m_jk唯一表示一个中心点，其中j∈{1，2，...，C}是中心点对应的类别，k∈{1，2，...，n_j}是中心点的编号，n_j表示第j类样本聚类后的中心点数目。对于每个聚类，把它的中心点称为其包含样本的目标中心点。一个样本的目标中心点，就是在度量学习过程中应该与其距离最近的中心点，目标中心点是在度量学习之前确定的并在学习过程中保持不变。为了表示目标中心点的信息，为每个训练样本x_i增加一个代表其目标中心点编号的标签t_i。显然，

就是x_i的目标中心点。

对于每个训练样本点(x_i，y_i)，L满足公式一的约束：对于任意一个不等于y_i的j，k为任意值时，

公式一 ${\left|\right|L(x_i-m_{y_i{t}_i})\left|\right|}_2^2+1<{\left|\right|L(x_i-m_{\mathrm{jk}})\left|\right|}_2^2$

即每个样本点与其目标中心点的距离和它与其它类别中心点的距离应至少保持一个单位间隔。维持一个单位间隔的目的是为了增加大间隔最近中心点分类的泛化能力。为了训练大间隔最近中心点模型，引入一个目标损失函数ε(L)，所述的目标损失函数ε(L)包括两项：第一项用来惩罚公式一的边界违背，第二项用来正则化线性变换矩阵L。其中，线性变换矩阵L使得目标损失函数ε(L)最小化：

$\mathrm{\&epsiv;}\left(L\right)=\underset{i}{\mathrm{\&Sigma;}}\underset{j\&NotEqual;y_i}{\mathrm{\&Sigma;}}\underset{k=1}{\overset{n_j}{\mathrm{\&Sigma;}}}{[1+{\left|\right|L(x_i-m_{y_i{t}_i})\left|\right|}_2^2-{\left|\right|L(x_i-m_{\mathrm{jk}})\left|\right|}_2^2]}_{+}+\mathrm{\&mu;}\left(\mathrm{tr}\left(L^{T}L\right)\right)$

$=\underset{i}{\mathrm{\&Sigma;}}\underset{j\&NotEqual;y_i}{\mathrm{\&Sigma;}}\underset{k=1}{\overset{n_j}{\mathrm{\&Sigma;}}}{[1+{(x_i-m_{y_i{t}_i})}^T{L}^{T}L(x_i-m_{y_i{t}_i})-{(x_i-m_{\mathrm{jk}})}^T{L}^{T}L(x_i-m_{\mathrm{jk}})]}_{+}+\mathrm{\&mu;}\left(\mathrm{tr}\left(L^{T}L\right)\right)$

其中，函数[z]₊＝max(z，0)表示铰链损失。当公式一中的不等式对任意样本x_i都满足时，所有铰链损失[z]₊的值都为0。此时，目标损失函数ε(L)达到最小值。

因为由线性变换矩阵L表示的目标损失函数ε(L)不是凸的，所以在使用梯度法求解时易陷入局部极小值。为了克服这一问题，引入一个新的矩阵变量M＝L^TL。可以看出矩阵M是半正定矩阵。通过用M替换L，可以把目标损失函数ε(L)表示为：

$\mathrm{\&epsiv;}\left(M\right)=\underset{i}{\mathrm{\&Sigma;}}\underset{j\&NotEqual;y_i}{\mathrm{\&Sigma;}}\underset{k=1}{\overset{n_j}{\mathrm{\&Sigma;}}}{[1+{(x_i-m_{y_i{t}_i})}^{T}M(x_i-m_{y_i{t}_i})-{(x_i-m_{\mathrm{jk}})}^{T}M(x_i-m_{\mathrm{jk}})]}_{+}+\mathrm{\&mu;}\left(\mathrm{tr}\left(M\right)\right)$

上式中的目标损失函数ε(L)，是关于矩阵M中元素的分段线性凸函数。为将所述的目标损失函数ε(L)形式化为半正定规划问题，对其进行标准化，定义非负松弛变量{ξ_ijk}来模拟上式中所有铰链损失[z]₊的影响。M是半正定矩阵，因此，将目标损失函数ε(M)的最小化形式化为凸半正定规划：

最小化： $\underset{i}{\mathrm{\&Sigma;}}\underset{j\&NotEqual;y_i}{\mathrm{\&Sigma;}}\underset{k=1}{\overset{n_j}{\mathrm{\&Sigma;}}}{\mathrm{\&xi;}}_{\mathrm{ijk}}+\mathrm{\&mu;}\left(\mathrm{tr}\left(M\right)\right)$

制约条件为：

$\left\{\begin{array}{c}{(x_i-m_{\mathrm{jk}})}^{T}M(x_i-m_{\mathrm{jk}})-{(x_i-m_{y_i{t}_i})}^{T}M(x_i-m_{y_i{t}_i})\&GreaterEqual;1-{\mathrm{\&xi;}}_{\mathrm{ijk}}\\ {\mathrm{\&xi;}}_{\mathrm{ijk}}\&GreaterEqual;0\\ M\underset{\&OverBar;}{f}0\end{array}\right.$

对于一些非线性多类决策问题，特征空间的一个全局线性变换可能不足以充分地改进最近中心点分类方法的预测精度。在这种情形下，对每一个样本类别分别学习一个线性变换。

多度量大间隔最近中心点分类模型，尝试学习C个线性变换矩阵L_j，理想情况下，每个样本点(x_i，y_i)，满足如下条件：

公式二 ${\left|\right|L_{y_i}(x_i-m_{y_i{t}_i})\left|\right|}_2^2+1<{\left|\right|L_j(x_i-m_{\mathrm{jk}})\left|\right|}_2^2$

其中，j为不等于y_i的从1到C的自然数，k为任意值，与公式一不同之处在于，公式二中样本与中心点的距离还依赖于中心点所对应的类别。为了获得满足条件的C个线性变换矩阵，定义目标损失函数ε(L₁，...，L_C)：

$\mathrm{\&epsiv;}(L_1,...,L_C)=\underset{i}{\mathrm{\&Sigma;}}\underset{j\&NotEqual;y_i}{\mathrm{\&Sigma;}}\underset{k=1}{\overset{n_j}{\mathrm{\&Sigma;}}}{[1+{\left|\right|L_{y_i}(x_i-m_{y_i{t}_i})\left|\right|}_2^2-{\left|\right|L_j(x_i-m_{\mathrm{jk}})\left|\right|}_2^2]}_{+}+\mathrm{\&mu;}\underset{j=1}{\overset{C}{\mathrm{\&Sigma;}}}\mathrm{tr}\left({L^T}_j{L}_j\right)$

最小化目标损失函数ε(L₁，...，L_C)，定义C个半正定矩阵M_j＝L_j ^TL_j，其中j为从1到C的自然数，定义松弛变量{ξ_ijk}，将目标损失函数ε(M_j)的最小化形式化为凸半正定规划：

最小化： $\underset{i}{\mathrm{\&Sigma;}}\underset{j\&NotEqual;y_i}{\mathrm{\&Sigma;}}\underset{k=1}{\overset{n_j}{\mathrm{\&Sigma;}}}{\mathrm{\&xi;}}_{\mathrm{ijk}}+\mathrm{\&mu;}\underset{j=1}{\overset{C}{\mathrm{\&Sigma;}}}\mathrm{tr}\left(M_j\right)$

制约条件为：

$\left\{\begin{array}{c}{(x_i-m_{\mathrm{jk}})}^T{M}_j(x_i-m_{\mathrm{jk}})-{(x_i-m_{y_i{t}_i})}^T{M}_{y_i}(x_i-m_{y_i{t}_i})\&GreaterEqual;1-{\mathrm{\&xi;}}_{\mathrm{ijk}}\\ {\mathrm{\&xi;}}_{\mathrm{ijk}}\&GreaterEqual;0\\ M_j\underset{\&OverBar;}{f}0j=1,...,C\end{array}\right.$

同理，C个线性变换矩阵可以用子梯度投影算法快速地求解，并且不存在局部极小值问题。

线性变换矩阵可以用快速的子梯度算法求解，具体过程如下：

在第t次迭代中，令M_t-1为迭代开始时的半正定矩阵，则此时的目标损失函数ε(M_j)为：

$\mathrm{\&epsiv;}\left(M_{t-1}\right)=\underset{i}{\mathrm{\&Sigma;}}\underset{j\&NotEqual;y_i}{\mathrm{\&Sigma;}}\underset{k=1}{\overset{n_j}{\mathrm{\&Sigma;}}}{[1+{(x_i-m_{y_i{t}_i})}^T{M}_{t-1}(x_i-m_{y_i{t}_i})-{(x_i-m_{\mathrm{jk}})}^T{M}_{t-1}(x_i-m_{\mathrm{jk}})]}_{+}+\mathrm{\&mu;}\left(\mathrm{tr}\left(M_{t-1}\right)\right)$

上式中，由于M_t-1是分段线性的。定义一个三元组集ψ^t，当i、j、k在ψ^t范围内触发铰链损失[z]₊，即ξ_ijk大于0时，可以计算目标损失函数ε(M_t-1)的梯度G_t：

$G_t=\underset{(i,j,k)\&Element;{\mathrm{\&psi;}}^t}{\mathrm{\&Sigma;}}[(x_i-m_{y_i{t}_i}){(x_i-m_{y_i{t}_i})}^T-(x_i-m_{\mathrm{jk}}){(x_i-m_{\mathrm{jk}})}^T]+\mathrm{\&mu;I}$

其中，I为单位矩阵。梯度G_t仅依赖于三元组集ψ^t。因此，连续两次迭代梯度的改变仅由ψ^t与ψ^t+1之差决定。因此，基于第t次迭代的梯度G_t快速地计算第t+1次迭代的梯度G_t+1：

$G_{t+1}=G_t+\underset{(i,j,k)\&Element;{\mathrm{\&psi;}}^{t+1}-{\mathrm{\&psi;}}^t}{\mathrm{\&Sigma;}}[(x_i-m_{y_i{t}_i}){(x_i-m_{y_i{t}_i})}^T-(x_i-m_{\mathrm{jk}}){(x_i-m_{\mathrm{jk}})}^T]$

$-\underset{(i,j,k)\&Element;{\mathrm{\&psi;}}^t-{\mathrm{\&psi;}}^{t+1}}{\mathrm{\&Sigma;}}[(x_i-m_{y_i{t}_i}){(x_i-m_{y_i{t}_i})}^T-(x_i-m_{\mathrm{jk}}){(x_i-m_{\mathrm{jk}})}^T]$

对于小的梯度步长，三元组集ψ^t在连续两次迭代中的改变很小。因此，上式中的梯度可以以极快的速度进行计算。

完成梯度G_t的计算后，目标损失函数ε(M_j)应沿着子梯度方向下降一步。即根据梯度步长α，参数矩阵M_t-1应作如下更新：

M′_t＝M_t-1-αG_t

为了保证更新后的参数矩阵是半正定矩阵，将M′_t对角化。用M′_t＝PΛP^T表示M′_t的特征分解，其中P是特征向量的正交矩阵，Λ是以对应的特征值为对角元素的对角矩阵。将Λ中所有负特征值都变为0，可以得到一个新的对角矩阵Λ⁺，Λ⁺＝max(Λ，0)。则M′_t到半正定锥投影为M_t，M_t＝PΛ⁺P^T。

在实际求解过程中，用子梯度投影算法直接对其进行最小化运算，当梯度步长足够小时，子梯度投影算法能够收敛到正确解。因此，大间隔最近中心点算法不会陷入局部极小值。根据收敛时的半正定矩阵M_t＝PΛ⁺P^T，可以获得线性变换矩阵L：L＝P^T(Λ⁺)^1/2。

具体实施方式二：本实施方式是对具体实施方式一中所述的一种基于大间隔最近中心点的蛋白质二级结构的工程预测方法的进一步说明，步骤三·三中所述的初始超参数μ取值范围为0、0.1、1、5、10或20，超参数μ在所述范围内的最优值是采用RS126非冗余数据集来快速确定的。

由于从步骤一所述的PDB数据库导出的PDB数据训练集包含相当多的蛋白质链，所述的PDB数据训练集的子梯度投影算法需要较长的时间才能收敛。因此，采用RS126非冗余数据集来快速地确定超参数μ，本实施方式中所述的超参数μ是被用于正则化线性变换矩阵。选取适合的超参数μ可以防止过学习，避免学习的模型过分拟合训练数据。

具体实施方式三：下面结合图3具体说明本实施方式。本实施方式是对具体实施方式一所述的一种基于大间隔最近中心点的蛋白质二级结构的工程预测方法的进一步限定，在步骤一中，构造非冗余蛋白质二级结构训练数据集是采用下列步骤实现的：

步骤一·一、基于PDB数据库中发布的由X-ray晶体衍射测定的PDB格式的蛋白质结构数据，应用DSSP程序将所述的PDB格式的蛋白质结构数据转化为DSSP格式的数据文件；

步骤一·二、基于DSSP格式的定义把DSSP格式的数据文件转化为FASTA格式的蛋白质序列数据文件。同时，将DSSP定义的8种二级结构归约为3类，其中，将H构象、G构象、I构象归约为H构象，即螺旋类构象，将B构象、E构象归约为E构象，即折叠类构象，将T构象、S构象、-构象归约为C构象，即卷曲类构象；

步骤一·三、调用序列一致性阈值为0.4的CD-HIT算法，将步骤一·二中得到的FASTA格式的蛋白质序列数据文件中冗余性蛋白质序列数据文件删除；

步骤一·四、调用Needleman-Wunsch算法对步骤一·三中所述的删除冗余性蛋白质序列数据文件后的蛋白质序列数据文件两两进行比对，对序列一致性大于30％的两个蛋白质链，只保留其中的一个，构造非冗余蛋白质二级结构训练数据集。

本实施方式中，所述的DSSP定义的8种二级结构是根据氢键模式来指定的，一般是用单一英文字母来描述。其中，H是alpha helix简称、G是3₁₀helix的简称、I是pi helix的简称、E是beta strand的简称、B是isolated beta-bridge的简称、T是turn的简称、S是bend的简称、一是rest的简称。

蛋白质是由20种氨基酸组成的，然而由于许多可能的氨基酸片段组成从不出现在蛋白质氨基酸的序列中，因此，PDB数据库中的蛋白质氨基酸序列含有高度的冗余性。而且，许多蛋白质链间的序列一致性可以达到80％以上。由于蛋白质的氨基酸序列决定了蛋白质二级结构序列，因此，蛋白质二级结构训练数据集应包含多样性的非冗余氨基酸序列。

具体实施方式四：下面结合图4具体说明本实施方式。本实施方式是对实施方式一所述的一种基于大间隔最近中心点的蛋白质二级结构的工程预测方法的进一步限定，在步骤二中，基于NCBI nr数据库，为目标蛋白质一级序列中的每个残基构造多序列比对特征向量是采用下列步骤实现的：

步骤二·一、调用pfilt程序，对步骤一中得到的NCBI nr数据库中的FASTA格式的蛋白质序列数据进行过滤；

步骤二·二、调用formatdb程序，将步骤二·一中得到的过滤后的蛋白质序列数据进行格式化，生成用于PSI-BLAST多序列比对的索引和数据文件；

步骤二·三、调用3次迭代的blastpgp程序，将目标蛋白质一级序列与步骤二·二中得到的格式化后的NCBI nr数据库中的数据文件进行PSI-BLAST多序列比对，生成一个二进制检查点文件；

步骤二·四、调用makemat程序，将二·三中得到的二进制检查点文件进行处理，获得一个N×20的矩阵，将所述的N×20的矩阵中所有元素除以1000，获得的矩阵中的行向量对应为目标蛋白质一级序列的残基，其中，N为目标蛋白质一级序列的长度；

步骤二·五、用长为13个残基的滑动窗口扫描目标蛋白质一级序列，将所述的窗口中每个残基在矩阵中对应的行向量依次结合，为目标蛋白质一级序列中的每个残基构造一个由260个元素构成的特征向量。

采用本实施方式所述的方法可以实现多序列比对特征的构造。

在步骤二·一中，调用pfilt程序对NCBI nr数据库中的FASTA格式的蛋白质序列数据进行过滤，能够删除低复杂性区域、跨膜区域和卷曲-卷曲片段。

在步骤二·五中，对于滑动窗口中没有对应残基的位置，用全零的元素表示其对应的行向量。滑动窗口中心残基的二级结构就是其对应的多序列比对特征的标签，因此，给定蛋白质一级序列中的每个残基都对应一个多序列比对特征向量。

具体实施方式五：本实施方式是对实施方式一所述的一种基于大间隔最近中心点的蛋白质二级结构的工程预测方法的进一步限定，在步骤三中还包括以下步骤：根据已知的蛋白质二级结构数据，采用评估准则评判步骤三中得到的蛋白质二级结构预测数据的预测可信度。

具体实施方式六：本实施方式是对实施方式五所述的一种基于大间隔最近中心点的蛋白质二级结构的工程预测方法的进一步限定，所述的评估准则是采用三态单残基精度计算的方式实现的，

三态单残基精度Q₃可表示为：

$Q_3=\frac{V}{U}\&times;100$

其中，V是待预测残基序列中被正确预测的残基数目，U是待预测残基序列的残基数目，

待预测序列属于H构象、E构象或C构象，单类型的三态单残基精度分别用Q_H、Q_E和Q_V表示，

$Q_i=\frac{V_i}{U_i}\&times;100$

其中，V_i是待预测序列中被正确预测的处于i构象的残基数目，U_i是待预测序列中处于i构象的残基数目，i属于H构象、E构象或C构象。

本发明采用三态单残基精度Q₃作为评价准则进行计算，蛋白质二级结构的工程预测方法的预测精度可以达到了80％以上。

具体实施方式七：本实施方式是对实施方式五所述的一种基于大间隔最近中心点的蛋白质二级结构的工程预测方法的进一步限定，所述的评估准则是采用Matthews相关系数计算的方式实现的，

Matthews相关系数为：

$R_i=\frac{p_i{n}_i-u_i{o}_i}{\sqrt{(p_i+u_i)(p_i+o_i)(n_i+u_i)(n_i+o_i)}}$

其中，p_i是实际处于i构象且被正确预测为i构象的氨基酸残基数目，n_i是实际不处于i构象且被正确预测为非i构象的氨基酸残基数目，u_i是实际处于i构象但未被预测为i构象的氨基酸残基数目，o_i是实际不处于i构象但被错误地预测为i构象的氨基酸残基数目，i属于H构象、E构象或C构象。

本实施方式中，对于任何一个给定的待评估预测结果，p_i与u_i之和、n_i与o_i之和均为常数，并且R_i的取值范围是从-1到1，R_i越大，预测的结果越好。当所有的构象都被正确预测时，u_i和o_i为0，R_i为1。

具体实施方式八：本实施方式是对实施方式五所述的一种基于大间隔最近中心点的蛋白质二级结构的工程预测方法的进一步限定，所述的评估准则是采用片段重叠度量的方式实现的，

其中，s₁为目标蛋白质的二级结构序列上的片段，s₂为预测蛋白质的二级结构序列上的片段，(s₁，s₂)为一个重叠的二级结构片段对，S(i)为处于i构象的重叠二级结构片段对的集合，S′(i)为处于i构象且不和任意s₂片段重叠的s₁片段的集合，i属于H构象、E构象或C构象，

${\mathrm{Sov}}_i=100\&times;\frac{1}{T\left(i\right)}\underset{S\left(i\right)}{\mathrm{\&Sigma;}}[\frac{\min \mathrm{ov}(s_1,s_2)+\mathrm{\&delta;}(s_1,s_2)}{\max \mathrm{ov}(s_1,s_2)}\&times;\mathrm{len}\left(s_1\right)]$

$T\left(i\right)=\underset{s\left(i\right)}{\mathrm{\&Sigma;}}\mathrm{len}\left(s_1\right)+\underset{s^{\&prime;}\left(i\right)}{\mathrm{\&Sigma;}}\mathrm{len}\left(s_1\right)$

δ(s₁，s₂)＝min{(maxov(s₁，s₂)-minov(s₁，s₂))；minov(s₁，s₂)；int(len(s₁)/2)；int(len(s₂)/2)}

其中，Sov_i为i构象的片段重叠度量，len(s₁)为片段s₁上的氨基酸残基数，minov(s₁，s2)为片段s₁与片段s₂的实际重叠长度，maxov(s₁，s₂)为片段s₁与片段s₂的总扩展长度，T(i)为归一化因子，

$\mathrm{Sov}=100\&times;[\frac{1}{T}\underset{i\&Element;\{H,E,C\}}{\mathrm{\&Sigma;}}\underset{S\left(i\right)}{\mathrm{\&Sigma;}}\frac{\min \mathrm{ov}(s_1,s_2)+\mathrm{\&delta;}(s_1,s_2)}{\max \mathrm{ov}(s_1,s_2)}\&times;\mathrm{len}\left(s_1\right)]$

其中，

Sov为总体三态的片段重叠度量。

具体实施方式六至具体实施方式八所述的三种蛋白质二级结构评估方法是目前国际上比较通用的评估准则，基于已知的蛋白质二级结构数据，采用评估准则可以让用户了解本发明的预测可信度。

Claims (8)

1.一种基于大间隔最近中心点的蛋白质二级结构的工程预测方法，其特征是：

采用下列步骤实现：

步骤一、下载发布的NCBI nr数据库和PDB格式的蛋白质结构数据，基于PDB格式的蛋白质结构数据构造非冗余蛋白质二级结构训练数据集；

步骤二、给定目标蛋白质一级序列数据，根据步骤一提供的NCBI nr数据库为目标蛋白质一级序列中的每个残基构造多序列比对特征向量；

步骤三、基于步骤二中构造的目标蛋白质序列的多序列比对特征向量，调用大间隔最近中心点算法，获得目标蛋白质的二级结构预测数据，

在步骤三中，所述的大间隔最近中心点算法是通过以下步骤实现的：

步骤三·一、基于步骤二中为非冗余蛋白质二级结构训练数据集中的所有残基构造的多序列比对特征向量，以残基对应的二级结构作为特征向量的标签构造大间隔最近中心点算法的训练样本集；

步骤三·二、基于步骤三·一构造的训练样本集，利用欧氏距离的K-means聚类算法确定各类样本的中心点，其中，螺旋类样本、卷曲类样本和折叠类样本对应的K值分别为3，3和2；

步骤三·三、基于步骤三·二确定的各类样本的中心点和给定的初始超参数μ，利用子梯度投影算法，通过最小化目标损失函数求解大间隔最近中心点模型的线性变换矩阵，其中，目标损失函数形式化为凸半定规划问题。

2.根据权利要求1所述的一种基于大间隔最近中心点的蛋白质二级结构的工程预测方法，其特征是：步骤三·三中所述的初始超参数μ取值范围为0、0.1、1、5、10或20，超参数μ在所述范围内的最优值是采用RS126非冗余数据集来快速确定的。

3.根据权利要求1所述的一种基于大间隔最近中心点的蛋白质二级结构的工程预测方法，其特征是：

在步骤一中，构造非冗余蛋白质二级结构训练数据集是采用下列步骤实现的：

步骤一·一、基于PDB数据库中发布的由X-ray晶体衍射测定的PDB格式的蛋白质结构数据，应用DSSP程序将所述的PDB格式的蛋白质结构数据转化为DSSP格式的数据文件；

步骤一·二、基于DSSP格式的定义把DSSP格式的数据文件转化为FASTA格式的蛋白质序列数据文件；同时，将DSSP定义的8种二级结构归约为3类，其中，将H构象、G构象、I构象归约为H构象，即螺旋类构象，将B构象、E构象归约为E构象，即折叠类构象，将T构象、S构象、-构象归约为C构象，即卷曲类构象；

步骤一·三、调用序列一致性阈值为0.4的CD-HIT算法，将步骤一·二中得到的FASTA格式的蛋白质序列数据文件中冗余性蛋白质序列数据文件删除；

步骤一·四、调用Needleman-Wunsch算法对步骤一·三中所述的删除冗余性蛋白质序列数据文件后的蛋白质序列数据文件两两进行比对，对序列一致性大于30％的两个蛋白质链，只保留其中的一个，构造非冗余蛋白质二级结构训练数据集。

4.根据权利要求1所述的一种基于大间隔最近中心点的蛋白质二级结构的工程预测方法，其特征是：

在步骤二中，基于NCBI nr数据库，为目标蛋白质一级序列中的每个残基构造多序列比对特征向量是采用下列步骤实现的：

步骤二·一、调用pfilt程序，对步骤一中得到的NCBI nr数据库中的FASTA格式的蛋白质序列数据进行过滤；

步骤二·二、调用formatdb程序，将步骤二·一中得到的过滤后的蛋白质序列数据进行格式化，生成用于PSI-BLAST多序列比对的索引和数据文件；

步骤二·三、调用3次迭代的blastpgp程序，将目标蛋白质一级序列与步骤二·二中得到的格式化后的NCBI nr数据库中的数据文件进行PSI-BLAST多序列比对，生成一个二进制检查点文件；

步骤二·四、调用makemat程序，将二·三中得到的二进制检查点文件进行处理，获得一个N×20的矩阵，将所述的N×20的矩阵中所有元素除以1000，获得的矩阵中的行向量对应为目标蛋白质一级序列的残基，其中，N为目标蛋白质一级序列的长度；

步骤二·五、用长为13个残基的滑动窗口扫描目标蛋白质一级序列，将所述的窗口中每个残基在矩阵中对应的行向量依次结合，为目标蛋白质一级序列中的每个残基构造一个由260个元素构成的特征向量。

5.根据权利要求1所述的一种基于大间隔最近中心点的蛋白质二级结构的工程预测方法，其特征是：在步骤三中还包括以下步骤：根据已知的蛋白质二级结构数据，采用评估准则评判步骤三中得到的蛋白质二级结构预测数据的预测可信度。

6.根据权利要求5所述的一种基于大间隔最近中心点的蛋白质二级结构的工程预测方法，其特征是：所述的评估准则是采用三态单残基精度计算的方式实现的，

三态单残基精度Q₃可表示为：

$Q_3=\frac{V}{U}\&times;100$

其中，V是待预测残基序列中被正确预测的残基数目，U是待预测残基序列的残基数目，

待预测序列属于H构象、E构象或C构象，单类型的三态单残基精度分别用Q_H、Q_E和Q_V表示，

$Q_i=\frac{V_i}{U_i}\&times;100$

其中，V_i是待预测序列中被正确预测的处于i构象的残基数目，U_i是待预测序列中处于i构象的残基数目，i属于H构象、E构象或C构象。

7.根据权利要求5所述的一种基于大间隔最近中心点的蛋白质二级结构的工程预测方法，其特征是：所述的评估准则是采用Matthews相关系数计算的方式实现的，

Matthews相关系数为：

$R_i=\frac{p_i{n}_i-u_i{o}_i}{\sqrt{(p_i+u_i)(p_i+o_i)(n_i+u_i)(n_i+o_i)}}$

其中，p_i是实际处于i构象且被正确预测为i构象的氨基酸残基数目，n_i是实际不处于i构象且被正确预测为非i构象的氨基酸残基数目，u_i是实际处于i构象但未被预测为i构象的氨基酸残基数目，o_i是实际不处于i构象但被错误地预测为i构象的氨基酸残基数目，i属于H构象、E构象或C构象。

8.根据权利要求5所述的一种基于大间隔最近中心点的蛋白质二级结构的工程预测方法，其特征是：所述的评估准则是采用片段重叠度量的方式实现的，

其中，s₁为目标蛋白质的二级结构序列上的片段，s₂为预测蛋白质的二级结构序列上的片段，(s₁，s₂)为一个重叠的二级结构片段对，S(i)为处于i构象的重叠二级结构片段对的集合，S′(i)为处于i构象且不和任意s₂片段重叠的s₁片段的集合，i属于H构象、E构象或C构象，

${\mathrm{Sov}}_i=100\&times;\frac{1}{T\left(i\right)}\underset{S\left(i\right)}{\mathrm{\&Sigma;}}[\frac{\min \mathrm{ov}(s_1,s_2)+\mathrm{\&delta;}(s_1,s_2)}{\max \mathrm{ov}(s_1,s_2)}\&times;\mathrm{len}\left(s_1\right)]$

$T\left(i\right)=\underset{s\left(i\right)}{\mathrm{\&Sigma;}}\mathrm{len}\left(s_1\right)+\underset{s^{\&prime;}\left(i\right)}{\mathrm{\&Sigma;}}\mathrm{len}\left(s_1\right)$

δ(s₁，s₂)＝min{(maxov(s₁，s₂)-minov(s₁，s₂))；minov(s₁，s₂)；int(len(s₁)/2)；int(len(s₂)/2)}

其中，Sov_i为i构象的片段重叠度量，len(s₁)为片段s₁上的氨基酸残基数，minov(s₁，s2)为片段s₁与片段s₂的实际重叠长度，maxov(s₁，s₂)为片段s₁与片段s₂的总扩展长度，T(i)为归一化因子，

$\mathrm{Sov}=100\&times;[\frac{1}{T}\underset{i\&Element;\{H,E,C\}}{\mathrm{\&Sigma;}}\underset{S\left(i\right)}{\mathrm{\&Sigma;}}\frac{\min \mathrm{ov}(s_1,s_2)+\mathrm{\&delta;}(s_1,s_2)}{\max \mathrm{ov}(s_1,s_2)}\&times;\mathrm{len}\left(s_1\right)]$

其中，Sov为总体三态的片段重叠度量。

Images