bio6

import numpy


seq1 = "GCATGCC"
seq2 = "GATTACA"

sim_matrix = {
    ('A','A'): +2,  ('A','G'): -1,  ('A','C'): -1,  ('A','T'): -1,
    ('G','A'): -1,  ('G','G'): +2,  ('G','C'): -1,  ('G','T'): -1,
    ('C','A'): -1,  ('C','G'): -1,  ('C','C'): +2,  ('C','T'): -1,
    ('T','A'): -1,  ('T','G'): -1,  ('T','C'): -1,  ('T','T'): +2
}

gap_penalty = -1

rows, cols = len(seq1)+1, len(seq2)+1

M = numpy.zeros((rows, cols), int)
P = numpy.zeros((rows, cols), int) #1 = diag, 2 = left, 3 = up
# dynamic programming
"""
for i in range(rows):
    M[i][0] = i * gap_penalty
for j in range(cols):
    M[0][j] = j * gap_penalty"""
for i in range(1, rows):
    for j in range(1, cols):
        diagonal = M[i-1][j-1] + sim_matrix[(seq1[i-1], seq2[j-1])]
        up   = M[i-1][j] + gap_penalty
        left = M[i][j-1] + gap_penalty

        M[i][j] = max(diagonal, up, left, 0)

        if M[i][j] == diagonal:
            P[i][j] = 1
        elif M[i][j] == left:
            P[i][j] = 2
        elif M[i][j] == up:
            P[i][j] =3

print (M)
print (P)
"""
# traceback

i, j = rows-1, cols-1
al_seq1 = al_seq2 = symbol = ""

while i > 0 and j > 0:
    score = M[i][j]
    diagonal = M[i-1][j-1]
    up       = M[i-1][j]
    left     = M[i][j-1]
    if score == diagonal + sim_matrix[(seq1[i-1], seq2[j-1])]:
        al_seq1 = seq1[i-1] + al_seq1
        al_seq2 = seq2[j-1] + al_seq2
        if seq1[i-1] == seq2[j-1]:
            symbol  = "|" + symbol
        else:
            symbol  = " " + symbol
        i -= 1
        j -= 1
    elif score == up + gap_penalty:
        al_seq1 = seq1[i-1] + al_seq1
        al_seq2 = "-" + al_seq2
        symbol  = " " + symbol
        i -= 1
    elif score == left + gap_penalty:
        al_seq1 = "-" + al_seq1
        al_seq2 = seq2[j-1] + al_seq2
        symbol  = " " + symbol
        j -= 1

while i > 0:
    al_seq1 = seq1[i-1] + al_seq1
    al_seq2 = "-" + al_seq2
    symbol  = " " + symbol
    i -= 1

while j > 0:
    al_seq1 = "-" + al_seq1
    al_seq2 = seq2[j-1] + al_seq2
    symbol  = " " + symbol
    j -= 1

print (al_seq1)
print (symbol)
print (al_seq2)"""