#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Jan 15 15:02:16 2019

@author: ariane.delrocq
"""

import Bio.SubsMat.MatrixInfo

mat = Bio.SubsMat.MatrixInfo.blosum62


def score(a, b, gap=-8, id=1, mut=-1):
    if not a or not b:
        return gap
    else:
        a = a.upper()
        b = b.upper()
        try:
            return mat[a,b]
        except KeyError:
            return mat[b,a]
    
#    elif a==b:
#        return id
#    else:
#        return mut

def align(x,y):
    n=len(x)
    m=len(y)
    mat = [[0 for _ in range(m+1)] for _ in range(n+1)]   # scores matrix
    path = [[0 for _ in range(m+1)] for _ in range(n+1)]   # last step for optimal score
    
    # first line and col
    for i in range(1, n+1):
        mat[i][0] = mat[i-1][0] + score('', x[i-1])
        path[i][0] = (i-1, 0)
    for j in range(1, m+1):
        mat[0][j] = mat[j-1][0] + score('', y[j-1])
        path[0][j] = (0, j-1)
    
    # fill table
    for i in range(1, n+1):
        for j in range(1, m+1):
            a = x[i-1]
            b = y[j-1]
            s1 = mat[i-1][j-1] + score(a,b)
            s2 = mat[i-1][j] + score('', a)
            s3 = mat[i][j-1] + score('', b)
            if s1 >= s2:
                if s1 >= s3:
                    mat[i][j] = s1
                    path[i][j] = (i-1, j-1)
                else:
                    mat[i][j] = s3
                    path[i][j] = (i, j-1)
            else:   # s2 > s1
                if s2 >= s3:
                    mat[i][j] = s2
                    path[i][j] = (i-1, j)
                else:
                    mat[i][j] = s3
                    path[i][j] = (i, j-1)
    
    bestScore = mat[n][m]
    alx = ""
    aly = ""
    i, j = n, m
    while i > 0 and j > 0:
        i2, j2 = path[i][j]
        if i==i2:
            alx+='-'
        else:
            alx+=x[i-1]
        if j==j2:
            aly+='-'
        else:
            aly+=y[j-1]
        i,j = i2, j2
    while i>0:
        aly+='-'
        alx+=x[i-1]
        i -=1
    while j>0:
        alx+='-'
        aly+=y[j-1]
        j -=1
    alx = alx[::-1]
    aly = aly[::-1]
    return bestScore, alx, aly

print(align("chat", "cat"))
print(align("chat", "cgat"))
print(align("chat", "at"))