CN101231677A - 基于序列谱中心和遗传优化处理的残基间的远程相互作用预测方法 - Google Patents

Abstract

本发明涉及基于序列谱中心和遗传优化处理的残基间的远程相互作用预测方法，每一个谱中心就是一个表示某类所有残基对的平均序列谱编码向量，所有残基间的空间距离小于8或者是从8到10的残基对组成一类，平均每类里的所有残基对序列谱向量就得到了这个类的谱中心。构建了一个基于遗传算法的多分类器来计算残基对的序列谱向量与序列谱中心之间的距离，并由此判断此残基对的两个残基是否处于接触状态。分析结果显示，44.48％的远程接触对位于序列谱中心的周围，分类器也得到了20.9％的正确率，残基对的序列上的距离至少为24个残基，并且选择length/5个最好的预测输出来计算预测器的性能指标。

Description

基于序列谱中心和遗传优化处理的残基间的远程相互作用预测方法

所属领域  本发明涉及生物信息学或者结构生物学领域，特别涉及基于序列谱中心和遗传优化处理的残基间的远程相互作用预测方法。

背景技术  在人类基因组计划执行以来，越来越多的蛋白质序列被高通量的测定出来。然后，由于X-ray晶体衍射技术以及核磁共振等技术来确定蛋白质的结构是非常困难而且极为耗时的，所以确定蛋白质空间结构的速度也就极其缓慢。然而，我们只有知道蛋白质的空间结构才能够进一步了解它们的功能并进行其它的分子设计。这样的话，由于蛋白质结构的确定速度远远低于其序列的测定速度，所以利用计算机技术来进行蛋白质结构的预测也就应运而生。

在目前现有的计算方法中，大致有三类方法来进行蛋白质结构的预测研究，它们是比较建模方法、反向折叠或者所谓的串线方法以及从头计算方法。不过到目前为止，它们仍然不能完全解决蛋白质结构预测这种极度复杂的问题。因此，我们有必要去寻找或者改进有效的计算方法来预测蛋白质的三维空间结构。在这些方法或技术路线中，蛋白质接触图谱的预测是一类相对简单的研究蛋白质3D结构的手段。

众所周知，蛋白质的接触图谱是蛋白质三级结构在二维平面上的一种表达形式，它可以被用来提高蛋白质的三级结构预测性能。我们知道，蛋白质残基间的非局部相互作用对于稳定蛋白质的自然状态是至关重要的。Gromiha曾经指出，残基间的远程相互作用在蛋白质的折叠和结构的稳定中发挥着重要的作用，预测出它们也就能够解析蛋白质的三维结构。如果我们能够知道蛋白质残基间的接触状况以及其二级结构的相关知识，我们就可以得出其三级结构的许多重要特征。更为重要的是，我们甚至可以把一个有错的接触图谱投射到其最接近的物理上允许的蛋白质构像上，从而复原其三维结构。最后一点需要指出的是，对于序列长度小于150个残基的蛋白质链，有结果显示至少具有50％正确率的残基间的接触图谱(阈值距离为8)才足以重建此蛋白链。

目前，许多的工作都集中在利用不同的计算方法来预测蛋白质的接触图谱。主要的方法或手段有：利用进化信息的、遗传规划集成了自组织映射、神经网络等等。然而，到目前为止，没有一种方法能够有效地预测不同类型的蛋白质的接触图谱或者长程相互作用。当考虑残基对的序列上的距离至少为6个残基的时候，PR0Fcon方法的预测精确度为30％，而其覆盖率为10％。虽然PR0Fcon方法获得了相对较高的预测精确度，它的覆盖率也仅为10％。当预测蛋白质残基间的远程相互作用的时候，Vullo的二段预测器可以得到19.8％的预测精度(残基对的序列上的距离至少为2 3个残基，并且选择length/5个最好的预测输出来计算预测器的性能指标，其中length为蛋白质链的序列长度)。综上所述，现有的接触图谱的预测方法还是不能满足其蛋白质三维结构的重建需要，那么我们就有必要发展新的思路和方法来预测蛋白质的远程相互作用。

发明内容  本发明的目的是分析残基对的序列谱中心，并在此基础上提供一种用于预测蛋白质残基间的远程相互作用的有效方法，并由此预测蛋白质的三维空间结构。

