Traceback in dynamic programming implementation of Needleman-Wunsch algorithm

Question

I almost have my needleman-wunsch implementation working but I am confused on how to handle the traceback on a specific case.

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.

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 ZZ to AAAC should be AAac. When I trace back by hand I end up with AAAc because I only visit the Ic matrix once. Here 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 AAAc. 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.

Note that the substitution matrix being used here is BLOSUM62. The recurrence relations are

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)

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.

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 > 0 or j > 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])

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.

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 > 1 and j > 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 >= 0 and s2_pos >= 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 >= 0 and gap is True:
            score -= self.extend*s1_pos
        elif s1_pos >= 0 and gap is False:
            score -= self.open+s1_pos*self.extend
        if s2_pos >= 0 and gap is True:
            score -= self.extend*s2_pos
        elif s2_pos >= 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'))

Blosum parser

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

Answer 1

def traceback_col_seq(self):
    """
    Traces back the column sequence to determine final global alignment.
    Recalculates the score using score_cell. 
    """
    i, j = self.maxI-1, self.maxJ-1
    self.traceback = list() #stores final sequence
    matrixDict = {0:'M',1:'Ic',2:'Ir'} #mapping between numeric value and matrix
    chars = self.col_seq[j-1] + self.row_seq[i-1] #store first characters
    thisM, thisC, thisR = self.score_cell(i,j, chars) 
    cell = max(thisM), max(thisC), max(thisR) #determine where to start
    prevMatrix = matrixDict[cell.index(max(cell))] #store where to go first
    while i > 0 or j > 0:
        #loop until the top left corner of the matrix is reached
        if prevMatrix == 'M':
            self.traceback.append(self.col_seq[j-1])
            i -= 1; j -= 1
            prevMatrix = matrixDict[thisM.index(max(thisM))]
        elif prevMatrix == 'Ic':
            self.traceback.append(self.col_seq[j-1].lower())
            j -= 1
            prevMatrix = matrixDict[thisC.index(max(thisC))]
        elif prevMatrix == 'Ir':
            self.traceback.append('-')
            i -= 1
            prevMatrix = matrixDict[thisR.index(max(thisR))]
        chars = self.col_seq[j-1] + self.row_seq[i-1] #store next characters
        thisM, thisC, thisR = self.score_cell(i,j, chars) #rescore next cell
    return ''.join(self.traceback[::-1])

Answer 2

As you alluded to with your comment about arrow color, the fundamental problem is that the second horizontal arrow is getting labeled as being a movement in M (black arrow) when it is really a movement in Ic (green arrow). This seems to be happening because the prevMatrix variable indicates the best matrix at (i-1, j-1), (i-1, j) or (i, j-1). This is the matrix in the previous cell from the forward perspective. Since you are doing trace*back* at this point, the cell that you are "coming from" in the trace is really (i, j) - the one that you are currently in. That determines the direction that you are moving in and thus whether you are aligning to a gap or not. Your best option is probably to have a variable that keeps track of the current matrix from loop iteration to loop iteration, and use this to determine what character to append to your alignment string. On each iteration you can update it after determining what cell you're going to next.

I think that making these changes will clarify the problem that you are running into with your re-written traceback_col_sequence function, but preliminarily it looks like you aren't taking into account the information that you have gained by tracing back. You're reusing score_cell(), which looks like it was designed to calculate scores going forward. When you are tracing back, you already know the ending, so you don't want to choose the maximum score for the previous cell. You want to choose the maximum score that could have lead to the location in the matrix that you are currently at. For instance, if you know that the matrix you are in in your traceback is M, then you know you didn't get there by extending a gap, so you shouldn't consider those options.