python

超轻量级php框架startmvc

使用python实现BLAST

更新时间:2020-05-22 04:54 作者:startmvc
最近在自学python,又用python实现了一下BLAST。这次更新了打分函数如下,空位罚分改为-5,

最近在自学python,又用python实现了一下BLAST。

这次更新了打分函数如下,空位罚分改为-5,但不区分gap open 和 gap extend。


''''' 
@author: JiuYu 
''' 
 
def score(a,b):#scoring function 
 score=0 
 lst=['AC','GT','CA','TG'] 
 if a==b: 
 score +=2 
 elif a+b in lst: 
 score += -5 
 else: 
 score += -7 
 return score 
 
def BLAST(seq1,seq2):#Basic Local Alignment Search Tool 
 l1 = len(seq1) 
 l2 = len(seq2) 
 GAP =-5 #-5 for any gap 
 scores =[] 
 point =[] 
 
 for j in range(l2+1): 
 if j == 0: 
 line1=[0] 
 line2=[0] 
 for i in range(1,l1+1): 
 line1.append(GAP*i) 
 line2.append(2) 
 else: 
 line1=[] 
 line2=[] 
 line1.append(GAP*j) 
 line2.append(3) 
 scores.append(line1) 
 point.append(line2) 
 
 #fill the blank of scores and point 
 for j in range(1,l2+1): 
 letter2 = seq2[j-1] 
 for i in range(1,l1+1): 
 letter1 = seq1[i-1] 
 diagonal_score = score(letter1, letter2) + scores[j-1][i-1] 
 left_score = GAP + scores[j][i-1] 
 up_score = GAP + scores[j-1][i] 
 max_score = max(diagonal_score, left_score, up_score) 
 scores[j].append(max_score) 
 
 if scores[j][i] == diagonal_score: 
 point[j].append(1) 
 elif scores[j][i] == left_score: 
 point[j].append(2) 
 else: 
 point[j].append(3) 
 
 #trace back 
 alignment1='' 
 alignment2='' 
 i = l2 
 j = l1 
 print 'scores =',scores[i][j] 
 while True: 
 if point[i][j] == 0: 
 break 
 elif point[i][j] == 1: 
 alignment1 += seq1[j-1] 
 alignment2 += seq2[i-1] 
 i -= 1 
 j -= 1 
 elif point[i][j] == 2: 
 alignment1 += seq1[j-1] 
 alignment2 += '-' 
 j -= 1 
 else: 
 alignment1 += '-' 
 alignment2 += seq2[i-1] 
 i -= 1 
 
 #reverse alignment 
 alignment1 = alignment1[::-1] 
 alignment2 = alignment2[::-1] 
 print 'The best alignment:' 
 print alignment1 
 print alignment2 
 
seq1=raw_input('Please input your first sequences:\n') 
seq2=raw_input('input second sequences:\n') 
BLAST(seq1, seq2) 

运行结果:

无疑python对字符串的处理更加强大,语言也更加简单,优雅。比如最后逆序输出alignment,java我是单独写了一个逆序函数,而python只用一个语句就可以完成相同任务。

以上就是本文的全部内容,希望对大家的学习有所帮助,也希望大家多多支持脚本之家。