hinagawa

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

Adblocker detected! Please consider disabling it...

We've detected AdBlock Plus or some other adblocking software preventing Pastebin.com from fully loading.

We don't have any obnoxious sound, or popup ads, we actively block these annoying types of ads!

Please add Pastebin.com to your ad blocker whitelist or disable your adblocking software.

×