CN104063632B

CN104063632B - 基于回归森林模型的蛋白质序列二硫键连接模式的预测方法

Info

Publication number: CN104063632B
Application number: CN201410303084.6A
Authority: CN
Inventors: 於东军; 李阳; 胡俊; 沈红斌; 杨静宇
Original assignee: Nanjing University of Science and Technology
Current assignee: Nanjing University of Science and Technology
Priority date: 2014-06-27
Filing date: 2014-06-27
Publication date: 2017-09-01
Anticipated expiration: 2034-06-27
Also published as: CN104063632A

Abstract

本发明公开一种基于回归森林模型的蛋白质序列二硫键连接模式的预测方法，包括以下步骤：步骤1、通过多视角特征提取和特征组合，得到蛋白质序列中每个半胱氨酸残基对的特征向量；步骤2、对于待预测的蛋白质序列信息和训练数据集合，生成待预测的蛋白质序列信息和训练数据集合中所有半胱氨酸残基对的特征向量，从而分别构成训练样本集合与待预测样本集合；步骤3、采用随机森林算法学习二硫键样本在特征空间中的分布规律，生成随机森林回归模型；步骤4、利用随机森林回归模型对待预测样本集合的特征向量进行预测，得到每个半胱氨酸残基对形成二硫键的倾向值，得分最高的二硫键连接模式即为最终预测得到的蛋白质序列中的二硫键连接模式。

Description

基于回归森林模型的蛋白质序列二硫键连接模式的预测方法

技术领域

本发明涉及生物信息学蛋白质序列中二硫键预测技术领域，具体而言涉及一种基于回归森林模型的蛋白质序列二硫键连接模式的预测方法。

背景技术

二硫键是最重要的蛋白质结构特性之一。二硫键是在蛋白质多肽链中两个半胱氨酸残基之间形成的主共价键，它们可以在肽链的链间或者链内形成。二硫键在蛋白质折叠方式以及稳定性方面有着非常重要的作用。因此，预测蛋白质里面的半胱氨酸残基组成二硫键的方式在预测蛋白质结构和功能上起着举足轻重的作用。

目前有很多种预测二硫键的方法，比如，DISULFIND法(A.Ceroni,A.Passerini,A.Vullo et al.,“DISULFIND:a disulfide bonding state and cysteine connectivityprediction server,”Nucleic Acids Research,vol.34,no.suppl2,pp.W177-W181,2006.)、Pair-Wise SVM法(C.-H.Tsai,B.-J.Chen,C.-h.Chan et al.,“Improvingdisulfide connectivity prediction with sequential distance between oxidizedcysteines,”Bioinformatics,vol.21,no.24,pp.4416-4419,2005.)、SS_SVR法(J.Song,Z.Yuan,H.Tan et al.,“Predicting disulfide connectivity from protein sequenceusing multiple sequence feature vectors and secondary structure,”Bioinformatics,vol.23,no.23,pp.3147-3154,2007.)、FS_SVR法(L.Zhu,J.Yang,J.N.Song et al.,“Improving the accuracy of predicting disulfide connectivityby feature selection,”J Comput Chem,vol.31,no.7,pp.1478-85,May,2010.)、DBCP法(H.-H.Lin,and L.-Y.Tseng,“DBCP:a web server for disulfide bondingconnectivity pattern prediction without the prior knowledge of the bondingstate of cysteines,”Nucleic acids research,vol.38,no.suppl2,pp.W503-W507,2010.)、DISLOCATE法(DISLOCATE+MIp+iCOV)(C.Savojardo,P.Fariselli,M.Alhamdooshet al.,“Improving the prediction of disulfide bonds in Eukaryotes withmachine learning methods and protein subcellular localization,”Bioinformatics,vol.27,no.16,pp.2224-30,Aug15,2011.)、DMC法(C.Savojardo,P.Fariselli,P.L.Martelli et al.,“Prediction of disulfide connectivity inproteins with machine-learning methods and correlated mutations,”BMCBioinformatics,vol.14,no.Suppl1,pp.S10,2013.)、DiANNA法(F.Ferrè,and P.Clote,“DiANNA1.1:an extension of the DiANNA web server for ternary cysteineclassification,”Nucleic Acids Research,vol.34,no.suppl2,pp.W182-W185,2006.)等等。

