1
\$\begingroup\$

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)
\$\endgroup\$
  • 1
    \$\begingroup\$ I'm afraid this question does not match what this site is about. Code Review is about improving existing, working code. Code Review is not the site to ask for help in fixing or changing what your code does. Once the code does what you want, we would love to help you do the same thing in a cleaner way! Please see our help center for more information. If the code in your answer does work, why not add the code from there into your question instead of the current code? Working code can be reviewed. \$\endgroup\$ – Mast Dec 11 '16 at 12:54
  • 1
    \$\begingroup\$ You posted some code in an answer that you say is up and working, I edited this code into your question. It is however a bit unclear right now what you want from us. Have you solved all your problems? Do you want us to help you improve the quality of your code? \$\endgroup\$ – Simon Forsberg Dec 11 '16 at 16:39
  • \$\begingroup\$ @SimonForsberg thank you for pointing that out! It would of course make me happy if you can help me out for optimisation. I wasn't aware this site was only for improving working code! Sorry for the trouble. \$\endgroup\$ – Iocust Dec 11 '16 at 20:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.