CN106529207B - 一种与核糖核酸结合的蛋白质的预测方法 - Google Patents

Abstract

本发明公开了一种与核糖核酸结合的蛋白质(RBP)的预测方法。首先根据氨基酸的特性，获取样本蛋白的候选特征总集；然后以候选特征之间的相关度和/或冗余度为标准，从候选特征总集中选取个最佳特征作为特征向量；根据样本蛋白的特征向量以及蛋白质特征，建立预测模型；最后根据待预测蛋白的特征向量，获得待预测蛋白的蛋白质特性的预测结果。本发明的候选特征总集涵盖了氨基酸的多种特性，全面考虑了影响蛋白质结合核糖核酸性能的多方面因素，准确度经验证超过90％，将现有技术的准确度提高了35％，马修相关系数为0.788，将现有技术提高了2倍，从而预测更为全面准确。

Description

一种与核糖核酸结合的蛋白质的预测方法

技术领域

本发明属于生物大分子间相互作用的预测领域，更具体地，涉及一种与核糖核酸结合的蛋白质(RBP)的预测方法。

背景技术

能够与核糖核酸(RNA)结合的蛋白质称为核糖核酸结合蛋白(RBP)，不能与核糖核酸结合的蛋白质称为非核糖核酸结合蛋白(non-RBP)。在生物体内，RBP通过和核糖核酸(RNA)相互作用形成复合物，在很多生物过程中起到重要作用，比如转录后基因调控、基因的可变剪切和翻译等，因此预测蛋白质是否为RBP很重要。

非专利文献(Zhao,H.,Y.Yang,and Y.Zhou,.RNA biology,2011.8(6):p.988-996)公开了一种高精度的核糖核酸结合蛋白质的预测方法(SPOT-seq)。该方法基于RBP的结构构造模板库，将每个靶蛋白质序列到模板库中搜索，把靶蛋白质序列和模板库中RBP的结构进行序列与结构的匹配，通过打分来预测靶蛋白是否是RBP。该方法在大多数情况下仅能对non-RBP作出正确预测，对于RBP的预测则常常出现错误。

发明内容

针对现有技术的以上缺陷或改进需求，本发明提供了一种与核糖核酸结合的蛋白质的预测方法，其目的在于选取合适的特征向量，由此解决现有技术对于RBP的预测的准确性低的技术问题。

为实现上述目的，按照本发明的一个方面，提供了一种与核糖核酸结合的蛋白质的预测方法，包括以下步骤：

S1.获得样本蛋白的候选特征总集

根据氨基酸的疏水性，获取第一候选特征的集合；根据氨基酸的极性，获得第二候选特征的集合；根据氨基酸的规范化的范德瓦尔斯体积，获取第三候选特征的集合；根据氨基酸的极化性，获取第四候选特征的集合；根据氨基酸的极化性，获取第五候选特征的集合；根据氨基酸的溶剂可及性，获取第六候选特征的集合；根据氨基酸的带电性和极性，获取第七候选特征的集合；根据氨基酸的位置特异性打分矩阵，获取第八候选特征的集合；

根据第一候选特征的集合至第八候选特征的集合的并集，获得样本蛋白的候选特征总集V_i＝{v_i,k}；其中，v_i,k表示第i个样本蛋白的第k个候选特征，i为1～ξ的任意整数，k为1～N₁的任意整数，ξ为样本蛋白的总数量，N₁为候选特征总集中候选特征的总数量，ξ≥5N₁；

S2.选取最佳特征

根据候选特征之间的相关度和/或冗余度，在所述第i个样本蛋白的候选特征总集中选取个候选特征作为第i个样本蛋白的最佳特征，其分别为并获得第i个样本蛋白的特征向量 其中，

S3.建立预测模型

根据样本蛋白的特征向量以及蛋白质特性，建立预测模型f(x_i)，并获取样品蛋白在所述预测模型中的评估指标σ，如果σ≥阈值η，进入S4，否则返回步骤S2；所述蛋白质特性为样本蛋白是否能与核糖核酸结合的特性；

S4.获得预测结果

获取待预测蛋白与样本蛋白对应的特征向量x，代入步骤S3中所建立的预测模型，获得待预测蛋白的蛋白质特性的预测结果f(x)。

优选地，所述步骤S2中选取候选特征的方法为最大相关最小冗余法。

优选地，所述步骤S3中建立预测模型的方法为支持向量机法、神经网络法、贝叶斯分类法或随机森林法。

作为进一步优选地，所述步骤S3中建立预测模型的方法为支持向量机法。

优选地，所述步骤S3中的评估指标为准确度、马修相关系数或受试者操作特性曲线下的面积。

作为进一步优选地，所述步骤S3中的评估指标为马修相关系数，所述阈值η为0.5。

总体而言，通过本发明所构思的以上技术方案与现有技术相比，具有以下有益效果：

1、候选特征总集涵盖了氨基酸的多种特性，全面考虑了影响蛋白质结合核糖核酸性能的多方面因素，从而预测更为全面准确；

2、以相关度和/或冗余度为标准，筛选获得最佳特征，减少了计算量的同时，避免了由于最佳特征的冗余而造成的计算误差；

3、以马修相关系数优选作为评估指标，比通过准确度评估更科学，建立的预测模型更为准确；

4、本发明方法预测获得的蛋白质特性的准确度经验证超过90％，将现有技术的准确度提高了35％，马修相关系数为0.788，将现有技术提高了2倍，从而更具有实用性和应用价值。

附图说明

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

图2为本发明的实施例1的预测方法的特征向量获取图；

图3为本发明的实施例1的预测方法的流程示意图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明的各个特性所涉及到的特征只要彼此之间未构成冲突就可以相互组合。

