使用python实现BLAST

最近在自学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只用一个语句就可以完成相同任务。

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