双序列局部比对的算法和全局比对算法的步骤相似,只是赋值规则稍有不同,每个单元格的分值由前面对角线的分值和引入一个空位得到的分值的最大值决定,但是分值不能为负值,如果为负值,则将其值设为0.
构建矩阵时先将第一列与第一行的值设为0,而在双序列全局比对中,此处的值为空位罚分的累加
除第一行和第一列以外每个位置的分值由一下四个值的最大值决定
1.S(i-1,j-1)+S(i,j) --匹配和不匹配
2.S(i,j-1)+gap -- 列引入空位
3.S(i-1,j)+gap -- 行引入空位
4.0
打分过程的目的是为了寻找比对的最大值,他代表了比对的结尾处,最大值不一定在最后一行和最后一列,这一点与全局比对不同
回溯过程从最大值开始,沿对角线向左直到碰到一位置(i,j)它的对角线上方,左方,上方的分值都小于匹配得分(在我的代码中是5时),停止寻找。
局部比对的比对区域较短,但相似性百分比较大
Python代码如下:
mport sys
import numpy
import string
## Score matrix for nucleotide alignment
NUC44=numpy.array([[5,-4,-4,-4,-2],\
[-4,5,-4,-4,-2],\
[-4,-4,5,-4,-2],\
[-4,-4,-4,5,-2],\
[-2,-2,-2,-2,-1]])
NBET='ATGCN'
## define the function for calculating score matrix and arrow matrix:
def scoreMat(NUC44,NBET,seq1,seq2,gap=-8):
len1,len2=len(seq1),len(seq2)
scorMat=numpy.zeros((len1+1,len2+1),int)
arrowMat=numpy.zeros((len1+1,len2+1),int)
# scorMat[0,:]=numpy.arange(len2+1)*gap
# scorMat[:,0]=numpy.arange(len1+1)*gap
arrowMat[0,:]=numpy.ones(len2+1)
arrowMat[1:,0]=numpy.zeros(len1)
for i in range(1,len1+1):
for j in range(1,len2+1):
s_mat=numpy.zeros(4)
s_mat[0]=scorMat[i-1,j]+gap
s_mat[1]=scorMat[i,j-1]+gap
n1,n2=NBET.index(seq1[i-1]),NBET.index(seq2[j-1])
s_mat[2]=scorMat[i-1,j-1]+NUC44[n1,n2]
scorMat[i,j]=numpy.max(s_mat)
arrowMat[i,j]=numpy.argmax(s_mat)
return scorMat,arrowMat
## obtain the alignment of the sequences
def DynAlign(scorMat,arrow,seq1,seq2