基本概念
双序列比对一般来说，是对两个DNA或蛋白质序列进行比较，从而找出两者之间最大的相似性匹配。主要是为了确定两个序列之间的相似性源自于同源性，按照一定的规律进行排序。
比对过程中，错配与突变相对应，而空位对应于插入或删除。该研究还可以拓展到现在热门的语言文本的研究中。
在生物信息处理中，我们希望找出两条序列S和T之间具有的某种相似性关系，这种寻找生物序列相似性关系的算法就是双序列比对算法。
我们通常利用两个序列之间的字符差异来测定序列之间的相似性，两条序列中相应位置的字符如果差异大，那么序列的相似性低，反之，序列的相似性就高——百度百科
全局比对是这将参与比对的两条序列里面的所有字符进行比对。但是目前来说，全局比对的运用范围并不广泛。
取代矩阵：当我们进行双序列比对时，需要对不同碱基的替换进行不同程度的分值设置。但是与此同时又不能想当然的进行设置。对于不同的研究目标或则对象来说，选择合适的替代矩阵至关重要，目前国际上比较常用的矩阵有两种PAM——point accepted mutation（点接受突变）和BLOSUM——blocks substitution matrix(块置换矩阵)。
替代矩阵中的数值，则为碱基或则氨基酸比对时的得分值，可以为正也可以为负。即match或mismatch的分数值。
线性罚分，即penalty＝g*d，其中g为空位长度，d为一个空位的罚分。
仿射型罚分，即penalty=d+(g-1)*e， 其中g为连续空位的长度，d为连续空位中第一个空位的罚分，e为连续空位中第2个及以后空位的罚分。
程序设计
**I.打分矩阵的设计**
在这里我们选取了BLOSUM62矩阵作为对角线的打分规则。BLOSUM矩阵在氨基酸序列比对常常作为常用打分矩阵。打分矩阵可以选择不同的打分矩阵，这里我们只用BLOSUM62矩阵进行研究。
**II.转换状态的设计**
在这里我们设计为X状态、Y状态、M状态。
X状态可以理解为Xi匹配到空位，Y状态可以理解为Yi匹配到空位，M状态可以理解为对角线，然后采用BLOSUM62中的氨基酸比对打分进行统计。
其中值得注意的是，X,Y状态与M状态有着很大的不同，且X,Y状态不能进行转换。
BLOSUM62打分矩阵为match和mismatch的分数值，即M态的对角线转换。X为水平转换、Y为竖直转换。 gap extending 为E代表X、Y自身状态的转换。gap opening为M态到非M态的转换。这里仿射空位罚分是把gap opening和gap extending都加起来了。
公式一中，我们使用X状态来代表。有以下两种情况
- 当X状态自身转换时，我们gap extending罚分即可。
- 当为M态转换而来时，根据上面的表达式。我们需要加入gap opening的空位罚分。即为第二个公式。
值得注意的时，我们在矩阵中可以用向下来直接表达为X状态。
公式二中，实则Y状态与X状态大同小异。在这里不再重复。
公式三中，也是这个设计核心所在。如果存在比对（不一定是相等）,我们可以根据给定的BLOSUM62进行打分。在矩阵中即表示为对角线的移动。这里M态的转换就分为三种情况
- M态自身转换，即左上角的状态加上BLOSUM比对的分值。对角线的转换为第一情况，也是主要情况。
当M状态对角线不存在时，那就只有两种，即为I,J中其中一个为0的情况。不存在对角线。我们这直接令两者相等。 - X状态即X[i][j]=M[i][j]
- Y状态即Y[i][j]=M[i][j]
代码实现
global S
global E
global MIN
global amino#氨基酸数组
global blosum #BLOSUM62
S = -11 #gap opening penalty 在这里表达M状态转换到非M状态的
E = -1 #gap extending penalty 在这里是表达X或则Y转换到其本身的罚分
MIN = -float("inf")#用来进行边值初始化
amino = ['A', 'R', 'N', 'D', 'C', 'Q', 'E', 'G', 'H', 'I', 'L', 'K', 'M', 'F', 'P', 'S', 'T', 'W', 'Y', 'V', 'B', 'Z', 'X', '*']
blosum = [
[ 4, -1, -2, -2, 0, -1, -1, 0, -2, -1, -1, -1, -1, -2, -1, 1, 0, -3, -2, 0, -2, -1, 0, -4],
[-1, 5, 0, -2, -3, 1, 0, -2, 0, -3, -2, 2, -1, -3, -2, -1, -1, -3, -2, -3, -1, 0, -1, -4],
[-2, 0, 6, 1, -3, 0, 0, 0, 1, -3, -3, 0, -2, -3, -2, 1, 0, -4, -2, -3, 3, 0, -1, -4],
[-2, -2, 1, 6, -3, 0, 2, -1, -1, -3, -4, -1, -3, -3, -1, 0, -1, -4, -3, -3, 4, 1, -1, -4],
[ 0, -3, -3, -3, 9, -3, -4, -3, -3, -1, -1, -3, -1, -2, -3, -1, -1, -2, -2, -1, -3, -3, -2, -4],
[-1, 1, 0, 0, -3, 5, 2, -2, 0, -3, -2, 1, 0, -3, -1, 0, -1, -2, -1, -2, 0, 3, -1, -4],
[-1, 0, 0, 2, -4, 2, 5, -2, 0, -3, -3, 1, -2, -3, -1, 0, -1, -3, -2, -2, 1, 4, -1, -4],
[ 0, -2, 0, -1, -3, -2, -2, 6, -2, -4, -4, -2, -3, -3, -2, 0, -2, -2, -3, -3, -1, -2, -1, -4],
[-2, 0, 1, -1, -3, 0, 0, -2, 8, -3, -3, -1, -2, -1, -2, -1, -2, -2, 2, -3, 0, 0, -1, -4],
[-1, -3, -3, -3, -1, -3, -3, -4, -3, 4, 2, -3, 1, 0, -3, -2, -1, -3, -1, 3, -3, -3, -1, -4],
[-1, -2, -3, -4, -1, -2, -3, -4, -3, 2, 4, -2, 2, 0, -3, -2, -1, -2, -1, 1, -4, -3, -1, -4],
[-1, 2, 0, -1, -3, 1, 1, -2, -1, -3, -2, 5, -1, -3, -1, 0, -1, -3, -2, -2, 0, 1, -1, -4],
[-1, -1, -2, -3, -1, 0, -2, -3, -2, 1, 2, -1, 5, 0, -2, -1, -1, -1, -1, 1, -3, -1, -1, -4],
[-2, -3, -3, -3, -2, -3, -3, -3, -1, 0, 0, -3, 0, 6, -4, -2, -2, 1, 3, -1, -3, -3, -1, -4],
[-1, -2, -2, -1, -3, -1, -1, -2, -2, -3, -3, -1, -2, -4, 7, -1, -1, -4, -3, -2, -2, -1, -2, -4],
[ 1, -1, 1, 0, -1, 0, 0, 0, -1, -2, -2, 0, -1, -2, -1, 4, 1, -3, -2, -2, 0, 0, 0, -4],
[ 0, -1, 0, -1, -1, -1, -1, -2, -2, -1, -1, -1, -1, -2, -1, 1, 5, -2, -2, 0, -1, -1, 0, -4],
[-3, -3, -4, -4, -2, -2, -3, -2, -2, -3, -2, -3, -1, 1, -4, -3, -2, 11, 2, -3, -4, -3, -2, -4],
[-2, -2, -2, -3, -2, -1, -2, -3, 2, -1, -1, -2, -1, 3, -3, -2, -2, 2, 7, -1, -3, -2, -1, -4],
[ 0, -3, -3, -3, -1, -2, -2, -3, -3, 3, 1, -2, 1, -1, -2, -2, 0, -3, -1, 4, -3, -2, -1, -4],
[-2, -1, 3, 4, -3, 0, 1, -1, 0, -3, -4, 0, -3, -3, -2, 0, -1, -4, -3, -3, 4, 1, -1, -4],
[-1, 0, 0, 1, -3, 3, 4, -2, 0, -3, -3, 1, -1, -3, -1, 0, -1, -3, -2, -2, 1, 4, -1, -4],
[ 0, -1, -1, -1, -2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -2, 0, 0, -2, -1, -1, -1, -1, -1, -4],
[-4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, 1],
]
#return match or mismatch score
def _match(s, t, i, j):#该函数用于对于双序列S、T中不同氨基酸的match或则mismatch打分
#读取t字符串中t[i-1]字符对于amino对于的顺序值
index1=amino.index(t[i-1])
#读取s字符串中t[j-1]字符对于amino对于的顺序值
index2 = amino.index(s[j-1])
#index2=____________________
return blosum[index1][index2]
#initializers for matrices
#对于X状态情况下初始化，gap extending=d+e*(n-1) 同等状态下的情况
def _init_x(i, j):
if i > 0 and j == 0:
return MIN
else:
if j > 0:
return S + E * j
else:
return 0
#对于Y状态情况下初始化，gap extending=d+e*(n-1) 同等状态下的情况
def _init_y(i, j):
if j > 0 and i == 0:
return MIN
else:
if i > 0:
return S + E * i
else:
return 0
#对M状态进行初始化
def _init_m(i, j):
if j == 0 and i == 0:
return 0
else:
if j == 0 or i == 0:
return MIN
else:
return 0
#生成距离矩阵
def distance_matrix(s, t):
dim_i = len(t) + 1
dim_j = len(s) + 1
#abuse list comprehensions to create matrices
#通过创建X、Y、M矩阵来实现创建三维空间
X = [[_init_x(i, j) for j in range(0, dim_j)] for i in range(0, dim_i)]
#print(X)
Y = [[_init_y(i, j) for j in range(0, dim_j)] for i in range(0, dim_i)]
M = [[_init_m(i, j) for j in range(0, dim_j)] for i in range(0, dim_i)]
# print(M)
for j in range(1, dim_j):
for i in range(1, dim_i):
#因为同状态之间不能相互转换（X不能转向Y)，所以X有两种情况存在。X代表竖直情况
# 第一种是X转向X，即X[i][j]=X[i][j-1]+e
# 第二种是M态转向，X[x][j-1]=M[i][j-1]+S+E
X[i][j] = max((S + E + M[i][j-1]), (E + X[i][j-1]))
#和上述X状态一致，实则X、Y为相同状态。只不过Y[i][j]中是另一个平面的情况。Y代表水平情况
#两者只是简单i,j的替换
# Y[i][j] = max(___________________, (E + Y[i-1][j]))
Y[i][j] = max((S + E + M[i-1][j]), (E + Y[i - 1][j]))
#M状态则存在三种情况
#一个是自身的转换，即M态到M态，这个时候就是对角线的转换。加上上述中的mismatch或则match积分
#第二种和第三种则是对于X、Y状态的转换。这种转换只能发生在当X、Y为边缘时，即为M本身的情况
#M[i][j] = max(_match(s, t, i, j) + M[i-1][j-1], ___________, ____________)
#M[i][j] = max(_match(s, t, i, j) + M[i - 1][j - 1],X[i-1][j-1]+_match(s, t, i, j),Y[i-1][j-1]+_match(s, t, i, j))
M[i][j] = max(_match(s, t, i, j) + M[i - 1][j - 1], X[i][j],
Y[i][j])
return [X, Y, M]
#回溯，得到双序列匹配的字符串
def backtrace(s, t, X, Y, M):
sequ1 = ''
sequ2 = ''
i = len(t)
j = len(s)
while (i>0 or j>0):
#当i，j不为0，即不在边界上时。且M状态是直接转换M态。说明对角线进行转换。
#蛋白质中的氨基酸序列中相匹配（不一定是相等）,但是为最优。
if (i>0 and j>0 and M[i][j] == M[i-1][j-1] + _match(s, t, i, j)):
sequ1 += s[j-1]
sequ2 += t[i-1]
i -= 1
j -= 1
#当j=0且i不为0时，说明只可以进行水平转移，说明在竖直状态下为空位
elif (i>0 and M[i][j] == Y[i][j]):
#sequ1 += ______
#当i>0,且j=0时，我们只对i进行相减.打分矩阵与Y矩阵值一致，说明为存在空位
sequ1 += '_'
sequ2 += t[i-1]
i -= 1
# 当i=0且j不为0时，说明只可以进行竖直转移，说明在水平状态下为空位
elif (j>0 and M[i][j] == X[i][j]):
#sequ1 += _______
#当j>0,且i=0时，我们只对j进行相减，打分矩阵与X矩阵值一致，说明为存在空位
sequ1 += s[j - 1]
sequ2 += '_'
j -= 1
#通过回溯后，倒序输出
sequ1r = ' '.join([sequ1[j] for j in range(-1, -(len(sequ1)+1), -1)])
sequ2r = ' '.join([sequ2[j] for j in range(-1, -(len(sequ2)+1), -1)])
return [sequ1r, sequ2r]
seq1 = input("Please input long sequence:")
seq2 = input("Please input short sequence:")
#需要将输入字符串转换为3个状态
[X, Y, M] = distance_matrix(seq1, seq2)
#print(X)
#print(M)
[str1, str2] = backtrace(seq1, seq2, X, Y, M)
score=M[len(seq2)][len(seq1)]
print("Alignment Score:"+ str(score))
print(str1)
print(str2)
实验结果
登录BLAST在线服务器[http://blast.ncbi.nlm.nih.gov/Blast.cgi]进行测验I.实验结果：
II.官方测验结果
使用指定网站
实验基本符合情况