本发明是这样实现的：

(1)蛋白质的训练和测试序列

我们利用PDB-REPRDB数据库来抽取蛋白质序列，此数据库是是基于版本号为2005_05_29的PDB数据集的。这些蛋白质由X-ray晶体学所测定结构的分辨率要小于2.0，并且任意两条蛋白质序列的相似性要小于25％。然后我们剔除没有文件的蛋白质链，最后我们就得到了480条蛋白质链。

为了验证我们所提出的方法，这480条蛋白质多肽链被分为两个子集，训练集拥有287条多肽链而测试集有193条多肽链。对于每一个测试的多肽链，我们的方法都是运行两次最后得到了平均的输出性能。

(2)蛋白质的序列谱

为了得到GaMC预测器的编码向量，按照前面所述我们从HSSP数据库抽取了残基的序列谱。我们已经知道，对每个在PDB数据库中的三维结构已知的蛋白质来说，HSSP数据库中都有一个唯一对应的文件，这个文件中包含了这个蛋白质的所有同源序列。HSSP数据库中的蛋白质序列谱包含了2 0种氨基酸出现在多序列比对中相应列的频率。当某一个位置上的某一氨基酸的值为100时，说明在这个位置上只出现这一种氨基酸，也就是说，这个位置非常保守。

由于我们要进行蛋白质残基间的远程相互作用，那么我们首先就要对每对残基i和j进行编码，并以每个残基为中心向蛋白质的N端和C端同时延伸以得到一个残基窗口，这样就得到两对以残基i和j为中心的残基窗口，在这里我们设定窗口的宽度为9(我们的实验结果显示，窗口宽度9可以得到最佳的性能)。在文献运用了第三个中心窗口并证实此连接残基i和j的中心窗口(其中心残基位置为int((i+j)/2))可以提高Contacts的预测性能。在这里，我们也运用了这个中心窗口，并取窗口的宽度为5并且得到了较好的预测结果。

为了构建用于残基对i和j间远程相互作用预测的输入向量，三个窗口共计包括了23个残基，结果训练和测试向量是一个(9+5+9)*1即460*1的列矩阵。

(3)蛋白质的Multi-Class Contact定义

一般的Contact Map的定义即是一个简单的二类样本的形式。由于在蛋白质的3D空间结构上，远程接触的残基对极少，所以为了避免在机器学习算法训练时出现的严重的正负输入样本的不平衡性，我们设计了一个新的接触图谱(Contact Map，CM)的定义，这个残基间的远程相互作用的多类接触的定义可如下式所示：

