# Traceback in sequence alignment with affine gap penalty (Needleman-Wunsch)

I am working on an implementation of the Needleman-Wunsch sequence alignment algorithm in python, and I've already implemented the one that uses a linear gap penalty equation for scoring, but now I'm trying to write one that uses an affine gap penalty equation.

I tried using as a reference the chapter that covers this part in the book Understanding Bioinformatics but to no avail. I turned to some easier to implement material: Link to the slides that explain the implementation

So far, I've managed to create the matrices needed for the alignment, but I have no idea how to create the traceback matrices. I've looked into -the- -two- answers on SO that covered what I was looking for, and still couldn't quite figure out how I had to create the three traceback matrices needed for the alignment.

Note that I want to return every possible optimal alignment and not just one of them, as most programs do.

Here is what I use to create the matrices M, X and Y. I made it so it keeps track of where the value came from, but from the answer -two- I found here, it seems that there should be only two possibilities for the matrices X and Y:

def create_affine_matrices(self, other_sequence, I, E):
NINF = float('-inf') # Necessary when using max() - used for impossible cases.

fromM, fromX, fromY = 1, 2, 3

my_seq_len = self.my_len()
other_seq_len = other_sequence.my_len()

M = [[(0, 0) for i in range(other_seq_len+1)] for j in range(my_seq_len+1)]
X = [[(0, 0) for i in range(other_seq_len+1)] for j in range(my_seq_len+1)]
Y = [[(0, 0) for i in range(other_seq_len+1)] for j in range(my_seq_len+1)]

for i in range(1, my_seq_len+1):
M[i][0] = (NINF, 0)
X[i][0] = (NINF, 0)
Y[i][0] = (I + i * E, fromY)

for i in range(1, other_seq_len+1):
M[0][i] = (NINF, 0)
X[0][i] = (I + i * E, fromX)
Y[0][i] = (NINF, 0)

for i in range(1, my_seq_len+1):
for j in range(1, other_seq_len+1):
sub_score = get_sub_score(self.seq[i - 1], other_sequence.seq[j - 1], self.sub_matrix)
M[i][j] = max((M[i-1][j-1][0]+sub_score, fromM), (X[i-1][j-1][0]+sub_score, fromX), (Y[i-1][j-1][0]+sub_score, fromY))

X[i][j] = max(
(I + E + M[i][j-1][0], fromM),
(E + X[i][j-1][0], fromX),
(I + E + Y[i][j-1][0], fromY))

Y[i][j] = max(
(I + E + M[i-1][j][0], fromM),
(I + E + X[i-1][j][0], fromX),
(E + Y[i-1][j][0], fromY))

return M, X, Y


Here is the code that I used for the traceback, up and working.

def traceback_affine(self, other_sequence, I, E , M , X , Y, seq1=[], seq2=[]):

NINF = -float('-inf') # Necessary when using max() - used for impossible cases.

fromM, fromX, fromY = 1, 2, 3

i = self.my_len()
j = other_sequence.my_len()

current = M

while (i > 0 or j > 0):

if current == M:
if debug:
print("Case M")
if current[i][j][1] == fromM:
if debug:
print("From M")
i -= 1
j -= 1
current = M
seq1.append(self.seq[i])
seq2.append(other_sequence.seq[j])

elif current[i][j][1] == fromX:
if debug:
print("From X")
i -= 1
j -= 1
current = X
seq1.append(self.seq[i])
seq2.append(other_sequence.seq[j])

elif current[i][j][1] == fromY:
if debug:
print("From Y")
i -= 1
j -= 1
current = Y
seq1.append(self.seq[i])
seq2.append(other_sequence.seq[j])

elif current == X:
if debug:
print("Case X")
if current[i][j][1] == fromM:
if debug:
print("From M")
j -= 1

current = M
seq1.append(self.seq[i])
seq2.append('-')

elif current[i][j][1] == fromX:
if debug:
print("From X")
j -= 1

current = X
seq1.append(self.seq[i])
seq2.append('-')

elif current[i][j][1] == fromY:
if debug:
print("From Y")
j -= 1

current = Y
seq1.append(self.seq[i])
seq2.append('-')

elif current == Y:
if debug:
print("Case Y")
if current[i][j][1] == fromM:
if debug:
print("From M")

i -= 1
current = M
seq1.append('-')
seq2.append(other_sequence.seq[j])

elif current[i][j][1] == fromX:
if debug:
print("From X")

i -= 1
current = X
seq1.append('-')
seq2.append(other_sequence.seq[j])
elif current[i][j][1] == fromY:
if debug:
print("From Y")

i -= 1
current = Y
seq1.append('-')
seq2.append(other_sequence.seq[j])

seq1.reverse()
seq2.reverse()

return ''.join(seq1), ''.join(seq2)