最近在自學python,又用python實現了一下BLAST。
這次更新了打分函數如下,空位罰分改為-5,但不區(qū)分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對字符串的處理更加強大,語言也更加簡單,優(yōu)雅。比如最后逆序輸出alignment,java我是單獨寫了一個逆序函數,而python只用一個語句就可以完成相同任務。
以上就是本文的全部內容,希望對大家的學習有所幫助,也希望大家多多支持腳本之家。
更多文章、技術交流、商務合作、聯系博主
微信掃碼或搜索:z360901061
微信掃一掃加我為好友
QQ號聯系: 360901061
您的支持是博主寫作最大的動力,如果您喜歡我的文章,感覺我的文章對您有幫助,請用微信掃描下面二維碼支持博主2元、5元、10元、20元等您想捐的金額吧,狠狠點擊下面給點支持吧,站長非常感激您!手機微信長按不能支付解決辦法:請將微信支付二維碼保存到相冊,切換到微信,然后點擊微信右上角掃一掃功能,選擇支付二維碼完成支付。
【本文對您有幫助就好】元