$S_{\mathrm{ij}}=\left\{\begin{array}{ccc}0& \mathrm{ifd}(i,j)\&le;d& \\ 1& d<d(i,j)\&le;d_1& \\ 2& d_1<d(i,j)\&le;d_2& ,|i-j|\&GreaterEqual;24\\ .& & \\ .& ...& \\ .& & \\ n-1& d(i,j)>d_{n-2}& \end{array}\right.---\left(1\right)$

在这个新的表达式里，我们仍然把那些在3D空间上距离小于或等于8(即d)并且在序列上的距离24个残基的残基对赋予0类，即这个残基对是接触的。为了对剩下的残基对归类，我们可以做相似的分类。我们以蛋白质多肽链PDB：1hg7A为例，此多肽链包含65个残基。考虑远程相互作用的话，就可以得到861个不重复的远程残基对，经过计算它们的空间距离，我们知道其中有82个残基对属于0类，即他们是远程接触的。这样的话，就剩余779个非远程接触对。按空间距离把这779个非远程接触对平均分配为4类，较近的分配为1类(即距离为d和d₁之间)，以此类推。在这里，我们把远程残基对归属为5类，即n＝5，这时我们可以得到较好的分类性能。结果，我们得到194个1类的远程残基对；2、3、4类含有195个远程残基对。这样的话，我们就可以把多类接触谱应用到一个多类分类器系统以得到比两类分类器更好的分类率。

(4)序列谱中心分析

为了确定残基对是不是处于远程接触的状态，即此残基对的样本向量是不是属于0类，我们拟计算此残基对的序列谱样本向量和每个序列谱的中心的距离。如果此序列谱样本向量与某个序列谱中心最近，那么就把此残基对归属于此序列谱中心所代表类别。按照上面多类Contact Map的定义，我们得到相应的n个序列谱中心。序列谱中心的定义在以前的文献中也有所涉及。而我们的谱中心的定义只是简单的计算每个远程接触类的平均序列谱样本向量。那么，序列谱中心Ci就可以给定为：

$C_i=\frac{1}{m_i}\underset{l=1}{\overset{m_i}{\mathrm{\&Sigma;}}}S\left(l\right)---\left(2\right)$

这里，S(*)指以残基对的序列所编码的一个列向量，其维数为460，m_i指远程残基类i的数量。

(5)基于遗传算法和柱状图编码的分类器

本专利的远程接触研究采用了基于遗传算法的优化器。首先，我们利用遗传算法在不损失分类率的基础上进行输入向量的维数压缩，然后分析残基间远程相互作用的分布。

我们知道，遗传算法(Genetic Algorithm，GA)是一种受生物进化启发的学习方法，是一个随机搜索和优化的技术，它所研究的问题是搜索候选假设空间并确定最佳假设，也就是GA基于一个适应度函数来建立问题可能解的种群并从中搜索到其中具有最小或最大适应度值的个体。为了搜寻到问题的最优解，GA从一系列的假设解(染色体)开始，并且按照一定的重复顺序对不同的但较佳的假设解进行进化。在每一代(重复次数)中，目标函数(适应度度量准则)决定了每一个假设解的适应程度，并且基于这些函数值，对当前的假设解集合(空间)进行遗传操作(选择、复制、交叉和变异)来产生下一代。通常情况下，即使针对的是一个搜索难度大的空间，遗传算法也能够快速地搜索到局部最优解。

下面我们就对GA算法最关键的染色体编码方式，染色体适应度函数的定义和遗传方式的操作等方面，来一一进行详细的描述。

1.染色体编码

为了开始我们的远程接触研究，首先我们来作一些假设。假设V是样本向量空间V＝{v₁，v2，…，v_m}的一个初始的离散随机变量，假设有N个样本向量X_n，n＝1，2，...N要依据它们的样本变量所构成的函数f(V_n)从而被分到K个类C₁，C₂，...C_K中。我们的目的是搜索一个最优的非线性变换W_s＝T(V_n)，使得样本向量被变换后更容易分类而且分类率最大，这时我们也可以得到相应的训练集的样本变量函数是f_n(W_s).为了达到我们的目的，GA就被用来在样本空间搜索得到可能的变换即问题的可能解。我们把样本的所有变量特征依次分成数个间隔

这里

是属于集合V的一个间隔。这样，每个可能的变换在GA中就被表达为一个字符串S^m，m＝1，2，...M，其中M是GA中所选的种群的大小也就是集合V的势大小。在遗传算法中染色体编码的方式有很多，最常用的是二进制编码和实数编码。在这里，染色体用一个三字符的字符串来编码，即一个三元字符组Z＝{a，b，c，}用来表达这些字符串，如附图1所示，字符串的长度即是分类器输入列向量的维数(即460位)。其中，字符‘a’表示它所代表的样本向量的相应位置的值从样本空间V中删除，而连续的字符‘b’和‘c’表示相应位置的值被合并到一起。

在这里，每一条染色体都对应了一个残基对样本所表达的列向量的变换。在这个变换过程中，如果在染色体中的连续字符是相同的话，那么在初始向量中的相应位置就被合到一起，而且它们的值相加成一个数值。染色体字符的操作规则如下：a).当染色体字符为a的时候，就删除此染色体字符所代表的特征变量；b).当染色体字符为b且有连续相同的b的时候，叠加所有字符b所代表的特征变量并删去被叠加的特征变量。如果相邻仅有一个b字符，则保留此字符所代表的特征变量；c).当染色体字符为c时，按b)规则处理。这里为了便于对变换过程的理解，我们假设研究一个只有8维的残基对样本，相应的每一条染色体就是一个由三字符组成的8位的字符串。然后，我们假设某个残基对样本向量是(v1，v2，v3，v4，v5，v6，v7，v8)，在当前的种群代数中假设有一条染色体为cbba bcca，则这条染色体所对应的残基对样本向量的特征变量的变换过程详见附图2所示。此图中，此染色体每个个体所代表的样本向量被删除或合并就得到了一个新的归一化的向量，此时样本向量的维数就有8维变换成了4维，这时新的样本向量就被用来进行我们的远程接触分类。