在本发明中，氨基酸的中文名称以及对应的字母简写分别如下：精氨酸(Arg,R)，天冬氨酸(Asp,D)，半胱氨酸(Cys,C)，谷氨酰胺Gln,Q)，谷氨酸(Glu,E)，组氨酸(His,H)，异亮氨酸(Ile,I)，甘氨酸(Gly,G)，天冬酰胺(Asn,N)，亮氨酸(Leu,L)，赖氨酸(Lys,K)，甲硫氨酸(Met,M)，苯丙氨酸(Phe,F)，脯氨酸(Pro,P)，丝氨酸(Ser,S)，苏氨酸(Thr,T)，色氨酸(Trp,W)，酪氨酸(Tyr,Y)，缬氨酸(Val,V)。

本发明提供了一种与核糖核酸结合的蛋白质的预测方法，包括以下步骤，如图1所示：

S1.根据氨基酸的特性，获取样本蛋白的候选特征总集V_i＝{v_i,k}；其中，v_i,_k表示第i个样本蛋白的第k个候选特征，i为1～N₁的任意整数，N₁为候选特征总集中候选特征的总数量；为了保证步骤S2中选取最佳特征的准确性，样本蛋白的总数量ξ≥5N₁，为了保证样本蛋白与待预测蛋白的真实情况接近，样本蛋白中通常RBP的数量要小于non-RBP；所述候选特征总集为第一特征集合至第八特征集合的并集；

其中，所述第一特征集合至第七特征集合分别对应的氨基酸的特性为：疏水性、极性、规范化的范德瓦尔斯体积、极化性、二级结构、溶剂可及性、带电性和极性；

所述第八特征集合对应的氨基酸的特性为位置特异性打分矩阵；

其中，获取第一特征集合至第七特征集合的方法具体为：

S11.根据氨基酸特性，将20种氨基酸分为K组，隶属于第1组至第K组的氨基酸分别以A₁～A_K表示，K为2、3或4；其中，由于疏水性，极性，规范化的范德瓦尔斯体积以及极化性为相对标准，可以根据该相对标准，将氨基酸平均分为2～4组；而二级结构，溶剂可及性，带电性和极性都为绝对标准，其中，二级结构即根据该氨基酸属于样本蛋白的螺旋结构、发卡结构或其它结构分为三组，溶剂可及性根据氨基酸的R基暴露于样本蛋白的表面或者埋藏于样本蛋白的内部分为两组，而带电性和极性则根据氨基酸为极性带负电、极性带正电、极性不带电以及非极性分为四组；

S12.以全局蛋白质序列描述符或者三联体作为编码方法，获取与样本蛋白的氨基酸特性对应的候选特征；

其中，以全局蛋白质序列描述符作为编码方法具体包括：A_i在氨基酸序列中的比例，第一个A_i在氨基酸序列中位置的比例，第25％个A_i在氨基酸序列中位置的比例，第50％个A_i在氨基酸序列中位置的比例，第75％个A_i在氨基酸序列中位置的比例，最后一个A_i在氨基酸序列中位置的比例，A_i和A_j无序组成的二联体在氨基酸序列中出现的概率；其中，i和j为1～K的任意整数，且i≠j；

以三联体作为编码方法具体为：A_iA_jA_n有序组成的三联体在氨基酸序列中出现的概率；i、j和n为1～K的任意整数；

S13.获取第八特征集合的方法具体为：

S131.利用样本蛋白的氨基酸序列搜索蛋白质序列数据库，获得位置特异性打分矩阵P_ij，其中，j为1～20，分别表示20种氨基酸，i表示样本蛋白含有的氨基酸个数；

S132.获得归一化后的位置特异性打分矩阵

S133.把位置特异性打分矩阵中每一列中处在氨基酸序列不同位置处的相同氨基酸对应的P_ij‘值相加，获得第八特征集合的候选特征P_mj”＝∑P_ij’(样本蛋白的氨基酸序列的第i个氨基酸为第m种氨基酸，m为1～20)，获得第八特征集合{P_mj”}，m、j为1～20；

S134.可将每种候选特征分别标准化，并整理为与样本蛋白一一对应的候选特征总集V_i＝{v_i,k}，也可直接获得候选特征总集V_i＝{v_i,k}；

S2.选取最佳特征

如果将所有候选特征都作为最佳特征进行预测模型的构建，则建立的预测模型中计算参数太多，从而容易增加计算量，同时由于候选特征之间可能有冗余，反而影响其计算准确度，而选取的最佳特征太少涵盖的氨基酸信息不全，从而影响预测准确度；因此需要从候选特征中选取个最佳特征S_n，n为的任意整数，且

选取方法可通过皮尔逊相关系数，残差分析以及最大相关最小冗余(MinimumRedundancy Maximum Relevance，mRMR)法，其中mRMR法由于同时考虑了候选特征之间的相关度以及冗余度，优选作为选取方法；

该方法具体包括以下子步骤：

S21.对所有样本蛋白的候选特征进行二值化处理，并获得二值化后的候选特征的概率；

令二值化的候选特征其中，通常设置为v_i,k的平均值，即令概率函数p(v_k1)为对应于v_i,k’＝1的v_i,k的概率，即令概率函数p(v_k2)为对应于v_i,k’＝－1的v_t,k的概率，则p(v_k1)+p(v_k2)＝1；

S22.令选取次数

其中，H(v_k)表示候选特征v_k包含的信息，I(v_k,v_n)表示有了候选特征v_k之后其它候选特征的信息的减少量，R₁(v_k)表示候选特征v_k与候选特征总集的平均互信息；

S23.获得候选特征总集的最大相关-最小冗余函数mRMR(v_k)＝R₁(v_k)，并选取mRMR(v_k)最大的v_k作为最佳特征集合中第一个最佳特征S₁；

