In [1]:

    

from numpy import zeros

def exampleCost(xc, yc):
    """ Cost function assigning 0 to match, 2 to transition, 4 to
        transversion, and 8 to a gap """
    if xc == yc: return 0 # match
    if xc == '-' or yc == '-': return 8 # gap
    minc, maxc = min(xc, yc), max(xc, yc)
    if minc == 'A' and maxc == 'G': return 2 # transition
    elif minc == 'C' and maxc == 'T': return 2 # transition
    return 4 # transversion

def globalAlignment(x, y, s):
    """ Calculate global alignment value of sequences x and y using
        dynamic programming.  Return global alignment value. """
    D = zeros((len(x)+1, len(y)+1), dtype=int)
    for j in range(1, len(y)+1):
        D[0, j] = D[0, j-1] + s('-', y[j-1])
    for i in range(1, len(x)+1):
        D[i, 0] = D[i-1, 0] + s(x[i-1], '-')
    for i in range(1, len(x)+1):
        for j in range(1, len(y)+1):
            D[i, j] = min(D[i-1, j-1] + s(x[i-1], y[j-1]), # diagonal
                          D[i-1, j  ] + s(x[i-1], '-'),    # vertical
                          D[i  , j-1] + s('-',    y[j-1])) # horizontal
    return D, D[len(x), len(y)]


    

In [2]:

    

D, globalAlignmentValue = globalAlignment('TACGTCAGC', 'TATGTCATGC', exampleCost)
globalAlignmentValue


    

    
Out[2]:





10

In [3]:

    

print(D)


    

    




[[ 0  8 16 24 32 40 48 56 64 72 80]
 [ 8  0  8 16 24 32 40 48 56 64 72]
 [16  8  0  8 16 24 32 40 48 56 64]
 [24 16  8  2 10 18 24 32 40 48 56]
 [32 24 16 10  2 10 18 26 34 40 48]
 [40 32 24 16 10  2 10 18 26 34 42]
 [48 40 32 24 18 10  2 10 18 26 34]
 [56 48 40 32 26 18 10  2 10 18 26]
 [64 56 48 40 32 26 18 10  6 10 18]
 [72 64 56 48 40 34 26 18 12 10 10]]