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]]