S24.令R₂(v_k,S_n)＝R₁(S_n)-R₁(S_n|v_k)…(7)；其中

S_n表示从候选特征中选取的最佳特征，n为的任意整数，R₁(S_n|v_k)表示最佳特征S_n对候选特征v_k的条件相关度，R₂(v_k,S_n)表示候选特征v_k对最佳特征S_n的冗余度；

令最大相关-最小冗余函数为根据方程(3)～(8)，选取mRMR(v_k)最大的v_k作为最佳特征集合中第个最佳特征

S25.判断最佳特征的数量是否已达到预先的设定值，是则进入下一步，否则返回上一步，直至获得所需数量的最佳特征组成的最佳特征集合{S_n}；令第i个样本蛋白的特征向量

S3.建立预测模型函数f(x_i)，为了覆盖所有已选特征，至少需要个样本蛋白，以保证最佳特征的相关系数的准确性，利用ξ个已知蛋白质特性的预测结果的样本蛋白迭代求解，即样本蛋白的总数量ξ也需满足建立最佳特征函数的方法包括支持向量机、神经网络法、贝叶斯分类法或随机森林法；

以神经网络法为例说明预测模型的建立过程如下：

S31.获取样本蛋白相应的最佳特征S_n(i)，并对于每个n，获取对应的聚类中心Z_nj ^*及其对应的聚类簇A_nj；其中，i为1～ξ的整数，ξ为样本蛋白的个数，j为1～k_n的整数，k_n为第n个最佳特征的聚类簇的个数；其具体方法为：

S311.由于S_n(i)为序号为i的样本蛋白对应的最佳特征，在n确定的情况下，可利用下式计算该数据点处的密度值D_i，

式中，i和j为小于等于ξ的自然数，作用半径r_a为0.3～0.5的常数。

S312.经过计算所有数据点处的密度值之后，选取密度值最大的数据点作为第一个聚类中心为Z₁ ^*，D₁ ^*记作第一个聚类中心的密度值。选定第一个聚类中心之后，令从而进行密度值修正；其中，β＝4/r_b ²，r_b是一个正值常数，定义密度值的显著减小邻域，为了避免取到过于靠近的聚类中心点，r_b必须大于r_a，通常定义r_b＝1.5r_a或r_b＝1.25r_a。

S313.然后，选取修正后的密度值中的最大值作为第二个聚类中心，利用和步骤S512中相同的方法进行密度修正，以此类推，经过k－1次修正之后，可以得到k个聚类中心；此时相应的密度值中的最大值为D_k ^*；

根据最大密度值D_k ^*，定义两个边界值其中，ε ^*D₁ ^*为下边界值，为上边界值，ε ^*为拒绝率，为接受率，一般定义 ε ^*＝0.15；

在密度修正中，分为三种情况：a.如果修正后的某数据点处的密度值D_k ^*＜ε ^*D₁ ^*；或且则该数据点处的密度值不被接受并且进入步骤S52；其中，δ_min为最大密度值D_k ^*对应的聚类中心Z_k ^*到之前所有聚类中心Z₁ ^*、Z₂ ^*、…、Z_k-2 ^*、Z_k-1 ^*的最小距离；

b.如果且则该点数据处的密度值被接受，令D_k ^*＝0，进入步骤S314；

c.如果修正后的密度值则该数据点处的密度值被接受，进入步骤S314；

S314.k＝k+1，返回步骤312.再进行新一轮的聚类。

每种最佳特征的数据点经过减法聚类后，获得了k个聚类中心Z₁ ^*、Z₂ ^*、…、Z_k ^*，以此进行聚类和模糊空间划分后能获得分别对应k个聚类中心的k个减法聚类簇A₁、A₂、…、A_k。

S32.利用所述步骤S51获得的聚类中心，进行自适应神经网络模型网络训练；

输入参数为最佳特征和蛋白质特性构成的数据集合其中，Y_n为针对第n个最佳特征的所有样本蛋白的数据集合，即Y_n＝[S_n(1),S_n(2),…,S_n(ξ)]；通常以1表示样本蛋白的蛋白质特性为RBP，以－1表示样本蛋白的蛋白质特性为non-RBP，即T为1和－1的值构成的集合；k_n为第n个最佳特征的聚类簇的个数；A_nj为第n个最佳特征的第j个减法聚类簇；对于任意最佳特征的数据点S_n(i)的组合，首先需判断其分别隶属于哪个减法聚类簇。记O^ε(ε＝1,2,3,4,5)表示第ε层的某个节点输出，每一层的功能及意义如下：

Layer1：模糊化层。每个节点代表一个语言变量值，对于每一个减法聚类簇A_nj所对应的聚类中心Z_nj ^*，利用公式(9)计算每个属于该减法聚类簇A_nj的最佳特征S_n(i)的隶属度函数，在该减法聚类簇对应多个最佳特征的数据点S_n(i)的情况下，则O_nj ¹取多个数据点对应的O_nj ¹的均值，没有数据点对应的情况下，则O_nj ¹＝1，其中，r_a表示步骤S511中的作用半径。

Layer2：计算规则使用度。令利用Layer1获得的O_nj ¹相乘，其中，每个j各自独立的表示1～k_n中的任意整数；每个节点代表一条模糊规则，利用公式(10)计算每条模糊规则的适用度：

Layer3：将Layer2获得的O_l ²归一化。利用公式(11)计算第l个节点的规则适用度与所有规则适用度之和的比值：

Layer4：由于对于每组最佳特征S_n(i)，其隶属于不同聚类中心的可能性共有种，因此根据如下模糊规则Rule1～Rulem，可对其输出值进行划分，模糊规则前件为模糊空间划分，模糊规则后件为一阶线性输出(规则中is的含义为最佳特征的数据点S_n(i)属于对应的减法聚类簇A_nj)：令向量模糊规则f_l＝P_l ^Tx_i+r₁，即：