2.适应度函数的定义

在定义了染色体编码以后，我们就开始描述蛋白质残基对样本向量的特征变换，并给出在变换过程中的适应度函数的定义以搜索得到最优的特征变换。我们所采用的残基对样本向量的特征变换，从其生物意义上，我们是为了更好地考虑去处多余的样本向量特征对蛋白质残基间的远程相互作用的影响；从计算角度来看，我们是为了追求一个较低的维数和一个较高的分类精度。显而易见，在维数降低到一定程度时，会与分类精度产生矛盾。如何在维数与分类精度中进行性能的平衡，已经被很多研究者研究，在这里就不加详述了。在这里，我们规定当重复80代以后如果适应度函数值仍稳定在一个值，我们就认为本次的样本向量的特征变换就趋于稳定并退出进行其他的优化计算。

在GA算法的每一代，按照所有染色体所对应的变换，我们将样本向量的特征或变量重新组合，并把这些变换后的蛋白质残基对的样本向量作为GaMC预测器的输入向量。在我们的实验中，使用所有的正样本和相同数目从负样本集中随机抽取的负样本来进行预测器的训练以平衡预测器的样本空间。然后数据集被随机的分成相等大小的两个子集，用这两个子集中的每一个进行训练，并在它的补集上进行测试。预测器的分类精度就可以用这试验的平均分类性能来进行定义。

这样的话，对于每一种变换T^m和每一个与字符串S^m应的间隔集

合我们就可以构造对于样本向量X_i的每一个间隔的样本变量函数。对于C_k类，我们有如下的基于T^m和

的初始样本变量函数定义：

在这里，|C_k|是类C_k的势。

给定这些初始的样本变量函数，我们就可以得到一个新的关于样本向量x的远程相互接触类的结构C_l，l＝1，2，...K：

$x\&Element;C_l\mathrm{ifd}(f,\stackrel{\&OverBar;}{f^k})\&le;d(f,\stackrel{\&OverBar;}{f^l}),l=\mathrm{1,2},...,K---\left(4\right)$

这里，f指未知远程接触类的样本向量x的样本变量函数，而d(·，·)指两个样本向量的变量函数的距离或相异值的一种度量。

我们知道，遗传算法的适应度函数

一般被用来度量基于字符串S^m的初始类的结构C1，C2，...Ck和转化后的类的结构C₁ ^m，C₂ ^m，C_K ^m的差异大小，GA的适应度函数定义如下：

${\&aleph;}^m=\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}|C_k\&cap;C_k^m|---\left(5\right)$

当两个类的结构正好相一致的时候，就可以得到的最大值；当它们之间的差异变大时，适应度函数值就相应的变小。

然而对于每一个远程接触类，在GA中可能不存在一个全局最优的样本变量的转换和间隔子集，所以我们把上述的GA方法扩展为一个多分类器系统并且对系统的每个类用一个特别的子分类器来训练和测试。扩展思想就是：对于每一个远程接触类，可以利用GA开发出K个子分类器，其目的就是每一个子分类器都搜索最优的变换和间隔子集。为了完成上述任务，归一化的样本变量函数作为第k个分类器的输入，这里v_t ^m，k和分别是此子分类器的最优变换和间隔子集。此外，由于一个子分类器所带来的错误可能被其他的子分类器所更正，所以融合所有这些独立的子分类器可以提高我们的分类率。

我们知道，第k个子分类器用来识别一个未知类别的样本是不是属于k类，因此我们可以定义它的相应的适应度函数

从而得到第k个子分类器的输出K_k。考虑此输出K_k，如果K_k等于1，那么此样本属于第k类；如果GA判定如果K_k等于0，那么此样本不属于第k类。这时多分类器系统的适应度函数定义为：

${\&aleph;}_k^m=\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}|C_k\&cap;C_k^m|+\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}|\stackrel{\&OverBar;}{C_k}\&cap;\stackrel{\&OverBar;}{C_k^m}|.---\left(6\right)$

这里，

和

分别表示C_k和C_k ^m的补集合。

