专利名称：基于平均场退火技术的蛋白质的立体结构比对方法
技术领域：
本发明属于生物信息工程技术领域。
背景技术：
蛋白质的结构比对(alignment)在生物信息学属于最重要的计算问题之一，它对蛋白质结构比对和分类，功能预测，主次结构的识别，多基因树的重构，遗传医疗，制药工程等起着关键作用。
相对于两列字符的比较，一对蛋白质结构的比对要困难得多，因为这需要最优地比对两列具有空间结构的刚体，且找出它们最贴近重叠的区域。从计算的复杂性而言，与现有的动态规划比对一对序列的多项式时间算法[S.B.Needleman and C.D.Wunsch，“A general methodapplicable to the search for similarities in the amino acid sequence of twoproteins”，Journal of Molecular Biology，48(1971)，pp.443-453]相比较，一双蛋白质的结构比对则属于NP困难问题。
至目前为止，已出现很多有关的算法，如基于距离的迭代动态规划(又称为LUND迭代法)[M.Gersten and M.Levitt，“Using iterative dynamic programming to obtain accuratepairwise and multiple alignments of protein structures”，in Proceedings of the FourthInternational Conference on Intelligent Systems in Molecular Biology(ISMB)，1996，pp.59-67]，Monte Carlo算法[S.H.Bryant and S.F.Altschul，“Statistics ofsequence-structure threading”，Current Opinion in Structure Biology，5(1995)，pp.236-244]，平均场近似方法[R.Blankenbecler，M.Ohlsson，C.Pettern and M.Ringner，“Matching protein structures with fuzzy alignments”，Proceedings of the NationalAcademy of Sciences，100(2003)，pp.11936-11940]和退火技术[L.Holm and C.Scander，“Protein structure comparison by alignment of distance matrices”，Journal ofMolecular Biology，233(1993)，pp.123-138]。特别是，基于距离的LUND迭代法，这种方法在退火过程中主要包括两个步骤首先，对一双被排列的蛋白质表示成相应的Potts形式，并估计相应的模糊分配变量；然后，用所得模糊分配来加权蛋白质原子坐标间的距离，并计算精确的旋转和平移变换，迭代上述两步直至收敛为止。这种方法的缺陷是在迭代过程中需要调整有关参数并需要事后处理。此外，尽管R.Blankenbecler，M.Ohlsson，C.Pettern andM.Ringner基于LUND方法，并结合平均场理论和退火技术，提出了平均场近似方法，但仍需要事后处理。本质上，这种方法仍是一种近似方法，因为在迭代过程中需要引进很多额外的因子以保证迭代过程的收敛。总之，还没有一种最优的方法根本解决一双蛋白质的结构比较问题。
发明内容
基于平均场的退火技术，我们提出了一种全新的且理论上更严格的蛋白质结构比对方法。更确切地说，我们把蛋白质的立体结构比对精确地归结为一非线性整数规划问题。基于平均场理论，进一步把此问题等价地转化为非线性连续型的最优化问题，然后利用工程上的退火技术，高效率地找出相应问题的最优结构比对解。在这一过程中，不需要任何附加的或人为的″软限制″条件，所获得的最优解是精确的，高效率且不需要事后处理。相应地，我们的发明，不仅适用于任何一组蛋白质的结构比对，而且可应用于蛋白质的结构分类，功能预测，主次结构的识别，多基因树的重构，制药工程等领域。
从技术上说，假定两列待比对的蛋白质长度分别是nx和ny，其中nx，ny均为正整数，即表示此两列蛋白质链分别含有nx个原子和ny个原子。注意到我们不考虑侧链，只考虑主链(即Cα链)且每一个蛋白质被认为是具有刚性的几何体，则蛋白质X中的第i个原子与蛋白质Y中的第j个原子相应的坐标在三维空间中分别为Xi＝(xi，1，xi，2，xi，3)，Yj＝(yj，1，yj，2，yj，3)，(i＝1，...，nx；j＝1，...，ny)。蛋白质坐标可由通用的Protein Data Bank数据库得到。定义蛋白质X中的第i个原子与蛋白质Y中的第j个原子间的距离平方为dij2=|Xi-Yj|2=Σk=13(xi,k-yj,k)2---(1)]]>三维空间中的刚性几何体变换由平移变换和旋转变换组成，它一般可表示为X^i=A+RXi---(2)]]>这里A是平移向量，R为旋转矩阵，Xi和 分别表示链X中变换前后的第i个原子。
如图一所示，我们定义整数变量sij来描述两个蛋白质原子间的相配 也即，当两个原子相配时，sij为1，否则为零。
注意到在图一中，右边的相配矩阵S被两条虚线分割成三部分，其中横向虚线以上的行向量表示s0j，纵向虚线左边列向量表示si0，余下部分则表示sij。图中左边的白色链条表示蛋白质链X，黑色链条表示蛋白质链Y。可见，蛋白质链X中的第1、2、6、7、8个原子(从左至右)分别与Y中的第1、2、4、5、8个原子相配。相应地有，右边相配矩阵S中的s11＝1，s22＝1，s64＝1，s75＝1，s88＝1。而蛋白质链X中的第3、4、5个原子以及Y中的第3、6、7、9个原子都相配于间隙，则有s30＝1，s40＝1，s50＝1，s03＝1，s06＝1，s07＝1，s09＝1。并且S中的其他元素为0。
故我们根据之前定义的原子间距离、平移旋转变换以及相配矩阵，引进误差函数E(S，A，R)，来描述两列蛋白质的相配状况
E(S,A,R)=Ec(S,A,R)+λΣi=1nxsi0-λ2Σi=2nxs(i-1)0si0+λΣj=1nys0j-λ2Σj=2nys0(j-1)s0j---(6)]]>这里λ是一惩罚参数，λ取正的实数值，即是0＜λ≤100，其中Ec(S,A,R)=Σi=1nxΣj=1nysij|A+RXi-Yj|2---(7)]]>并且 项表示蛋白质X中所有相配于间隙的原子，而-λ2Σi=2nxs(i-1)0si0]]>项则描述了相配于间隙的情形连续出现的状态，即蛋白质X中相邻的两个原子均相配于间隙。同理地， 项以及-λ2Σj=2nys0(j-1)s0j]]>项也有类似的含义。
显然，E(S，A，R)是一个关于S、A和R的函数，我们的目标是在满足各种制约条件下，找出平移向量A和旋转矩阵R，使变换后得到的误差函数E最小。也就是说，只要各种制约条件得到满足并且误差函数E最小，所对应的两个蛋白质间的结构比对解则为最佳。这样得到的一双蛋白质的原子，它们之间相差的总距离也为最小。
虽然只要求解以上问题，就可得到一双最佳蛋白质的结构比对。但是，由于相配项sij是整数，以上问题是一个非线性整数规划问题，还没有很有效的解法。因此我们采用平均场的退火算法，该算法的出发点是基于物理中固体物质的退火过程与一般组合优化问题之间的相似性。在某一初温下，伴随温度参数的不断下降，在解空间中寻找目标函数的全局最优解。这是解决此类问题一个有效的求解算法，为了应用该技术，我们引进取值在0和1之间的实数值变量vij来代替整数变量sij。则误差函数变为E(V,A,R)=EC(V,A,R)+λΣi=1nxvi0-λ2Σi=2nxv(i-1)0vi0+λΣj=1nyv0j-λ2Σj=2nyv0(j-1)v0j---(8)]]>其中EC(V,A,R)=Σi=1nxΣj=1nyVij|A+RXi-Yj|2---(9)]]>此外，我们再引进自由能量函数F(V,A,R)=E(V,A,R)+TΣi=0nxΣj=0nyvijlnvij---(10)]]>这里T是退火过程中的温度。现在，我们把蛋白质的结构比对问题转化为非线性连续型的最优化问题。也即是，在满足各种制约条件下，找出平移向量A和旋转矩阵R，使自由能量函数F最小。
进一步，根据最优化问题的Karush-Kuhn-Tucker(KKT)条件，我们通过引入辅助变量uij和wj，将所要求解的相配变量vij表达成代数形式。这样，我们把一个非常复杂的蛋白质结构比对优化问题转换为简单的连续变量的代数问题。只要退火温度T充分小，这一代数问题的解就趋近于整数0或1，也即是得到收敛于最佳蛋白质的结构比对解。
本发明的特征在于它依次含有以下步骤(1)对计算机进行初始化● 设定以下各种初始值退火温度T，退火温度递降系数ε，惩罚函数λ，取正的实数值，0＜λ≤100，平移向量A为3阶零向量，旋转矩阵R为3×3单位矩阵，实数值相配变量vij(0)＝1/max(nx+1，ny+1)，0≤vij≤1，0≤i≤nx，0≤j≤ny，其中nx、ny分别为只计主链时，两列待比对的蛋白质X、Y长度，nx、ny均为正整数，i、j分别为X中的第i个原子和Y中的第j个原子，nx、ny也分别是X、Y中蛋白质原子的数目；vij的收敛阈值0.001；l是迭代次数，初始值为0；● 规范化数据从通用的Protein Data Bank数据库中导入待排蛋白质链{Xi0}i＝1nx、{Yj0}j＝1ny，则得到蛋白质X中的第i个原子与蛋白质Y中的第j个原子在三维空间中的相应坐标为Xi＝(xi，1，xi，2，xi，3)，Yj＝(yj，1，yj，2，yj，3)，把上述两个蛋白质X、Y平移到几何中心后Xi1=Xi0-cc,Yj1=Yj0-cc,]]>其中cc=(Σi=1nxXi0+Σj=1nyYj0)/(nx+ny),]]>再把Xi1、Yj1单位化，得到Xi=Xi1/fac,Yj=Yj1/fac,]]>其中fac=max1≤i≤nx,1≤j≤ny|Xi1-YJ1|;]]>(2)把要求解的相配变量vij用辅助变量uij、wj表达成代数形式，计算机用下式求解、更新uijuij(l)＝-|A+RXi-Yj|2/T(l-1)，1≤i≤nx，1≤j≤ny，u0j(l)=(-λ+λv02(l-1)/2)/T(l-1)j=1(-λ+λ(v0(j-1)(l-1)+v0(j+1)(l-1))/2)/T(l-1)2≤(-λ+λv0(ny-1)(l-1)/2)/T(l-1)j=nyj<ny,]]>
ui0(l)=(-λ+λv20(l-1)/2)/T(l-1)i-1(-λ+λ(V(i-1)0(l-1)+v(i+1)0(l-1))/2)/T(l-1)2≤(-λ+λv(nx-1)0(l-1)/2)/T(l-1)i=nxi<nx,]]>计算机通过迭代法对下式求解另一变量wjwj(l)=e-1+u0j(l)+Σi=1nx(euij(l)/Σj^=1ny(euij^(l)/wj^(l)))]]>其中，l为迭代到第l次的序号， 表示从1到ny的整数；(3)计算机用下式求解相配变量vij(l)=(euij(l)/wj(l))/Σj^=0ny(euij^(l)/wj^(l))]]>V0j(l)=e-1+u0j(l)/wj(l)]]>(4)根据下述收敛判断准则判别vij是否收敛(1/nxny)Σi=1nxΣj=1ny|vij(l)-vij(l-1)|≤0.001]]>(5)若不收敛，回到步骤(2)，进行下一次迭代，若收敛，则执行下一步骤；(6)计算机根据奇异值分解法用下面式子计算旋转矩阵R和平移向量AR=V0U0T]]>A=Σj=1ny(Σi=1nxvij/vtot)Yj-RΣi=1nx(Σj=1nyvij/vtot)Xi]]>其中vtot=Σi=1nxΣj=1nyvij,]]>且矩阵V0、U0满足下式 即以上等式的右边是左边的SVD分解，p表示从1到nx的整数，q表示从1到ny的整数；(7)计算机根据下式计算新的坐标X^i=A+RX]]>(8)降低退火温度T(l)＝T(l-1)×ε(9)判断相配变量vij是否趋近于整数即Σi=1nxΣj=1nyvij2/nx∈[0.99,1]]]>(10)若相配变量不趋近于整数，则返回步骤(2)，进行再次迭代；若相配变量趋近于整数，则进行下一步；
(11)对vij四舍五入，并令sij＝vij，按下式计算蛋白质的相配数m和平均距离rmsm=Σi=1nxΣj=1nysij,]]>rms=fac1mΣi=1nxΣj=1nysij|Xi-Yj|2]]>我们的发明，不仅能高效率地找出一双蛋白质相应的最优结构比对解，而且所获得的最优解是精确的，不需要事后处理。除了适用于任何一双蛋白质的结构比对，还可应用于蛋白质的结构分类，功能预测，主次结构的识别，多基因树的重构，制药工程等领域。
表1为部分数值模拟结果，并列举了CE法和Lund法的结果，以便比较。所有的蛋白质数据来自Protein Data Bank(PDB)。很显然，我们的发明优于CE法和Lund法。
表1蛋白质组相应得最优解构排定解的比较(rms两个蛋白质间的距离，m两个蛋白质相配原子个数，δ＝λ/2)
还原酶双氢原子还原酶由该表可知，在Reductases系和Globins系的蛋白质中，我们的算法都取得了与其他算法一致的结果；而对于10组公认的“困难”结构，我们的算法不仅在一些组中得到与其他算法相当接近的结果，如1FXIα-1UBQ，并且在更多的组中得到明显优于其它算法的结果，如1CRL-1EDE和2SIM-1NSBα。
对图说明
图1.两个蛋白质原子间的相配矩阵S的示意图1aX、Y两个蛋白质的相对位置图；1bS矩阵图。
图2.本发明的程序流程框图。
图3.本发明的一个实施图3a蛋白质组1DHFα和8DFR在原坐标系中；
3b本发明所得最优结构排定解。
具体实施例方式本方法在普通PC机(硬件)平台和Windows软件系统中即可运行。
结合图3，下面介绍本方法的各个步骤1.初始化(a)设定各种初始值取T(0)＝20，ε＝0.8，A是3阶零向量，R为3×3单位矩阵，迭代次数l初值取0；再令λ＝9，vij(0)＝1/max(nx+1，ny+1)，其中0≤i≤nx，0≤j≤ny，vij的收敛阈值0.001；(b)规范化初始数据首先从数据库中导入蛋白质链的初始数据{Xi0}i＝1nx}j＝1ny，然后把待排的蛋白质链平移到它们的共同几何中心ccXi1=Xi0-cc,Yj1=Yj0=Yj0-cc,]]>其中cc=(Σi=1nxXi0+Σj=1nyYj0)/(nx+ny)]]>并使不同蛋白质链的原子间的最大距离为一个单位长Xi=Xi1/fac,Yj=Yj1/fac,]]>其中fac=max1≤i≤nx,1≤j≤ny|Xi1-Yj1|;]]>2.退火过程我们根据退火算法的基本思想，在给定的初温下，不断降低温度参数进行迭代，从而在解空间中寻找自由能量函数F的全局最优解。
● 第一步计算uij(l)，uij(l)=-1T(l-1)∂E(V(l-1))∂vij---(11)]]>其中l表示第l次迭代，我们整理公式(11)得以下表达式uij(l)＝-|A+RXi-Yj|2/T(l-1)，1≤i≤nx，1≤j≤nyu0j(l)=(-λ+λv02(l-1)/2)/T(l-1)j=1(-λ+λ(v0(j-1)+v0(j+1)(l-1))/2)/T(l-1)2≤j<ny(-λ+λv0(ny-1)/2)/T(l-1)j=ny]]>ui0(l)=(-λ+λv20(l-1)/2)/T(l-1)i=1(-λ+λ(v(i-1)0(l-1)+v(i+1)0(l-1)/2)/T(l-1)i≤i<nx(-λ+λv(nx-1)0(l-1)/2)/T(l-1)i=nx]]>● 第二步用迭代法求解wj(l)
wj(l)=Gj(w1(l),...,wny(l))---(12)]]>其中Gj(w1(l),...,wny(l))=e-1+uDj(l)+Σi=1nx(ewij(l)/Σj^=1ny(euij^(l)/wj^(l)))]]>● 第三步根据前两步结果，将uij(l)和wj(l)代入下式求解vij(l)，vij(l)=(euij(l)/wj(l))/Σj^=0ny(euij(l)/wj^(l))----(13)]]>v0j(l)=e-1+u0j(l)/wj(l)---(14)]]>● 若变量V收敛，即是(1/nxny)Σi=1nxΣj=1ny|vij(l)-vij(l-1)|≤0.001]]>则往下继续至第四步，否则的话，转至第一步● 第四步(a)用SVD法计算平移向量A和旋转矩阵R；R=V0U0T,]]>A=Σj=1ny(Σi=1nxvij/vtot)Yj-RΣi=1nx(Σj=1nyvij/vtot)Xi,]]>其中vtot=Σi=1nxΣj=inyvij,]]>且矩阵V0、U0满足下式 即以上等式的右边是左边的SVD分解；(b)计算新的坐标X^i=A+RXi;]]>(c)降低温度T(l)＝T(l-1)×ε；● 若相配矩阵的元素趋近于整数，即是Σi=1nxΣj=1nyvij2/nx∈
]]>则结束退火过程，否则，转至第一步3.输出根据退火过程得到的变量，接下来整理数据，准备输出● 对vij四舍五入，并让sij＝vij，这样即可得到一双蛋白质的最佳结构比对解。
● 计算一双蛋白质链的相配数m及它们间的平均距离(rms)m=Σi=1nxΣj=1nysij,]]>rms=fac1mΣi=1nxΣj=1nysij|Xi-Yj|2]]>
权利要求
1.基于平均场退火技术的蛋白质的立体结构比对方法，其特征在于，它依次含有以下步骤(1)对计算机进行初始化●设定以下各种初始值退火温度T，退火温度递降系数ε，惩罚函数λ，取正的实数值，0＜λ≤100，平移向量A为3阶零向量，旋转矩阵R为3×3单位矩阵，实数值相配变量vij(0)＝1/max(nx+1，ny+1)，0≤vij≤1，0≤i≤nx，0≤j≤ny，其中nx、ny分别为只计主链时，两列待比对的蛋白质X、Y长度，nx、ny均为正整数，i、j分别为X中的第i个原子和Y中的第j个原子，nx、ny也分别是X、Y中蛋白质原子的数目；vij的收敛阈值0.001；l是迭代次数，初始值为0；●规范化数据从通用的Protein Data Bank数据库中导入待排蛋白质链{Xi0}j＝1nx、{Yj0}j＝1ny，则得到蛋白质X中的第i个原子与蛋白质Y中的第j个原子在三维空间中的相应坐标为Xi＝(xi，1，xi，2，xi，3)，Yj＝(yj，1，yj，2，yj，3)，把上述两个蛋白质X、Y平移到几何中心后Xi1=Xi0-cc,Yj1=Yj0-cc,]]>其中cc=(Σi=1nxXi0+Σj=1nyYj0)/(nx+ny),]]>再把Xi1、Yj1单位化，得到Xi=Xi1/fac,Yj=Yj1/fac,]]>其中fac=max1≤i≤nx,1≤j≤ny|Xi1-Yj1|;]]>(2)把要求解的相配变量vij用辅助变量uij、wj表达成代数形式，计算机用下式求解、更新uijuij(l)＝-|A+RXi-Yj|2/T(l-1)，1≤i≤nx，1≤j≤ny，u0j(l)=(-λ+λv02(l-1)/2)/T(l-1)j=1(-λ+λ(v0(j-1)(l-1)+v0(j+1)(l-1))/2)/T(l-1)2≤j<ny(-λ+λv0(ny-1)(l-1)/2)/T(l-1)j=ny,]]>ui0(l)=(-λ+λv20(l-1)/2)/T(l-1)i=1(-λ+λ(v(j-1)0(l-1)+v(j+1)0(l-1))/2)/T(l-1)2≤i<nx(-λ+λv(nx-1)0(l-1)/2)/T(l-1)i=nx,]]>计算机通过迭代法对下式求解另一变量wjwj(l)=e-1+u0j(l)+Σi=1nx(euij(l)/Σj^=1ny(euij^(l)/wj^(l)))]]>其中，l为迭代到第l次的序号， 表示从1到ny的整数；(3)计算机用下式求解相配变量vij(l)=(euij(l)/wj(l))/Σj^=0ny(euij^(l)/wj^(l))]]>v0j(l)=e-1+u0j(l)/wj(l)]]>(4)根据下述收敛判断准则判别vij是否收敛(1/nxny)Σi=1nxΣj=1ny|vij(l)-vij(l-1)|≤0.001]]>(5)若不收敛，回到步骤(2)，进行下一次迭代，若收敛，则执行下一步骤；(6)计算机根据奇异值分解法用下面式子计算旋转矩阵R和平移向量AR=V0U0T]]>A=Σj=1ny(Σi=1nxvij/vtot)Yj-RΣi=1nx(Σj=1nyvij/vtot)Xi]]>其中vtot=Σi=1nxΣj=1nyvij,]]>且矩阵V0、U0满足下式 即以上等式的右边是左边的SVD分解，p表示从1到nx的整数，q表示从1到ny的整数；(7)计算机根据下式计算新的坐标X^i=A+RXi]]>(8)降低退火温度T(l)＝T(l-1)×ε(9)判断相配变量vij是否趋近于整数即Σi=1nxΣj=1nyvij2/nx∈
]]>(10)若相配变量不趋近于整数，则返回步骤(2)，进行再次迭代；若相配变量趋近于整数，则进行下一步；(11)对vij四舍五入，并令sij＝vij，按下式计算蛋白质的相配数m和平均距离rmsm=Σi=1nxΣj=1nysij,rms=fac1mΣi=1nxΣj=1nysij|Xi-Yj|2.]]>
全文摘要
基于平均场退火技术的蛋白质的立体结构比对方法属于生物信息技术领域，其特征在于它把蛋白质的立体比对精确地归结为一非线性整数规划问题。再基于平均场理论，把问题等价为非线性连续型的最优化问题，然后利用工程上的退火技术，高效率地找出相应问题的最优结构比对解。它不需要附加的或人为的“软限制”条件，而得到了精确的最优解，不需要事后处理。它还可用于蛋白质的结构分类，功能预测，主次结构的识别，多基因树的重构以及制药工程等领域。
文档编号G06F19/00GK1670754SQ20041006884
公开日2005年9月21日 申请日期2004年7月9日 优先权日2004年7月9日
发明者陈洛南, 周天寿, 唐云 申请人:清华大学, 陈洛南