利用公式(13)计算每条模糊规则的输出；如果有多个最佳特征的数据点对应规则f_l时，则令f_l为多个数据点对应的f_l的均值，如果没有最佳特征的数据点S_n(i)对应规则f_l时，则计f_l＝0：

式中，是从第3层传来的归一化后的适用度，是该节点的参数集。

Layer5：总输出。利用公式(14)计算所有传感信号之和作为总输出：

以蛋白质特性作为ANFIS的总输出O⁵，采用最小二乘法更新模糊规则后件的一阶线性方程f_l的系数集然后将Layer4层的系数集设置为定值，根据总输出O⁵与蛋白质特性的误差值，从Layer4～Layer1采用梯度下降法更新模糊规则前件的聚类中心及其作用半径{Z_nj ^*,r_a}；然后将Layer1层的{Z_nj ^*,r_a}设置为定值，以S_n(i)作为输入值，计算总输出O⁵与蛋白质特性的误差值，并与上一次获得的误差值进行比较；如果误差值减少，则自适应神经网络模型尚未获得最佳系数集，则继续从Layer1～Layer4采用最小二乘法更新模糊规则后件的一阶线性方程的系数集并重复以上步骤，直至误差值不再减小为止；最后获得所述自适应神经网络模型的系数则该自适应神经网络模型为该预测模型；当目标蛋白对应的输出O⁵≥0，且该目标蛋白为RBP时，则认为预测准确；或者当目标蛋白对应的输出O⁵＜0，且该目标蛋白为non-RBP时，则认为预测准确；否则认为预测错误；

此外，由于支持向量机SVM有严格的理论和数学基础，如避免了神经网络模型中的经验成分；2)SVM基于结构风险最小化原则，保证具有良好的泛化能力；3)可保证算法的全局最优性；4)通过引用核函数，将输入空间中的非线性问题映射到高维特征空间，在高维空间构造线性函数判别；以支持向量机为例，构建预测模型的过程如下：

S31.将样本蛋白的最佳特征整理为特征向量 且y_i为样本蛋白的标签，通常以y_i＝1表示第i个样本蛋白为RBP，以y_i＝－1表示第i个样本蛋白为non-RBP；

令样品蛋白的特征向量为维空间上的点，该维空间上具有一个最优分类超平面，用于把样品蛋白根据其对应的特征向量分为RBP以及non-RBP；设置该最优分类超平面的法向量ω＝[θ₁,θ₂,…,θ_n,…,θ_ξ]，最优分类超平面的偏移量为b，将数据分为两类的预测模型函数则为f(x_i)＝sgn[ω^TΦ(x_i)+b]，当f(x_i)≥0时，样品蛋白为RBP，否则为non-RBP；

以下列目标函数组为优化目标

y_i[ω^TΦ(x_i)+b]－1+ε_i≥0…(16)；其中，C为错分的惩罚参数，ε_n为松弛项，C>0，ε_i≥0

利用拉格朗日乘数法，引入拉格朗日乘子α_i构造函数(15)的拉格朗日函数沃尔夫对偶型：

引入核函数将特征向量映射至高维空间，核函数k(xi,x_j)＝Φ(xi)^TΦ(x_j)，i、j为1～ξ的任意整数，Φ(x_i)为x_i到高维空间的非线性映射函数；核函数可选用径向基核函数、线性核函数以及皮尔逊VII广义核函数等；

例如，径向基核函数的公式为||·||表示欧式距离，γ表示核参数，γ>0；

则公式(17)可进一步变化为则需满足公式(15)～(18)需：尽可能多的α_i满足0≤α_i≤C，且 利用优化算法(如序贯最小化优化算法)对公式(19)求解，直至所有α_i满足沃尔夫对偶函数的卡罗需-库恩-塔克条件。

S32.将样本蛋白大致平均地分为N₂个子集，子集的数量越多，验证越准确，然而计算时间也会相应受到影响，这里我们以子集的数量N₂＝10为例进行说明：

S321.对log₂C以及log₂γ各自设定初步的搜索范围；例如可将log₂C的搜索范围设置为[－6,20]，log₂γ的搜索范围设置为[－20,6]，在该搜索范围内，各自均匀地取5～20个数据点，如果各自都取14个数据点，则C以及γ的数据组合个数为14×14；令搜索次数ρ＝1；

S322.以十个子集中的一个作为测试集，九个作为训练集，进行交叉验证，并获得交叉验证识别率最高的C以及γ的数据组合，具体包括以下子步骤：

S322a.令C以及γ的数据组合的序号λ＝1；

S322b.针对第一对C以及γ的数据组合，利用九个子集作为训练集，获得预测模型函数的参数ω和b，将剩余的一个子集作为测试集，测试预测模型函数的准确率Pre₁；

S322c.更换作为测试集的子集，并重复步骤S322a，获得针对不同测试集的准确率Pre₂～Pre₁₀；并获得针对第λ对C以及γ的数据组合的平均准确率

S322d.如果所有C以及γ的数据组合都已交叉验证完毕，则选取平均准确率最高的Pre_λ对应的C以及γ的数据组合，作为下一次搜索的C以及γ的范围中心，否则λ＝1+λ，返回步骤S322b；

S323.搜索次数ρ＝ρ+1，以新的log₂C以及log₂γ作为各自搜索范围的中心，缩小搜索范围再依次取5～20个数据点，例如，log₂C可以[log₂C－2,log₂C,log₂C+2]为搜索范围，以0.5为搜索间隔，取9个数据点；搜索次数ρ越大，搜索间隔设置得越小；

