devbox
笔触狂放9
这个屌丝很懒，什么也没留下！
关注作者
热门标签
热门文章
当前位置:   article > 正文

1.序列比较算法（全局序列比对及局部序列比对的python实现）_全局比对函数和局部比对函数

作者：笔触狂放9 | 2024-08-20 13:34:02

赞

踩

全局比对函数和局部比对函数

1.序列比较算法（全局序列比对及局部序列比对的python实现）

前言

阶段性地完成了DNA序列比较算法，还有很多不足和需要完善的地方有待日后改进。

算法思想介绍

一个很详细完整的算法介绍

双序列全局比对及算法
Needleman-Wunsch 算法：动态规划法
输入值：两条序列、替换记分矩阵以确定不同字母间的相似度得分，以及空位罚分
全局比对计算公式

双序列局部比对算法
局部比对的计算公式在全局比对的基础上增加了第四个元素“0”。得分矩阵初始值仍是0，但第一行和第一列与全局比对不同，全是0。
局部比对计算公式

实现功能及实现方法

  1. 使用已定义的记分矩阵
    替换记分矩阵
  2. 设置需比对的序列、gap大小及要进行的比对选择
  3. 打印最高得分的序列比对结果
    方法：倒序查找最终路进行序列比对
  4. 打印得分矩阵及比对路径
    方法：使用递归和栈记录最终路径

运行结果演示

  • 双序列全局比对

    输入序列：atcggtac；aatcgta在这里插入图片描述
    输入序列：atcggt；aacg
    全局比对

  • 双序列局部比对

输入序列：atccga；tcga在这里插入图片描述
输入序列：acgtc；cg
局部比对

源代码

#!/usr/bin/env python
# -*- coding:utf-8 -*-

from numpy import *
import copy
from matplotlib import pyplot as plt
from pandas import *

#创建全局比对得分矩阵
def GlobalScoreMatrix(m,n,w,replace,s,path,senquence1,senquence2,gap):
    for i in range(m):
        for j in range(n):
            #判断s（0,0）这一特殊情况
            if i == 0 and j == 0:
                s[i][j] = 0
            elif i-1 < 0:#判断第一行的特殊情况
                s[i][j] = s[i][j - 1] + gap
                path[i,j,0] = 1
            elif j-1 < 0:#判断第一列的特殊情况
                s[i][j] = s[i - 1][j] + gap
                path[i,j,1] = 1
            else:
                #获取最大值
                s[i][j] = max(s[i - 1][j - 1] + w[replace[senquence2[i - 1]]][replace[senquence1[j - 1]]],
                              s[i - 1][j] + gap, s[i][j - 1] + gap)
                #记录最大值来的方向
                if s[i - 1][j - 1] + w[replace[senquence2[i - 1]]][replace[senquence1[j - 1]]] == s[i][j]:
                    path[i,j,2] = 1
                if s[i - 1][j] + gap == s[i][j]:
                    path[i,j,1] = 1
                if s[i][j - 1] + gap == s[i][j]:
                    path[i,j,0] = 1

#创建局部比对得分矩阵
def LocalScoreMatrix(m,n,w,replace,s,path,senquence1,senquence2,gap):
    for i in range(1,m):
        #局部矩阵第一行及第一列均为0，不需要再初始化
        for j in range(1,n):
            #获取最大值,与全局比对不同之处在于选取最大值时将0加入选择
            s[i][j] = max(0,s[i - 1][j - 1] + w[replace[senquence2[i - 1]]][replace[senquence1[j - 1]]],
                          s[i - 1][j] + gap, s[i][j - 1] + gap)
            #记录最大值来的方向，若最大值为0则不必记录
            if s[i - 1][j - 1] + w[replace[senquence2[i - 1]]][replace[senquence1[j - 1]]] == s[i][j]:
                path[i,j,2] = 1
            if s[i - 1][j] + gap == s[i][j]:
                path[i,j,1] = 1
            if s[i][j - 1] + gap == s[i][j]:
                path[i,j,0] = 1

