title: How do I decide which way to backtrack in the Smith–Waterman algorithm?
body: <p>I am trying to implement local sequence alignment in Python using the <a href="http://en.wikipedia.org/wiki/Smith%E2%80%93Waterman_algorithm" rel="nofollow">Smith–Waterman algorithm</a>.</p>
  
  <p>Here's what I have so far. It gets as far as building the <a href="http://en.wikipedia.org/wiki/Similarity_matrix" rel="nofollow">similarity matrix</a>:</p>
  
  <pre><code>import sys, string
  from numpy import *
  
  f1=open(sys.argv[1], 'r')
  seq1=f1.readline()
  f1.close()
  seq1=string.strip(seq1)
  
  f2=open(sys.argv[2], 'r')
  seq2=f2.readline()
  f2.close()
  seq2=string.strip(seq2)
  
  a,b =len(seq1),len(seq2)
  
  penalty=-1;
  point=2;
  
  #generation of matrix for local alignment
  p=zeros((a+1,b+1))
  
  # table calculation and matrix generation
  for i in range(1,a+1):
      for j in range(1,b+1):
          vertical_score =p[i-1][j]+penalty;
          horizontal_score= p[i][j-1]+penalty;
          if seq1[i-1]==seq2[j-1]:
              diagonal_score =p[i-1][j-1]+point;
          else:
              diagonal_score = p[i-1][j-1]+penalty;
          p[i][j]=max(0,vertical_score,horizontal_score,diagonal_score);
  
  print p
  </code></pre>
  
  <p>For example, with the two sequences:</p>
  
  <pre><code>agcacact
  acacacta
  </code></pre>
  
  <p>my code outputs the similarity matrix:</p>
  
  <pre><code>[[  0.   0.   0.   0.   0.   0.   0.   0.   0.]
   [  0.   2.   1.   2.   1.   2.   1.   0.   2.]
   [  0.   1.   1.   1.   1.   1.   1.   0.   1.]
   [  0.   0.   3.   2.   3.   2.   3.   2.   1.]
   [  0.   2.   2.   5.   4.   5.   4.   3.   4.]
   [  0.   1.   4.   4.   7.   6.   7.   6.   5.]
   [  0.   2.   3.   6.   6.   9.   8.   7.   8.]
   [  0.   1.   4.   5.   8.   8.  11.  10.   9.]
   [  0.   0.   3.   4.   7.   7.  10.  13.  12.]]
  </code></pre>
  
  <p>Now I am stuck on the next step of the algorithm, the backtracking to construct the best alignment.</p>
  
  <p><a href="http://en.wikipedia.org/wiki/Smith%E2%80%93Waterman_algorithm#Example" rel="nofollow">Wikipedia says</a>,</p>
  
  <blockquote>
    <p>To obtain the optimum local alignment, we start with the highest value in the matrix (<em>i</em>, <em>j</em>). Then, we go backwards to one of positions (<em>i</em> − 1, <em>j</em>), (<em>i</em>, <em>j</em> − 1), and (<em>i</em> − 1, <em>j</em> − 1) depending on the direction of movement used to construct the matrix. We keep the process until we reach a matrix cell with zero value, or the value in position (0, 0).</p>
  </blockquote>
  
  <p>I am having trouble determining which position to backtrack to. What does Wikipedia mean by "depending on the direction of movement used to construct the matrix" mean and how would I implement this in Python?</p>
  
  <p>Finally I did this </p>
  
  <pre><code>import sys, string
  from numpy import*
  import re
  
  # read first sequence
  
  fasta_sequence1=open(sys.argv[1], 'r')
  
  seq1=""
  for i in fasta_sequence1:
      if i.startswith("&gt;"):
          pass
      else:
          seq1 = seq1 + i.strip()   
  fasta_sequence1.close()
  
  fasta_sequence2=open(sys.argv[2], 'r')
  
  seq2 = ""
  for i in fasta_sequence2:
      if i.startswith('&gt;'):
          pass
      else:
          seq2 = seq2+ i.strip()
  fasta_sequence2.close()
  
  a,b =len(seq1),len(seq2)
  
  
  penalty=-1;
  point=2;
  
  #generation of matrix for local alignment 
  
  p=zeros((a+1,b+1))
  
  #intialization of max score
  max_score=0;
  #pointer to store the traceback path
  
  pointer=zeros((a+1,b+1))
  
  # table calculation and matrix generation
  for i in range(1,a+1):
      for j in range(1,b+1):
          vertical_score =p[i-1][j]+penalty;
          horizontal_score= p[i][j-1]+penalty;
          if seq1[i-1]==seq2[j-1]:
              diagonal_score =p[i-1][j-1]+point;
          else:
              diagonal_score = p[i-1][j-1]+penalty;
  
  for i in range(1,a+1):
      for j in range(1,b+1):            
          p[i][j]=max(0,vertical_score,horizontal_score,diagonal_score);
  
  for i in range(1,a+1):
      for j in range(1,b+1):
          if p[i][j]==0:
              pointer[i][j]=0; #0 means end of the path
          if p[i][j]==vertical_score:
              pointer[i][j]=1; #1 means trace up
          if p[i][j]==horizontal_score:
              pointer[i][j]=2; #2 means trace left
          if p[i][j]==diagonal_score:
              pointer[i][j]=3; #3 means trace diagonal
          if p[i][j]&gt;=max_score:
              maximum_i=i;
              maximum_j=j;
              max_score=p[i][j];
  
  
  #for i in range(1,a+1):
     # for j in range(1,b+1):
          #if p[i][j]&gt;max_score:
              #max_score=p[i][j]
  
  print max_score
  
  # traceback of all the pointers 
  
  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;
  
  
  
  align1=align1[::-1]; #reverse sequence 1
  align2=align2[::-1]; #reverse sequence 2
  
  #output_file = open(sys.argv[3],"w")
  #output_file.write(align1)
  
  #output_file.write(align2)
  
  print align1
  print align2
  </code></pre>
  
  <p>But i think there could be a better and more efficient way of doing this</p>
  
  <pre><code>output_file = open(sys.argv[3],"w")
  output_file.write(align1)
  output_file.write(align2)
  </code></pre>
  
  <p>the result shows like </p>
  
  <pre><code>agcacacta-cacact
  </code></pre>
  
  <p>on the contraryc for :</p>
  
  <pre><code>print align1
  print align2
  </code></pre>
  
  <p>shows correct output:</p>
  
  <pre><code>agcacact
  a-cacact
  </code></pre>
  
  <p>How can i get the above output in the file writer . Thanks</p>