S324.重复步骤S322.获得对应第ρ次搜索的平均准确率Pre_λ，如果该平均准确率优于第ρ－1次搜索，则返回步骤S223；否则以第ρ次搜索的平均准确率Pre_λ对应的C、γ、ω和b，作为预测模型f(x_i)＝sgn[ω^TΦ(x_i)+b]的最终参数；

S4.获得利用预测模型进行样本蛋白的评估指标；所述评估指标包括准确度，马修相关系数，或受试者操作特性曲线下的面积；其中，准确度是指样本蛋白的RBP中被正确预测的比例；马修相关系数

其中，真阳性TP为阳性样本(样本蛋白的RBP)中被正确预测的个数，假阴性FN为阳性样本中被错误预测的个数，真阴性TN为阴性样本(样本蛋白的non-RBP)中被正确预测的个数，假阳性FP为阴性样本中被错误预测的个数；

受试者操作特性曲线下的面积是以灵敏度(真阳性率)为纵坐标，特异性(真阳性率)为横坐标绘制的曲线与x轴围成的面积；其中，灵敏度为阳性样本中被正确预测的比例；特异性为阴性样本中被正确预测的比例。

如果利用准确度(阳性样本和阴性样本中被正确预测的比例)来作为评价标准，但该评价标准往往在阳性数据和阴性数据大致相等的时候才准确；当两者有较大差异时，例如大部分为阴性样本时，即使模型建立有误差，使得全部预测结果都为阴性样本，也能获得较高的准确度，这对于判断预测模型是否成功建立将造成影响，因此本发明优选采用马修相关系数作为预测模型的判断标准；

将所述步骤S3中所用的样本蛋白的最佳特征再次代入步骤S3中建立的预测模型，获得样本蛋白的评估指示；当以马修相关系数σ作为评估指标时，如果马修相关系数σ≥阈值η(一般设置为50％)，证实该预测模型的预测精度满足需求，进入S5；否则证明选取的最佳特征的数量可能偏少，返回步骤S2重新选取最佳特征；

S5.获得与样本蛋白所对应的待预测蛋白的最佳特征，并作为输入参数S_n代入步骤S3中所建立的预测模型，获得待预测蛋白的预测结果；例如，当预测模型的函数为支持向量机的f(x_i)＝sgn[ω^TΦ(x_i)+b]时，当f(x_i)≥0时，样品蛋白为RBP，否则为non-RBP。

实施例1

实施例1的预测方法包括特征向量的选取、预测模型的构建以及待预测目标蛋白的预测；

S1.图2为发明实施例1中选取特征向量流程图，分别包括第一特征向量至第八特征向量；分别对应蛋白质所含氨基酸的疏水性，蛋白质所含氨基酸的极性，蛋白质所含氨基酸的规范化的范德瓦尔斯体积，蛋白质所含氨基酸的极化性，蛋白质的二级结构，蛋白质的溶剂可及性；蛋白质所含氨基酸的侧链的带电性和极性；以及蛋白质的进化信息；

S11.其中，第一特征向量至第六特征向量采用的是全局蛋白质序列描述符的编码方法而获得；

这里以规范化的范德瓦尔斯体积为例说明第一特征向量至第四特征向量的获取方法：

根据20个氨基酸的规范化的范德瓦尔斯体积的数值从小到大，将氨基酸划分为第1类{G,A,S,C,T,P,D}、第2类{N,V,E,Q,I,L}、以及第3类{M,H,K,F,R,Y,W}；

获取第i类的氨基酸个数占蛋白质的氨基酸序列的总长度的比例和该类中氨基酸的分布情况(考虑第i类中第一个、第25％个、第50％个、第75％个和最后一个氨基酸在整个氨基酸序列中的位置)，和不同类氨基酸的过渡情况(相邻的氨基酸属于不同类的数目比例)，i为1～3的任意整数。

