2007-09-25

Smith-Waterman Algorithm in Python (for amino acid sequences)



This code is for amino acid sequences rather than nucleic acid sequences. But you can modify it slightly, such as changing the scoring matrix, to make it work for nucleic acid sequences.

Please note that I wrote this code when I first started programming in Python. Yes, the code is ugly. I will update later when I have time. And, I still do not like OOP.

#This software is a free software. Thus, it is licensed under GNU General Public License version 3 or later.
#Python implementation to Smith-Waterman Algorithm for Homework 1 of Bioinformatics class.
#Forrest Bao, Sept. 26, 2007 <http://fsbao.net> <forrest.bao aT gmail.com>

import sys, string
from numpy import *
from matplotlib import *

#read the first sequence
f1=open(sys.argv[1], 'r')
seq1=f1.readline()
seq1=string.strip(seq1)

#read the second sequence
f2=open(sys.argv[2], 'r')
seq2=f2.readline()
seq2=string.strip(seq2)

#read in BLOSUM62 matrix
f3=open(sys.argv[3],'r')
BLOSUM62=[]
for line in f3.readlines():
  BLOSUM62.append(map(int, line.split()))

#define similar amino acids
similarAA=['ST','TS','SP','PS','SA','AS','SG','GS','TP','PT','TA','AT','TG','GT','PA','AP','PG','GP','AG','GA','ND','DN','NE','EN','NQ','QN','DE','ED','DQ',
'QD','EQ','QE','HR','RH','HK','KH','RK','KR','MI','IM','ML','LM','MV','VM','IL','LI','IV','VI','LV','VL','FY','YF','FW','WF','YW','WY'];

m,n =  len(seq1),len(seq2) #length of two sequences

penalty=-4;   #define the gap penalty

#generate DP table and traceback path pointer matrix
score=zeros((m+1,n+1))   #the DP table
pointer=zeros((m+1,n+1))  #to store the traceback path

P=0;

def match_score(alpha,beta,BLOSUM62): #the function to find match/dismatch score from BLOSUM62 by letters of AAs
  alphabet={}
  alphabet["A"] = 0
  alphabet["R"] = 1
  alphabet["N"] = 2
  alphabet["D"] = 3
  alphabet["C"] = 4
  alphabet["Q"] = 5
  alphabet["E"] = 6
  alphabet["G"] = 7
  alphabet["H"] = 8
  alphabet["I"] = 9
  alphabet["L"] = 10
  alphabet["K"] = 11
  alphabet["M"] = 12
  alphabet["F"] = 13
  alphabet["P"] = 14
  alphabet["S"] = 15
  alphabet["T"] = 16
  alphabet["W"] = 17
  alphabet["Y"] = 18
  alphabet["V"] = 19
  alphabet["B"] = 20
  alphabet["Z"] = 21
  alphabet["X"] = 22
  lut_x=alphabet[alpha]
  lut_y=alphabet[beta]
  return BLOSUM62[lut_x][lut_y]

max_score=P;  #initial maximum score in DP table

#calculate DP table and mark pointers
for i in range(1,m+1):
  for j in range(1,n+1):
    score_up=score[i-1][j]+penalty;
    score_down=score[i][j-1]+penalty;
    score_diagonal=score[i-1][j-1]+match_score(seq1[i-1],seq2[j-1],BLOSUM62);
    score[i][j]=max(0,score_up,score_down,score_diagonal);
    if score[i][j]==0:
      pointer[i][j]=0; #0 means end of the path
    if score[i][j]==score_up:
      pointer[i][j]=1; #1 means trace up
    if score[i][j]==score_down:
      pointer[i][j]=2; #2 means trace left
    if score[i][j]==score_diagonal:
      pointer[i][j]=3; #3 means trace diagonal
    if score[i][j]>=max_score:
      max_i=i;
      max_j=j;
      max_score=score[i][j];
#END of DP table


align1,align2='',''; #initial sequences

i,j=max_i,max_j; #indices of path starting point

#traceback, follow pointers
while pointer[i][j]!=0:
  if pointer[i][j]==3:
    align1=align1+seq1[i-1];
    align2=align2+seq2[j-1];
    i=i-1;
    j=j-1;
  elif pointer[i][j]==2:
    align1=align1+'-';
    align2=align2+seq2[j-1];
    j=j-1;
  elif pointer[i][j]==1:
    align1=align1+seq1[i-1];
    align2=align2+'-';
    i=i-1;
#END of traceback

align1=align1[::-1]; #reverse sequence 1
align2=align2[::-1]; #reverse sequence 2

i,j=0,0;

#calcuate identity, similarity, score and aligned sequeces
def result(align1,align2,BLOSUM62):
  symbol='';
  found=0;
  score=0;
  penalty=-4;
  identity,similarity=0,0;
  for i in range(0,len(align1)):
    if align1[i]==align2[i]:     #if two AAs are the same, then output the letter
      symbol=symbol+align1[i];
      identity=identity+1;
      similarity=similarity+1;
      score=score+match_score(align1[i],align2[i],BLOSUM62);
    elif align1[i]!=align2[i] and align1[i]!='-' and align2[i]!='-': #if they are not identical and none of them is gap
      score=score+match_score(align1[i],align2[i],BLOSUM62);
      for j in range(0,len(similarAA)-1):
        if align1[i]+align2[i]==similarAA[j]: #search whether these two AAs form a pair in similar AA table
          found=1;
      if found==1:
        symbol=symbol+':';   #if they are similar AA, output ':'
        similarity=similarity+1;
      if found==0:
        symbol=symbol+' ';  #o/w, output a space
      found=0;
    elif align1[i]=='-' or align2[i]=='-':   #if one of them is a gap, output a space
      symbol=symbol+' ';
      score=score+penalty;

  identity=float(identity)/len(align1)*100;
  similarity=float(similarity)/len(align2)*100;

  print 'Identity =', "%3.3f" % identity, 'percent';
  print 'Similarity =', "%3.3f" % similarity, 'percent';
  print 'Score =', score;
  print align1
  print symbol
  print align2
  #END of function result

  result(align1,align2,BLOSUM62)

