modig_align.py - dev - Tsitic/bioinf_align

modig_align.py
133 строки3 KB
Tsitic
upload files
28 ноя 2025, 04:04
28 ноя 2025, 04:0455a13bf
#Параметры модели

delta = 0.2
tau = 0.1
eps = 0.1

a_MM = 1 - 2*delta - tau
a_MX = delta
a_MY = delta

a_XX = eps
a_YY = eps

a_XM = 1 - eps - tau
a_YM = 1 - eps - tau

def p_nucl(a, b):
    pair = a + b
    if pair in ["TT","CC","AA","GG"]:
        return 0.5
    if pair in ["CT","AG","GA","TC"]:
        return 0.05
    if pair in ["AT","GC","TA","CG"]:
        return 0.3
    if pair in ["GT","AC","TG","CA"]:
        return 0.15

q = 0.25


X = "TGTATC"
Y = "TAG"

n = len(X)
m = len(Y)


DP = [[None]*(m) for _ in range(n)]

#первая клетка
DP[0][0] = {
    "M": (1.0, None, None),
    "X": (0.0, None, None),
    "Y": (0.0, None, None),
}

#Первый столбец
for i in range(1, n):
    prev_M = DP[i-1][0]["M"][0] * a_MX
    prev_X = DP[i-1][0]["X"][0] * a_XX

    if prev_M >= prev_X:
        val =q* prev_M
        DP[i][0] = {
            "M": (0.0, None, None),
            "X": (val, (i-1,0), "M"),
            "Y": (0.0, None, None),
        }
    else:
        val = q * prev_X
        DP[i][0] = {
            "M": (0.0, None, None),
            "X": (val, (i-1,0), "X"),
            "Y": (0.0, None, None),
        }



#первая строка
for j in range(1, m):

    prev_M = DP[0][j-1]["M"][0] * a_MY
    prev_Y = DP[0][j-1]["Y"][0] * a_YY

    if prev_M >= prev_Y:
        val = q * prev_M
        DP[0][j] = {
            "M": (0.0, None, None),
            "X": (0.0, None, None),
            "Y": (val, (0,j-1), "M"),
        }
    else:
        val =q* prev_Y
        DP[0][j] = {
            "M": (0.0, None, None),
            "X": (0.0, None, None),
            "Y": (val, (0,j-1), "Y"),
        }



for i in range(1, n):
    for j in range(1, m):

        #M
        cand = [
            (DP[i-1][j-1]["M"][0] * a_MM, (i-1,j-1), "M"),
            (DP[i-1][j-1]["X"][0] * a_XM, (i-1,j-1), "X"),
            (DP[i-1][j-1]["Y"][0] * a_YM, (i-1,j-1), "Y"),
        ]
        best = max(cand, key=lambda x: x[0])
        Mval = p_nucl(X[i-1], Y[j-1]) * best[0]
        #X
        cand = [
            (DP[i-1][j]["M"][0] * a_MX, (i-1,j), "M"),
            (DP[i-1][j]["X"][0] * a_XX, (i-1,j), "X"),
        ]
        best2 = max(cand, key=lambda x: x[0])
        Xval =q * best2[0]
        #Y
        cand = [
            (DP[i][j-1]["M"][0] * a_MY, (i,j-1), "M"),
            (DP[i][j-1]["Y"][0] * a_YY, (i,j-1), "Y"),
        ]
        best3 = max(cand, key=lambda x: x[0])
        Yval = q* best3[0]
        #rezult
        DP[i][j] = {
            "M": (Mval, best[1], best[2]),
            "X": (Xval, best2[1], best2[2]),
            "Y": (Yval, best3[1], best3[2]),
        }



for i in range(n):
    for j in range(m):
        cell = DP[i][j]
        print(f"[{i},{j}]  "
              f"M={cell['M'][0]} from {cell['M'][1]}: {cell['M'][2]} | "
              f"X={cell['X'][0]} from {cell['X'][1]}: {cell['X'][2]} | "
              f"Y={cell['Y'][0]} from {cell['Y'][1]}: {cell['Y'][2]}")
    print()