论文部分内容阅读
摘要：序列比对是生物信息处理中非常重要的一类方法，基本的序列比对算法是基于动态规划思想提出的。本文提出了一种基于动态规划思想的全局双序列比对优化算法(Optimized Global Pairwise Sequence Alignment based on the idea of Dynamic Programming) OGPSADP，在保持基本动态规划敏感性的前提下，GOPSA方法计算替换矩阵时只需存储当前相邻两列的元素，同时引用checkpoint技术以减少计算迭代次数，有效降低了时间复杂度和空间复杂度。
关键词：生物信息学；序列全局比对；动态规划；替换矩阵
中图分类号：TP301文献标识码：A文章编号：1009-3044(2007)06-11594-03
1 引言
生物信息学是生物学的一个分支，它采用信息科学、计算机科学、生物数学、比较生物学等学科的观点和方法对生命的现象及其组成分子(核酸、蛋白质等)进行研究，主要研究生命中的本质和规律，包括物质组成、结构功能、生命体的能量和信息的交换传递等。
序列比对是生物学中最基本、最重要的方法，是生物学计算的核心。序列比对又叫序列联配，其意义在于从核酸、氨基酸的层次中分析序列的相似性，推测其结构功能及进化上的联系，是基因识别、分子进化、生命起源研究的基础[1]。最常见的比对是蛋白质序列之间或核酸序列之间的两两比对，通过比较两个序列之间的相似性区域，寻找二者可能的分子进化关系。进一步的比对是将多个蛋白质序列或核酸序列同时进行比较，寻找这些有进化关系的序列之间共同的保守区域、位点和剖面信息，从而探索导致它们产生相同功能的序列模式。
2 序列比对
序列同源(homology)指的是序列来自相同的祖先，意味着这些序列具有相同的进化历史，而序列的相似性(similarity)指的是两序列在某参数条件下的相像，它可以用相同残基的百分比或是其他的方法来表示。序列之间的相似度是可以量化的参数，而序列是否同源需要有进化事实的验证，显著的相似性通常意味着同源。
序列比对是运用某种特定的数学模型或算法，找出两个或多个序列之间的最大匹配碱基或残基数，比对算法的结果在很大程度上反映了序列之间的相似性程度以及它们的生物学特征[2]。序列比对根据同时进行比对的序列数目多少可分为双序列比对(pair-wise sequence alignment)和多序列比对(multiple sequence alinment)。序列比对从比对范围考虑也可分为全局比对(global alignment)和局部比对(local alignment)，全局比对考虑序列的全局相似性，局部比对考虑序列片断之间的相似性。如下所示。
全局比对：
LGPSSKQTGKGS-SRIWDN
LN-ITKSAGKGAIMRLGDA
局部比对：
在实际应用中，用全局比对方法企图找出只有局部相似性的两个序列之间的关系显然是徒劳的；而用局部比对得到的局部相似性结果则同样不能说明这两个序列的三维结构或折叠方式是否相同。从这个意义上讲，局部相似性搜索比整体相似性比对更加灵敏，也更具有生物学意义[3]。
3 动态规划思想
动态规划[4]解决序列比对问题的基本思想：使用迭代法计算出两个序列的相似分值，并存入一个得分矩阵中，根据得分矩阵回溯寻找最优的比对序列。
全局比对中的动态规划：
假设存在两条序列As=(a1a2…an)和Bs=(b1b2…bn)，如图1所示，将As和Bs分别作为横坐标和纵坐标放置，组成一个路径矩阵，即得分矩阵，矩阵元素(i,j)值为比对的得分值。在得分矩阵中到达位置为(i,j)的某一个元素有三种可能的路径：通过位置i-1,j-1的对角方向，没有空位罚分；通过列j的垂直方向，通过行i的水平方向，空位罚分的值取决于插入空格的个数。
得分矩阵的元素值通过公式(1)进行迭代计算，得分矩阵的元素值的迭代方式如图2所示。
图2 迭代逻辑图
其中Si,j=s(a1a2…an, b1b2…bn,)是到达序列As中第i位字符与序列Bs中第j位字符的比对得分值；S(ai,bi)是根据替换矩阵得到的ai和bi匹配得分；wx和wy是As和Bs序列的i和j位置前长度为x和y的空位罚分，Sij是动态决定序列As和Bs的子序列的最优得分，矩阵的每个位置都具有得分时，序列整体的最优比对得分将是矩阵最后一行一列的位置的得分。要决定两条序列的最优比对，需要用到回溯技术，即沿着产生最高序列比对得分的路径向相反方向移动，记录相应位置的序列残基字母，这样就可以得到最优的序列的比对方式。
计算得分矩阵的算法描述如下：
步骤一：
步骤二：
For (i=|As|; j=|Bs|; i>0 && j>0) do
If S[i,j]=S[i-1,j-1]+σ(a[i],b[j]) then i--, j--
else if S[i,j]=S[i-1,j]+σ(a[i],_) then i--, insert(“_”,Bs,j)
else if S[i,j]=S[i,j-1]+σ(-,T[j]) then j--, insert(“_”,As,i)
以图1为例，序列As=“acgctg”，Bs=“catgt”，记分函数：σ(x,x)=+2,σ(x,y)=σ(x,-)=σ(-,y)=-1，x和y表示字符变量，|As|、|Bs|表示序列As和Bs的长度。首先要计算矩阵S，元素位置为(0,0)的S[0,0]赋值0，然后根据上述动态规划算法步骤1连续和循环计算矩阵的各个元素。然后，从矩阵的右下角元素开始根据动态规划算法步骤2回溯寻找最优路径。经过动态规划算法得到如图1所示得分矩阵，图中箭头表示回溯路径。
比对结果1：
A：acgctg-
B：-c-atgt
比对结果2：
A：acgctg-
B：-ca-tgt
比对结果3：
A：-acgctg
B：catg-t-
从上述的算法描述和举例中可以看出，若m=|As|，n=|Bs|，基本动态规划算法的时间和空间复杂度为O(m×n)。随着目前研究的序列长度的增加，利用基本动态规划方法解决双序列比对问题无法满足实际的需要。
4 OGPSADP算法
OGPSADP算法的基本思想是：在U矩阵的计算过程中只存储前一列和当前列，从而节省存储空间，增加checkpoint点的个数，可以减少迭代次数；且在U矩阵计算过程中同时记录元素的来源关系，最佳比对路径的获得不需要回溯，减少运算时间。
(1)基于得分矩阵S计算替换矩阵U，矩阵元素U[ab,d]表示得分矩阵S中沿着序列As在对角线ab方向上值为d最大距离值。U矩阵的行与S矩阵的对角线相对应的对角线相对应，列与S矩阵中“等高”的固定计分值对应。ab=a-b对应M对角线，d对应S矩阵中的得分值，U[ab,d]的值为a，该值由S矩阵中ab、ab+1和ab-1上值为d-1的元素所决定。计算过程中，只存储相邻两列的元素，同时记录矩阵元素的数据来源关系，以便于不用回溯就可以获得最佳比对结果。
U[ab,d]=max a s.t.D[a,b]=d where ab=a-b
=-infinity if no such a exists
U[0,0]=max a s.t. As[1…a]=Bs[1…b]
U[ab,d]=-infinity, if |ab|>d
{outer loop,iterated until U[|As|-|Bs|,d]=|As|}
U[ab,d]=max(U[ab+1,d-insertCost],
U[ab,d-mismatchCost]+1,
U[ab-1,d-deleteCost]+1)
{inner Loop,extends diagonal on a run of matches}
While (As[U[ab,d]+1]=Bs[U[ab,d]-ab+1])
U[ab,d]+=1
(2)将替换矩阵U均匀地分割成k2个小块，对每一块从上到下，和从下到上两个方向向中间行计算，根据计算的结果在中间行得到回溯路径经过该行的点，这个点成为一个checkpoint点，checkpoint点将矩阵分为左上角和右下角的两个矩阵部分，然后接着对这两个字矩阵进行同样的处理，直到最后分割的子矩阵只剩下一行或列，这样得到的所有checkpoint点构成最优化路径。
对于基本动态规划算法Needleman-Wunsch、OGPSADP算法进行了对比试验。从美国国家生物信息技术中心的核算序列数据库Genbank中选取两对长度大约为1000、5000的DNA序列，所用参数：单一的记分矩阵，线性空位罚分，试验结果如下表1。
表1 算法空间需求和运行时间比较
从两组试验结果可看出，OGPSADP算法的运行时间和存储空间需求小于基本动态规划算法。
5 结束语
目前在生物信息学领域中用于两两序列全局比对的算法研究并不多，许多可用的比对软件是基于局部比对算法或者直接对局部比对结果进一步加工得到全局比对结果，但是得到的结果并不理想。OGPSADP算法能够达到序列大片段连续对齐，并且对碱基突变、插入和缺失的频率非常小的序列的全局比对具有较满意的计算效率和精度，能有效降低时间和空间复杂度。经过实验发现，将算法用于序列同源性研究和进化分析时，效果并不理想。作者下一步的工作是将此双序列比对算法推广到多序列比对中。
参考文献：
[1] Martin Tompa. Biological Sequence Analysis. Technical Report of Washington University. Department of Computer Sciences, 2000.
[2] T K Attwood, D J Parry-Smith. 罗静初，等译. 生物信息学概论[M]. 北京:北京大学出版社,2001:141-145.
[3] 张敏. 生物序列比对算法研究现状与展望[J]. 大连大学学报,2004,25(4):75-78.
[4] W.J.Kent, A.M.Zahler. Conservation, regulation, synteny, and introns in a large-scale C.briggsae-C.elegans genome alignment based on segment-to-segment, comparison. Proc. Natl. Acad. Sci. USA, 1996,93(22):12098-12103.
本文中所涉及到的图表、注解、公式等内容请以PDF格式阅读原文。
关键词：生物信息学；序列全局比对；动态规划；替换矩阵
中图分类号：TP301文献标识码：A文章编号：1009-3044(2007)06-11594-03
1 引言
生物信息学是生物学的一个分支，它采用信息科学、计算机科学、生物数学、比较生物学等学科的观点和方法对生命的现象及其组成分子(核酸、蛋白质等)进行研究，主要研究生命中的本质和规律，包括物质组成、结构功能、生命体的能量和信息的交换传递等。
序列比对是生物学中最基本、最重要的方法，是生物学计算的核心。序列比对又叫序列联配，其意义在于从核酸、氨基酸的层次中分析序列的相似性，推测其结构功能及进化上的联系，是基因识别、分子进化、生命起源研究的基础[1]。最常见的比对是蛋白质序列之间或核酸序列之间的两两比对，通过比较两个序列之间的相似性区域，寻找二者可能的分子进化关系。进一步的比对是将多个蛋白质序列或核酸序列同时进行比较，寻找这些有进化关系的序列之间共同的保守区域、位点和剖面信息，从而探索导致它们产生相同功能的序列模式。
2 序列比对
序列同源(homology)指的是序列来自相同的祖先，意味着这些序列具有相同的进化历史，而序列的相似性(similarity)指的是两序列在某参数条件下的相像，它可以用相同残基的百分比或是其他的方法来表示。序列之间的相似度是可以量化的参数，而序列是否同源需要有进化事实的验证，显著的相似性通常意味着同源。
序列比对是运用某种特定的数学模型或算法，找出两个或多个序列之间的最大匹配碱基或残基数，比对算法的结果在很大程度上反映了序列之间的相似性程度以及它们的生物学特征[2]。序列比对根据同时进行比对的序列数目多少可分为双序列比对(pair-wise sequence alignment)和多序列比对(multiple sequence alinment)。序列比对从比对范围考虑也可分为全局比对(global alignment)和局部比对(local alignment)，全局比对考虑序列的全局相似性，局部比对考虑序列片断之间的相似性。如下所示。
全局比对：
LGPSSKQTGKGS-SRIWDN
LN-ITKSAGKGAIMRLGDA
局部比对：
在实际应用中，用全局比对方法企图找出只有局部相似性的两个序列之间的关系显然是徒劳的；而用局部比对得到的局部相似性结果则同样不能说明这两个序列的三维结构或折叠方式是否相同。从这个意义上讲，局部相似性搜索比整体相似性比对更加灵敏，也更具有生物学意义[3]。
3 动态规划思想
动态规划[4]解决序列比对问题的基本思想：使用迭代法计算出两个序列的相似分值，并存入一个得分矩阵中，根据得分矩阵回溯寻找最优的比对序列。
全局比对中的动态规划：
假设存在两条序列As=(a1a2…an)和Bs=(b1b2…bn)，如图1所示，将As和Bs分别作为横坐标和纵坐标放置，组成一个路径矩阵，即得分矩阵，矩阵元素(i,j)值为比对的得分值。在得分矩阵中到达位置为(i,j)的某一个元素有三种可能的路径：通过位置i-1,j-1的对角方向，没有空位罚分；通过列j的垂直方向，通过行i的水平方向，空位罚分的值取决于插入空格的个数。
得分矩阵的元素值通过公式(1)进行迭代计算，得分矩阵的元素值的迭代方式如图2所示。
图2 迭代逻辑图
其中Si,j=s(a1a2…an, b1b2…bn,)是到达序列As中第i位字符与序列Bs中第j位字符的比对得分值；S(ai,bi)是根据替换矩阵得到的ai和bi匹配得分；wx和wy是As和Bs序列的i和j位置前长度为x和y的空位罚分，Sij是动态决定序列As和Bs的子序列的最优得分，矩阵的每个位置都具有得分时，序列整体的最优比对得分将是矩阵最后一行一列的位置的得分。要决定两条序列的最优比对，需要用到回溯技术，即沿着产生最高序列比对得分的路径向相反方向移动，记录相应位置的序列残基字母，这样就可以得到最优的序列的比对方式。
计算得分矩阵的算法描述如下：
步骤一：
步骤二：
For (i=|As|; j=|Bs|; i>0 && j>0) do
If S[i,j]=S[i-1,j-1]+σ(a[i],b[j]) then i--, j--
else if S[i,j]=S[i-1,j]+σ(a[i],_) then i--, insert(“_”,Bs,j)
else if S[i,j]=S[i,j-1]+σ(-,T[j]) then j--, insert(“_”,As,i)
以图1为例，序列As=“acgctg”，Bs=“catgt”，记分函数：σ(x,x)=+2,σ(x,y)=σ(x,-)=σ(-,y)=-1，x和y表示字符变量，|As|、|Bs|表示序列As和Bs的长度。首先要计算矩阵S，元素位置为(0,0)的S[0,0]赋值0，然后根据上述动态规划算法步骤1连续和循环计算矩阵的各个元素。然后，从矩阵的右下角元素开始根据动态规划算法步骤2回溯寻找最优路径。经过动态规划算法得到如图1所示得分矩阵，图中箭头表示回溯路径。
比对结果1：
A：acgctg-
B：-c-atgt
比对结果2：
A：acgctg-
B：-ca-tgt
比对结果3：
A：-acgctg
B：catg-t-
从上述的算法描述和举例中可以看出，若m=|As|，n=|Bs|，基本动态规划算法的时间和空间复杂度为O(m×n)。随着目前研究的序列长度的增加，利用基本动态规划方法解决双序列比对问题无法满足实际的需要。
4 OGPSADP算法
OGPSADP算法的基本思想是：在U矩阵的计算过程中只存储前一列和当前列，从而节省存储空间，增加checkpoint点的个数，可以减少迭代次数；且在U矩阵计算过程中同时记录元素的来源关系，最佳比对路径的获得不需要回溯，减少运算时间。
(1)基于得分矩阵S计算替换矩阵U，矩阵元素U[ab,d]表示得分矩阵S中沿着序列As在对角线ab方向上值为d最大距离值。U矩阵的行与S矩阵的对角线相对应的对角线相对应，列与S矩阵中“等高”的固定计分值对应。ab=a-b对应M对角线，d对应S矩阵中的得分值，U[ab,d]的值为a，该值由S矩阵中ab、ab+1和ab-1上值为d-1的元素所决定。计算过程中，只存储相邻两列的元素，同时记录矩阵元素的数据来源关系，以便于不用回溯就可以获得最佳比对结果。
U[ab,d]=max a s.t.D[a,b]=d where ab=a-b
=-infinity if no such a exists
U[0,0]=max a s.t. As[1…a]=Bs[1…b]
U[ab,d]=-infinity, if |ab|>d
{outer loop,iterated until U[|As|-|Bs|,d]=|As|}
U[ab,d]=max(U[ab+1,d-insertCost],
U[ab,d-mismatchCost]+1,
U[ab-1,d-deleteCost]+1)
{inner Loop,extends diagonal on a run of matches}
While (As[U[ab,d]+1]=Bs[U[ab,d]-ab+1])
U[ab,d]+=1
(2)将替换矩阵U均匀地分割成k2个小块，对每一块从上到下，和从下到上两个方向向中间行计算，根据计算的结果在中间行得到回溯路径经过该行的点，这个点成为一个checkpoint点，checkpoint点将矩阵分为左上角和右下角的两个矩阵部分，然后接着对这两个字矩阵进行同样的处理，直到最后分割的子矩阵只剩下一行或列，这样得到的所有checkpoint点构成最优化路径。
对于基本动态规划算法Needleman-Wunsch、OGPSADP算法进行了对比试验。从美国国家生物信息技术中心的核算序列数据库Genbank中选取两对长度大约为1000、5000的DNA序列，所用参数：单一的记分矩阵，线性空位罚分，试验结果如下表1。
表1 算法空间需求和运行时间比较
从两组试验结果可看出，OGPSADP算法的运行时间和存储空间需求小于基本动态规划算法。
5 结束语
目前在生物信息学领域中用于两两序列全局比对的算法研究并不多，许多可用的比对软件是基于局部比对算法或者直接对局部比对结果进一步加工得到全局比对结果，但是得到的结果并不理想。OGPSADP算法能够达到序列大片段连续对齐，并且对碱基突变、插入和缺失的频率非常小的序列的全局比对具有较满意的计算效率和精度，能有效降低时间和空间复杂度。经过实验发现，将算法用于序列同源性研究和进化分析时，效果并不理想。作者下一步的工作是将此双序列比对算法推广到多序列比对中。
参考文献：
[1] Martin Tompa. Biological Sequence Analysis. Technical Report of Washington University. Department of Computer Sciences, 2000.
[2] T K Attwood, D J Parry-Smith. 罗静初，等译. 生物信息学概论[M]. 北京:北京大学出版社,2001:141-145.
[3] 张敏. 生物序列比对算法研究现状与展望[J]. 大连大学学报,2004,25(4):75-78.
[4] W.J.Kent, A.M.Zahler. Conservation, regulation, synteny, and introns in a large-scale C.briggsae-C.elegans genome alignment based on segment-to-segment, comparison. Proc. Natl. Acad. Sci. USA, 1996,93(22):12098-12103.
本文中所涉及到的图表、注解、公式等内容请以PDF格式阅读原文。