利用python处理dna序列_详解基于python的全局与局部序列比对的实现(DNA)
程序能实现什么
a.完成gap值的自定义输入以及两条需比对序列的输入
b.完成得分矩阵的计算及输出
c.输出序列比对结果
d.使用matplotlib对得分矩阵路径的绘制
一、实现步骤
1.用户输入步骤
a.输入自定义的gap值
b.输入需要比对的碱基序列1（A,T,C,G）换行表示输入完成
b.输入需要比对的碱基序列2（A,T,C,G）换行表示输入完成
输入（示例）：
2.代码实现步骤
1.获取到用户输入的gap，s以及t
2.调用构建得分矩阵函数，得到得分矩阵以及方向矩阵
3.将得到的得分矩阵及方向矩阵作为参数传到回溯函数中开始回溯得到路径，路径存储使用的是全局变量，存的仍然是方向而不是坐标位置减少存储开销，根据全局变量中存储的方向将比对结果输出。
4.根据全局变量中存储的方向使用matplotlib画出路径
全局比对代码如下：
import matplotlib.pyplot as plt
import numpy as np
#定义全局变量列表finalList存储最后回溯的路径 finalOrder1,finalOrder2存储最后的序列 finalRoad用于存储方向路径用于最后画图
def createList():
global finalList
global finalOrder1
global finalOrder2
global finalRoad
finalList = []
finalOrder1 = []
finalOrder2 = []
finalRoad = []
#创建A G C T 对应的键值对，方便查找计分矩阵中对应的得分
def createDic():
dic = {'A':0,'G':1,'C':2,'T':3}
return dic
#构建计分矩阵
# A G C T
def createGrade():
grade = np.matrix([[10,-1,-3,-4],
[-1,7,-5,-3],
[-3,-5,9,0],
[-4,-3,0,8]])
return grade
#计算两个字符的相似度得分函数
def getGrade(a,b):
dic = createDic() # 碱基字典 方便查找计分矩阵
grade = createGrade() # 打分矩阵grade
return grade[dic[a],dic[b]]
#构建得分矩阵函数 参数为要比较序列、自定义的gap值
def createMark(s,t,gap):
a = len(s) #获取序列长度a,b
b = len(t)
mark = np.zeros((a+1,b+1)) #初始化全零得分矩阵
direction = np.zeros((a+1,b+1,3)) #direction矩阵用来存储得分矩阵中得分来自的方向 第一个表示左方 第二个表示左上 第三个表示上方 1表示能往哪个方向去
#由于得分可能会来自多个方向，所以使用三维矩阵存储
direction[0][0] = -1 #确定回溯时的结束条件 即能够走到方向矩阵的值为-1
mark[0,:] = np.fromfunction(lambda x, y: gap * (x + y), (1, b + 1), dtype=int) #根据gap值将得分矩阵第一行计算出
mark[:,0] = np.fromfunction(lambda x, y: gap * (x + y), (1, a + 1), dtype=int) #根据gap值将得分矩阵第一列计算出
for i in range(1,b+1):
direction[0,i,0] = 1
for i in range(1, a + 1):
direction[i, 0, 2] = 1
for i in range(1,a+1):
for j in range(1,b+1):
threeMark = [mark[i][j-1],mark[i-1][j-1],mark[i-1][j]] #threeMark表示现在所要计算得分的位置的左边 左上 上边的得分
threeGrade = [gap,getGrade(s[i-1],t[j-1]),gap] #threeGrade表示经过需要计算得左边 左上 上边的空位以及相似度得分
finalGrade = np.add(threeMark,threeGrade) #finalGrade表示最终来自三个方向上的得分
mark[i][j] = max(finalGrade) #选取三个方向上的最大得分存入得分矩阵
#可能该位置的得分可以由多个方向得来，所以进行判断并循环赋值
for k in range(0,len([y for y,x in enumerate(finalGrade) if x == max(finalGrade)])):
directionList = [y for y,x in enumerate(finalGrade) if x == max(finalGrade)]
direction[i][j][directionList[k]] = 1
return mark,direction
#回溯函数 参数分别为 得分矩阵 方向矩阵 现在所处得分矩阵的位置 以及两个序列
def remount(mark,direction,i,j,s,t):
if direction[i][j][0] == 1 :
if direction[i][j-1][0] == -1: #如果该位置指向左边 先判断其左边是否是零点
finalList.append(0) #如果是 将该路径存入路径列表
finalList.reverse() #将列表反过来得到从零点开始的路径
index1 = 0 #记录现在所匹配序列s的位置 因为两个字符串可能是不一样长的
index2 = 0 #记录现在所匹配序列t的位置
for k in finalList:
if k == 0 :
finalOrder1.append("-")
finalOrder2.append(t[index2])
index2 += 1
if k == 1 :
finalOrder1.append(s[index1])
finalOrder2.append(t[index2])
index1 += 1
index2 += 1
if k == 2 :
finalOrder1.append(s[index1])
finalOrder2.append("-")
index1 += 1
finalList.reverse() # 将原来反转的路径再返回来
finalRoad.append(np.array(finalList)) # 将此次的路径添加到最终路径记录用于最后画图
finalList.pop() #输出后将当前方向弹出 并回溯
return
else :
finalList.append(0) #如果不是零点 则将该路径加入路径矩阵，继续往下走
remount(mark,direction,i,j-1,s,t)
finalList.pop() #该方向走完后将这个方向弹出 继续下一轮判断 下面两个大的判断同理
if direction[i][j][1] == 1 :
if direction[i-1][j-1][0] == -1:
finalList.append(1)
finalList.reverse() # 将列表反过来得到从零点开始的路径
index1 = 0 # 记录现在所匹配序列s的位置 因为两个字符串可能是不一样长的
index2 = 0 # 记录现在所匹配序列t的位置
for k in finalList:
if k == 0 :
finalOrder1.append("-")
finalOrder2.append(t[index2])
index2 += 1
if k == 1 :
finalOrder1.append(s[index1])
finalOrder2.append(t[index2])
index1 += 1
index2 += 1
if k == 2 :
finalOrder1.append(s[index1])
finalOrder2.append("-")
index1 += 1
finalList.reverse() # 将原来反转的路径再返回来
finalRoad.append(np.array(finalList)) # 将此次的路径添加到最终路径记录用于最后画图
finalList.pop()
return
else :
finalList.append(1)
remount(mark,direction,i-1,j-1,s,t)
finalList.pop()
if direction[i][j][2] == 1 :
if direction[i-1][j][0] == -1:
finalList.append(2)
finalList.reverse() # 将列表反过来得到从零点开始的路径
index1 = 0 # 记录现在所匹配序列s的位置 因为两个字符串可能是不一样长的
index2 = 0 # 记录现在所匹配序列t的位置
for k in finalList:
if k == 0 :
finalOrder1.append("-")
finalOrder2.append(t[index2])
index2 += 1
if k == 1 :
finalOrder1.append(s[index1])
finalOrder2.append(t[index2])
index1 += 1
index2 += 1
if k == 2 :
finalOrder1.append(s[index1])
finalOrder2.append("-")
index1 += 1
finalList.reverse() # 将原来反转的路径再返回来
finalRoad.append(np.array(finalList)) # 将此次的路径添加到最终路径记录用于最后画图
finalList.pop()
return
else :
finalList.append(2)
remount(mark,direction,i-1,j,s,t)
finalList.pop()
#画箭头函数
def arrow(ax,sX,sY,aX,aY):
ax.arrow(sX,sY,aX,aY,length_includes_head=True, head_width=0.15, head_length=0.25, fc='w', ec='b')
#画图函数
def drawArrow(mark, direction, a, b, s, t):
#a是s的长度为4 b是t的长度为6
fig = plt.figure()
ax = fig.add_subplot(111)
val_ls = range(a+2)
scale_ls = range(b+2)
index_ls = []
index_lsy = []
for i in range(a):
if i == 0:
index_lsy.append('#')
index_lsy.append(s[a-i-1])
index_lsy.append('0')
for i in range(b):
if i == 0:
index_ls.append('#')
index_ls.append('0')
index_ls.append(t[i])
plt.xticks(scale_ls, index_ls) #设置坐标字
plt.yticks(val_ls, index_lsy)
for k in range(1,a+2):
y = [k for i in range(0,b+1)]
x = [x for x in range(1,b+2)]
ax.scatter(x, y, c='y')
for i in range(1,a+2):
for j in range(1,b+2):
ax.text(j,a+2-i,int(mark[i-1][j-1]))
lX = b+1
lY = 1
for n in range(0,len(finalRoad)):
for m in (finalRoad[n]):
if m == 0:
arrow(ax,lX,lY,-1,0)
lX = lX - 1
elif m == 1:
arrow(ax,lX,lY,-1,1)
lX = lX - 1
lY = lY + 1
elif m == 2:
arrow(ax, lX, lY, 0, 1)
lY = lY + 1
lX = b + 1
lY = 1
ax.set_xlim(0, b + 2) # 设置图形的范围，默认为[0,1]
ax.set_ylim(0, a + 2) # 设置图形的范围，默认为[0,1]
ax.set_aspect('equal') # x轴和y轴等比例
plt.show()
plt.tight_layout()
if __name__ == '__main__':
createList()
print("Please enter gap:")
gap = int(input()) #获取gap值 转换为整型 tip：刚开始就是因为这里没有进行类型导致后面的计算部分报错
print("Please enter sequence 1:")
s = input() #获取用户输入的第一条序列
print("Please enter sequence 2:")
t = input() #获取用户输入的第二条序列
a = len(s) #获取s的长度
b = len(t) #获取t的长度
mark,direction = createMark(s,t,gap)
print("The scoring matrix is as follows：") #输出得分矩阵
print(mark)
remount(mark,direction,a,b,s,t) #调用回溯函数
c = a if a > b else b #判断有多少种比对结果得到最终比对序列的长度
total = int(len(finalOrder1)/c)
for i in range(1,total+1): #循环输出比对结果
k = str(i)
print("Sequence alignment results "+k+" is：")
print(finalOrder1[(i-1)*c:i*c])
print(finalOrder2[(i-1)*c:i*c])
drawArrow(mark, direction, a, b, s, t)
局部比对代码如下
import matplotlib.pyplot as plt
import numpy as np
import operator
#在局部比对中 回溯结束的条件是方向矩阵中该位置的值全为0
#定义全局变量列表finalList存储最后回溯的路径 finalOrder1,finalOrder2存储最后的序列
def createList():
global finalList
global finalOrder1
global finalOrder2
global finalRoad
finalList = []
finalOrder1 = []
finalOrder2 = []
finalRoad = []
#创建A G C T 对应的键值对，方便查找计分矩阵中对应的得分
def createDic():
dic = {'A':0,'G':1,'C':2,'T':3}
return dic
#构建计分矩阵
# A G C T
def createGrade():
grade = np.matrix([[10,-1,-3,-4],
[-1,7,-5,-3],
[-3,-5,9,0],
[-4,-3,0,8]])
return grade
#计算两个字符的相似度得分函数
def getGrade(a,b):
dic = createDic() # 碱基字典 方便查找计分矩阵
grade = createGrade() # 打分矩阵grade
return grade[dic[a],dic[b]]
#构建得分矩阵函数 参数为要比较序列、自定义的gap值
def createMark(s,t,gap):
a = len(s) #获取序列长度a,b
b = len(t)
mark = np.zeros((a+1,b+1)) #初始化全零得分矩阵
direction = np.zeros((a+1,b+1,3)) #direction矩阵用来存储得分矩阵中得分来自的方向 第一个表示左方 第二个表示左上 第三个表示上方 1表示能往哪个方向去
#由于得分可能会来自多个方向，所以使用三维矩阵存
for i in range(1,a+1):
for j in range(1,b+1):
threeMark = [mark[i][j-1],mark[i-1][j-1],mark[i-1][j]] #threeMark表示现在所要计算得分的位置的左边 左上 上边的得分
threeGrade = [gap,getGrade(s[i-1],t[j-1]),gap] #threeGrade表示经过需要计算得左边 左上 上边的空位以及相似度得分
finalGrade = np.add(threeMark,threeGrade) #finalGrade表示最终来自三个方向上的得分
if max(finalGrade) >= 0: #如果该最大值是大于0的则 选取三个方向上的最大得分存入得分矩阵 否则不对矩阵进行修改
mark[i][j] = max(finalGrade)
for k in range(0,len([y for y,x in enumerate(finalGrade) if x == max(finalGrade)])): #可能该位置的得分可以由多个方向得来，所以进行判断并循环赋值
directionList = [y for y,x in enumerate(finalGrade) if x == max(finalGrade)]
direction[i][j][directionList[k]] = 1
return mark,direction
#回溯函数 参数分别为 得分矩阵 方向矩阵 现在所处得分矩阵的位置 以及两个序列
def remount(mark,direction,i,j,s,t):
if direction[i][j][0] == 1 :
if all(direction[i][j-1] == [0,0,0]): #如果该位置指向左边 先判断其左边是否是零点
finalList.append(0) #如果是 将该路径存入路径列表
finalList.reverse() #将列表反过来得到从零点开始的路径
index1 = i #记录现在所匹配序列s的位置 因为两个字符串可能是不一样长的
index2 = j-1 #记录现在所匹配序列t的位置
for k in finalList:
if k == 0 :
finalOrder1.append("-")
finalOrder2.append(t[index2])
index2 += 1
if k == 1 :
finalOrder1.append(s[index1])
finalOrder2.append(t[index2])
index1 += 1
index2 += 1
if k == 2 :
finalOrder1.append(s[index1])
finalOrder2.append("-")
index1 += 1
finalList.reverse()
finalRoad.append(np.array(finalList)) # 将此次的路径添加到最终路径记录用于最后画图
finalList.pop() #输出后将当前方向弹出 并回溯
return
else :
finalList.append(0) #如果不是零点 则将该路径加入路径矩阵，继续往下走
remount(mark,direction,i,j-1,s,t)
finalList.pop() #该方向走完后将这个方向弹出 继续下一轮判断 下面两个大的判断同理
if direction[i][j][1] == 1 :
if all(direction[i-1][j-1] == [0,0,0]):
finalList.append(1)
finalList.reverse() # 将列表反过来得到从零点开始的路径
index1 = i-1 # 记录现在所匹配序列s的位置 因为两个字符串可能是不一样长的
index2 = j-1 # 记录现在所匹配序列t的位置
for k in finalList:
if k == 0 :
finalOrder1.append("-")
finalOrder2.append(t[index2])
index2 += 1
if k == 1 :
finalOrder1.append(s[index1])
finalOrder2.append(t[index2])
index1 += 1
index2 += 1
if k == 2 :
finalOrder1.append(s[index1])
finalOrder2.append("-")
index1 += 1
finalList.reverse()
finalRoad.append(np.array(finalList)) # 将此次的路径添加到最终路径记录用于最后画图
finalList.pop()
return
else :
finalList.append(1)
remount(mark,direction,i-1,j-1,s,t)
finalList.pop()
if direction[i][j][2] == 1 :
if all(direction[i-1][j] == [0,0,0]):
finalList.append(2)
finalList.reverse() # 将列表反过来得到从零点开始的路径
index1 = i-1 # 记录现在所匹配序列s的位置 因为两个字符串可能是不一样长的
index2 = j # 记录现在所匹配序列t的位置
for k in finalList:
if k == 0 :
finalOrder1.append("-")
finalOrder2.append(t[index2])
index2 += 1
if k == 1 :
finalOrder1.append(s[index1])
finalOrder2.append(t[index2])
index1 += 1
index2 += 1
if k == 2 :
finalOrder1.append(s[index1])
finalOrder2.append("-")
index1 += 1
finalList.reverse()
finalRoad.append(np.array(finalList)) # 将此次的路径添加到最终路径记录用于最后画图
finalList.pop()
return
else :
finalList.append(2)
remount(mark,direction,i-1,j,s,t)
finalList.pop()
#画箭头函数
def arrow(ax,sX,sY,aX,aY):
ax.arrow(sX,sY,aX,aY,length_includes_head=True, head_width=0.15, head_length=0.25, fc='w', ec='b')
#画图函数
def drawArrow(mark, direction, a, b, s, t,mx,my):
#a是s的长度为4 b是t的长度为6
fig = plt.figure()
ax = fig.add_subplot(111)
val_ls = range(a+2)
scale_ls = range(b+2)
index_ls = []
index_lsy = []
for i in range(a):
if i == 0:
index_lsy.append('#')
index_lsy.append(s[a-i-1])
index_lsy.append('0')
for i in range(b):
if i == 0:
index_ls.append('#')
index_ls.append('0')
index_ls.append(t[i])
plt.xticks(scale_ls, index_ls) #设置坐标字
plt.yticks(val_ls, index_lsy)
for k in range(1,a+2):
y = [k for i in range(0,b+1)]
x = [x for x in range(1,b+2)]
ax.scatter(x, y, c='y')
for i in range(1,a+2):
for j in range(1,b+2):
ax.text(j,a+2-i,int(mark[i-1][j-1]))
lX = my + 1
lY = a - mx + 1
for n in range(0,len(finalRoad)):
for m in (finalRoad[n]):
if m == 0:
arrow(ax,lX,lY,-1,0)
lX = lX - 1
elif m == 1:
arrow(ax,lX,lY,-1,1)
lX = lX - 1
lY = lY + 1
elif m == 2:
arrow(ax, lX, lY, 0, 1)
lY = lY + 1
lX = b + 1
lY = 1
ax.set_xlim(0, b + 2) # 设置图形的范围，默认为[0,1]
ax.set_ylim(0, a + 2) # 设置图形的范围，默认为[0,1]
ax.set_aspect('equal') # x轴和y轴等比例
plt.show()
plt.tight_layout()
if __name__ == '__main__':
createList()
print("Please enter gap:")
gap = int(input()) #获取gap值 转换为整型 tip：刚开始就是因为这里没有进行类型导致后面的计算部分报错
print("Please enter sequence 1:")
s = input() #获取用户输入的第一条序列
print("Please enter sequence 2:")
t = input() #获取用户输入的第二条序列
a = len(s) #获取s的长度
b = len(t) #获取t的长度
mark,direction = createMark(s,t,gap)
print("The scoring matrix is as follows：") #输出得分矩阵
print(mark)
maxDirection = np.argmax(mark) #获取最大值的位置
i = int(maxDirection/(b+1))
j = int(maxDirection - i*(b+1))
remount(mark,direction,i,j,s,t) #调用回溯函数
print(finalOrder1)
print(finalOrder2)
drawArrow(mark, direction, a, b, s, t, i, j)
二、实验结果截图
1.全局比对
2.局部比对
到此这篇关于详解基于python的全局与局部序列比对的实现(DNA)的文章就介绍到这了,更多相关python全局与局部序列比对内容请搜索脚本之家以前的文章或继续浏览下面的相关文章希望大家以后多多支持脚本之家！
总结
本次实验使用动态规划对全局序列比对进行了实现，自己卡的最久的地方是回溯以及画图的时候。刚开始在实现回溯的过程中，老是找不准回溯的条件以及将所有的路径都记录下来的方法，最后是使用的方向矩阵，也就是重新定义一个与得分矩阵等大的矩阵（但是这个矩阵是三维），存放的是每个位置能够回溯的方向，第一个数值表示左边，第二个表示左上，第三个表示上方，为0时表示当前方向不能回溯，没有路径，为1时表示能回溯，当该位置的所有能走的方向都走完时即可返回。将所有路径记录下来的方法是定义全局变量，当有路径能够走到终点时便将这条路径存放入该全局变量中。
绘图的时候使用的是matplotlib中的散点图，然后将每个点的得分以注释的形式标记在该点的右上角，并用箭头将路径绘出。不得不说的是，这个图确实太丑了，我学识浅薄，也没想到能画出这个图的更好的方法，还希望老师指点。
总的来说这次实验经历的时间还比较长，主要是因为python也没有很熟悉，很多函数也是查了才知道，然后可视化更是了解的少，所以画出来的图出奇的丑，还有回溯的时候也是脑子转不过弯来，所以要学习的东西还有很多，需要更加努力。
本次实验还能够有所改进的地方是：
1.把两个比对算法结合，让用户能够选择使用哪种比对方式。
2.作出一个更好看的界面，增加用户体验感。
3.把图画的更美观。
（老丁已阅，USC的同学们谨慎借鉴）
利用python处理dna序列_详解基于python的全局与局部序列比对的实现(DNA)相关推荐
- python selenium爬虫_详解基于python +Selenium的爬虫
详解基于python +Selenium的爬虫 一.背景 1. Selenium Selenium 是一个用于web应用程序自动化测试的工具,直接运行在浏览器当中,支持chrome.firefox等主 ...
- python编写数据库连接工具_详解使用Python写一个向数据库填充数据的小工具(推荐)...
一. 背景 公司又要做一个新项目,是一个合作型项目,我们公司出web展示服务,合作伙伴线下提供展示数据. 而且本次项目是数据统计展示为主要功能,并没有研发对应的数据接入接口,所有展示数据源均来自数据库 ...
- 用python写管理系统局域网_详解用python -m http.server搭一个简易的本地局域网
工作时同事间几mb小文件的传输,一般使用QQ或者微信就足够了,但当传输文件几百MB或者几十G时,这种方法的效率就显得不足了.本篇就是简单说明一个python小功能,让大家能利用python方便的搭建一 ...
- python解释器在哪里_详解查看Python解释器路径的两种方式
进入python的安装目录, 查看python解释器 进入bin目录 # ls python(看一下是否有python解释器版本) # pwd (查看当前目录) 复制当前目录即可 1. 通过脚本查看 ...
- python登录系统文件版_详解用python实现基本的学生管理系统(文件存储版)（python3）...
这个是升级的版本,比较进阶一点的,相对与之前的文件管理系统,数据只是存储到了内存中,而不是存储到硬盘上,我们想让文件存储到硬盘上的话,一个是存储到文件里,一个是存储到数据库中,存储到数据库中的版本会后 ...
- python 写入excel 日期_详解：Python实现大数据收集至excel的思路大牛分享（建议收藏）...
一.在工程目录中新建一个excel文件 二.使用python脚本程序将目标excel文件中的列头写入,本文省略该部分的code展示,可自行网上查询 三.以下code内容为:实现从接口获取到的数据值写入 ...
- python微博评论爬虫_详解用python写网络爬虫-爬取新浪微博评论 基于Python的新浪微博爬虫研究...
怎样爬取新浪微博的评论信息 针对八爪鱼在微博的应用上,除了用户信息之外还包括话题内容方面的采集,目前绝大多数企业均在微博设有官方微博,八爪鱼可以协助企业快速及时的抓取与企业产品相关联的话题信息,规则市 ...
- python生成随机数方法_详解用python生成随机数的几种方法
今天学习了用python生成仿真数据的一些基本方法和技巧,写成博客和大家分享一下. 本篇博客主要讲解如何从给定参数的的正态分布/均匀分布中生成随机数以及如何以给定概率从数字列表抽取某数字或从区间列表的 ...
- python numpy 子数组_详解：Python 取numpy数组的某几行某几列方法（含对与错示例）...
前言: 今天为大家带来的内容是Python 取numpy数组的某几行某几列方法,文中通过示例代码介绍的非常详细,对大家的学习或者工作具有一定的参考学习价值,要是喜欢的话记得点赞转发收藏不迷路哦!!! ...
最新文章
- Uva 10537 过路费
- Linux minicom 软件退出方法
- Linux Mint 18.2 Sonya 将于 6 月上旬发布公测版本
- 网络安全分析 | 用OpenFEA定位WebShell木马后门
- 2016 Multi-University Training Contests
- Mac上好用的web前端开发工具
- web 打印时 表格多页 导致缺少线段 看上去页面不完整的处理方案
- 动作捕捉用于索并联机构中的理论验证
- 洞态IAST Agent正式开源
- 对比汇新云和猪八戒？
- php后台管理修改密码,重置网站后台管理员密码
- 两台计算机如何共享文档,两台电脑怎么共享文件夹？
- SpringCloud整合Sa-Token登录认证+Gateway网关拦截
- 小程序＋动易SF制作随手拍实例全景式操作(3)
- 哔哩哔哩2020校园招聘前端笔试题（卷一）
- 学习笔记——神经网络压缩
- 使用python实现图像对比度增强
- activiti流程图查看
- dll转换成.a /dll转换lib
- 儿子的读后感：有爱，有趣——读《大地之歌》有感