前述方法大体可以分为三类：(Ⅰ)二硫键成键状态预测；(Ⅱ)已知半胱氨酸是否形成二硫键的信息来预测二硫键的组成模式；(Ⅲ)对上述两者都进行预测。最近，利用机器学习的方法预测蛋白质二硫键模式受到越来越多人的关注，从实验的结果来看，运用先进的机器学习算法可以很大程度上提高预测性能，选取有效的特征和利用强大的机器学习算法成为可以显著影响影响二硫键连接模式预测性能的两个重要方面。综上所述，现有技术中虽然已经在预测二硫键模式方面取得了一定的发展，不过在预测精度方面仍有提升的空间。

发明内容

本发明目的在于提供一种基于回归森林模型的蛋白质序列二硫键连接模式的预测方法，提高预测性能。

为达成上述目的，本发明所采用的的技术方案如下：

一种基于回归森林模型的蛋白质序列二硫键连接模式的预测方法，其实现包括以下步骤：

步骤1、特征提取，基于输入的蛋白质序列信息，进行多视角特征提取和特征组合，即：

分别使用PSI-BLAST算法抽取蛋白质序列的进化信息，利用PSIPRED算法抽取蛋白质序列的二级结构信息，再利用滑动窗口与特征串行组合方式从前述进化信息和二级结构信息中提取每个半胱氨酸残基对的多视角特征；

然后计算蛋白质序列中半胱氨酸残基对在蛋白质序列中的距离特征信息(DOC)，利用PSI-BLAST算法抽取蛋白质序列中每个半胱氨酸残基对的协变异信息(CM)，利用MODELLER同源建模软件获取半胱氨酸残基对的三维结构距离信息(PDTCR)；

最后将上述从进化信息和二级结构信息中提取每个半胱氨酸残基对的多视角特征与前述距离特征信息、协变异信息、三维结构距离信息进行特征组合，得到蛋白质序列中每个半胱氨酸残基对最终的特征向量；

步骤2、对于待预测的蛋白质序列信息和训练数据集合，采用上述步骤1的方法，生成待预测的蛋白质序列信息和训练数据集合中所有半胱氨酸残基对的特征向量，从而构成训练样本集合与待预测样本集合；

步骤3、特征学习

采用随机森林算法(Random Forests,RF)学习二硫键样本在特征空间中的分布规律，生成随机森林回归模型，其中，模型的决策树个数(nTree)采用大于等于300的整数，每次分裂时随机选择的特征个数(mTry)由前述步骤1所得的特征向量的维度决定；

步骤4、二硫键连接模式预测

利用前述步骤3的随机森林回归模型对待预测样本集合中的特征向量进行预测，得到每个半胱氨酸残基对形成二硫键的倾向值，记为ppcp；

令p为所有可能的二硫键连接模式的数量，则第i个二硫键连接模式的得分S_i用下述公式得到：

其中B代表第i个二硫键连接模式中二硫键的个数，ppcp_j是第i个二硫键连接模式中第j个半胱氨酸残基对的倾向值：

如果第i个二硫键的连接模式的得分最高，则该第i个模式就是最终预测得到的结果：

即：该得分最高的二硫键连接模式i^*为最终预测得到的蛋白质序列中的二硫键连接模式。

由以上本发明的技术方案可知，本发明的有益效果在于：

1、提高预测速度与精度：PDTCR特征直接计算了两个半胱氨酸之间的距离，而这个距离值有跟二硫键的成键有着直接的关系，所以可以显著提高预测精度；

2、在随机森林中运用的随机策略已经证明要优于许多其他的分类器，并且随机森林还有抗过拟合的特性，随机森林的训练和预测过程相比支持向量回归更加快速，训练时参数获取也更加容易，预测精度也有一定的提高。

附图说明

图1为本发明一实施方式基于回归森林模型的蛋白质序列二硫键连接模式的预测方法的实现流程示意图。

具体实施方式

为了更了解本发明的技术内容，特举具体实施例并配合所附图式说明如下。

如图1所示，根据本发明的较优实施例，一种基于回归森林模型的蛋白质序列二硫键连接模式的预测方法，其实现包括以下步骤：

步骤3、特征学习

步骤4、二硫键连接模式预测

