Python设计双序列全局比对的程序——生物信息pythonlaixiangwei的博客-

双序列比对一般来说，是对两个DNA或蛋白质序列进行比较，从而找出两者之间最大的相似性匹配。主要是为了确定两个序列之间的相似性源自于同源性，按照一定的规律进行排序。
比对过程中，错配与突变相对应，而空位对应于插入或删除。该研究还可以拓展到现在热门的语言文本的研究中。

在生物信息处理中，我们希望找出两条序列S和T之间具有的某种相似性关系，这种寻找生物序列相似性关系的算法就是双序列比对算法。

我们通常利用两个序列之间的字符差异来测定序列之间的相似性，两条序列中相应位置的字符如果差异大，那么序列的相似性低，反之，序列的相似性就高——百度百科

全局比对是这将参与比对的两条序列里面的所有字符进行比对。但是目前来说，全局比对的运用范围并不广泛。

取代矩阵：当我们进行双序列比对时，需要对不同碱基的替换进行不同程度的分值设置。但是与此同时又不能想当然的进行设置。对于不同的研究目标或则对象来说，选择合适的替代矩阵至关重要，目前国际上比较常用的矩阵有两种PAM——point accepted mutation（点接受突变）和BLOSUM——blocks substitution matrix(块置换矩阵)。
替代矩阵中的数值，则为碱基或则氨基酸比对时的得分值，可以为正也可以为负。即match或mismatch的分数值。

空位罚分：补偿插入和缺失对序列相似性的影响。当双序列中不存在比对情况时，我们进行空位罚分。分为gap opening 和 gap extending。gap opening为两个或则两个以上的连续空位罚分，gap extending则为出现第一次的空位罚分。

线性罚分，即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状态来代表。有以下两种情况

1. 当X状态自身转换时，我们gap extending罚分即可。
2. 当为M态转换而来时，根据上面的表达式。我们需要加入gap opening的空位罚分。即为第二个公式。
   值得注意的时，我们在矩阵中可以用向下来直接表达为X状态。

公式二中，实则Y状态与X状态大同小异。在这里不再重复。

公式三中，也是这个设计核心所在。如果存在比对（不一定是相等）,我们可以根据给定的BLOSUM62进行打分。在矩阵中即表示为对角线的移动。这里M态的转换就分为三种情况

1. M态自身转换，即左上角的状态加上BLOSUM比对的分值。对角线的转换为第一情况，也是主要情况。
   当M状态对角线不存在时，那就只有两种，即为I,J中其中一个为0的情况。不存在对角线。我们这直接令两者相等。
2. X状态即X[i][j]=M[i][j]
3. 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) 

I.实验结果：

II.官方测验结果
使用指定网站

实验基本符合情况