使用python實現(xiàn)BLAST的方法

這篇文章將為大家詳細(xì)講解有關(guān)使用python實現(xiàn)BLAST的方法，小編覺得挺實用的，因此分享給大家做個參考，希望大家閱讀完這篇文章后可以有所收獲。

這次更新了打分函數(shù)如下，空位罰分改為-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)

運行結(jié)果：

無疑python對字符串的處理更加強(qiáng)大，語言也更加簡單，優(yōu)雅。比如最后逆序輸出alignment，java我是單獨寫了一個逆序函數(shù)，而python只用一個語句就可以完成相同任務(wù)。

關(guān)于“使用python實現(xiàn)BLAST的方法”這篇文章就分享到這里了，希望以上內(nèi)容可以對大家有一定的幫助，使各位可以學(xué)到更多知識，如果覺得文章不錯，請把它分享出去讓更多的人看到。