Aligning sequences

Finding homology when the sequence has changed

Last time

We searched a genome for exact matches to sequences of interest.

CGATGCATGC
  ATG  <- match!

Every method we built — linear search, \(k\)-mer dictionary, suffix array — depends on an exact character-for-character match.

But sequences evolve

Two copies of the same gene in two individuals or species are often not identical

ACGTCA
ACTGA

Are these related? Which positions correspond to which?

The moment evolution edits a sequence, exact matching stops working.

Why align sequences?

Homology

Two sequences are homologous if they are descended from a common ancestral sequence.

This is a statement about history, not about similarity.

  • Homology is not a matter of degree
    • There is no such thing as “80% homologous”
  • Two sequences either share ancestry or not
  • Similarity is the evidence. Homology is the hypothesis

An alignment is a hypothesis

Aligning two sequences means proposing, for every column, which nucleotide residues descend from the same ancestral residue.

A C G T C A
A C - T G A
  • Columns 1, 2, 4, 6: Homologous positions, unchanged
  • Column 5: Homologous position, a substitution happened
  • Column 3: An insertion or deletion happened in one lineage

An alignment is a hypothesis about a series of evolutionary events

What is an “indel”?

A C G T C A
A C - T G A
  • Column 3: An insertion or deletion happened in one lineage

We don’t know if G was inserted in the top sequence or deleted from the bottom sequence.

So, we call it an indel

Why we need alignments

  1. Evolutionary inference. We cannot learn about evolutionary history without assuming which positions correspond. Every molecular evolutionary analysis starts from an alignment!

  2. Locating a sequence. E.g., where in the genome did this read come from? Exact search fails on sequencing errors and biological variation.

  3. Transferring knowledge. A new gene of unknown function is compared to genes of known function. The inference is, similar sequence, shared ancestry, probably shared function.

Three choices to make when aligning

Choice 1: How many sequences?

Pairwise

Two sequences at a time.

ACGTCA
AC-TGA

Tractable and exactly solvable.

Multiple

Many sequences at once.

ACGTCA
AC-TGA
ACGT-A
AGGTCA

What we often want for analysis, but requires heuristics

Exact multiple alignment of \(N\) sequences means filling an \(N\)-dimensional table. This is infeasible beyond a handful of sequences.

Choice 2: How much of each sequence?

Global

Align the sequences end to end.

Assumes homology runs the full length of both.

Local

Find the best-matching region and ignore the rest.

Assumes homology is confined to a subsequence.

This is a biological assumption, not a technical preference.

E.g., Aligning two homologs of the same short gene? Use global.

E.g., Aligning a conserved domain in two otherwise unrelated proteins? Use local.

Choice 3: Optimal or approximate?

  • Optimal: Guaranteed to find the highest-scoring alignment under your scoring scheme
  • Approximate: Fast, usually accurate, but no optimality guarantee

Recall from our algorithms lecture: BLAST1 is a heuristic. It finds exact \(k\)-mer seeds, then extends only around those.

BLAST can miss an optimal alignment. It is used billions of times a year because an approximate answer today is better than an exact answer in 9 years.

Today, we will focus on

  • Pairwise, not multiple
  • Optimal, not approximate
  • Global (Needleman-Wunsch) and local (Smith-Waterman)

Both are dynamic programming algorithms, the strategy we used to compute Fibonacci numbers

  • Solve each subproblem once, store the answer, and reuse it

Scoring alignments

Which alignment is better?

A C G T C A
A C - T G A
A C G T C A
A C T G A -

Both are valid hypotheses. We need a way to choose.

An objective function: a rule that assigns a score to any alignment, so “best” means something specific.

A simple scoring scheme

Score each column independently, then add them up:

Column Score
Match \(+1\)
Mismatch \(-1\)
Gap \(-2\) 
A  C  G  T  C  A
A  C  -  T  G  A
+1 +1 -2 +1 -1 +1 = 1
A  C  G  T  C  A
A  C  T  G  A  -
+1 +1 -1 -1 -1 -2 = -3

The left alignment wins, under this scoring scheme.

The scoring scheme is a model

  • Gaps cost more than mismatches because indels tend to be rarer than substitutions
  • Changing the numbers can change the answer
  • The “optimal” alignment is only optimal with respect to your scoring scheme

The algorithm is guaranteed to find a best-scoring alignment. It is not guaranteed to find the historically correct one.

Can we try all alignments?

The number of possible global alignments of two sequences of length \(n\):

\(n\) Alignments
5 1,683
10 \(\approx 8 \times 10^{6}\)
100 \(> 10^{75}\)

Exhaustive search becomes impossible quickly!

We need to find the best alignment without looking at all of them

Needleman-Wunsch

Their insight

The best alignment of ACGTCA and ACTGA ends in one of three ways:

  • A aligned with A (match or mismatch)
  • A aligned with a gap
  • a gap aligned with A

Whichever it is, the rest of the alignment must itself be the best alignment of the two shorter sequences that remain.

A larger problem made of smaller versions of itself, solved once and reused

The recurrence1

\[ F(i,j) = \max \begin{cases} F(i-1,\,j-1) + s(A_i, B_j) & \text{(match/mismatch)}\\ F(i-1,\,j) + g & \text{(gap in Seq }B)\\ F(i,\,j-1) + g & \text{(gap in Seq }A) \end{cases} \]

\(F(i,j)\) is the score of the best alignment of the first \(i\) residues of Sequence \(A\) with the first \(j\) residues of Sequence \(B\).

Every cell depends only on cells above and to the left, so we can fill the table in one pass.

Setting up the table

\(A\) = ACGTCA down the rows, \(B\) = ACTGA across the columns.

– A C T G A
– 0 -2 -4 -6 -8 -10
A -2
C -4
G -6
T -8
C -10
A -12

Row 0 and Column 0 are alignments against nothing but gaps, so they are just multiples of the gap penalty.

Filling the first cell

\(F(1,1)\): best alignment of A with A

– A C T G A
– 0 -2 -4 -6 -8 -10
A -2 ?
C -4
G -6
T -8
C -10
A -12
From Calculation Score
Diagonal (A with A, match) 0 + 1 1
Up (gap in \(B\)) -2 + -2 -4
Left (gap in \(A\)) -2 + -2 -4

Take the maximum: \(F(1,1) = 1\)

Filling a cell where a gap winds

\(F(3,2)\): best alignment of ACG with AC

– A C T G A
– 0 -2 -4
A -2 1 -1
C -4 -1 2
G -6 -3 ?
T
C
A
From Calculation Score
Diagonal (G with C, mismatch) -1 + -1 -2
Up (gap in \(B\)) 2 + -2 0
Left (gap in \(A\)) -3 + -2 -5

\(F(3,2) = 0\), from above: the G is aligned to a gap.

This one cell is where the indel in our final alignment gets decided.

The completed table

– A C T G A
– 0 -2 -4 -6 -8 -10
A -2 1 -1 -3 -5 -7
C -4 -1 2 0 -2 -4
G -6 -3 0 1 1 -1
T -8 -5 -2 1 0 0
C -10 -7 -4 -1 0 -1
A -12 -9 -6 -3 -2 1

The bottom-right cell is the score of the best global alignment: 1.

Every cell was computed from cells already filled. No cell was computed twice.

The table is NOT an alignment

– A C T G A
– 0 -2 -4 -6 -8 -10
A -2 1 -1 -3 -5 -7
C -4 -1 2 0 -2 -4
G -6 -3 0 1 1 -1
T -8 -5 -2 1 0 0
C -10 -7 -4 -1 0 -1
A -12 -9 -6 -3 -2 1

We know the best score (1), but we do not know the alignment that achieves it.

Traceback to get the alignment

– A C T G A
– 0 -2 -4 -6 -8 -10
A -2 1 -1 -3 -5 -7
C -4 -1 2 0 -2 -4
G -6 -3 0 1 1 -1
T -8 -5 -2 1 0 0
C -10 -7 -4 -1 0 -1
A -12 -9 -6 -3 -2 1

Start at the bottom-right and walk backward, asking at each cell: which of the three directions did the score come from?

  • A diagonal step consumes one nucleotide from each sequence.
  • A vertical step consumes one from \(A\) and adds a gap to \(B\).
  • A horizontal step consumes one from \(B\) and adds a gap to \(A\).

The (an) optimal alignment

Seq A: A  C  G  T  C  A
Seq B: A  C  -  T  G  A
Score: +1 +1 -2 +1 -1 +1 = 1
  • The score in the bottom-right corner of the scoring matrix and the score of the alignment agree
  • Ties are possible: Several different alignments can share the best score
  • The algorithm returns one of them, which one depends on how you break ties

What did it cost? (Complexity)

  • Time: one pass over an \((n{+}1) \times (m{+}1)\) table, so \(O(nm)\)
  • Memory: the whole table, also \(O(nm)\)

Two 100,000-bp sequences: \(10^{10}\) cells, about 40 GB at 4 bytes each.

Optimal alignment is doable for genes, but not genomes.

This is why read mappers index instead of aligning each read to the genome.

When global fails

Two sequences sharing one region

A = TTTTACGTCA
B = ACGTCAGGGG

They share ACGTCA, but it sits at opposite ends.

What does the Needleman-Wunsch algo do?

Needleman-Wunsch on these sequences

T T T T A C G T C A
A C G T C A G G G G