作为可选的实施方式，前述步骤1的实现过程中，其中的蛋白质序列的进化信息(即位置特异性得分矩阵)、蛋白质序列的二级结构信息、半胱氨酸对在序列中的距离特征(DOC)、半胱氨酸残基对的协变异信息(CM)、以及半胱氨酸残基对的三维结构距离信息(PDTCR)，其具体的提取和组合方式如下。

A：蛋白质序列的进化信息，即位置特异性得分矩阵

对于某一给定的含有l个氨基酸残基的特定蛋白质P，其原始位置特异性得分矩阵PSSM首先由PSI-BLAST程序用默认的E-value参数从Swiss-Prot数据库(该数据库为现有技术中的公共数据库)经过三次迭代得到的，该矩阵为一个n行20列矩阵，表达如下：

其中：A、R、...、V表示20种常见氨基酸残基，O_i,j表示蛋白质的第i个氨基酸残基在进化过程中突变成20种常见氨基酸残基中的第j种氨基酸残基的可能性。

然后，对原始位置特异性得分矩阵中的每一个元素x，利用下述函数公式(2)归一化。

归一化后的位置特异性得分矩阵，表达如下：

再次，对于蛋白质序列P中的一个半胱氨酸C_i(i为半胱氨酸残基在序列中的位置)，以位置特异性得分矩阵中的第i行为中心，使用一个宽度为W的窗口，该窗口内所有元素构成第i个氨基酸残基的PSSM特征矩阵，表达如下：

最后，将上述特征矩阵(4)按行优先的组合方式形成维数为20×W的特征向量，取滑动窗口大小为13，即W＝13，则半胱氨酸残基C_i的遗传信息可以用下面特征向量表示：

对于一个半胱氨酸残基对，分别计算两个半胱氨酸残基的遗传信息特征向量(即前述表达5所示的特征向量)，并串行组合，即可得到一个半胱氨酸残基对的遗传信息特征向量。

B、蛋白质序列的二级结构特征(PSS)

对于一个含有l个氨基酸残基的蛋白质P，首先通过PSIPRED算法生成蛋白质序列P的二级结构矩阵，该矩阵大小为l×3，表达如下：

其中，C、H、E表示蛋白质的三种二级结构(coil、helix、strand)，S_i,1表示蛋白质P中半胱氨酸残基C_i的二级结构是coil的概率，S_i,2表示蛋白质中半胱氨酸残基C_i的二级结构是helix的概率，S_i,3表示蛋白质中半胱氨酸残基C_i的二级结构是strand的概率；然后，使用上述同样的宽度为W的滑动窗口与按行优先的组合方式，一个半胱氨酸残基可以由一个维数为3×W的特征向量表示：

对于一个半胱氨酸残基对，分别计算两个半胱氨酸残基的二级结构信息特征向量(即前述表达7所示的特征向量)，并串行组合，即可得到一个半胱氨酸残基对的二级结构信息特征向量。

C、半胱氨酸对在序列中的距离特征(DOC)

根据Tsai等人首先提出的半胱氨酸对的距离特征(DOC)研究，本实施例中采用下述公式来表示一个半胱氨酸对的距离特征：

DOC(i,j)＝|i-j|

其中，i和j分别是两个半胱氨酸残基在蛋白质序列中的位置。

然后，再使用DOC_log作为归一化方式，得到距离特征(DOC)的最终表达：

D、协变异特征(CM)

如果半胱氨酸残基对存在某种较强的联系，那么这个半胱氨酸残基对就可能独立于进化过程，鉴于此，借鉴现有技术中Rubinstein和Fiser等人通过协变异分析进行二硫键预测研究，本实施例中采用其所提出的打分模式，对每个半胱氨酸残基对进行打分，这样每个半胱氨酸残基对都获得一个介于[0,1]的得分。

关于协变异的具体打分模式，如Rubinstein和Fiser所著“Predicting disulfidebond connectivity in proteins by correlated mutations analysis[J].Bioinformatics,2008,24(4):498-504.”

E、半胱氨酸残基对三维结构距离特征(PDTCR)

使用MODELLER同源建模软件从蛋白质序列得到该蛋白质的三维结构；然后再利用得到的结构中的坐标值计算出两个半胱氨酸的距离。半胱氨酸的间距越小，它们形成一个二硫键的可能性就越大。

