In [1]:

    

# Using an edit-distance-like dynamic programming formulation, we can
# look for approximate occurrences of p in t.

import sys
import numpy

# Assume x is the string labeling rows of the matrix and y is the
# string labeling the columns

def trace(D, x, y):
    ''' Backtrace edit-distance matrix D for strings x and y '''
    i, j = len(x), len(y)
    xscript = []
    while i > 0:
        diag, vert, horz = sys.maxsize, sys.maxsize, sys.maxsize
        delt = None
        if i > 0 and j > 0:
            delt = 0 if x[i-1] == y[j-1] else 1
            diag = D[i-1, j-1] + delt
        if i > 0:
            vert = D[i-1, j] + 1
        if j > 0:
            horz = D[i, j-1] + 1
        if diag <= vert and diag <= horz:
            # diagonal was best
            xscript.append('R' if delt == 1 else 'M')
            i -= 1; j -= 1
        elif vert <= horz:
            # vertical was best; this is an insertion in x w/r/t y
            xscript.append('I')
            i -= 1
        else:
            # horizontal was best
            xscript.append('D')
            j -= 1
    # j = offset of the first (leftmost) character of t involved in the
    # alignment
    return j, (''.join(xscript))[::-1] # reverse and string-ize

def kEditDp(p, t):
    ''' Find and return the alignment of p to a substring of t with the
        fewest edits.  We return the edit distance, the offset of the
        substring aligned to, and the edit transcript.  If multiple
        alignments tie for best, we report the leftmost. '''
    D = numpy.zeros((len(p)+1, len(t)+1), dtype=int)
    # Note: First row gets zeros.  First column initialized as usual.
    D[1:, 0] = range(1, len(p)+1)
    for i in range(1, len(p)+1):
        for j in range(1, len(t)+1):
            delt = 1 if p[i-1] != t[j-1] else 0
            D[i, j] = min(D[i-1, j-1] + delt, D[i-1, j] + 1, D[i, j-1] + 1)
    # Find minimum edit distance in last row
    mnJ, mn = None, len(p) + len(t)
    for j in range(len(t)+1):
        if D[len(p), j] < mn:
            mnJ, mn = j, D[len(p), j]
    # Backtrace; note: stops as soon as it gets to first row
    off, xcript = trace(D, p, t[:mnJ])
    # Return edit distance, offset into T, edit transcript
    return mn, off, xcript, D


    

In [2]:

    

p, t = 'TACGTCAGC', 'AACCCTATGTCATGCCTTGGA'
mn, off, xscript, D = kEditDp(p, t)
mn, off, xscript


    

    
Out[2]:





(2, 5, 'MMRMMMMDMM')

In [3]:

    

D


    

    
Out[3]:





array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 0, 0, 1, 1, 1],
       [2, 1, 1, 2, 2, 2, 1, 0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 1, 1, 1, 2, 1],
       [3, 2, 2, 1, 2, 2, 2, 1, 1, 2, 2, 1, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2],
       [4, 3, 3, 2, 2, 3, 3, 2, 2, 1, 2, 2, 2, 3, 2, 2, 2, 3, 3, 2, 2, 3],
       [5, 4, 4, 3, 3, 3, 3, 3, 2, 2, 1, 2, 3, 2, 3, 3, 3, 2, 3, 3, 3, 3],
       [6, 5, 5, 4, 3, 3, 4, 4, 3, 3, 2, 1, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4],
       [7, 6, 5, 5, 4, 4, 4, 4, 4, 4, 3, 2, 1, 2, 3, 4, 4, 4, 4, 4, 5, 4],
       [8, 7, 6, 6, 5, 5, 5, 5, 5, 4, 4, 3, 2, 2, 2, 3, 4, 5, 5, 4, 4, 5],
       [9, 8, 7, 6, 6, 5, 6, 6, 6, 5, 5, 4, 3, 3, 3, 2, 3, 4, 5, 5, 5, 5]])

In [ ]: