Needleman-Wunsch Global Alignment, in depth: recurrence, traceback, affine gaps, variants and scaling

The scoring model

An alignment of strings a and b writes them in two rows of equal length, each column holding two characters, or one character and a gap symbol. A column with two gaps is not allowed. The score is the sum over columns: s(x, y) for a column pairing x with y, and a gap penalty for gap columns. Global means every character of both strings appears, so the alignment spans both sequences completely.

The substitution score s encodes biology or domain knowledge. For DNA, a simple model uses one score for a match and another for a mismatch. For proteins, matrices such as BLOSUM62 or PAM250 give log-odds scores. Positive values mean the pair appears in related proteins more often than chance predicts, and negative values mean it appears less often. Two different amino acids with similar chemistry can score positively. In a text-diff application you might score case-insensitive matches slightly below exact matches.

The gap model is the other half. A linear penalty charges g per gap character. An affine penalty charges an opening cost plus an extension cost per additional character. That fits reality better, because one insertion event of ten bases is far more likely than ten separate insertions. Conventions differ on whether a gap of length k costs open + k times extend or open + (k - 1) times extend. BLAST and EMBOSS both document their own convention, so read the docs before copying numbers between tools. Common protein defaults are BLOSUM62 with gap existence 11 and extension 1 in BLAST, and EBLOSUM62 with gap open 10 and extend 0.5 in EMBOSS needle.

The recurrence

Let F[i][j] be the best score for aligning the first i characters of a with the first j characters of b. The last column of such an alignment is one of three things: a_i paired with b_j, a_i against a gap, or a gap against b_j. Removing that column leaves an optimal alignment of shorter prefixes, which is the optimal substructure that makes dynamic programming work. With a linear gap g (a negative number):

F[0][0] = 0
F[i][0] = i * g                 # a's first i characters against gaps
F[0][j] = j * g
F[i][j] = max( F[i-1][j-1] + s(a[i], b[j]),   # diagonal: pair the characters
               F[i-1][j]   + g,               # up: a[i] against a gap
               F[i][j-1]   + g )              # left: gap against b[j]
answer  = F[n][m]

The first row and column are what make the alignment global. They charge for leading gaps, so neither sequence can start part of the way in for free. The table has (n + 1) times (m + 1) cells, each computed in constant time, so time and space are both O(nm). Each cell depends only on its upper, left and upper-left neighbours, so every anti-diagonal can be computed in parallel. That fact drives the SIMD and GPU implementations discussed later.

To recover the alignment, not just its score, start at F[n][m] and walk backwards. At each cell, find which predecessor produced the value, emit the corresponding column, and move there. Either store a pointer per cell or recompute the three candidates during traceback, which costs no extra memory.

Worked example

Align a = GATTACA with b = GTCA, scoring match +1, mismatch -1 and gap -1. Rows are prefixes of a, columns are prefixes of b, and bold cells mark the traceback path.

Check a few cells by hand. F[1][1] pairs G with G: 0 + 1 = 1. F[2][1] is the best of the diagonal F[1][0] + s(A, G) = -1 - 1 = -2, up F[1][1] - 1 = 0 and left F[2][0] - 1 = -3, so 0. The bottom-right cell is 1, the optimal global score. The traceback, preferring diagonal, then up, then left, produces:

GATTACA
G--T-CA      score: 1 - 1 - 1 + 1 - 1 + 1 + 1 = 1

Four matches and three gaps. Other alignments also score 1; for example, G-T--CA pairs b's T with a's first T instead of the second. Both are optimal, and which one you get depends only on tie-breaking order. That matters operationally: two correct implementations can emit different alignments with identical scores, so tests should compare scores and validate alignments, not diff alignment strings.

Code: linear and affine gaps

The linear-gap version, with traceback, is a direct transcription of the recurrence.

def needleman_wunsch(a, b, match=1, mismatch=-1, gap=-1):
    n, m = len(a), len(b)
    F = [[0] * (m + 1) for _ in range(n + 1)]
    for i in range(1, n + 1):
        F[i][0] = i * gap
    for j in range(1, m + 1):
        F[0][j] = j * gap
    s = lambda x, y: match if x == y else mismatch
    for i in range(1, n + 1):
        for j in range(1, m + 1):
            F[i][j] = max(F[i-1][j-1] + s(a[i-1], b[j-1]),
                          F[i-1][j] + gap, F[i][j-1] + gap)
    top, bot, i, j = [], [], n, m                  # traceback: diagonal, up, left
    while i > 0 or j > 0:
        if i and j and F[i][j] == F[i-1][j-1] + s(a[i-1], b[j-1]):
            top.append(a[i-1]); bot.append(b[j-1]); i -= 1; j -= 1
        elif i and F[i][j] == F[i-1][j] + gap:
            top.append(a[i-1]); bot.append("-"); i -= 1
        else:
            top.append("-"); bot.append(b[j-1]); j -= 1
    return F[n][m], "".join(reversed(top)), "".join(reversed(bot))

print(needleman_wunsch("GATTACA", "GTCA"))   # (1, 'GATTACA', 'G--T-CA')

Affine gaps cannot be handled by one table, because the cost of extending a gap depends on whether the previous column was already a gap. Osamu Gotoh's 1982 formulation keeps three tables, one per possible state of the last column, and stays O(nm). The figure shows the transitions.

NEG = float("-inf")

def gotoh_score(a, b, s, gap_open, gap_extend):
    """Affine gaps: a gap of length k scores gap_open + (k - 1) * gap_extend."""
    n, m = len(a), len(b)
    M = [[NEG] * (m + 1) for _ in range(n + 1)]   # ends with a[i] paired with b[j]
    X = [[NEG] * (m + 1) for _ in range(n + 1)]   # ends with a[i] against a gap
    Y = [[NEG] * (m + 1) for _ in range(n + 1)]   # ends with a gap against b[j]
    M[0][0] = 0
    for i in range(1, n + 1):
        X[i][0] = gap_open + (i - 1) * gap_extend
    for j in range(1, m + 1):
        Y[0][j] = gap_open + (j - 1) * gap_extend
    for i in range(1, n + 1):
        for j in range(1, m + 1):
            M[i][j] = s(a[i-1], b[j-1]) + max(M[i-1][j-1], X[i-1][j-1], Y[i-1][j-1])
            X[i][j] = max(M[i-1][j] + gap_open, X[i-1][j] + gap_extend, Y[i-1][j] + gap_open)
            Y[i][j] = max(M[i][j-1] + gap_open, Y[i][j-1] + gap_extend, X[i][j-1] + gap_open)
    return max(M[n][m], X[n][m], Y[n][m])

sub = lambda x, y: 2 if x == y else -1
print(gotoh_score("AAAGGGTTT", "AAATTT", sub, -5, -1))   # 5: one gap of 3 costs -7

This version was checked against exhaustive enumeration of every alignment on 200 random pairs of short strings. Do the same for your own variant before trusting it. Traceback with affine gaps must track which of the three tables you are in, not only the cell, a common source of alignments whose recomputed score does not match the reported score. In production, prefer a maintained library. Biopython's Bio.Align.PairwiseAligner supports mode = "global", open_gap_score, extend_gap_score and a substitution_matrix loaded with substitution_matrices.load("BLOSUM62"). Parasail and EMBOSS needle cover SIMD-accelerated and command-line use.

Global, semi-global and local

Three variants differ only in boundary conditions and where you read the answer. Global is the recurrence above. Semi-global (also called glocal or overlap alignment) leaves end gaps free: initialise the first row and/or column to zero and take the maximum over the last row and/or column. Use it to place a short read inside a long reference, or to align overlapping fragments, where global alignment would wrongly penalise the unaligned overhang. Local alignment, Smith-Waterman, adds a zero option to every cell's max and takes the maximum anywhere in the table. It finds the best-matching substrings and ignores unrelated flanks.

Choosing the wrong variant is the most common modelling mistake. Global alignment of a gene against a whole chromosome produces an absurd alignment full of gap penalties. Local alignment of two full-length homologous proteins can trim divergent but real ends. Ask whether both sequences are expected to be related from end to end. If yes, use global. If one contains the other, use semi-global. If only parts are related, use local.

Scaling: memory, banding and parallelism

O(nm) memory is the first wall. Two 30,000-character sequences need 9 times 10 to the 8 cells. At 4 bytes each that is 3.6 GB, and affine gaps triple it. Three techniques push the wall back.

• Score only, two rows. The score needs just the previous row, so memory drops to O(m). This is enough for ranking or filtering.
• Linear-space alignment. Hirschberg's divide and conquer finds the optimal crossing point of the middle row by running the forward pass on the top half and the reverse pass on the bottom half, then recursing. It recovers the full alignment in O(n + m) space for roughly twice the time. Myers and Miller extended the idea to affine gaps.
• Banding. If the sequences are similar, the optimal path stays near the main diagonal. Compute only cells within k of it, for O(k times max(n, m)) work. The result is exact only if the true path stays inside the band, so check whether the path touches the band edge and widen the band if it does.

For throughput, the anti-diagonal independence enables parallelism. Farrar's striped SIMD layout, used by Parasail and others, computes many cells per instruction with narrow 8- or 16-bit integers, and GPU aligners process thousands of pairs at once, one pair or one tile per thread block. Narrow integers bring a specific bug: scores saturate or overflow on long, similar sequences. Libraries detect this and recompute at wider precision; hand-written kernels often do not.

Failure modes

Failure modes seen in real pipelines:

• Gap-convention mismatch. Copying open 10, extend 1 between tools with different definitions changes every gap's cost by one extension. Scores will not reproduce across tools.
• Unknown characters. N in DNA, X or B in proteins, lowercase soft-masked bases. Decide their scores explicitly; a matrix lookup that throws or silently returns 0 changes results.
• Comparing scores across lengths or matrices. Raw scores grow with length and depend on the matrix. Use normalised identity or statistical significance for ranking, not raw scores.
• Memory blowups. A single unexpectedly long input in a batch job can exhaust memory. Guard on n times m before allocating, and route large pairs to a linear-space or banded path.
• Assuming a unique alignment. Ties are normal. Downstream code that counts gap positions should be robust to equally optimal alternatives, or apply a documented tie rule such as leftmost gap placement.

Trade-offs and relatives

Needleman-Wunsch is part of a family. Longest common subsequence is the special case with match 1, mismatch forbidden and gap 0. Edit distance is the minimising mirror image with unit costs, and longest common substring is the contiguous cousin of local alignment. For the general design method behind all of them, see dynamic programming.

What to do next

1. Decide global, semi-global or local from whether the sequences are related end to end, before choosing parameters.
2. Pick a substitution matrix and gap scheme, and write down the tool's gap-length convention next to the numbers.
3. Implement the linear-gap version above once by hand, then verify a library's scores against it on small random inputs.
4. Test affine-gap code against brute-force enumeration and recompute the score of every returned alignment.
5. Add a size guard on n times m, and route long pairs to banded or linear-space alignment.
6. Compare implementations by score, never by alignment string, because ties are legitimate.