最终，每个半胱氨酸残基对的特征向量将由上面提到的五个类型的特征组成，是位置特异性得分矩阵PSSM(520维，520＝2×260)、二级结构特征PSS(78维，78＝2×39)、距离特征DOC(1维)、协变异特征CM(1维)、三维结构距离特征PDTCR(1维)，最终得到的特征向量，其维度为520+78+1+1+1＝601维。

如图1所示，对于待预测的蛋白质序列信息和训练数据集合，均采用上述步骤1的方法(特征提取和特征串行组合)，生成待预测的蛋白质序列信息和训练数据集合中所有半胱氨酸残基对的特征向量，从而分别构成训练样本集合与待预测样本集合。

然后，在步骤3的特征学习中，采用随机森林算法(Random Forests,RF)学习二硫键样本在特征空间中的分布规律，生成随机森林回归模型，其中，模型的决策树个数(nTree)采用大于等于300的整数，每次分裂时随机选择的特征个数(mTry)由前述步骤1所得的特征向量的维度决定。

本实施例中，前述决策树个数(nTree)优选500。每次分裂时随机选择的特征个数(mTry)取值为其中601即前述特征向量的特征维度。

接下来在步骤4中，采用前述步骤3的随机森林回归模型对待预测样本集合中的特征向量进行预测，得到每个半胱氨酸残基对形成二硫键的倾向值，记为ppcp；

如图1所示，以图1中的一对半胱氨酸残基对作为示例，该两个半胱氨酸的位置分别是3和n-1，如前述步骤1所述的，获取此半胱氨酸残基对的PSSM、PSS、CM、DOC还有PDTCR特征后，进行特征组合可得到一个601维的一个特征向量。

再用已经训练完毕的随机森林回归模型预测这个特征向量，得到该个半胱氨酸残基对形成二硫键的倾向值。

如前述步骤4的处理，取得分最高的二硫键的连接模式记为最终的预测结果。

如图1所示，其中所有可能的二硫键模式包括三种：[C1-C2,C3-C4]，[C1-C3,C2-C4]，[C1-C4,C2-C3]，按照图1所示，其得分别为：

S₁＝0.3+0.4＝0.7；

S₂＝0.1+0.2＝0.3；

S₃＝0.0+0.5＝0.5。

由此可见，第一个模式的得分最高，即为最终所预测的二硫键连接模式。

本实施例中，将PDTCR特征与PSSM、PSS、CM和DOC特征相结合，在SP39数据集和PDBCYS-R数据集的基础上，进行了对比试验：只用不包含PDTCR特征的四个特征(即PSSM、PSS、CM和DOC特征)通过交叉验证，同样地，在前述两个数据集(SP39数据集和PDBCYS-R数据集)上用四个原始特征(即PSSM、PSS、CM和DOC特征)加上PDTCR特征再进行交叉验证，然后得到预测性能的结果对比。

其中，前述SP39数据集如FARISELLI P,CASADIO R.“Prediction of disulfideconnectivity in proteins”[J].(Bioinformatics,2001,17(10):957-64).

前述PDBCYS-R数据集是根据SAVOJARDO C,FARISELLI P,ALHAMDOOSH M,et al.“Improving the prediction of disulfide bonds in Eukaryotes with machinelearning methods and protein subcellular localization”([J].Bioinformatics,2011,27(16):2224-30.)所提出的PDBCYS数据集，本实施例中，基于此PDBCYS数据集把二硫键个数小于2或者大于5的序列都移除，剩下的263条蛋白质序列形成了PDBCYS-R数据集。

由于现有的预测方法多数是基于上述两个数据集，因此本实施例中，在SP39数据集和PDBCYS-R两个数据集上，分别采用4轮交叉验证和20轮交叉验证，表4.1显示了分别使用随机森林RF(本发明提出的方法)和支持向量回归SVR在有无PDTCR 特征的情况下的实验结果。

表4.1—分别用支持向量回归和随机森林在有无PDTCR特征的情况下的实验结果

^awithout-PDTCR:四个传统特征,PSSM,PSS,CM,和DOC；

^bwith-PDTCR:四个传统特征加上PDTCR特征；

^c通过才用4轮交叉验证得到的结果。

