smith-waterman-affine-gap-gotoh

Answer

The fix is a correct stateful traceback for the affine-gap (Gotoh) local alignment. The DP itself (recurrences below) was fine; the subtle bug was in traceback:

The bug. A naive traceback that compares H[i][j] == E[i][j] and then re-evaluates H vs E vs F at (i, j-1) is wrong. When H[i][j] == E[i][j] was produced by extending a gap (E[i][j] = E[i][j-1] − gap_extend), the previous cell must be followed through the E-chain itself, not through the plain H argmax — H[i][j-1] may have come from a diagonal, and jumping to it silently re-roots the alignment at a different source, changing its score (verified: it returned 6 for an alignment worth 1).

The fix. Track which matrix's value you are tracing (state ∈ {H, E, F}). In E/F state, decide how the gap value was produced: E[i][j] == E[i][j-1] − gap_extend means the gap run continues; otherwise it was opened (E[i][j] == H[i][j-1] − gap_open − gap_extend) and you return to the plain H chain.

"""Gotoh local alignment with affine gaps. O(n*m) time, O(n*m) memory.

Scoring: a gap of length L costs gap_open + L * gap_extend.
Recurrences:
  E[i][j] = max(E[i][j-1] - gap_extend, H[i][j-1] - gap_open - gap_extend)  # gap in seq1
  F[i][j] = max(F[i-1][j] - gap_extend, H[i-1][j] - gap_open - gap_extend)  # gap in seq2
  H[i][j] = max(0, H[i-1][j-1] + (match|mismatch), E[i][j], F[i][j])
"""
from array import array

_NEG = -(10 ** 12)

def gotoh_local(seq1, seq2, match, mismatch, gap_open, gap_extend):
    """Return (score, aligned_seq1, aligned_seq2)."""
    n, m = len(seq1), len(seq2)
    H = [array('q', [0]) * (m + 1) for _ in range(n + 1)]        # best local score
    E = [array('q', [_NEG]) * (m + 1) for _ in range(n + 1)]     # ends with gap in seq1
    F = [array('q', [_NEG]) * (m + 1) for _ in range(n + 1)]     # ends with gap in seq2

    best = 0
    bi = bj = 0
    go = gap_open + gap_extend

    for i in range(1, n + 1):
        Hprev, Eprev, Fprev = H[i-1], E[i-1], F[i-1]
        Hcur, Ecur, Fcur = H[i], E[i], F[i]
        s1 = seq1[i-1]
        for j in range(1, m + 1):
            e = Ecur[j-1] - gap_extend
            o = Hcur[j-1] - go
            Ecur[j] = e if e > o else o
            f = Fprev[j] - gap_extend
            o = Hprev[j] - go
            Fcur[j] = f if f > o else o
            v = Hprev[j-1] + (match if s1 == seq2[j-1] else mismatch)
            if Ecur[j] > v: v = Ecur[j]
            if Fcur[j] > v: v = Fcur[j]
            if v < 0: v = 0
            Hcur[j] = v
            if v > best:
                best = v; bi, bj = i, j

    if best == 0:
        return (0, "", "")

    # stateful traceback: follow the matrix that actually produced each value
    aln1, aln2 = [], []
    i, j, state = bi, bj, 'H'
    while i > 0 and j > 0 and H[i][j] > 0:
        if state == 'E':                                   # gap in seq1
            aln1.append('-'); aln2.append(seq2[j-1]); j -= 1
            if not (E[i][j] == E[i][j-1] - gap_extend):    # opened here?
                state = 'H'
        elif state == 'F':                                 # gap in seq2
            aln1.append(seq1[i-1]); aln2.append('-'); i -= 1
            if not (F[i][j] == F[i-1][j] - gap_extend):
                state = 'H'
        elif H[i][j] == E[i][j]:
            state = 'E'
        elif H[i][j] == F[i][j]:
            state = 'F'
        else:
            aln1.append(seq1[i-1]); aln2.append(seq2[j-1])
            i -= 1; j -= 1

    return (best, ''.join(reversed(aln1)), ''.join(reversed(aln2)))

Notes: - Full H, E, F matrices are kept for traceback; using array('q') keeps 2000×2000 at ~96 MB instead of ~450 MB with Python lists. - The DP is strictly O(n·m); 2000×2000 runs in ~1.2–1.9 s. - Empty/no-positive-score inputs return (0, "", "").