本实施例中，包括由40个氨基酸组成的蛋白质，其具有如SEQ NO:1所示的氨基酸序列。根据规范化的范德瓦尔斯体积分类，上述序列可标记为3333132222313123223233313312131222111133。因此，其包含的第1类氨基酸有11个，第2类氨基酸有12个，第3类氨基酸有17个，故第i类的氨基酸个数占序列总长度的比例分别为11/40＝0.275,12/40＝0.300,17/40＝0.425；而第1类和第2类氨基酸之间的过渡有5次，第1类和第3类之间的过渡有11次，第2类和第3类氨基酸之间的过渡有7次，故不同类氨基酸的过渡为5/39＝0.128,11/39＝0.282,7/39＝0.179；第一类中第一个氨基酸处在序列5位置处，第25％个氨基酸处在序列位置14处，第50％个氨基酸处在序列位置29处，第75％个氨基酸处在序列位置36处，最后一个氨基酸处在序列位置38处，故第一类中氨基酸的分布情况为5/40＝0.125,14/40＝0.350,29/40＝0.725,36/40＝0.900,38/40＝0.950；同样，第2类氨基酸的分布为7/40＝0.175,9/40＝0.225，17/40＝0.425，28/40＝0.700，34/40＝0.850；第三类氨基酸的分布为1/40＝0.025,6/40＝0.150,19/40＝0.475,25/40＝0.625,40/40＝1.00。因此，SEQ NO:1的规范化的范德瓦尔斯体积所对应的第三特征向量集合(v43～v63分别为(0.275,0.300,0.425,0.128,0.282,0.179,0.125,0.350,0.725,0.900,0.950,0.175，0.225,0.425,0.700,0.850,0.025,0.150,0.475,0.625,1.000)。

同样，氨基酸也可按疏水性、极性以及极化性的大小划分为三类，从而获得相应的第一特征向量(v1～v21)、第二特征向量(v22～v42)以及第四特征向量(v64～v84)。

其中，第五特征向量集合和第六特征向量集合分别对应蛋白质的二级结构和溶剂可及性。蛋白质的二级结构，是指蛋白质的多肽链借助氢键形成的有规则的局部结构；蛋白质的溶剂可及性，是指蛋白质的残基是暴露在蛋白质表面，还是被埋藏在蛋白质内部。通过利用蛋白质的序列相似性与基于序列的结构相似性的结合，来预测蛋白质的二级结构和溶剂可及性。

所述利用蛋白质的序列相似性是指利用与目标蛋白相似的其他蛋白质序列计算出该目标蛋白的进化信息的序列图谱；所述利用基于序列的结构相似性是指如果目标蛋白的氨基酸序列有部分与蛋白质结构数据库PDB中有结构的蛋白质的序列相似，就可用PDB中序列的注释信息去注释该目标蛋白。在本实施例中采用SSPro程序预测蛋白质的二级结构，其精度为92％左右；采用ACCPro程序预测蛋白质的溶剂可及性，其精度为90％左右。其中，二级结构预测结果H表示二级结构为螺旋(helix)，E表示发卡结构(strand)，C表示其他二级结构(others)；溶剂可及性预测结果e表示残基暴露在表面，-表示残基被埋藏；

根据上述三类二级结构，采用全局蛋白质序列描述符，获得第五描述向量(v85～v105)；

对于溶剂可及性，预测结果是两个互补的状态，故只考虑一种状态即可，本实施例中选择编码暴露的残基e。同样根据全局蛋白质序列描述符，针对e这一种溶剂可及性状态，可得到一个7维的特征向量，即第六特征向量(v106～v112)，分别表示e的比例，e向-过渡的比例，排序第一的e在蛋白质中所占的顺序比例，排序为第25％的e在蛋白质中所占的顺序比例，排序为第50％的e在蛋白质中所占的顺序比例，排序为第75％的e在蛋白质中所占的顺序比例，排序最后的e在蛋白质中所占的顺序比例。

S12.第七特征向量则对于侧链的带电性和极性这一特性，采用三联体的编码方式。20个氨基酸根据侧链的带电性和极性的特性被分成4类，第1类-极性带负电{D、E}、第2类-极性带正电{H、R、K}、第3类-极性不带电{C、G、N、Q、S、T、Y}、第4类-非极性{A、F、I、L、M、P、V、W}。三联体的编码方式用来表达每个氨基酸和它邻近氨基酸的特征，一个三联体指的就是三个相邻氨基酸构成的单元，不考虑同一类氨基酸的差异性，认为是相同的。根据以上的分类，遍历整个蛋白质序列，每计算一次向后移一个氨基酸，统计序列中每种三联体出现的频率，并将其归一化。

同样，对SEQ NO:1所示的蛋白质序列，给出侧链带电性和极性这一特性的编码过程。SEQ NO:1所示的序列可标记为2242344343413344133122432434124441443142。序列中出现的三联体频率的最大值为3，最小值为0。如434这类三联体在序列中出现3次，对应的出现概率为(3-0)/3＝1，122三联体出现1次，对应的出现概率为(1-0)/3＝0.333，111三联体出现0次，对应的概率为(0-0)/3＝0。最终的向量包括了每种三联体类型出现的概率，共64种三联体类型，所以得到的是一个64维的向量(v113～v176)，具体值对应(0,0,0,0,0,0.333,0,0.333,0,0,0.666,0,0,0.333,0,0.333,0,0,0,0,0,0,0,0.666,0,0,0,0.333,0,0.333,0.666,0.333,0,0.333,0,0.333,0,0,0,0.333,0.333,0,0,0.333,0.666,0,0.333,0.666,0,0.333,0.666,0.333,0,0,0.333,0,0.333,0.333,0,1,0.666,0,0.666,0.333)。

S13.第八特征向量对应蛋白质的进化信息，即蛋白质序列的位置特异性打分矩阵。蛋白质序列的位置特异性打分矩阵是蛋白质序列上每个位置的氨基酸被替换成上述矩阵中每列位置上的氨基酸的对数似然值；

在本实例中，采用位置特异迭代的基础局部比对搜索工具(PSI-BLAST)进行迭代的多序列比对，获取目标蛋白的位置特异性打分矩阵，所述位置特异性打分矩阵的获得方法如下：

S131.通过将目标蛋白的氨基酸序列与指定蛋白质序列数据库中的蛋白质序列做比较，从该蛋白质序列数据库查找与上述目标蛋白的氨基酸序列具有最优局部比对结果的序列或片段，进行多序列比对，获得目标蛋白的初始的位置特异性打分矩阵；

S132.进行第一次迭代，用S131.步骤产生的目标蛋白的初始的位置特异性打分矩阵去搜索上述指定蛋白质序列数据库，进行矩阵与蛋白质序列数据库中的序列的比对，获得目标蛋白的第一次迭代的位置特异性打分矩阵；

S133.进行第二次迭代，用S132.步骤产生的目标蛋白的第一次迭代的位置特异性打分矩阵去搜索上述指定蛋白质序列数据库，进行矩阵与数据库中蛋白质序列的比对，获得第二次迭代的位置特异性打分矩阵；

S134.进行第三次迭代，用S133.步骤产生的目标蛋白的第二次迭代的位置特异性打分矩阵去搜索上述指定蛋白质序列数据库，进行矩阵与数据库中蛋白质序列的比对，获得所述目标蛋白的位置特异性打分矩阵。

其中，上述指定蛋白质序列数据库指的是序列一致性为90％的蛋白质非冗余序列数据库(即由实验验证的，人工注释的非冗余的蛋白质序列)，本实施例中非冗余序列数据库的数据来源于6个数据库(GenBank CDS translations+RefSeq+PDB+SwissProt+PIR+PRF)；

对上述位置特异性打分矩阵，P_ij代表矩阵中的值，表示目标蛋白序列上第i位置的氨基酸被替换成矩阵中第j列位置上的氨基酸的对数似然值，其中i表示目标蛋白上的氨基酸的序号，为1～N的任意整数，N表示目标蛋白上的氨基酸的总数，j表示矩阵中列的序号；

令把上述位置特异性打分矩阵中的每个值归一化，获得归一化后的位置特异性打分矩阵；

对上述归一化后的位置特异性打分矩阵，把每一列中处在序列不同位置处的相同氨基酸对应的值相加，获得第八特征向量(v177～v576)。

S2.目标特征向量的选择

把按照上述S1.编码的8种属性整合成一个综合特征向量集合，向量维度为576的向量(v1,v2,v3,...,v576)。对于RBP序列，在所述综合特征向量前添加“+1”标签；对于no-RBP序列，在所述综合特征向量前添加“-1”标签。为了减少计算量和特征之间可能存在的冗余性，我们以特征之间最小的冗余性和最大的相关性(mRMR法)为标准，选择了300个特征向量作为目标特征向量；

S3.预测模型的构建

图3是按照本发明的预测方法的流程示意图，以下结合图3进行详细说明。

S31.训练集的样本蛋白的选取

在本实施例中，从UniProt蛋白质数据库中选取RBP序列，序列一致性阈值为25％，利用聚类工具对这些检索出的RBP序列去冗余，得到2780个RBP序列，从而作为训练集的阳性样本。同时，从蛋白质结构数据库PDB出发，提取出X射线分辨率0.0～3.0A,序列一致性小于或等于25％，序列长度50～10000氨基酸的由X射线解出的蛋白质链，从中去除与RNA结合或PDB记录中包含“核糖体的”、“未知功能”、“RNA”、“核蛋白质”的蛋白质链，最终获得7093个不与核糖核酸结合的非冗余蛋白质序列(non-RBP)作为训练集的阴性样本；

S32.支持向量机法构建预测模型

支持向量机的核心思想是通过用内积定义的核函数k(x_i,x_i)进行非线性映射，经非线性映射将输入数据的空间变换射到一个高维空间，在该高维空间中寻找最优分类超平面，使不同类别的数据线性可分，该方法基于结构风险最小化原则，具有良好的泛化能力，保证算法的全局最优性的优点。在本实施例中，采用了LIBSVM软件进行这一步骤。

S321.最优分类超平面

最优分类超平面需满足两个条件：

1)分类间隔最大；