^d通过采用20轮交叉验证得到的结果。

其中：

Q_C用于衡量在所有的二硫键中可以正确预测出来的二硫键的百分比，它的定义如下：

其中N_C是正确预测出来的二硫键个数，T_C是测试数据集中所有的二硫键个数。

Q_P用于衡量所有蛋白质中可以正确预测出来的蛋白质的百分比。只有一个蛋白质中所有的二硫键都被预测出来并且预测出来的二硫键数量要跟真是的数量相同，这样才可以判断这个蛋白质被正确预测出来。Q_P的定义如下：

其中N_P是正确预测出来的蛋白质的数量，而T_P是数据集中所有蛋白质的数量。

表4.1中的数据表明，结合新的PDTCR特征，预测性能可以得到明显的提升。在二硫键个数(B＝2、3、4、5)不同的情况下，Q_P和Q_C的值在用了PDTCR特征之后相比不用新特征几乎都有提高。在SP39数据集跟PDBCYS-R数据集上无论用支持向量回归亦或随机森林，结果都是如此。我们以总的Q_P和Q_C为例，在SP39数据集上加上了PDTCR特征之后，用SVR方法两者分别提高了3.6％跟2.8％，用RF方法两者分别提高了1.6％和1.9％。而在PDBCYS-R数据集上，用SVR方法两者分别提高了8.0％和7.3％，用RF则分别提高了6.1％和4.4％。表2中的数据足以说明们新加特征的有效性。预测性能得到很大提升的原因是PDTCR特征直接计算了两个半胱氨酸之间的距离，而这个距离值有跟二硫键的成键有着直接的关系。明显，在同样的特征下，使用随机森林比使用SVR可以获得更高预测的精度。

而且，在训练随机森林模型时，随机森林还进行了特征选择，从而降低了维度并消除冗余，提高预测性能同时减少预测时间。

虽然本发明已以较佳实施例揭露如上，然其并非用以限定本发明。本发明所属技术领域中具有通常知识者，在不脱离本发明的精神和范围内，当可作各种的更动与润饰。因此，本发明的保护范围当视权利要求书所界定者为准。

Claims

1.一种基于回归森林模型的蛋白质序列二硫键连接模式的预测方法，其特征在于，包括以下步骤：

步骤1、特征提取，即基于输入的蛋白质序列信息，进行多视角特征提取和特征组合，得到蛋白质序列中每个半胱氨酸残基对最终的特征向量；

步骤2、对于待预测的蛋白质序列信息和训练数据集合，采用上述步骤1的方法，生成待预测的蛋白质序列信息和训练数据集合中所有半胱氨酸残基对的特征向量，从而分别构成训练样本集合与待预测样本集合；

步骤3、特征学习

采用随机森林算法学习二硫键样本在特征空间中的分布规律，生成随机森林回归模型，其中，模型的决策树个数采用大于等于300的整数，每次分裂时随机选择的特征个数由前述步骤1所得的特征向量的维度决定；

步骤4、二硫键连接模式预测

令q为所有可能的二硫键连接模式的数量，则第i个二硫键连接模式的得分S_i用下述公式得到：

<mrow> <msub> <mi>S</mi> <mi>i</mi> </msub> <mo>=</mo> <msubsup> <mi>&Sigma;</mi> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>B</mi> </msubsup> <msub> <mi>ppcp</mi> <mi>j</mi> </msub> <mo>,</mo> <mn>1</mn> <mo>&le;</mo> <mi>i</mi> <mo>&le;</mo> <mi>q</mi> </mrow>

即：该得分最高的二硫键连接模式i^*为最终预测得到的蛋白质序列中的二硫键连接模式；

其中，所述步骤1中特征提取的具体实现包括以下步骤：

