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

权利要求书

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

采用下列步骤实现：

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

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

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

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

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

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

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

所述的大间隔最近中心点，学习一个线性变换矩阵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，y_i)，L满足公式一的约束；对于任意一个不等于y_i的j，k为任 意值时，

公式一||L(xi-myiti)||22+1<||L(xi-mjk)||22]]>

即每个样本点与其目标中心点的距离和它与其它类别中心点的距离应至少保持一个单 位间隔，目标损失函数ε(L)包括两项：第一项用来惩罚公式一的边界违背，第二项用来正 则化线性变换矩阵L，其中，线性变换矩阵L使得目标损失函数ε(L)最小化：

ϵ(L)=ΣiΣj≠yiΣk=1nj[1+||L(xi-myiti)||22-||L(xi-mjk)||22]++μ(tr(LTL))]]>

=ΣiΣj≠yiΣk=1nj[1+(xi-myiti)TLTL(xi-myiti)-(xi-mjk)TLTL(xi-mjk)]++μ(tr(LTL))]]>

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

引入一个矩阵变量M＝L^TL，M是半正定矩阵，通过用M替换L，可以把目标损失函 数ε(L)表示为：

ϵ(M)=ΣiΣj≠yiΣk=1nj[1+(xi-myiti)TM(xi-myiti)-(xi-mjk)TM(xi-mjk)]++μ(tr(M))]]>

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

最小化：ΣiΣj≠yiΣk=1njξijk+μ(tr(M))]]>

制约条件为：

(xi-mjk)TM(xi-mjk)-(xi-myiti)TM(xi-myiti)≥1-ξijkξijk≥0M≥0]]>

对每一个样本类别分别学习一个线性变换；

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

公式二||Lyi(xi-myiti)||22+1<||Lj(xi-mjk)||22]]>

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

ϵ(L1,...,LC)=ΣiΣj≠yiΣk=1nj[1+||Lyi(xi-myiti)||22-||Lj(xi-mjk)||22]++μΣj=1Ctr(LjTLj)]]>

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

最小化：ΣiΣj≠yiΣk=1njξijk+μΣj=1Ctr(Mj)]]>

制约条件为：

(xi-mjk)TMj(xi-mjk)-(xi-myiti)TMyi(xi-myiti)≥1-ξijkξijk≥0Mj≥0j=1,...,C]]>

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

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

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

ϵ(Mt-1)=ΣiΣj≠yiΣk=1nj[1+(xi-myiti)TMt-1(xi-myiti)-](xi-mjk)TMt-1(xi-mjk)]+·+μ(tr(Mt-1))]]>

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

Gt=Σ(i,j,k)∈ψt[(xi-myiti)(xi-myiti)T-(xi-mjk)(xi-mjk)T]+μI]]>

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

Gt+1=Gt+Σ(i,j,k)∈ψt+1-ψt[(xi-myiti)(xi-myiti)T-(xi-mjk)(xi-mjk)T]]]>

-Σ(i,j,k)∈ψt-ψt+1[(xi-myiti)(xi-myiti)T-(xi-mjk)(xi-mjk)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。