任选一种编程语言,设计一个双序列全局比对的程序。要求:
1) 输入两条蛋白质序列,输出比对结果例如:
Alignment Score: 12345
E E E E E K K K K K A A A A A F F F
E E E E E - - - - - B B B B B F F F
2) 使用BLOSUM62矩阵来对序列中氨基酸符号的match和mismatch打分。
3) 使用动态规划的方法(Needleman-Wunsch算法)。
4) 计算空位(gap)的罚分时,使用仿射空位罚分策略,gap opening的分数为-11,gap extension的分数为-1。
5) 测试程序时,使用Horse Hemoglobin subunit theta-1(马的theta-1亚基血红蛋白)序列和Bullfrog Hemoglobin subunit alpha-type(牛蛙的alpha型亚基血红蛋白)序列。
6) 将程序输出结果与BLAST软件的进行比较。
登录BLAST在线服务器http://blast.ncbi.nlm.nih.gov/Blast.cgi。

代码如下:

global S
global E
global MIN
global amino
global blosum
S   = -11
E   = -1
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):index1=amino.index(t[i-1])index2=amino.index(s[j-1])return blosum[index1][index2]#initializers for matrices
def _init_x(i, j):if i > 0 and j == 0:return MINelse:if j > 0:return S + E * jelse:return 0def _init_y(i, j):if j > 0 and i == 0:return MINelse:if i > 0:return S + E * ielse:return 0def _init_m(i, j):if j == 0 and i == 0:return 0else:if j == 0 or i == 0:return MINelse:return 0def distance_matrix(s, t):dim_i = len(t) + 1dim_j = len(s) + 1#abuse list comprehensions to create matricesX = [[_init_x(i, j) for j in range(0, dim_j)] for i in range(0, dim_i)]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)]for j in range(1, dim_j):for i in range(1, dim_i):X[i][j] = max((S + E + M[i][j-1]), (E + X[i][j-1]))Y[i][j] = max((S + E + M[i-1][j]), (E + Y[i-1][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):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 -= 1j -= 1elif (i>0 and M[i][j] == Y[i][j]):sequ1 += '_'sequ2 += t[i-1]i -= 1elif (j>0 and M[i][j] == X[i][j]):sequ1 += s[j-1]sequ2 += '_'j -= 1sequ1r = ' '.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:")[X, Y, M] = distance_matrix(seq1, seq2)
[str1, str2] = backtrace(seq1, seq2, X, Y, M)score=M[len(seq2)][len(seq1)]
print("Alignment Score:"+ str(score))
print(str1)
print(str2)

生物信息学:任选一种编程语言,设计一个双序列全局比对的程序相关推荐

  1. 假如一种编程语言代表一个国家,哎哟,有意思了!

    "如果每种编程语言代表一个国家" Dimage Sapelkin 这是多么有趣的一个问题啊! ◆  ◆  ◆ C-俄罗斯 一切都必须以一种倒退的方式进行,但一切皆有可能,而且有很多 ...

  2. 设计一个基于GUI的扑克程序

    2019独角兽企业重金招聘Python工程师标准>>> 在本课程教材扑克牌代码的基础上,设计一个基于GUI的扑克程序 a) 可以显示 52 张扑克牌,包括洗牌,发牌在内(2) b) ...

  3. c语言双序列全局比对,基于动态规划进行双序列全局比对

    说明 核酸序列打分算法脚本,基于动态规划进行双序列全局比对,得到两条DNA序列的相似度并打分,但程序还有一些问题,在匹配长序列的时候还不够完善. 环境 Linux.Python3.6 实例 comma ...

  4. Pyqt5入门--用qtdesigner设计一个计算屏幕PPI小程序(qtdesigner/pyuic/pyinstaller/python)

    本文利用python中的pyqt5包,设计一个计算PPI小程序的界面,再利用pyuic将界面的ui文件转为py文件.再新建一个py文件继承界面py文件中类,并定义每一个按钮对应的函数,完成后利用pyi ...

  5. HDLBits 系列(15) 如何设计一个双边沿采样的电路?

    目录 背景 原题复现 审题 我的设计1 我的设计2 背景 曾经专门写过这个话题,可是今天在练习HDLBits时候,又发现了这个问题,但是以前的思路我已经忘了,不得不回顾. FPGA中如何实现双边沿采样 ...

  6. vb设计一个由计算机,计算机VB程序的设计第一章.ppt

    Visual Basic程 序 设 计 ;1.初期的程序设计 高运行效率.少占用内存为目标2.结构化程序设计程序的可读性.可维护性为目标 程序=算法+数据结构 的面向过程的程序设计3.面向对象的程序设 ...

  7. 用jk触发器设计一个011序列检测器的设计分析过程

    心得体会:经过此次设计,加深了对时序逻辑电路的理解,当要求对一个连续的一串信号进行输入输出处理时可以用到有记忆存储.反馈功能的jk触发器或者d触发器.

  8. 用74194设计一个00011101序列信号发生器

    序列信号发生器可由移位寄存器(74194)和反馈逻辑电路构成,如下图所示: 一.列出所要产生的序列(周期为8,最右边信号先输出)与移位寄存器状态的关系 如上图,序列下面的线段所对应的数码表示移位寄存器 ...

  9. 用c语言设计一个统计字符个数的程序,「第6篇」「C程序上机题」「统计输入的字符个数思路与实现」...

    一.统计输入的字符个数 同学们在学习C语言课程中,经常会遇到一道题,就是要求你写一个C程序,这个C程序能够读取你从键盘上输入的字符,并且统计其中的字符个数,最后输出总的字符个数并且把这些输入的字符再输 ...

最新文章

  1. python 3.5-python3.5
  2. LOJ#2087 国王饮水记
  3. 纪事本 乱码_纪事地图和Yahoo Cloud服务基准
  4. Spring4.x集成xfire1.26 问题汇总
  5. oracle批量替换保留字,oracle保留字大全
  6. devops测试_使用DevOps管道自动执行用户验收测试
  7. php fpm 内存增加,不断增加php-fpm的内存使用量?
  8. 实现简单的Java内存缓存
  9. 利用vue-gird-layout 制作可定制桌面 (一)
  10. Qt_解决Qt5.15 + Xcode12iOS端qmake不可用的问题
  11. 变色龙引导安装黑苹果 遇到的问题的解决办法
  12. 协方差矩阵的计算方法
  13. sql服务器字段顺序怎么修改,你可能不知道SQL Server索引列的升序和降序带来的性能问题...
  14. 值得反复研读的表连接之CARTESIAN JOIN方式
  15. 编译php8,Centos编译体验PHP8 Alpha 2
  16. eclipse如何去掉无用的validation、优化eclipse
  17. 网易企业邮箱用Python发邮件
  18. 【基础代码】python 一些常用的基础代码
  19. 区块链开发者观点:来自 MYKEY 的胖哥 Ricky
  20. Threes_位置变动

热门文章

  1. Ophone平台2D游戏引擎实现——物理引擎(一)(二)
  2. 如何在自己的网站上添加HTML悬浮音乐播放器?
  3. 基于单片机的智能窗控制系统设计(电路+流程)
  4. Spring自动装配 autowire
  5. 根据标签进行群发php,根据标签进行群发
  6. netstat命令实战详解
  7. 好用实用稳定的API接口【物流快递篇】
  8. 从信息包围到信息追踪,你在网络中还剩下多少自由?
  9. 叫我修吧全新升级 品牌售后保障更全
  10. TMS320F28335 CAN通信