公式(6)的第一项计算分类器输出的正确的正样本的总数，也即是属于k类实际上也被分类器分到k类的样本的数量；而其第二项计算分类器输出的正确的负样本的总数，也即是不属于k类实际上也没有被分类器分到k类的样本的数量。

基于公式(4)～(6)，我们可以把每一个序列谱输入样本赋予一个相应的序列谱中心所对应的类，即一个序列谱输入样本可能是属于类0，也可能是类1、2、3或4的。也就是说，如果一个序列谱样本向量S和序列谱中心C_k之间的距离要小于S与其他序列谱中心C_l(这里，l＝1，...K，and l≠k)间的距离的话，那么S就被赋予一个序列谱中心对应的类k。即，对应公式(4)，我们有，

$x\&Element;C_k\mathrm{ifd}(f\left(s\right),\stackrel{\&OverBar;}{{f\left(s\right)}^k})\&le;d(f\left(s\right),\stackrel{\&OverBar;}{{f\left(s\right)}^l}),l=\mathrm{1,2},...,K$

同样的，在这里f(s)指未知远程接触类的样本向量s的变量函数，而d(·，·)指两个样本向量的变量函数的距离或相异值的一种度量。

3.GA参数的选择

在我们所使用的遗传算法中，GA的选择操作使用的是轮盘赌(roulette wheel)方式。使用这种方式来决定每个个体的选择份数，能够使最优秀的个体获得最多的生存繁殖机会。在选择操作后，我们使用了一个简单的交叉变异方式，即使用单点交叉方式，交叉概率定为0.95，变异概率定为0.01。

本发明的创新之处在于：

1、对残基对的序列谱中心的定义和分析。

本方法给出了一个全新的残基接触对的定义，并基于此定义分析了每个类的序列谱中心，以这些序列谱中心为基础构建了一个全新的预测器。多类残基接触对的定义可以降低数据集的正负不平衡性，从而增加分类的正确率。

2、基于遗传算法和序列谱中心的远程残基接触对的预测方法。

本方法首先把数据经过遗传算法的优化处理后，得到了数据维数的压缩，降低了计算的复杂性和数据处理的复杂性。然后根据定义的序列谱中心得到一个距离分类器，距离这些序列谱中心最近的序列谱(或者样本变量函数)就被赋予这个类，从而达到了分类的目的。

附图说明  下面结合附图对本发明作进一步的说明。

图1是染色体的编码方式。

图2是蛋白质残基对样本向量特征的变换过程。初始归一化向量的每个柱状图表示一个特征或变量，而柱状图的高低与变量的大小相一致。为了便于计算，在变量变换以后的样本向量也要被归一化。

图3是输入向量转换图，其中(a)子图表示一个初始输入样本向量，(b)子图表示类1分类器转换后的的输入样本向量。为了使输入样本向量的变量范围相等，这两种样本向量都经过了归一化的处理。

图4是蛋白质的残基间的远程接触对的数量与相应的蛋白质链序列长度的关系图。其中，方框表示训练蛋白质链而‘*’表示测试蛋白质链。

图5是表示远程接触对的氨基酸组成。每个柱状图表示在各个子集合中对应的氨基酸比率。

图6是GA分类器的精度。对应5种CATH类型，图中显示了对应的5组蛋白质残基队的远程接触分类精度，并且每组精度值包括了6个精度，也就是对于选择最佳2L、L、L/2、L/5、L/10、和L/20个数量的输出分类的精度值。其中，每个柱状图的字符“I”表示标准偏差。

图7是蛋白质链1hg7A的接触图谱和预测的远程接触对。(a)图表示其接触图谱和考虑2L输出数量的预测的接触图谱比较图；(b)，(c)和(d)分别描述了类似的当考虑L，L/2，L/5输出数量的预测的比较图。其中，每个子图的两个斜线表示了序列上残基对的两个残基相隔24个残基的分割线，而左上角和右下角的点分别表示实际的远程接触对和预测的接触对。

表1是GA的选择参数。

表2是GA预测器的平均精确度。

具体实施方式

采用PDB-REPRDB数据库里抽取的数据进行训练和测试。把获得的480条蛋白质多肽链被分为两个子集，训练集拥有287条多肽链而测试集有193条多肽链。对于每一个测试的多肽链，我们的方法都是运行两次最后得到了平均的输出性能。