#寻找全局序列比对路径
def FindGlobalPath(i,j,path,OnePath,LastGlobalPath):
    #递归终止条件：回到原点（0，0）
    if i == 0 and j == 0:
        OnePath.append((i, j))
        #将OnePath进行深拷贝再加入至最终路径列表LastGlobalPath中
        LastGlobalPath.append(copy.deepcopy(OnePath))
        #将该点出栈
        OnePath.pop()
    else:
        for k in range(3):
            #判断每个点来的方向
            if path[i,j,k] == 1:
                #下标0处记录从左来的方向
                if k == 0:
                    #将该点入栈
                    OnePath.append((i,j))
                    #进行递归记录下一个点
                    FindGlobalPath(i,j - 1,path,OnePath,LastGlobalPath)
                    #递归返回后将该点出栈，记录另一方向
                    OnePath.pop()
                #下标1处记录从上来的方向
                elif k == 1:
                    OnePath.append((i, j))
                    FindGlobalPath(i - 1, j, path,OnePath,LastGlobalPath)
                    OnePath.pop()
                #下标2处记录从左上来的方向
                else:
                    OnePath.append((i, j))
                    FindGlobalPath(i - 1, j - 1, path,OnePath,LastGlobalPath)
                    OnePath.pop()

# 寻找局部序列比对路径
def FindLocalPath(i, j, path, OnePath, LastLocalPath):
    #递归终止条件：某个没有记录方向的点
    if all(path[i][j] == [0, 0, 0]):  
        OnePath.append((i, j))
        # 将OnePath进行深拷贝再加入至最终路径列表LastLocalPath中
        LastLocalPath.append(copy.deepcopy(OnePath))
        # 将该点出栈
        OnePath.pop()
    else:
        for k in range(3):
            # 判断每个点来的方向
            if path[i, j, k] == 1:
                # 下标0处记录从左来的方向
                if k == 0:
                    # 将该点入栈
                    OnePath.append((i, j))
                    # 进行递归记录下一个点
                    FindLocalPath(i, j - 1, path, OnePath, LastLocalPath)
                    # 递归返回后将该点出栈，记录另一方向
                    OnePath.pop()
                # 下标1处记录从上来的方向
                elif k == 1:
                    OnePath.append((i, j))
                    FindLocalPath(i - 1, j, path, OnePath, LastLocalPath)
                    OnePath.pop()
                # 下标2处记录从左上来的方向
                else:
                    OnePath.append((i, j))
                    FindLocalPath(i - 1, j - 1, path, OnePath, LastLocalPath)
                    OnePath.pop()

#输出比对后的序列
def ShowContrastResult(LastPath,senquence1,senquence2):
    #依次输出每条路径
    for k,aPath in enumerate(LastPath):
        rowS = ''
        colS = ''
        #每条路径倒序遍历
        for i in range(len(aPath) -1,0,-1):
            #方向从左边来
            if aPath[i][0] == aPath[i - 1][0]:
                rowS += senquence1[aPath[i - 1][1] - 1]
                colS += '-'
            #方向从上面来
            elif aPath[i][1] == aPath[i - 1][1]:
                colS += senquence2[aPath[i - 1][0] -1]
                rowS += '-'
            #方向从左上来
            else:
                rowS += senquence1[aPath[i - 1][1] - 1]
                colS += senquence2[aPath[i - 1][0] - 1]
        #依次输出每条路的序列比对结果
        print("======比对结果",k+1,"======")
        print("序列1:",rowS)
        print("序列2:",colS)

# 判断是否为最终路径中的点
def judgePath(point, LastPath):
    for aPath in LastPath:
        if point in aPath:
            return True
    return False