Score: \(-6\). No gaps, and the shared region is missed!

  • Lining up the two copies of ACGTCA would require eight gaps: \(-16\), plus \(+6\) for the matches, so \(-10\)
  • A bad global alignment at \(-6\) beats the truth at \(-10\)

The algorithm did not fail. We asked the wrong question! We assumed homology runs the full length.

Smith-Waterman1 algorithm

Three changes to Needleman-Wunsch:

  1. Add zero to the options. Any cell that would go negative becomes 0. A bad alignment is abandoned rather than carried forward.
  2. Start the traceback at the highest cell anywhere, not the bottom-right corner.
  3. Stop the traceback at a zero, not at the top-left corner.

The local alignment table

– A C G T C A G G G G
– 0 0 0 0 0 0 0 0 0 0 0
T 0 0 0 0 1 0 0 0 0 0 0
T 0 0 0 0 1 0 0 0 0 0 0
T 0 0 0 0 1 0 0 0 0 0 0
T 0 0 0 0 1 0 0 0 0 0 0
A 0 1 0 0 0 0 1 0 0 0 0
C 0 0 2 0 0 1 0 0 0 0 0
G 0 0 0 3 1 0 0 1 1 1 1
T 0 0 0 1 4 2 0 0 0 0 0
C 0 0 1 0 2 5 3 1 0 0 0
A 0 1 0 0 0 3 6 4 2 0 0

\[ F(i,j) = \max \begin{cases} 0 \\ F(i-1,\,j-1) + s(A_i, B_j)\\ F(i-1,\,j) + g\\ F(i,\,j-1) + g \end{cases} \]

Mostly zeros, with one diagonal ridge climbing to 6.

The local answer

A C G T C A
A C G T C A

Score 6, positions 5–10 of \(A\) against positions 1–6 of \(B\).

Use global when

You believe the sequences are homologous along their whole length

Use local when

You believe homology is confined to a region, or the sequences differ a lot in length

Our examples were simple, but generalizing is easy

Change the scoring, not the algorithm

Nothing in either algorithm “knows” what a nucleotide is.

The only place nucleotides appear is \(s(A_i, B_j)\).

  • Proteins: Replace \(+1/-1\) with a 20×20 matrix such as BLOSUM621, where a leucine/isoleucine substitution costs far less than a leucine/aspartate substitution
  • DNA, more nuanced: E.g., penalize transversions more than transitions
  • Same table, same traceback

Two more refinements

  • Affine gaps: A large cost to open a gap and small cost to extend it, because one 10-base indel is far more likely than ten 1-base indels
  • Multiple alignment: Exact is infeasible, so tools align the most similar pair first and progressively add sequences
    • A greedy heuristic built on pairwise alignment

Every multiple sequence alignment tool is an approximation built upon the algorithm we just walked through.

Summary

  • An alignment is a hypothesis of homology, column by column
  • A scoring scheme is a model of evolution that gives us an optimality criterion
  • Dynamic programming finds the optimal alignment without enumerating the \(10^{75}\) alternatives
  • Global (Needleman–Wunsch) assumes end-to-end homology; Local (Smith–Waterman) does not
  • Optimal means optimal under your assumptions, not evolutionarily correct

In the exercise

You will build a global aligner from scratch:

  1. score_alignment: Score a pair of aligned sequences
  2. get_dotplot: Visualize the homology before any aligning
  3. fill_scoring_matrix: The Needleman-Wunsch table
  4. traceback_alignment: Turn the table into an alignment

Watch for

The dotplot of the two sequences from the “when global fails” slide. You will be able to see why the global algorithm made the choice it did.

Exercise

  1. Log into molab
  2. Open the exercise in molab
  3. Click “Fork” button

References

Altschul, Stephen F., Warren Gish, Webb Miller, Eugene W. Myers, and David J. Lipman. 1990. “Basic Local Alignment Search Tool.” Journal of Molecular Biology 215 (3): 403–10. https://doi.org/10.1016/S0022-2836(05)80360-2.
Henikoff, Steven, and Jorja G Henikoff. 1992. “Amino Acid Substitution Matrices from Protein Blocks.” Proceedings of the National Academy of Sciences 89 (22): 10915–19. https://doi.org/10.1073/pnas.89.22.10915.
Needleman, Saul B., and Christian D. Wunsch. 1970. “A General Method Applicable to the Search for Similarities in the Amino Acid Sequence of Two Proteins.” Journal of Molecular Biology 48 (3): 443–53. https://doi.org/10.1016/0022-2836(70)90057-4.
Smith, T. F., and M. S. Waterman. 1981. “Identification of Common Molecular Subsequences.” Journal of Molecular Biology 147 (1): 195–97. https://doi.org/10.1016/0022-2836(81)90087-5.

Acknowledgements

I used Anthropic’s Claude (Opus 5 model) to generate a draft of these slides.