本专利里，我们分析了氨基酸残基的序列谱中心并提出了一个新的思路来解决蛋白质接触图谱的问题。我们的出发点是基于此表示平均序列谱的序列谱中心的，每一个谱中心就是一个表示某类所有残基对的平均序列谱编码向量，比如说所有残基间的空间距离小于8或者是从8到10的残基对组成一类，平均每类里的所有残基对序列谱向量就得到了这个类的谱中心。然后，我们深入的探讨我们的多分类器能否构建和学习蛋白质的序列信息与其相应的接触图谱间的联系。下面，我们就简单介绍一下本专利的内容。首先，我们对残基对进行序列谱的编码；其次，我们计算序列谱中心。在这些准备工作完成以后，我们构建了一个基于遗传算法的多分类器(GaMC)来计算残基对的序列谱向量与序列谱中心之间的距离，并由此判断此残基对的两个残基之间是否处于接触状态。我们利用GA对每一个子分类器来搜索输入样本向量的最优变换形式。结果，再用分类器融合方法来合并多个子分类器的输出。我们的分析结果显示，44.48％的远程接触对位于序列谱中心的周围，我们的分类器也得到了20.9％的正确率(残基对的序列上的距离至少为24个残基，并且选择length/5个最好的预测输出来计算预测器的性能指标)。

Claims (1)

1.一种基于序列谱中心和遗传优化处理的残基间的远程相互作用预测方法，其特征是：

a、蛋白质的训练和测试序列

首先利用数据库来抽取蛋白质序列，此数据库是是基于版本号为2005_05_29的PDB数据集的，这些蛋白质由X-ray晶体学所测定结构的分辨率要小于2.0，并且任意两条蛋白质序列的相似性要小于25％，然后我们剔除没有文件的蛋白质链，得到480条蛋白质链；

480条蛋白质多肽链被分为两个子集，训练集拥有287条多肽链而测试集有193条多肽链，对于每一个测试的多肽链，方法都是运行两次最后得到了平均的输出性能；

b、蛋白质的序列谱

GaMC预测器的编码向量从HSSP数据库抽取了残基的序列谱，对每个在PDB数据库中的三维结构已知的蛋白质来说，HSSP数据库中都有一个唯一对应的文件，这个文件中包含了这个蛋白质的所有同源序列，HSSP数据库中的蛋白质序列谱包含了20种氨基酸出现在多序列比对中相应列的频率，当某一个位置上的某一氨基酸的值为100时，说明在这个位置上只出现这一种氨基酸；

对每对残基i和j进行编码，并以每个残基为中心向蛋白质的N端和C端同时延伸以得到一个残基窗口，这样就得到两对以残基i和j为中心的残基窗口，设定窗口的宽度为5或9；

为了构建用于残基对i和j间远程相互作用预测的输入向量，三个窗口共计包括了23个残基，结果训练和测试向量是一个(9+5+9)*1即460*1的列矩阵；

c、蛋白质的Multi-Class Contact定义

设计一个新的接触图谱的定义，这个残基间的远程相互作用的多类接触的定义可如下式所示：

$S_{\mathrm{ij}}=\left\{\begin{array}{ccc}0& \mathrm{ifd}(i,j)\&le;d& \\ 1& d<d(i,j)\&le;d_1& \\ 2& d_1<d(i,j)\&le;d_2& ,|i-j|\&GreaterEqual;24\\ .& & \\ .& ...& \\ .& & \\ n-1& d(i,j)>d_{n-2}& \end{array}\right.---\left(1\right)$

d、序列谱中心分析

拟计算此残基对的序列谱样本向量和每个序列谱的中心的距离，如果此序列谱样本向量与某个序列谱中心最近，那么就把此残基对归属于此序列谱中心所代表类别，按照上面多类Contact Map的定义，得到相应的n个序列谱中心；

e、基于遗传算法和柱状图编码的分类器

本专利的远程接触研究采用了基于遗传算法的优化器，利用遗传算法在不损失分类率的基础上进行输入向量的维数压缩，然后分析残基间远程相互作用的分布；

GA算法最关键的染色体编码方式，染色体适应度函数的定义和遗传方式的操作包括GA参数的选择、染色体编码和适应度函数的定义。

Images