warming up your workspace

Sequence alignment from scratch, the algorithm behind every DNA comparison

"These two sequences are 92 percent identical." You see that number in every genomics paper and every BLAST result. But two DNA strings are rarely the same length, and mutations insert and delete letters, so before you can compare them you have to line them up, deciding where to slide one against the other and where to open gaps. That lining-up is not guesswork. It is an optimization problem with one provably best answer, and the algorithm that finds it, Needleman-Wunsch, is one of the cleanest pieces of dynamic programming there is.

The one idea

To align two sequences optimally, you do not try all the ways to line them up, there are exponentially many. Instead you build a table where cell (i, j) holds the best possible score for aligning the first i letters of one sequence against the first j letters of the other. Each cell has only three ways to be reached: match the two current letters, or insert a gap in one sequence, or in the other. You take the best of those three. Because every cell is built from cells you already solved, the whole table fills in one sweep, and the bottom-right corner holds the best score for the full alignment.

The trick that makes it work: a big alignment's best score is built out of smaller alignments' best scores. Solve the small ones once, reuse them.

Build the table

Score a match +1, a mismatch -1, and a gap -1. Fill the grid, then walk backward from the corner to recover which choices were made.

def needleman_wunsch(a, b, match=1, mismatch=-1, gap=-1):
    n, m = len(a), len(b)
    H = [[0]*(m+1) for _ in range(n+1)]
    for i in range(1, n+1): H[i][0] = i*gap          # aligning a's prefix to gaps
    for j in range(1, m+1): H[0][j] = j*gap
    for i in range(1, n+1):
        for j in range(1, m+1):
            s = match if a[i-1] == b[j-1] else mismatch
            H[i][j] = max(H[i-1][j-1] + s,            # align a[i] with b[j]
                          H[i-1][j] + gap,            # gap in b
                          H[i][j-1] + gap)            # gap in a
    return H

Three details that matter:

  • The first row and column are not zeros. Aligning the first three letters of a sequence against nothing means three gaps, so the edges count -1, -2, -3. Get this boundary wrong and every alignment drifts.
  • Each cell asks one question three ways: does the best alignment ending here come from a diagonal step (both letters used), a step down (a gap), or a step right (a gap)? The max picks the winner and, implicitly, the decision.
  • The score numbers are policy, not law. Change the match, mismatch, and gap values and you change what "best" means. Real aligners use a substitution matrix that knows some mismatches are more likely than others.

Recover the alignment

The table gives the score. To get the actual lined-up strings you trace back from the bottom-right corner, at each step asking which of the three neighbors produced this cell, and emitting a letter or a gap accordingly.

score, x, y = needleman_wunsch_align("GATTACA", "GCATGCU")
print("score:", score)
print(x); print(y)

Running the full version on the textbook pair prints:

score: 0
G-ATTACA
GCA-TGCU

Read the two lines together: the algorithm inserted a gap after the G in the first sequence and after the A in the second, so the shared letters, the G, the A, the T, the middle C, stack up. That is the single highest-scoring way to align these two strings out of all possibilities, found without enumerating any of them. A quick sanity check confirms the scoring: two identical five-letter strings score exactly 5, one point per matched letter.

Where this shows up

This is the ancestor of every alignment tool in biology. Smith-Waterman changes one rule, never let a score go below zero, and suddenly you get local alignment, the best matching region rather than the whole length, which is what you want when a short gene sits inside a long genome. BLAST, the tool that runs billions of these comparisons against genome databases, is a fast heuristic wrapped around this same core. The dynamic-programming pattern itself, a table of subproblems built once and reused, is the same one behind edit distance in your spell checker and diff in your version control.

If you want to build Smith-Waterman, substitution matrices, and the scoring that real aligners use, that is the path the bioinformatics track on IWTLP takes from here, one sequence at a time.