BioInfa 9.2

import math
import sys


def affine_align(x, y, p1, p2, g, s):
    """
    Scoring parameters:
    p1 = match
    p2 = mismatch
    g = gap open penalty
    e = gap extension penalty
    """
    M = [[0] * (len(x) + 1) for i in range(len(y) + 1)]
    Ix = [[0] * (len(x) + 1) for i in range(len(y) + 1)]
    Iy = [[0] * (len(x) + 1) for i in range(len(y) + 1)]
    for i in range(1, len(y) + 1):
        M[i][0] = -math.inf
    for j in range(1, len(x) + 1):
        M[0][j] = -math.inf
    for i in range(0, len(y) + 1):
        Ix[i][0] = -math.inf
    for j in range(1, len(x) + 1):
        Ix[0][j] = -g if Ix[0][j - 1] == -math.inf else Ix[0][j - 1] - s
    for j in range(0, len(x) + 1):
        Iy[0][j] = -math.inf
    for i in range(1, len(y) + 1):
        Iy[i][0] = -g if Iy[i - 1][0] == -math.inf else Iy[i - 1][0] - s
    for i in range(1, len(y) + 1):
        for j in range(1, len(x) + 1):
            M[i][j] = max(M[i - 1][j - 1] + delta(x[j - 1], y[i - 1], p1, p2),
                          Ix[i - 1][j - 1] + delta(x[j - 1], y[i - 1], p1, p2),
                          Iy[i - 1][j - 1] + delta(x[j - 1], y[i - 1], p1, p2))
            Ix[i][j] = max(M[i][j - 1] - g,
                           Ix[i][j - 1] - s)
            Iy[i][j] = max(M[i - 1][j] - g,
                           Iy[i - 1][j] - s)
    x_ret = ""
    y_ret = ""
    i = len(y)
    j = len(x)
    align_scores = (M[i][j], Iy[i][j], Ix[i][j])
    matrix_idx = align_scores.index(max(align_scores))
    # matrix_key will track the current matrix through the traceback
    matrix_key = ["M", "Iy", "Ix"][matrix_idx]
    while i > 0 and j > 0:
        if matrix_key == "M":
            if M[i][j] == M[i - 1][j - 1] + p1 or M[i][j] == M[i - 1][j - 1] - p2:
                x_ret = x[j - 1] + x_ret
                y_ret = y[i - 1] + y_ret
                i -= 1
                j -= 1
                matrix_key = "M"
            elif M[i][j] == Iy[i - 1][j - 1] + p1 or M[i][j] == Iy[i - 1][j - 1] - p2:
                x_ret = x[j - 1] + x_ret
                y_ret = y[i - 1] + y_ret
                i -= 1
                j -= 1
                matrix_key = "Iy"
            elif M[i][j] == Ix[i - 1][j - 1] + p1 or M[i][j] == Ix[i - 1][j - 1] - p2:
                x_ret = x[j - 1] + x_ret
                y_ret = y[i - 1] + y_ret
                i -= 1
                j -= 1
                matrix_key = "Ix"
        elif matrix_key == "Iy":
            if Iy[i][j] == M[i - 1][j] - g:
                x_ret = "_" + x_ret
                y_ret = y[i - 1] + y_ret
                i -= 1
                matrix_key = "M"
            elif Iy[i][j] == Iy[i - 1][j] - s:
                x_ret = "_" + x_ret
                y_ret = y[i - 1] + y_ret
                i -= 1
                matrix_key = "Iy"
        elif matrix_key == "Ix":
            if Ix[i][j] == M[i][j - 1] - g:
                x_ret = x[j - 1] + x_ret
                y_ret = "_" + y_ret
                j -= 1
                matrix_key = "M"
            elif Ix[i][j] == Ix[i][j - 1] - s:
                x_ret = x[j - 1] + x_ret
                y_ret = "_" + y_ret
                j -= 1
                matrix_key = "Ix"
    if i > 0:
        x_ret = ("_" * i) + x_ret
        y_ret = y[0:i] + y_ret
    if j > 0:
        x_ret = x[0:j] + x_ret
        y_ret = ("_" * j) + y_ret
    print("x", x_ret, "\ny", y_ret)
    return (x_ret, y_ret)


def delta(char1, char2, p1, p2):
    if char1 == char2:
        return p1
    else:
        return -p2


def find_sum(match, mismatch, open_gap, ext_gap, first, second):
    new_first = list(first)
    new_second = list(second)
    l = 0  # count of matches
    f = 0  # count of mismatches
    u = 0  # count of open
    r = 0  # count of all _
    if new_first[0] == "_":
        u += 1
    if new_second[0] == "_":
        u += 1
    for i in range(0, len(new_first) - 1):
        if new_first[i] == new_second[i] and (new_first[i] != "_" or new_second[i] != "_"):
            l += 1
        elif new_first[i] != "_" and new_second[i] != "_" and new_first[i] != new_second[i]:
            f += 1
        if new_first[i] == "_" or new_second[i] == "_":
            r += 1
        if i != 0 and (new_first[i - 1] != "_" and new_first[i] == "_") or (new_second[i-1] != "_" and new_second[i] == "_"):
            u += 1
    print("Matches=", l, "Mismatches=", f, "Open=", u, "All=", r, "Function=",
          l * match - f * mismatch - u * open_gap - r * ext_gap)

if __name__ == "__main__":
    match = int(input("Enter match bonus\n"))
    mismatch = int(input("Enter mismatch fine\n"))
    open_gap = int(input("Enter open gap penalty\n"))
    ext_gap = int(input("Enter gap extension penalty\n"))
    x = input("Enter first DNA\n")
    y = input("Enter second DNA\n")
    x_ret, y_ret = affine_align(x, y, match, mismatch, open_gap, ext_gap)
    find_sum(match, mismatch, open_gap, ext_gap, x_ret, y_ret)