import numpy as np seq1 = "ACTAC" seq2 = "AACTGATGA" m = len(seq1) n = len(seq2) penal = np.zeros((m+1, n+1)) trace = np.zeros((m+1, n+1)) for j in range(0, n+1): penal[0][j] = j for i in range(0, m+1): penal[i][0] = i for j in range(1, n+1): trace[0][j] = 1 for i in range(1, m+1): trace[i][0] = 0 for i in range(1, m+1): for j in range(1, n+1): if seq1[i-1] == seq2[j-1]: pij = 0 else: pij = 1 x = [penal[i-1][j] + 1, penal[i][j-1] + 1, penal[i-1][j-1] + pij] penal_min = np.amin(x) index_min = np.argmin(x) penal[i][j] = penal_min trace[i][j] = index_min result = [[], []] i = m j = n while i > 0 or j > 0: if trace[i][j] == 0: result[0].append(seq1[i-1]) result[1].append("-") i -= 1 elif trace[i][j] == 1: result[0].append("-") result[1].append(seq2[j-1]) j -= 1 else: result[0].append(seq1[i-1]) result[1].append(seq2[j-1]) i -= 1 j -= 1 for k in range(len(result[0])-1, -1, -1): print(result[0][k], " <-> ", result[1][k])