Modeling sequence evolution

Revealing hidden evolution with Markov chains

An aside about terms

Term Refers to Example
Site A position; in an alignment, a column of homologous positions “3rd codon positions evolve faster than other sites”
Nucleotide The building block of DNA: sugar + phosphate + base “nucleotide substitution”
Base The part of a nucleotide that makes it A, C, G, or T (the letter we write) “base frequencies”, “base pair”
Residue A building block linked into a chain: Nucleotides or amino acids “amino acid residues” 
              Site 3
                ↓
Species 1:  A C G T
Species 2:  A C A T

The site persists; a mutation changes the nucleotide/base that occupies it

So far

We have measured how different two sequences are by counting differences

  • Hamming distance in the complexity exercise
  • Pairwise distances to build a tree of SARS-CoV-2 variants

But are we measuring the thing we care about?

The number of differences between two sequences is not the same as the amount of evolution that separates them.

Hidden substitutions

Three sites in an ancestor and a descendant:

Site History Substitutions Differences we see
1 A → G 1 1
2 A → G → T 2 1
3 A → G → A 2 0

5 substitutions happened, but we only observe 2 differences

Multiple hits: A later change can hide an earlier one

Differences saturate

Figure 1

To account for hidden changes, we need a model

A model of sequence evolution describes the probabilities of substitution over time

With a model, we can ask:

  • How many differences should we expect after time \(t\)?
  • Given the differences we see, how much time most likely passed (or evolution happened)?

This is our first probabilistic model of the semester

The Jukes-Cantor model of sequence evolution

The simplest possible model

Chop time into many tiny steps. In each step, every base follows the same rule1:

  1. With a small probability (the rate), the base is substituted
  2. If it is substituted, it becomes one of the other three bases, each equally likely

That’s the entire model!

You will write a Python function to do this in the exercise

Measuring time

  • If each base has a 1% chance of changing per step, then after 100 steps, we expect 1 substitution per site (on average)
  • We measure time (\(t\)) in expected substitutions per site
  • This is how branch lengths are measured on most phylogenetic trees

\(t = 1\) does not mean every site changed once. Some changed zero times, some twice, some more.

The model as a table of rates

Rows: The base now. Columns: The base it changes to.

A C G T
A — \(\mu\) \(\mu\) \(\mu\)
C \(\mu\) — \(\mu\) \(\mu\)
G \(\mu\) \(\mu\) — \(\mu\)
T \(\mu\) \(\mu\) \(\mu\) —

Every possible change happens at the same rate (\(\mu\))

A Markov chain

  • What happens to a base next depends only on what it is now
  • Not on what it used to be
  • Not on how long it has been this base

A random process with this memoryless property is a Markov chain

Nearly every model of sequence evolution used in biology is a Markov chain

The model forgets where it started

Start 1,000 sites as A and let them evolve:

  • At first, almost all are still A
  • Over time, A is lost and C, G, T are gained
  • Eventually, about 1/4 of the sites are each base, and it stays that way

This is the model’s stationary distribution

What the model predicts

The probability that a site differs from its ancestor after time \(t\):

\[ P(\text{different}) = \frac{3}{4}\left(1 - e^{-\frac{4}{3}t}\right) \]

  • \(t = 0\): \(e^{0} = 1\), so \(P(\text{different}) = 0\)
  • \(t\) large: \(e^{-\frac{4}{3}t} \to 0\), so \(P(\text{different}) \to \frac{3}{4}\)

This is the “differences we see” curve from the saturation plot

Turning it around

We observe \(p\), the proportion of sites that differ. Solve for \(t\) (and replace with \(d\) for distance):

\[ d = -\frac{3}{4}\ln\left(1 - \frac{4}{3}p\right) \]

\(p\) (what we see) 0.01 0.10 0.30 0.50 0.70 \(\geq 0.75\)
\(d\) (Jukes-Cantor distance) 0.010 0.107 0.383 0.824 2.031 \(\infty\)
  • Small \(p\): Barely any correction needed
  • Large \(p\): The correction is large, and very sensitive to noise
  • \(p \geq 0.75\): Differences are saturated; no estimate is possible

Models are assumptions

What Jukes-Cantor assumes

  • All four bases are equally frequent
  • All substitutions are equally likely
  • Every site evolves independently and at the same rate
  • The process is the same at all times and on all lineages

Real sequences break these assumptions

Two violations of the assumptions

Transitions are more common

  • Purines: A ↔︎ G
  • Pyrimidines: C ↔︎ T
  • Transitions (within a type) usually outnumber transversions (between types), even though there are twice as many possible transversions

Sites differ in rate

  • In protein-coding genes, many 3rd codon position changes are synonymous
  • So 3rd positions usually evolve much faster than 1st and 2nd positions

You will see both in a real yeast gene in the exercise

A family of models

Same machinery, more entries in the table:

Model Adds Free parameters
JC691 Nothing: all changes equal, all bases equal 0
K802 Transitions and transversions have different rates (\(\kappa\)) 1
HKY853 Unequal base frequencies 4
GTR4 Every pair of bases has its own rate 8

Plus: Let rates vary among sites (e.g., “+Γ”5)

Each model is a special case of the next

Where does the Jukes-Cantor corrected distance come from?

Likelihood: A coin-flipping example

Flip a coin 10 times: 7 heads. If \(P(\text{heads}) = \theta\):

\[ L(\theta) = \theta^{7}(1 - \theta)^{3} \]

\(\theta\) 0.3 0.5 0.7 0.9
\(L(\theta)\) 0.000075 0.00098 0.0022 0.00048
  • The probability of the data, as a function of the parameter, is the likelihood
  • The best value (\(\theta = 0.7\)) is the maximum-likelihood estimate

Sites are coin flips

Each site “comes up different” with probability \(P(\text{different})\), which depends on \(t\)

\(n\) sites compared, \(k\) differ:

\[ L(t) = P(\text{different})^{k} \times \big(1 - P(\text{different})\big)^{n - k} \]

The value of \(t\) that maximizes \(L(t)\) is the Jukes-Cantor distance

The Jukes-Cantor distance formula we used is a maximum likelihood estimate!

Multiplying many probabilities

Thousands of sites means multiplying thousands of probabilities (numbers less than 1)

  • The likelihood of a 1,900-base gene can be \(\approx 10^{-361}\)
  • The smallest number Python can store is \(\approx 5 \times 10^{-324}\)
  • So the answer silently becomes 0.0!

The fix: \(\ln(a \times b) = \ln a + \ln b\)

\[ \ln L(t) = k \ln P(\text{different}) + (n - k) \ln\big(1 - P(\text{different})\big) \]

Why learn likelihood if we have a formula?

  • For K80, HKY, GTR, and for trees, there is usually no formula
  • But we can always calculate the likelihood, and search for the parameter values that maximize it
  • Likelihood also lets us compare how well different models (or trees) explain the same data

Coming up: the likelihood of a tree1

Summary

  • Counting differences underestimates evolution: Multiple hits, and differences saturate at 0.75
  • The Jukes-Cantor model is a Markov chain: Each step depends only on the current base
  • The Jukes-Cantor distance corrects for hidden changes, but gets noisy near saturation
  • Real genes break the model’s assumptions, so there is a family of richer models
  • The Jukes-Cantor distance is a maximum likelihood estimate

In the exercise

  1. Write the model: get_different_base and evolve_one_step
  2. Watch 1,000 sites that start as A evolve over time
  3. Watch the differences between an ancestor and its descendant saturate
  4. jc_distance: Correct for the changes you cannot see
  5. Compare two species of yeast, codon position by codon position
  6. log_likelihood: Find where the corrected distance comes from

Watch for

How good a simple model can be versus counting differences

Exercise

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

References

Felsenstein, Joseph. 1981. “Evolutionary Trees from DNA Sequences: A Maximum Likelihood Approach.” Journal of Molecular Evolution 17: 368–76.
Hasegawa, Masami, Hirohisa Kishino, and Taka-aki Yano. 1985. “Dating of the Human-Ape Splitting by a Molecular Clock of Mitochondrial DNA.” Journal of Molecular Evolution 22 (2): 160–74.
Jukes, Thomas. H., and Charles. R. Cantor. 1969. “Evolution of Protein Molecules.” Chap. 24 in Mammalian Protein Metabolism, edited by H. N. Munro, III. Academic Press.
Kimura, Motoo. 1980. “A Simple Method for Estimating Evolutionary Rates of Base Substitutions Through Comparative Studies of Nucleotide Sequences.” Journal of Molecular Evolution 16: 111–20. https://doi.org/10.1007/BF01731581.
Tavaré, Simon. 1986. “Some Probabilistic and Statistical Problems in the Analysis of DNA Sequences.” In Some Mathematical Questions in Biology: DNA Sequence Analysis, edited by Robert M. Miura. American Mathematical Society.
Yang, Ziheng. 1993. “Maximum Likelihood Estimation of Phylogeny from DNA Sequences When Substitution Rates Differ over Sites.” Molecular Biology and Evolution 10 (6): 1396–401.
Yang, Ziheng. 1994. “Maximum Likelihood Phylogenetic Estimation from DNA Sequences with Variable Rates over Sites: Approximate Methods.” Journal of Molecular Evolution 39 (3): 306–14. https://doi.org/10.1007/BF00160154.

Acknowledgements

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