Traceback in Smith-Wateman algorithm with affine gap penalty

Question

I'm trying to implement the Smith-Waterman algorithm for local sequence alignment using the affine gap penalty function. I think I understand how to initiate and compute the matrices required for calculating alignment scores, but am clueless as to how to then traceback to find the alignment. To generate the 3 matrices required I have the following code

for j in range(1, len2):
    for i in range(1, len1):
        fxOpen = F[i][j-1] + gap
        xExtend = Ix[i][j-1] + extend
        Ix[i][j] = max(fxOpen, xExtend)

        fyOpen = F[i-1][j] + gap
        yExtend = Iy[i-1][j] + extend
        Iy[i][j] = max(fyOpen, yExtend)

        matchScore = (F[i-1][j-1]  + simMatrixDict[seq1[i-1]+seq2[j-1]])
        xScore = Ix[i-1][j-1] + simMatrixDict[seq1[i-1]+seq2[j-1]]
        yScore = Iy[i-1][j-1] + simMatrixDict[seq1[i-1]+seq2[j-1]]
        F[i][j] = max(0, matchScore, xScore, yScore)

I am unsure if I need a single matrix for traceback, or just 1? Any clarification on how to go about tracing back from the max score in F would be much appreciated.

Answer 1

The important thing to remember about traceback in Smith-Waterman is that the matrix a value is in determines the direction that you move. So, if you are in Fyou're moving diagonally, if you're in Ix, you're moving horizontally, and if you're in Iy, you're moving vertically. This means that all you need to store in the pointer matrix is the matrix that you arrived at a square from. The matrix that you are coming from, not the one that you are going to, determines the dirction which to go.

For example:

Say you are at F[5][5]:

• If pointer matrix says to go to Ix, go to Ix[4][4]
• If pointer matrix says to go to Iy, go to Iy[4][4]
• If pointer matrix says to go to F, go to F[4][4]

Whereas if you are at Ix[5][5]:

• If pointer matrix says to go to Ix, go to Ix[4][5]
• If pointer matrix says to go to F, go to F[4][5]

Or if you are at Iy[5][5]:

• If pointer matrix says to go to Iy, go to Iy[5][4]
• If pointer matrix says to go to F, go to F[5][4]

Assuming that the first index is the x coordinate and the second is the y coordinate.

Continue tracing back until you reach a cell with a maximum value of 0.

Building the pointer matrix: You need one pointer matrix each for F, Ix, and Iy. These matrices only need to indicate which matrix a value came from, because that tells you which direction you were moving in. So, as you're running through the dynamic programming phase of the algorithm, you should also be building up the pointer matrices. Every time you store a new maximum value in a cell in F, Ix, or Iy, you should update the corresponding matrix to indicate where it came from. If, for instance, the highest value you can have in F[5][5] comes by aligning the two next bases when you're in F[4][4], the Fpointer[5][5] should be set to F, because you got there from the F matrix.