The core alignment algorithms and their tradeoffs

In BIO-3.1 you learned what an alignment means: it is a diff between two biological sequences, where each column tells a small evolutionary story. Two equal residues in a column is a match, two different residues is a mismatch (a substitution, one of the mutation types from S7.1), and a residue paired with a dash is a gap (an insertion or deletion). BIO-3.1 gave you the vocabulary. This lesson gives you the machine that actually computes the best diff, and, just as important, teaches you that best is a word with strings attached.

Why you cannot just try every alignment

Start with the obvious idea and watch it explode. To align two sequences you decide, at every position, whether to line residues up or slip in a gap. For two sequences of length n, the count of distinct gapped alignments grows faster than 2 to the n. Two sequences of length 30, tiny by biology's standards, already admit more alignments than there are seconds in the age of the universe. You cannot enumerate them. You need an algorithm that finds the best one without looking at all of them, and that is exactly what dynamic programming was built for.

Dynamic programming: solve one column, reuse it forever

The trick that tames the exponential blowup is a single observation about structure. Ask for the best alignment of the whole of sequence a against the whole of sequence b. Now look only at the very last column of that best alignment. It can be exactly one of three things:

• a[last] aligned against b[last] (a diagonal step, a match or mismatch),
• a[last] aligned against a gap (a[last] was inserted, or b's residue was deleted),
• a gap aligned against b[last] (the mirror case).

Whichever of the three it is, everything before that last column must itself be the best alignment of the two shorter prefixes. If it were not, you could swap in a better prefix alignment and beat the supposed best, a contradiction. That property, the best answer is built from best answers to smaller subproblems, is called optimal substructure, and it is what lets you fill a table once and never recompute.

Lay the two sequences on the edges of a grid. Let matrix[i][j] hold the score of the best alignment of the first i residues of a against the first j residues of b. Row i indexes a, column j indexes b, and matrix[0][0] is the empty-against-empty origin. The three cases above become three neighbors you already computed, and the recurrence writes itself:

matrix[i][j] = max of three candidates. The diagonal candidate is matrix[i-1][j-1] plus the substitution score for aligning a[i-1] with b[j-1]. The up candidate is matrix[i-1][j] plus the gap penalty (consume a residue, leave a gap in b). The left candidate is matrix[i][j-1] plus the gap penalty (consume a b residue, leave a gap in a). Seed the borders first: matrix[i][0] is i times the gap penalty and matrix[0][j] is j times the gap penalty, because a prefix aligned to nothing is all gaps.

That is the entire idea. Each cell costs three additions and a max, and there are (m+1) times (n+1) cells, so the whole table fills in O(m times n) time. The code below is a faithful Python mirror of the real fill in src/lib/bio/alignment.ts, which the interactive you are about to use runs under the hood.

def needleman_wunsch(a, b, match=1, mismatch=-1, gap=-2):
    m, n = len(a), len(b)
    # matrix has (m+1) rows and (n+1) cols. Row i indexes a, col j indexes b.
    matrix = [[0] * (n + 1) for _ in range(m + 1)]

    # Seed the borders. A prefix of length k aligned to nothing costs k * gap.
    for i in range(1, m + 1):
        matrix[i][0] = i * gap
    for j in range(1, n + 1):
        matrix[0][j] = j * gap

    # Fill every cell from the three neighbors already computed.
    for i in range(1, m + 1):
        for j in range(1, n + 1):
            s = match if a[i - 1] == b[j - 1] else mismatch
            diag = matrix[i - 1][j - 1] + s   # align a[i-1] with b[j-1]
            up = matrix[i - 1][j] + gap        # a gap in b (consume a[i-1])
            left = matrix[i][j - 1] + gap      # a gap in a (consume b[j-1])
            matrix[i][j] = max(diag, up, left)

    # The global score sits in the bottom-right corner.
    return matrix[m][n], matrix

The score is now in the bottom-right corner, but a score is not an alignment. To recover the actual lined-up sequences you run traceback: start at the corner and, at each cell, ask which of the three candidates produced its value, step to that neighbor, and record a match, a mismatch, or a gap accordingly. Walk back to the origin and reverse what you collected. The path you traced is the optimal alignment, read out column by column.

Global versus local: two algorithms, one recurrence

Everything above describes Needleman-Wunsch, a global alignment. Global means end to end: it forces the entire a to line up against the entire b, spending gaps wherever it must to reach the far corner. That is the right tool when your two sequences are believed to be homologous along their whole length and are roughly the same size, for example the same gene sequenced in two species.

Now change one question. What if a shared motif is buried inside two otherwise unrelated sequences, or you are hunting a short fragment inside a long one? Forcing a global end-to-end alignment would drown the real signal under a sea of gaps at the ends. You want the best matching subregion, not the best full-length march. That is local alignment, and Smith-Waterman gets it with two small edits to the same recurrence:

• Add a fourth candidate, zero, to the max. Now matrix[i][j] equals max of (zero, diagonal, up, left). A stretch that has gone badly negative resets to zero instead of dragging the score down, so a fresh high-scoring region can start anywhere.
• Read the answer from the largest cell anywhere in the grid, not the corner, and trace back from there until you hit a cell whose value is zero. That trims the alignment down to just the good part.

Same table, same O(m times n) cost, a different reading of it. The interactive below runs both from the same code path.

Drive the grid yourself

Here is the flagship. Type in two sequences, toggle between global (Needleman-Wunsch) and local (Smith-Waterman), and set the three scores. Watch the full DP grid fill in and the highlighted traceback thread through it to the reconstructed alignment. Do this deliberately: start with the defaults, then drop the gap penalty from -2 to -1 and watch the optimal path shift toward using more gaps. Then push it to -4 and watch gaps become too expensive to afford. You are not reading about parameter sensitivity, you are causing it.

matrix[i][j] = max(
  matrix[i-1][j-1] + s(a[i-1], b[j-1]),   // diagonal: align the two residues
  matrix[i-1][j]   + gap,                 // up: a gap in sequence b
  matrix[i][j-1]   + gap                  // left: a gap in sequence a
)
s(x, y) = match when x == y, else mismatch

// A thin sketch of needlemanWunsch (global). Local alignment is the same
// loop with a 0 floor and a traceback that starts at the largest cell.
function fill(a: string, b: string, match: number, mismatch: number, gap: number) {
  const m = a.length;
  const n = b.length;
  const matrix: number[][] = [];
  for (let i = 0; i <= m; i++) matrix.push(new Array(n + 1).fill(0));

  // Seed the borders: a run of k gaps costs k * gap. (Local leaves these at 0.)
  for (let i = 1; i <= m; i++) matrix[i][0] = i * gap;
  for (let j = 1; j <= n; j++) matrix[0][j] = j * gap;

  // Fill. Row i indexes a, column j indexes b.
  for (let i = 1; i <= m; i++) {
    for (let j = 1; j <= n; j++) {
      const s = a[i - 1] === b[j - 1] ? match : mismatch;
      const diag = matrix[i - 1][j - 1] + s;   // align a[i-1] with b[j-1]
      const up = matrix[i - 1][j] + gap;        // a gap in b
      const left = matrix[i][j - 1] + gap;      // a gap in a
      matrix[i][j] = Math.max(diag, up, left);
      // Smith-Waterman (local) instead: Math.max(0, diag, up, left)
    }
  }

  // The global score is the bottom-right corner.
  return matrix[m][n];
}

Scoring is where the biology actually lives

The algorithm is fixed. The scores are the science, and everything interesting hides in them.

The simplest scheme, and the one the widget defaults to, is match, mismatch, and gap: a flat reward for equal residues, a flat penalty for different ones, and a flat penalty per gap column. For DNA that crude scheme is often fine. For proteins it throws away real information, because not all substitutions are equally likely or equally damaging. Swapping leucine for isoleucine (two similar residues from S4.2, both greasy and nearly identical in shape, since they are branched-chain isomers) happens constantly in evolution and barely perturbs a protein. Swapping either one for charged aspartate is rare and disruptive. A flat mismatch penalty cannot tell those apart.

Substitution matrices fix this by pricing every one of the possible residue-against-residue pairings from real data. The two families you will meet are PAM and BLOSUM. Both are log-odds matrices: each entry is roughly the logarithm of (how often this substitution is actually observed in aligned homologous proteins) divided by (how often it would occur by pure chance). A positive score means the substitution shows up more than chance expects, a signal of a tolerated, conservative swap. A negative score means it shows up less than chance, a signal of a disruptive one. The scores are empirical, read off from thousands of already-aligned proteins, not derived from any law.

• PAM (Point Accepted Mutation, from Margaret Dayhoff) is built from a model of accepted mutations in closely related proteins, then extrapolated to longer evolutionary distances by matrix multiplication. Higher PAM number means more divergence: PAM250 is for distant relatives, PAM30 for close ones.
• BLOSUM (BLOcks Substitution Matrix, from the Henikoffs) is built directly from blocks of already-conserved protein regions, with no extrapolation. Higher BLOSUM number means less divergence: BLOSUM80 for similar sequences, BLOSUM45 for distant ones. BLOSUM62 is the workhorse default in protein BLAST.

There is one more modeling knob, and the frozen library deliberately keeps it simple so you can see the seam. Gaps in real biology come from insertion or deletion events, and a single event can insert or delete a whole run of residues at once (recall that a frameshift from S7.1 is one indel, not many). So a five-residue gap should be cheaper than five separate one-residue gaps, because it is more likely one event than five. Production aligners model this with an affine gap penalty: a larger gap-open cost charged once, plus a smaller gap-extend cost per additional residue. The alignment library behind this widget uses a linear gap penalty instead, one flat charge per gap column, which slightly over-penalizes long indels. That is an honest simplification worth naming: linear gaps are clearer for learning the recurrence, affine gaps are what you want for real work.

Why exact dynamic programming hits a wall

The recurrence is beautiful and it does not scale. O(m times n) is wonderful for two genes: a 1000 by 1000 grid is a million cells, done in a blink. Now aim it at modern data. A single sequencing run produces hundreds of millions of short reads, each maybe 100 to 150 bases (BIO-4 covers how they are made), and you want to place each read against the roughly 3.2 billion bases of a human genome. One read against the whole genome by full DP is 100 times 3.2 billion, about 320 billion cells. Multiply by hundreds of millions of reads and you are past 10 to the 19th cell computations for a single experiment, and the table alone would need more memory than exists on any machine.

Exact O(m times n) dynamic programming is the right answer for pairs of genes and the wrong answer for the genome-scale problem. It is not that the algorithm is slow per cell, it is that the number of cells grows as the product of two enormous numbers. This is the gap that motivates the next lesson.

Where this leaves you

Alignment is dynamic programming with a biologically motivated cost function. One recurrence, matrix[i][j] as the max of a diagonal step (score the substitution), an up step, and a left step (charge a gap), fills an O(m times n) grid, and a traceback reads the alignment back out. Needleman-Wunsch reads it end to end for whole-length homologs, Smith-Waterman adds a zero floor to find the best shared subregion. The algorithm is the easy part. The scores, match and mismatch, the gap model, the empirical PAM or BLOSUM matrix, are modeling choices, so optimal always means optimal-for-these-parameters and never optimal-as-truth. And because the grid grows as the product of the two lengths, exact DP simply cannot be pointed at a genome. That wall is what BIO-3.3 is built to climb.