Code: Semi-global pairwise alignment

def gap_penalty():
    return -2.0

def match_score(letterA,letterB):
    if letterA == '-' or letterB == '-':
        return gap_penalty()
    elif letterA == letterB:
        return 3.0
    else:
        return -1.0

# Initiating dynamic programming matrices, S and trace
def initiate_semiglobal_dp(m,n):
    S = np.zeros((m,n))
    trace = np.zeros((m,n,2))
    S[0,0] = 0.
    trace[0,0,:] = (0.,0.)
    for i in range(1,m):
        S[i,0] = 0
        trace[i,0,:] =(-1,0)
    for j in range(1,n):
        S[0,j] = 0
        trace[0,j,:] =(0,-1)
    return S,trace

def format_semiglobal_alignment(seqA,seqB,S,trace):
    outA,outB = "",""
    m,n = S.shape[0]-1, S.shape[1]-1
    i,j = S[:,n].argmax(), S[m,:].argmax()
    if S[i,n] > S[m,j]:
        print("Best score: " + str(S[i,n]))
        j = n
        # point the alignmnént from (m,n) to here
        for ix in range(i+1,m+1):
            print(ix,n)
            trace[ix,n,:] = (-1,0)
    else:
        print("Best score: " + str(S[m,j]))
        i = m    
        # point the alignmnént from (m,n) to here
        for ix in range(j+1,n+1):
            print(m,ix)
            trace[m,ix,:] = (0,-1)
    i,j = m,n
    while min(trace[i,j])<0:
        di,dj = trace[i,j]
        i += int(di)
        j += int(dj)
        if di == 0:
            outA = "-" + outA
        else:
            outA = seqA[i] + outA
        if dj == 0:
            outB = "-" + outB
        else:
            outB = seqB[j] + outB
    return outA,outB

def semiglobal_align(seqA,seqB,print_dynamic_matrix = False):
    # Initiating variables
    m, n = len(seqA)+1, len(seqB)+1
    S,trace = initiate_semiglobal_dp(m,n)
    # Fill in the rest of the dynamic programming matrix
    for i in range(1,m):
        for j in range(1,n):
            match = S[i-1][j-1] + match_score(seqA[i-1],seqB[j-1])
            delete = S[i-1,j] + match_score(seqA[i-1],'-') 
            insert = S[i,j-1] + match_score('-',seqB[j-1]) 
            S[i,j] = max(match, delete, insert)
            if match >= max(delete,insert):
                trace[i,j,:] = (-1,-1.)
            elif delete >= insert:
                trace[i,j,:] = (-1.,0.)
            else:
                trace[i,j,:] = (0.,-1.)
    if print_dynamic_matrix:
        print_dynamic(seqA,seqB,S)
    return format_semiglobal_alignment(seqA,seqB,S,trace)

seqA,seqB = semiglobal_align("CTATCTATCTCGCTATCCAA","CTACGCTATTTCA",True)
print_alignment(seqA,seqB)