devbox
爱喝兽奶帝天荒
这个屌丝很懒，什么也没留下！
关注作者
热门标签
热门文章
当前位置:   article > 正文

Python设计双序列全局比对的程序——生物信息_编写代码找出blosum62矩阵中得分最高的氨基酸配对

作者：爱喝兽奶帝天荒 | 2024-08-20 13:38:53

赞

踩

编写代码找出blosum62矩阵中得分最高的氨基酸配对

基本概念

双序列比对一般来说，是对两个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)

  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41
  • 42
  • 43
  • 44
  • 45
  • 46
  • 47
  • 48
  • 49
  • 50
  • 51
  • 52
  • 53
  • 54
  • 55
  • 56
  • 57
  • 58
  • 59
  • 60
  • 61
  • 62
  • 63
  • 64
  • 65
  • 66
  • 67
  • 68
  • 69
  • 70
  • 71
  • 72
  • 73
  • 74
  • 75
  • 76
  • 77
  • 78
  • 79
  • 80
  • 81
  • 82
  • 83
  • 84
  • 85
  • 86
  • 87
  • 88
  • 89
  • 90
  • 91
  • 92
  • 93
  • 94
  • 95
  • 96
  • 97
  • 98
  • 99
  • 100
  • 101
  • 102
  • 103
  • 104
  • 105
  • 106
  • 107
  • 108
  • 109
  • 110
  • 111
  • 112
  • 113
  • 114
  • 115
  • 116
  • 117
  • 118
  • 119
  • 120
  • 121
  • 122
  • 123
  • 124
  • 125
  • 126
  • 127
  • 128
  • 129
  • 130
  • 131
  • 132
  • 133
  • 134
  • 135
  • 136
  • 137
  • 138
  • 139
  • 140
  • 141
  • 142
  • 143
  • 144
  • 145
  • 146
  • 147
  • 148
  • 149
  • 150
  • 151
  • 152
  • 153
  • 154

实验结果

登录BLAST在线服务器[http://blast.ncbi.nlm.nih.gov/Blast.cgi]进行测验

I.实验结果：

输入氨基酸序列用来判断

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

BLAST官网测试

实验基本符合情况

声明：本文内容由网友自发贡献，不代表【wpsshop博客】立场，版权归原作者所有，本站不承担相应法律责任。如您发现有侵权的内容，请联系我们。转载请注明出处：https://www.wpsshop.cn/w/爱喝兽奶帝天荒/article/detail/1007438
推荐阅读
相关标签

Copyright © 2003-2013 www.wpsshop.cn 版权所有，并保留所有权利。

  

闽ICP备14008679号