最近在自学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只用一个语句就可以完成相同任务。
以上就是本文的全部内容，希望对大家的学习有所帮助，也希望大家多多支持脚本之家。