2)对所有样本尽可能多地分类正确；

以2维的情况示例：上述最优分类超平面对应于最优分类线，所述最优分类线是指该分类线不仅能将两类样本尽可能地分开，且满足分类距离最大；所述分类间隔是指位于分类线两侧且与分类线平行的两条直线的距离，上述两条直线分别是指过第一和第二类样本中离分类线最近的样本且平行于分类线的直线；

S322.最优分类超平面的求解

求解最优分类超平面可转化为如下约束问题：

设训练集样本特征向量为x_ki∈Rⁿ，属于两类，分类向量y∈R^m，且y_k∈{1,-1}，k＝1,2,…m。其中x_ki表示训练集中第k条蛋白质的第i个特征向量，对应于上述综合特征向量集合；y_k表示训练集中第k条蛋白质的标签，即“+1”或“-1”；Rⁿ表示由训练集样本蛋白质的n维特征向量的全体构成的向量空间(在本实施例中，n＝300)；R^m表示m维向量的全体构成的向量空间，m为训练集的样本蛋白的数量。

满足上述的条件1，分类间隔最大即：

取最小值：

满足上述的条件2，对所有样本尽可能多地分类正确，即尽可能多的x_k满足：

ε_k≥0，k＝1,2,…m

其中，ω为最优分类超平面的法向量，b为最优分类超平面的偏移量，ε_k为松弛项，

在本实例中，采用径向基函数作为核函数，径向基函数如下：

其中，为映射函数，γ为核参数，C为错分的惩罚参数，C>0，γ>0；

以上述条件1和条件2为标准，得到最优解，即用于分类的决策函数：

故在此分类向量机中，有两个参数C和γ需要学习优化，以用于模型的构建。

S323.参数C和γ的优化

本实施例采用十倍交叉验证的方式来优化参数C和γ的值。所谓十倍交叉验证，是指整个训练集被随机分成10等份，每次选择其中9份作为训练集，剩余1份作为测试集，重复10次，以保证每一份都分别作为测试集一次。这样的方式，能够有效避免过度训练，参数的优化过程如下：

1)设定一个初步的搜索空间范围，log₂C∈{-6,-4,…,20}，以2为间隔；log₂γ∈{6,4,…,-20}，以-2为间隔，共196(14x 14)对(C,γ)；

2)针对上述搜索空间，对于每对(C,γ),对训练集做十倍交叉验证，计算十个子测试集的平均准确度，并重新将(C,γ)值设置为十倍交叉验证准确度最高的那对(C,γ)值；

3)分别以步骤2)选出的log₂C和log₂γ为中心，设置搜索空间为log₂C∈{log₂C－2,log₂C,log₂C+2}，log₂γ∈{log₂γ－2,log₂γ,log₂γ+2}，以0.5为间隔为新的搜索空间，共81对(C,γ)”返回步骤2)；

4)比较步骤3)获得的新的十倍交叉验证的精确度与上一次搜索获得的精确度的差值，若差值小于或等于0，优化结束；否则重新返回步骤2)；在本实施例中，c＝185363.800047，γ＝0.000690533966002。

S4.获取目标靶蛋白的目标特征向量，并代入所述步骤S3.中建立的预测模型，获得目标靶蛋白的预测结果；当决策函数f(x)结果为+1表示该蛋白是RBP,-1表示该蛋白是non-RBP。