分别使用PSI-BLAST算法抽取蛋白质序列的进化信息，利用PSIPRED算法抽取蛋白质序列的二级结构信息，再利用滑动窗口与特征串行组合方式从前述进化信息和二级结构信息中提取每个半胱氨酸残基对的多视角特征；然后计算蛋白质序列中半胱氨酸残基对在蛋白质序列中的距离特征信息，利用PSI-BLAST算法抽取蛋白质序列中每个半胱氨酸残基对的协变异信息，利用MODELLER同源建模软件获取半胱氨酸残基对的三维结构距离信息；最后将上述从进化信息和二级结构信息中提取每个半胱氨酸残基对的多视角特征与前述距离特征信息、协变异信息、三维结构距离信息进行特征组合，得到蛋白质序列中每个半胱氨酸残基对最终的特征向量；具体地，其中的蛋白质序列的进化信息即位置特异性得分矩阵、蛋白质序列的二级结构信息、半胱氨酸对在序列中的距离特征、半胱氨酸残基对的协变异信息以及半胱氨酸残基对的三维结构距离信息，其具体的提取和组合方式如下：

A：蛋白质序列的进化信息，即位置特异性得分矩阵

对于某一给定的含有l个氨基酸残基的特定蛋白质P，其原始位置特异性得分矩阵PSSM首先由PSI-BLAST程序用默认的E-value参数从Swiss-Prot数据库经过三次迭代得到的，该矩阵为一个n行20列矩阵，表达如下：

其中：A、R、...、V表示20种常见氨基酸残基，O_i,j表示蛋白质的第i个氨基酸残基在进化过程中突变成20种常见氨基酸残基中的第j种氨基酸残基的可能性；

然后，对原始位置特异性得分矩阵中的每一个元素x，利用下述函数公式(2)归一化

归一化后的位置特异性得分矩阵，表达如下：

再次，对于蛋白质序列P中的一个半胱氨酸C_i，i为半胱氨酸残基在序列中的位置，以位置特异性得分矩阵中的第i行为中心，使用一个宽度为W的窗口，该窗口内所有元素构成第i个氨基酸残基的PSSM特征矩阵，表达如下：

对于一个半胱氨酸残基对，分别计算两个半胱氨酸残基的遗传信息特征向量即前述表达5所示的特征向量，并串行组合，即可得到一个半胱氨酸残基对的遗传信息特征向量；

B、蛋白质序列的二级结构特征

其中，C、H、E表示蛋白质的三种二级结构即coil、helix、strand，S_i,1表示蛋白质P中半胱氨酸残基C_i的二级结构是coil的概率，S_i,2表示蛋白质中半胱氨酸残基C_i的二级结构是helix的概率，S_i,3表示蛋白质中半胱氨酸残基C_i的二级结构是strand的概率；然后，使用上述同样的宽度为W的滑动窗口与按行优先的组合方式，一个半胱氨酸残基可以由一个维数为3×W的特征向量表示：

对于一个半胱氨酸残基对，分别计算两个半胱氨酸残基的二级结构信息特征向量即前述表达7所示的特征向量，并串行组合，即可得到一个半胱氨酸残基对的二级结构信息特征向量；

C、半胱氨酸对在序列中的距离特征

采用下述公式来表示一个半胱氨酸对的距离特征：

DOC(i,j)＝|i-j|

其中，i和j分别是两个半胱氨酸残基在蛋白质序列中的位置；

然后，再使用DOC_log作为归一化方式，得到距离特征的最终表达：

D、协变异特征

如果半胱氨酸残基对存在某种较强的联系，那么这个半胱氨酸残基对就可能独立于进化过程，采用打分模式，对每个半胱氨酸残基对进行打分，这样每个半胱氨酸残基对都获得一个介于[0,1]的得分；

E、半胱氨酸残基对三维结构距离特征

使用MODELLER同源建模软件从蛋白质序列得到该蛋白质的三维结构；然后再利用得到的结构中的坐标值计算出两个半胱氨酸的距离；半胱氨酸的间距越小，它们形成一个二硫键的可能性就越大；

最终，每个半胱氨酸残基对的特征向量将由上面提到的五个类型的特征组成，是520维的位置特异性得分矩阵PSSM、78维的二级结构特征PSS、1维的距离特征DOC、1维的协变异特征CM、1维的三维结构距离特征PDTCR，最终得到的特征向量，其维度为520+78+1+1+1＝601维。

2.根据权利要求1所述的基于回归森林模型的蛋白质序列二硫键连接模式的预测方法，其特征在于，前述随机森林回归模型的每次分裂时随机选择的特征个数为24。

3.根据权利要求1所述的基于回归森林模型的蛋白质序列二硫键连接模式的预测方法，其特征在于，前述随机森林回归模型的决策树个数取值为500。