The execution is like this:
forrest@rainbow:/forrest/work/BioInfomatics/SW$ python SW.py Q9GMA3 Q9M276 BLOSUM62
Identity = 35.294 percent
Similarity = 52.941 percent
Score = 118
KRKKRRHRTVFTAHQLEELEKAF-SEAHY-P-D-VY-AREMLAMKTELPEDRIQVWFQNRRAKWR-KR-EKRWGGSSVMAEYG-L
K KK      F:  Q:  LE  F SE::  P   V  ARE L::  : P   : :WFQN:RA:W: K  EK     :: A:Y  L
KMKKSNNQKRFSEEQIKSLELIFESETRLEPRKKVQVARE-LGL--Q-PR-QVAIWFQNKRARWKTKQLEKEY--NTLRANYNNL

Update Jun. 20, 2010: Several users have written to me regarding the format of my BLOSUM62. It is as follows:
  4 -1 -2 -2  0 -1 -1  0 -2 -1 -1 -1 -1 -2 -1  1  0 -3 -2  0 -2 -1  0 -4 
 -1  5  0 -2 -3  1  0 -2  0 -3 -2  2 -1 -3 -2 -1 -1 -3 -2 -3 -1  0 -1 -4 
 -2  0  6  1 -3  0  0  0  1 -3 -3  0 -2 -3 -2  1  0 -4 -2 -3  3  0 -1 -4 
 -2 -2  1  6 -3  0  2 -1 -1 -3 -4 -1 -3 -3 -1  0 -1 -4 -3 -3  4  1 -1 -4 
  0 -3 -3 -3  9 -3 -4 -3 -3 -1 -1 -3 -1 -2 -3 -1 -1 -2 -2 -1 -3 -3 -2 -4 
 -1  1  0  0 -3  5  2 -2  0 -3 -2  1  0 -3 -1  0 -1 -2 -1 -2  0  3 -1 -4 
 -1  0  0  2 -4  2  5 -2  0 -3 -3  1 -2 -3 -1  0 -1 -3 -2 -2  1  4 -1 -4 
  0 -2  0 -1 -3 -2 -2  6 -2 -4 -4 -2 -3 -3 -2  0 -2 -2 -3 -3 -1 -2 -1 -4 
 -2  0  1 -1 -3  0  0 -2  8 -3 -3 -1 -2 -1 -2 -1 -2 -2  2 -3  0  0 -1 -4 
 -1 -3 -3 -3 -1 -3 -3 -4 -3  4  2 -3  1  0 -3 -2 -1 -3 -1  3 -3 -3 -1 -4 
 -1 -2 -3 -4 -1 -2 -3 -4 -3  2  4 -2  2  0 -3 -2 -1 -2 -1  1 -4 -3 -1 -4 
 -1  2  0 -1 -3  1  1 -2 -1 -3 -2  5 -1 -3 -1  0 -1 -3 -2 -2  0  1 -1 -4 
 -1 -1 -2 -3 -1  0 -2 -3 -2  1  2 -1  5  0 -2 -1 -1 -1 -1  1 -3 -1 -1 -4 
 -2 -3 -3 -3 -2 -3 -3 -3 -1  0  0 -3  0  6 -4 -2 -2  1  3 -1 -3 -3 -1 -4 
 -1 -2 -2 -1 -3 -1 -1 -2 -2 -3 -3 -1 -2 -4  7 -1 -1 -4 -3 -2 -2 -1 -2 -4 
  1 -1  1  0 -1  0  0  0 -1 -2 -2  0 -1 -2 -1  4  1 -3 -2 -2  0  0  0 -4 
  0 -1  0 -1 -1 -1 -1 -2 -2 -1 -1 -1 -1 -2 -1  1  5 -2 -2  0 -1 -1  0 -4 
 -3 -3 -4 -4 -2 -2 -3 -2 -2 -3 -2 -3 -1  1 -4 -3 -2 11  2 -3 -4 -3 -2 -4 
 -2 -2 -2 -3 -2 -1 -2 -3  2 -1 -1 -2 -1  3 -3 -2 -2  2  7 -1 -3 -2 -1 -4 
  0 -3 -3 -3 -1 -2 -2 -3 -3  3  1 -2  1 -1 -2 -2  0 -3 -1  4 -3 -2 -1 -4 
 -2 -1  3  4 -3  0  1 -1  0 -3 -4  0 -3 -3 -2  0 -1 -4 -3 -3  4  1 -1 -4 
 -1  0  0  1 -3  3  4 -2  0 -3 -3  1 -1 -3 -1  0 -1 -3 -2 -2  1  4 -1 -4 
  0 -1 -1 -1 -2 -1 -1 -1 -1 -1 -1 -1 -1 -1 -2  0  0 -2 -1 -1 -1 -1 -1 -4 
 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4  1 

No comments: