title: Traceback in dynamic programming implementation of Needleman-Wunsch algorithm
body: <p>I almost have my needleman-wunsch implementation working but I am confused on how to handle the traceback on a specific case.</p>
  
  <p>The idea is that in order to re construct the sequence (the longest path) we re-calculate to determine the matrix the score came from. The edge case I am having a problem with is when the bottom right score is not in the match matrix, but is in the insert column matrix (meaning that the resulting traced back sequence should have a insert.</p>
  
  <p>These sequences are being recorded in the a2m format, where inserts in the sequence are recorded as a lower case character. So in the final output the alignment of <code>ZZ</code> to <code>AAAC</code> should be <code>AAac</code>. When I trace back by hand I end up with <code>AAAc</code> because I only visit the Ic matrix once. <a href="http://imgur.com/CveocHz">Here</a> is a picture of my whiteboard. As you can see, I have three black arrows and one green arrow, which is why my traceback gives me <code>AAAc</code>. Should I be counting the first cell, then stopping at position 1,1? I am not sure how I would change the way I implemented this to do so.</p>
  
  <p>Note that the substitution matrix being used here is BLOSUM62. The recurrence relations are</p>
  
  <pre><code>M(i,j) = max(Ic(i-1,j-1)+subst, M(i-1,j-1)+subst, Ir(i-1,j-1)+subst)
  Ic(i,j) = max(Ic(i,j-1)-extend, M(i,j-1)-open, Ir(i,j-1)-double)
  Ir(i,j) = max(Ic(i-1,j)-double, M(i-1,j)-open, Ir(i-1,j)-extend)
  </code></pre>
  
  <p>EDIT: Here is the traceback_col_seq function re-written to be cleaner. Note that score_cell now returns thisM,thisC,thisR instead of the max of these. This version scores the alignment as AaAc, still having the same problem, and now with another problem of why does it go into Ic again at 1,2. This version is much more legible however.</p>
  
  <pre><code>def traceback_col_seq(self):
      i, j = self.maxI-1, self.maxJ-1
      self.traceback = list()
      matrixDict = {0:'M',1:'Ic',2:'Ir',3:'M',4:'Ic',5:'Ir',6:'M',7:'Ic',8:'Ir'}
      while i &gt; 0 or j &gt; 0:
          chars = self.col_seq[j-1] + self.row_seq[i-1]
          thisM, thisC, thisR = self.score_cell(i, j, chars)
          cell = thisM + thisC + thisR
          prevMatrix = matrixDict[cell.index(max(cell))]
          print(cell, prevMatrix,i,j)
          if prevMatrix == 'M':
              i -= 1; j -= 1
              self.traceback.append(self.col_seq[j])
          elif prevMatrix == 'Ic':
              j -= 1
              self.traceback.append(self.col_seq[j].lower())
          elif prevMatrix == 'Ir':
              i -= 1
              self.traceback.append('-')
      return ''.join(self.traceback[::-1])
  </code></pre>
  
  <p>Here is the python class that generates the dynamic programming matrix and traces back the alignment. There is also a score function used to check if the alignment produced is correct.</p>
  
  <pre><code>class global_aligner():
      def __init__(self, subst, open=12, extend=1, double=3):
          self.extend, self.open, self.double, self.subst = extend, open, double, subst
      def __call__(self, row_seq, col_seq):
          #add alphabet error checking?
          score_align(row_seq, col_seq)
          return traceback_col_seq()
      def init_array(self):
          """initialize three numpy arrays, set values in 1st column and row"""
          self.M = zeros((self.maxI, self.maxJ), float)
          self.Ic = zeros((self.maxI, self.maxJ), float)
          self.Ir = zeros((self.maxI, self.maxJ), float)
          for i in xrange(1,self.maxI):
              self.M[i][0], self.Ic[i][0], self.Ir[i][0] = \
                      -float('inf'), -float('inf'), -(self.open+self.extend*(i-1))
          for j in xrange(1,self.maxJ):
              self.M[0][j], self.Ir[0][j], self.Ic[0][j] = \
                      -float('inf'), -float('inf'), -(self.open+self.extend*(j-1))
          self.Ic[0][0] = self.Ir[0][0] = -float('inf')
      def score_cell(self, i, j, chars):
          """score a matrix cell based on the 9 total neighbors (3 each direction)"""
          thisM = [self.M[i-1][j-1]+self.subst[chars], self.Ic[i-1][j-1]+ \
                  self.subst[chars], self.Ir[i-1][j-1]+self.subst[chars]]
          thisC = [self.M[i][j-1]-self.open, self.Ic[i][j-1]-self.extend, \
                          self.Ir[i][j-1]-self.double]
          thisR = [self.M[i-1][j]-self.open, self.Ic[i-1][j]-self.double, \
                  self.Ir[i-1][j]-self.extend]
          return max(thisM), max(thisC), max(thisR)
      def score_align(self, row_seq, col_seq):
          """build dynamic programming matrices to align two sequences"""
          self.row_seq, self.col_seq = list(row_seq), list(col_seq)
          self.maxI, self.maxJ = len(self.row_seq)+1, len(self.col_seq)+1
          self.init_array() #initialize arrays
          for i in xrange(1, self.maxI):
              row_char = self.row_seq[i-1]
              for j in xrange(1, self.maxJ):
                  chars = row_char+self.col_seq[j-1]
                  self.M[i][j], self.Ic[i][j], self.Ir[i][j] = self.score_cell(i, j, chars)
      def traceback_col_seq(self):
          """trace back column sequence in matrices in a2m format"""
          i, j = self.maxI-1, self.maxJ-1
          self.traceback = list()
          #find which matrix to start in
          #THIS IS WHERE THE PROBLEM LIES I THINK
          cell = (self.M[i][j], self.Ic[i][j], self.Ir[i][j])
          prevMatrix = cell.index(max(cell))
          while i &gt; 1 and j &gt; 1:
              if prevMatrix == 0: #M matrix
                  i -= 1; j -= 1 #step up diagonally
                  prevChars = self.row_seq[i-1]+self.col_seq[j-1]
                  diag = self.score_cell(i, j, prevChars) #re-score diagonal cell
                  prevMatrix = diag.index(max(diag)) #determine which matrix that was
                  self.traceback.append(self.col_seq[j])
              elif prevMatrix == 1: #Ic matrix
                  j -= 1 
                  prevChars = self.row_seq[i-1]+self.col_seq[j-1]
                  left = self.score_cell(i, j, prevChars)
                  prevMatrix = left.index(max(left))
                  self.traceback.append(self.col_seq[j].lower())
              elif prevMatrix == 2: #Ir matrix
                  i -= 1
                  prevChars = self.row_seq[i-1]+self.col_seq[j-1]
                  up = self.score_cell(i, j, prevChars)
                  prevMatrix = up.index(max(up))
                  self.traceback.append('-')
          for j in xrange(j,0,-1): #hit top of matrix before ending, add chars
              self.traceback.append(self.col_seq[j-1])
          for i in xrange(i,0,-1): #hit left of matrix before ending, add gaps
              self.traceback.append('-')
          print(''.join(self.row[::-1]))
          return ''.join(self.traceback[::-1])
      def score_a2m(self, s1, s2):
          """scores an a2m alignment of two sequences. I believe this function correctly
          scores alignments and is used to test my alignments. The value produced by this
          function should be the same as the largest value in the bottom right of the three
          matrices"""
          s1, s2 = list(s1.strip('.')), list(s2.strip('.'))
          s1_pos, s2_pos = len(s1)-1, len(s2)-1
          score, gap = 0, False
          while s1_pos &gt;= 0 and s2_pos &gt;= 0:
              if s1[s1_pos].islower() and gap is False:
                  score -= self.open; s1_pos -= 1; gap = True
              elif s1[s1_pos].islower() and gap is True:
                  score -= self.extend; s1_pos -= 1
              elif s2[s2_pos].islower() and gap is False:
                  score -= self.open; s2_pos -= 1; gap = True
              elif s2[s2_pos].islower() and gap is True:
                  score -= self.extend; s2_pos -= 1
              elif s1[s1_pos] == '-' and gap is False:
                  score -= self.open; s1_pos -= 1; s2_pos -= 1; gap = True
              elif s1[s1_pos] == '-' and gap is True:
                  score -= self.extend; s1_pos -= 1; s2_pos -= 1
              elif s2[s2_pos] == '-' and gap is False:
                  score -= self.open; s1_pos -= 1; s2_pos -= 1; gap = True
              elif s2[s2_pos] == '-' and gap is True:
                  score -= self.extend; s1_pos -= 1; s2_pos -= 1
              elif gap is True:
                  score += self.subst[s1[s1_pos].upper() + s2[s2_pos].upper()]
                  s1_pos -= 1; s2_pos -= 1; gap = False
              else:
                  score += self.subst[s1[s1_pos].upper() + s2[s2_pos].upper()]
                  s1_pos -= 1; s2_pos -= 1
          if s1_pos &gt;= 0 and gap is True:
              score -= self.extend*s1_pos
          elif s1_pos &gt;= 0 and gap is False:
              score -= self.open+s1_pos*self.extend
          if s2_pos &gt;= 0 and gap is True:
              score -= self.extend*s2_pos
          elif s2_pos &gt;= 0 and gap is False:
              score -= self.open+s2_pos*self.extend
          return score
  
  
  test = global_aligner(blosumMatrix)
  s1,s2 = 'ZZ','AAAC'
  test.score_align(s1, s2)
  align = test.traceback_col_seq()
  print('This score: ', test.score_a2m(s1,align)
  print('Correct score: ', test.score_a2m(s1,'AAac'))
  </code></pre>
  
  <p>Blosum parser</p>
  
  <pre><code>def parse_blosum(blosumFile):
      blosumMatrix, commentFlag = dict(), False
      for line in blosumFile:
          if not line.startswith('#') and not commentFlag:
              alphabet = line.rstrip().split()
              commentFlag = True
          elif commentFlag:
              line = line.rstrip().split()
              thisChar, line = line[0], line[1:]
              for x in xrange(len(line)):
                  alphaChar, thisValue = alphabet[x], line[x]
                  blosumMatrix[thisChar+alphaChar] = int(thisValue)
      return blosumMatrix
  </code></pre>