# 画出路径图
def ShowPaths(senquence1, senquence2, LastPath):
    s1 = "0" + senquence1
    s2 = "0" + senquence2
    # 列索引
    col = list(s1)
    # 行索引
    row = list(s2)
    # 设置画布大小
    fig = plt.figure(figsize=(5, 5))
    ax = fig.add_subplot(111, frameon=True, xticks=[], yticks=[], )
    the_table = plt.table(cellText=s, rowLabels=row, colLabels=col, rowLoc='right',
                          loc='center', cellLoc='bottom right', bbox=[0, 0, 1, 1])
    # 设置表格文本字体大小
    the_table.set_fontsize(8)
    # 画出每个点的路径图
    for i in range(m):
        for j in range(n):
            for k in range(3):
                if path[i, j, k] == 1:  # 画出记录的方向
                    # 下标0处记录从左来的方向
                    if k == 0:
                        if judgePath((i, j), LastPath):  # 若某点在在最终路径中
                            # 画出红色箭头
                            plt.annotate('', xy=(j / n, (2 * m - 2 * i - 1) / (2 * (m + 1))),
                                         xytext=((2 * j + 1) / (2 * n), (2 * m - 2 * i - 1) / (2 * (m + 1))),
                                         arrowprops=dict(arrowstyle="->", color='r', connectionstyle="arc3"))
                        else:
                            # 未在最终路径中则画出黑色箭头
                            plt.annotate('', xy=(j / n, (2 * m - 2 * i - 1) / (2 * (m + 1))),
                                         xytext=((2 * j + 1) / (2 * n), (2 * m - 2 * i - 1) / (2 * (m + 1))),
                                         arrowprops=dict(arrowstyle="->", connectionstyle="arc3"))
                    # 下标1处记录从上来的方向
                    elif k == 1:
                        if judgePath((i, j), LastPath):
                            plt.annotate('', xy=((2 * j + 1) / (2 * n), (2 * m - 2 * i - 1) / (2 * (m + 1))),
                                         xytext=((2 * j + 1) / (2 * n), (m - i) / (m + 1)),
                                         arrowprops=dict(arrowstyle="<-", color='r', connectionstyle="arc3"))
                        else:
                            plt.annotate('', xy=((2 * j + 1) / (2 * n), (2 * m - 2 * i - 1) / (2 * (m + 1))),
                                         xytext=((2 * j + 1) / (2 * n), (m - i) / (m + 1)),
                                         arrowprops=dict(arrowstyle="<-", connectionstyle="arc3"))
                    # 下标1处记录从上来的方向
                    elif k == 2:
                        if judgePath((i, j), LastPath):
                            plt.annotate('', xy=((2 * j + 1) / (2 * n), (2 * m - 2 * i - 1) / (2 * (m + 1))),
                                         xytext=(j / n, (m - i) / (m + 1)),
                                         arrowprops=dict(arrowstyle="<-", color='r', connectionstyle="arc3"))
                        else:
                            plt.annotate('', xy=((2 * j + 1) / (2 * n), (2 * m - 2 * i - 1) / (2 * (m + 1))),
                                         xytext=(j / n, (m - i) / (m + 1)),
                                         arrowprops=dict(arrowstyle="<-", connectionstyle="arc3"))
    plt.show()


#将碱基转换为集合下标
replace = {'A':0,'G':1,'C':2,'T':3}
#构造替换计分矩阵
w = [[10,-1,-3,-4],[-1,7,-5,-3],[-3,-5,9,0],[-4,-3,0,8]]
#定义需要比对的序列
senquence1 = input("请输入序列1：").upper()
senquence2 = input("请输入序列2：").upper()
#定义输入的gap
gap = int(input("请输入gap："))
choise = int(input("请选择要进行的序列比对(1-全局序列比对  2-局部序列比对) : "))
# 获取序列的长度
m = len(senquence2) + 1
n = len(senquence1) + 1
#构建m*n全0矩阵
s = zeros((m,n))
#记录每个点的方向，下标0处存储从左来的方向，下标1处存储从上来的方向，下标2处存储从左上来的方向
#初始值均为0，若存在从某方向上来则将其对应下标的值置为1
path = zeros((m,n,3))
#记录每条路径
OnePath = []
#记录所有全局序列比对路径
LastGlobalPath = []
#记录所有局部序列比对路径
LastLocalPath = []

