Omicswell

BioAlgo

Needleman–Wunsch, step by step

Global alignment fills a table one cell at a time, then reads the best alignment back from the corner. Fill one yourself below.

Two DNA sequences that share an ancestor rarely line up letter for letter. A base may have changed, or a short stretch may have been inserted in one sequence and lost from the other. Global alignment asks for the best way to line the two sequences up from the first letter to the last, with gaps allowed. In 1970 Saul Needleman and Christian Wunsch described a method that finds that best alignment without trying every possibility. Its central idea, dynamic programming, still runs inside many alignment tools.

Scoring an alignment

Write one sequence above the other, column by column. A column with two equal letters is a match, a column with two different letters is a mismatch, and a column with a letter against a dash is a gap. Each column earns a score, and the score of the alignment is the sum. This note uses match +1, mismatch −1 and gap −2, so a gap costs more than a mismatch.

xGATTCAyGA-TCAscore+1+1−2+1+1+1= 3
match +1mismatch −1gap −2
x = GATTCA over y = GATCA with one gap. Five matches and one gap give a score of 3.

Here x is GATTCA and y is GATCA. The second sequence is one letter shorter, so every end-to-end alignment needs at least one gap. The alignment above has five matches and one gap, for 1 + 1 − 2 + 1 + 1 + 1 = 3. Is there a better one? Listing every alignment to find out is hopeless, because their number grows exponentially with the length of the sequences.

From alignments to a table

Instead of whole alignments, work with prefixes. Let F(i, j) be the best score for aligning the first i letters of x with the first j letters of y. The table has |x| + 1 rows and |y| + 1 columns, with x down the side and y across the top. The answer is the corner cell, F(|x|, |y|).

Row 0 and column 0 align a prefix against nothing, so every letter there faces a gap. That gives F(i, 0) = i · g and F(0, j) = j · g, and with g = −2 the borders read 0, −2, −4, −6 and so on.

Three ways into a cell

Look at the last column of the best alignment of the first i letters of x and the first j letters of y. It can only be one of three things, and each one leads back to a smaller cell that is already filled.

  • xiyjDiagonal, from (i−1, j−1). Pair xi with yj and add s(xi, yj), the match or mismatch score.
  • xi–Down, from (i−1, j). Put xi against a gap and add g.
  • –yjRight, from (i, j−1). Put yj against a gap and add g.
F(i, j) = maxF(i−1, j−1) + s(xi, yj)F(i−1, j) + gF(i, j−1) + g

The cell keeps the largest of the three candidates. Each candidate builds on the best score of shorter prefixes, so the largest one is the best score for this pair of prefixes as well. Filling the table row by row, from left to right, guarantees that the three neighbours of every cell are ready when the cell is reached.

Fill the table yourself

Step through the example below. Each step fills one cell and shows its three candidates, and the small arrows mark every move that reaches the cell's value. Try your own short sequences or other scores too. A larger gap penalty, for example, pushes the alignment towards mismatches instead of gaps.

0
G1
A2
T3
C4
A5
0
0
−2
−4
−6
−8
−10
G1
−2
A2
−4
T3
−6
T4
−8
C5
−10
A6
−12
Borders

Row 0 and column 0 align a prefix against nothing, so every step along a border adds g = −2. Press Next to fill cell (1, 1), or Play to watch the whole table.

  • cell being filled
  • diagonal source
  • source above
  • source on the left
  • optimal path, after traceback
Arrows in a cell mark every move that reaches its value: ↖ diagonal, ↑ from above, ← from the left. Letters other than A, C, G and T are ignored, and each sequence takes up to 9 letters.

Reading the alignment back

The corner holds the score, 3, but not the alignment itself. To recover it, start at the corner and walk back to (0, 0), each time stepping to a neighbour whose move produced the value of the current cell. A diagonal step writes a pair of letters, a step up writes a letter of x against a gap, and a step left writes a letter of y against a gap. The letters come out from the end towards the start, so the two rows are reversed at the end.

Sometimes two moves reach the same value, and then both lead to optimal alignments. In the example, cell (4, 3) can be reached diagonally from (3, 2) and from above, from (3, 3), both with the value 1. The walk splits there. One branch gives GA-TCA, the other GAT-CA, so the gap sits opposite one T or the other. Both alignments score 3, and the scores alone give no reason to prefer either. A tool that reports a single alignment breaks such ties with a fixed rule.

What it costs

  • Time: each of the (|x| + 1) · (|y| + 1) cells is computed once from three neighbours, so filling the table takes O(|x| · |y|) time.
  • Memory: the whole table is kept for the traceback, so memory is O(|x| · |y|) too. Two sequences of 10,000 bases need about 100 million cells, which is why long alignments use a linear-space variant.
  • Traceback: one optimal path has at most |x| + |y| steps, so reading one alignment back is cheap compared with filling the table.

Where this leads

  • Proteins: match and mismatch are replaced by a substitution matrix such as BLOSUM62, where each pair of amino acids has its own score.
  • Gap runs: an insertion of several bases is usually one event, so an affine penalty charges more to open a gap than to extend it. That needs three tables instead of one.
  • Local alignment: Smith–Waterman keeps the same three moves but lets a cell fall back to zero, so it finds the strongest shared stretch instead of aligning the full length.

Related: Two ways to align two sequences

All notes