# BioInfa 9.2

Mar 1st, 2021
1,354
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. import math
2. import sys
3.
4.
5. def affine_align(x, y, p1, p2, g, s):
6.     """
7.    Scoring parameters:
8.    p1 = match
9.    p2 = mismatch
10.    g = gap open penalty
11.    e = gap extension penalty
12.    """
13.     M = [[0] * (len(x) + 1) for i in range(len(y) + 1)]
14.     Ix = [[0] * (len(x) + 1) for i in range(len(y) + 1)]
15.     Iy = [[0] * (len(x) + 1) for i in range(len(y) + 1)]
16.     for i in range(1, len(y) + 1):
17.         M[i][0] = -math.inf
18.     for j in range(1, len(x) + 1):
19.         M[0][j] = -math.inf
20.     for i in range(0, len(y) + 1):
21.         Ix[i][0] = -math.inf
22.     for j in range(1, len(x) + 1):
23.         Ix[0][j] = -g if Ix[0][j - 1] == -math.inf else Ix[0][j - 1] - s
24.     for j in range(0, len(x) + 1):
25.         Iy[0][j] = -math.inf
26.     for i in range(1, len(y) + 1):
27.         Iy[i][0] = -g if Iy[i - 1][0] == -math.inf else Iy[i - 1][0] - s
28.     for i in range(1, len(y) + 1):
29.         for j in range(1, len(x) + 1):
30.             M[i][j] = max(M[i - 1][j - 1] + delta(x[j - 1], y[i - 1], p1, p2),
31.                           Ix[i - 1][j - 1] + delta(x[j - 1], y[i - 1], p1, p2),
32.                           Iy[i - 1][j - 1] + delta(x[j - 1], y[i - 1], p1, p2))
33.             Ix[i][j] = max(M[i][j - 1] - g,
34.                            Ix[i][j - 1] - s)
35.             Iy[i][j] = max(M[i - 1][j] - g,
36.                            Iy[i - 1][j] - s)
37.     x_ret = ""
38.     y_ret = ""
39.     i = len(y)
40.     j = len(x)
41.     align_scores = (M[i][j], Iy[i][j], Ix[i][j])
42.     matrix_idx = align_scores.index(max(align_scores))
43.     # matrix_key will track the current matrix through the traceback
44.     matrix_key = ["M", "Iy", "Ix"][matrix_idx]
45.     while i > 0 and j > 0:
46.         if matrix_key == "M":
47.             if M[i][j] == M[i - 1][j - 1] + p1 or M[i][j] == M[i - 1][j - 1] - p2:
48.                 x_ret = x[j - 1] + x_ret
49.                 y_ret = y[i - 1] + y_ret
50.                 i -= 1
51.                 j -= 1
52.                 matrix_key = "M"
53.             elif M[i][j] == Iy[i - 1][j - 1] + p1 or M[i][j] == Iy[i - 1][j - 1] - p2:
54.                 x_ret = x[j - 1] + x_ret
55.                 y_ret = y[i - 1] + y_ret
56.                 i -= 1
57.                 j -= 1
58.                 matrix_key = "Iy"
59.             elif M[i][j] == Ix[i - 1][j - 1] + p1 or M[i][j] == Ix[i - 1][j - 1] - p2:
60.                 x_ret = x[j - 1] + x_ret
61.                 y_ret = y[i - 1] + y_ret
62.                 i -= 1
63.                 j -= 1
64.                 matrix_key = "Ix"
65.         elif matrix_key == "Iy":
66.             if Iy[i][j] == M[i - 1][j] - g:
67.                 x_ret = "_" + x_ret
68.                 y_ret = y[i - 1] + y_ret
69.                 i -= 1
70.                 matrix_key = "M"
71.             elif Iy[i][j] == Iy[i - 1][j] - s:
72.                 x_ret = "_" + x_ret
73.                 y_ret = y[i - 1] + y_ret
74.                 i -= 1
75.                 matrix_key = "Iy"
76.         elif matrix_key == "Ix":
77.             if Ix[i][j] == M[i][j - 1] - g:
78.                 x_ret = x[j - 1] + x_ret
79.                 y_ret = "_" + y_ret
80.                 j -= 1
81.                 matrix_key = "M"
82.             elif Ix[i][j] == Ix[i][j - 1] - s:
83.                 x_ret = x[j - 1] + x_ret
84.                 y_ret = "_" + y_ret
85.                 j -= 1
86.                 matrix_key = "Ix"
87.     if i > 0:
88.         x_ret = ("_" * i) + x_ret
89.         y_ret = y[0:i] + y_ret
90.     if j > 0:
91.         x_ret = x[0:j] + x_ret
92.         y_ret = ("_" * j) + y_ret
93.     print("x", x_ret, "\ny", y_ret)
94.     return (x_ret, y_ret)
95.
96.
97. def delta(char1, char2, p1, p2):
98.     if char1 == char2:
99.         return p1
100.     else:
101.         return -p2
102.
103.
104. def find_sum(match, mismatch, open_gap, ext_gap, first, second):
105.     new_first = list(first)
106.     new_second = list(second)
107.     l = 0  # count of matches
108.     f = 0  # count of mismatches
109.     u = 0  # count of open
110.     r = 0  # count of all _
111.     if new_first[0] == "_":
112.         u += 1
113.     if new_second[0] == "_":
114.         u += 1
115.     for i in range(0, len(new_first) - 1):
116.         if new_first[i] == new_second[i] and (new_first[i] != "_" or new_second[i] != "_"):
117.             l += 1
118.         elif new_first[i] != "_" and new_second[i] != "_" and new_first[i] != new_second[i]:
119.             f += 1
120.         if new_first[i] == "_" or new_second[i] == "_":
121.             r += 1
122.         if i != 0 and (new_first[i - 1] != "_" and new_first[i] == "_") or (new_second[i-1] != "_" and new_second[i] == "_"):
123.             u += 1
124.     print("Matches=", l, "Mismatches=", f, "Open=", u, "All=", r, "Function=",
125.           l * match - f * mismatch - u * open_gap - r * ext_gap)
126.
127. if __name__ == "__main__":
128.     match = int(input("Enter match bonus\n"))
129.     mismatch = int(input("Enter mismatch fine\n"))
130.     open_gap = int(input("Enter open gap penalty\n"))
131.     ext_gap = int(input("Enter gap extension penalty\n"))
132.     x = input("Enter first DNA\n")
133.     y = input("Enter second DNA\n")
134.     x_ret, y_ret = affine_align(x, y, match, mismatch, open_gap, ext_gap)
135.     find_sum(match, mismatch, open_gap, ext_gap, x_ret, y_ret)
136.
RAW Paste Data