本实施例采用了更全面的评估指标，包括灵敏度、特异性、准确度、马修相关系数，受试者操作特性曲线下的面积来对预测结果进行评估。所述评估指标的定义如下：

灵敏度是指阳性数据(即RBP)中被正确预测为阳性的比例；

特异性是指阴性数据(即non-RBP)中被正确预测为阴性的比例；

准确度是指阳性和阴性数据中被正确预测的比例；

马修相关系数

其中，真阳性TP指的是阳性数据中被正确预测为阳性的个数；真阴性TN是阴性数据中被正确预测为阴性的个数；假阳性FP是阴性数据中被预测为阳性的个数；假阴性FN是阳性数据被预测为阴性的个数；

受试者操作特性曲线下的面积是以灵敏度为纵坐标，特异性为横坐标绘制的曲线，曲线与x轴围成的面积是一种评估预测方法结果的更客观的指标。

根据上述评估指标，在步骤S3中所用的训练集样本蛋白质的2780个RBP和7093个non-RBP上进行十倍交叉验证。本实施例的预测结果具有高灵敏度和特异性，取得的结果如下(均以每个指标的平均值表示)：灵敏度为83.07％，特异性为96.00％，准确度为92.36％，马修相关系数为0.808，受试者操作特性曲线下的面积为0.975。

实施例2

按照上述S4的步骤，把人类蛋白质数据集作为目标蛋白进行测试，所述人类蛋白质数据集包括967个RBP和579个非冗余non-RBP。预测结果如下：967个RBP中的84％被正确预测，597个non-RBP中的97％被正确预测，马修相关系数为0.788。

实施例3

以所述的相同步骤重复实施例2，区别在于，所述目标蛋白为酵母的蛋白质数据集，相应的马修相关系数为0.729。

实施例4

以所述的相同步骤重复实施例2，区别在于，所述目标蛋白为拟南芥的蛋白质数据集，相应的马修相关系数为0.537。

对比例1

以所述的相同步骤重复实施例1，区别在于，步骤S2中选取的特征的数量为100个；在所述步骤S4中相同的2780个RBP和7093个non-RBP上进行十倍交叉验证，相应的预测结果为：马修相关系数为0.777，受试者操作特性曲线下的面积为0.965。

对比例2

以所述的相同步骤重复实施例1，区别在于，步骤S2中选取的特征的数量为200个；在所述步骤S4中相同的2780个RBP和7093个non-RBP上进行十倍交叉验证，相应的预测结果为：马修相关系数为0.795，受试者操作特性曲线下的面积为0.970。

对比例3

以所述的相同步骤重复实施例1，区别在于，不经过步骤S2而直接进入步骤S3，所述步骤S3中的c＝46340.9500118，γ＝0.000345266983001；在所述步骤S4中相同的2780个RBP和7093个non-RBP上进行十倍交叉验证，马修相关系数为0.814，受试者操作特性曲线下的面积为0.975。

对比例4

采用SPOT-seq方法对实施例2中所用的人类蛋白质数据集进行相同的预测，结果表明，967个RBP中只有35％的被正确预测；597个non-RBP中有94％的被正确预测，准确度为57％，马修相关系数为0.330，与实施例2相比，仅为实施例2的一半不到。可知，本发明的评估指标与现有技术相比，有了显著提升。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (6)

1.一种与核糖核酸结合的蛋白质的预测方法，其特征在于，包括以下步骤：

S1.根据氨基酸的疏水性，获取第一候选特征的集合；根据氨基酸的极性，获得第二候选特征的集合；根据氨基酸的规范化的范德瓦尔斯体积，获取第三候选特征的集合；根据氨基酸的极化性，获取第四候选特征的集合；根据蛋白质的二级结构，获取第五候选特征的集合；根据氨基酸的溶剂可及性，获取第六候选特征的集合；根据氨基酸的带电性和极性，获取第七候选特征的集合；根据氨基酸的位置特异性打分矩阵，获取第八候选特征的集合；

根据第一候选特征的集合至第八候选特征的集合的并集，获得样本蛋白的候选特征总集V_i＝{v_i,k}；其中，v_i,k表示第i个样本蛋白的第k个候选特征，i为1～ξ的任意整数，k为1～N₁的任意整数，ξ为样本蛋白的总数量，N₁为候选特征总集中候选特征的总数量，ξ≥5N₁；

S2.根据候选特征之间的相关度和/或冗余度，在所述第i个样本蛋白的候选特征总集中选取个候选特征作为第i个样本蛋白的最佳特征，其分别为并获得第i个样本蛋白的特征向量 其中，

S3.根据样本蛋白的特征向量以及蛋白质特性，建立预测模型f(x_i)，并获取样品蛋白在所述预测模型中的评估指标σ，如果σ≥阈值η，进入S4，否则返回步骤S2；

S4.获取待预测蛋白与样本蛋白对应的特征向量x，代入所述步骤S3中所建立的预测模型，获得待预测蛋白的蛋白质特性的预测结果f(x)，从而判断该待预测蛋白是属于核糖核酸结合蛋白，还是属于非核糖核酸结合蛋白。

2.如权利要求1所述的预测方法，其特征在于，所述步骤S2中选取候选特征的方法为最大相关最小冗余法。

3.如权利要求1所述的预测方法，其特征在于，所述步骤S3中建立预测模型的方法为支持向量机法、神经网络法、贝叶斯分类法或随机森林法。

4.如权利要求3所述的预测方法，其特征在于，所述步骤S3中建立预测模型的方法为支持向量机法。

5.如权利要求1所述的预测方法，其特征在于，所述步骤S3中的评估指标为准确度、马修相关系数或受试者操作特性曲线下的面积。

6.如权利要求5所述的预测方法，其特征在于，所述步骤S3中的评估指标为马修相关系数，所述阈值η为0.5。