if choise == 1:#进行全局序列比对
    #构建得分矩阵
    GlobalScoreMatrix(m,n,w,replace,s,path,senquence1,senquence2,gap)
    FindGlobalPath(m-1,n-1,path,OnePath,LastGlobalPath)
    ShowContrastResult(LastGlobalPath,senquence1,senquence2)
    ShowPaths(senquence1, senquence2, LastGlobalPath)
elif choise == 2:#进行局部序列比对
    #构建得分矩阵
    LocalScoreMatrix(m, n, w, replace, s, path, senquence1, senquence2, gap)
    row = where(s == np.max(s))[0]
    # 获取得分矩阵中最大值的列索引
    col = where(s == np.max(s))[1]
    for i in range(len(row)):#依次寻找不同局部比对的比对路径并输出比对结果
        FindLocalPath(row[i], col[i], path, OnePath, LastLocalPath)
        ShowContrastResult(LastLocalPath, senquence1, senquence2)
        ShowPaths(senquence1, senquence2, LastLocalPath)
else:
    print("无效选择！")

  • 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
  • 155
  • 156
  • 157
  • 158
  • 159
  • 160
  • 161
  • 162
  • 163
  • 164
  • 165
  • 166
  • 167
  • 168
  • 169
  • 170
  • 171
  • 172
  • 173
  • 174
  • 175
  • 176
  • 177
  • 178
  • 179
  • 180
  • 181
  • 182
  • 183
  • 184
  • 185
  • 186
  • 187
  • 188
  • 189
  • 190
  • 191
  • 192
  • 193
  • 194
  • 195
  • 196
  • 197
  • 198
  • 199
  • 200
  • 201
  • 202
  • 203
  • 204
  • 205
  • 206
  • 207
  • 208
  • 209
  • 210
  • 211
  • 212
  • 213
  • 214
  • 215
  • 216
  • 217
  • 218
  • 219
  • 220
  • 221
  • 222
  • 223
  • 224
  • 225
  • 226
  • 227
  • 228
  • 229
  • 230
  • 231
  • 232
  • 233
  • 234
  • 235
  • 236
  • 237
  • 238
  • 239
  • 240
  • 241
  • 242
  • 243

遇到的问题及总结

  1. 思维误区： 最初在存储最终路径的问题上，认为在每次路径递归至初始结点时才将其入栈，导致遇到分岔路则无法记录前一条路的完整路径
    经过高人指点后解决：每次到达一个结点就将其入栈，递归结束后将其出栈

  2. 思维误区2：最初以采用存储每个点的前一点坐标为存储方式，导致所有得分矩阵只能存储一条路径
    再次经过高人指点后解决：改存储方式为存储每个点计算得来的方向

  3. 递归终止条件存储最终路径 -涉及问题：list的赋值、浅拷贝、深拷贝
    Python List的赋值方法

  4. 画出路径矩阵表格 -涉及问题：表格与画布大小不一致且导致无法确定箭头坐标
    解决方式：使用bbox参数

  5. 获取得分矩阵中最大值的索引 - 涉及问题：获取where结果的索引值

  6. 局部比对递归终止条件 - 涉及问题：列表比较 any() all()

总结
这次的作业因为拖延很晚才开始，遇到的问题也绝不仅仅只有以上几个，最深刻的还是最开始的思维误区，直接卡到崩溃，但其他问题也都耗费了或多或少的时间才解决，更有无数小问题调试了无数遍才得以做出这个半成品。其实还有很多待完善的地方，比如输入的字符判断，比如一致度和相似度的计算，比如图形化界面的编写，写出来的代码也还有很多不成熟的地方，但是因为时间问题，暂时就到这了，有时间再改进咯。

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

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

  

闽ICP备14008679号