1.序列比较算法（全局序列比对及局部序列比对的python实现）
- 前言
-
- 算法思想介绍
- 实现功能及实现方法
- 运行结果演示
- 源代码
- 遇到的问题及总结
前言
阶段性地完成了DNA序列比较算法，还有很多不足和需要完善的地方有待日后改进。
算法思想介绍
一个很详细完整的算法介绍
双序列全局比对及算法
Needleman-Wunsch 算法：动态规划法
输入值：两条序列、替换记分矩阵以确定不同字母间的相似度得分，以及空位罚分
双序列局部比对算法
局部比对的计算公式在全局比对的基础上增加了第四个元素“0”。得分矩阵初始值仍是0，但第一行和第一列与全局比对不同，全是0。
实现功能及实现方法
- 使用已定义的记分矩阵
- 设置需比对的序列、gap大小及要进行的比对选择
- 打印最高得分的序列比对结果
方法：倒序查找最终路进行序列比对 - 打印得分矩阵及比对路径
方法：使用递归和栈记录最终路径
运行结果演示
-
双序列全局比对
输入序列：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("无效选择！")
遇到的问题及总结
-
思维误区： 最初在存储最终路径的问题上，认为在每次路径递归至初始结点时才将其入栈，导致遇到分岔路则无法记录前一条路的完整路径
经过高人指点后解决：每次到达一个结点就将其入栈，递归结束后将其出栈 -
思维误区2：最初以采用存储每个点的前一点坐标为存储方式，导致所有得分矩阵只能存储一条路径
再次经过高人指点后解决：改存储方式为存储每个点计算得来的方向 -
递归终止条件存储最终路径 -涉及问题：list的赋值、浅拷贝、深拷贝
Python List的赋值方法 -
画出路径矩阵表格 -涉及问题：表格与画布大小不一致且导致无法确定箭头坐标
解决方式：使用bbox参数 -
获取得分矩阵中最大值的索引 - 涉及问题：获取where结果的索引值
-
局部比对递归终止条件 - 涉及问题：列表比较 any() all()
总结
这次的作业因为拖延很晚才开始，遇到的问题也绝不仅仅只有以上几个，最深刻的还是最开始的思维误区，直接卡到崩溃，但其他问题也都耗费了或多或少的时间才解决，更有无数小问题调试了无数遍才得以做出这个半成品。其实还有很多待完善的地方，比如输入的字符判断，比如一致度和相似度的计算，比如图形化界面的编写，写出来的代码也还有很多不成熟的地方，但是因为时间问题，暂时就到这了，有时间再改进咯。
