### 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 *

f1=open(sys.argv[1], 'r')
seq1=string.strip(seq1)

f2=open(sys.argv[2], 'r')
seq2=string.strip(seq2)

f3=open(sys.argv[3],'r')
BLOSUM62=[]
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

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