Phylogenetics

The likelihood of a tree, and the ancestors within it

Last time

Any model of sequence evolution is a rate matrix \(Q\), and \(P(t) = e^{Qt}\) gives us its transition probabilities

With \(P(t)\), we calculated the likelihood of the divergence between two sequences

But most questions in biology involve many sequences, related by a tree

Trees

Tree terms

Term Meaning
Tip nodes The sequences we have
Internal nodes Their ancestors (no data!)
Branch Line connecting nodes; a lineage evolving through time
Branch length Expected substitutions per site (like \(t\))
Root The ancestor of everything
Topology The branching pattern
Clade An ancestor and all its descendants

Representing trees in computers

Newick format: Parentheses hold descendants; numbers after colons are branch lengths

((Human:0.1,Chimp:0.1)Primates:0.3,(Mouse:0.2,Rat:0.2)Rodents:0.2)Root;
Index Name Children Length
0 Human 0.1
1 Chimp 0.1
2 Primates 0, 1 0.3
3 Mouse 0.2
4 Rat 0.2
5 Rodents 3, 4 0.2
6 Root 2, 5

In the exercise, we represent a tree as a list of nodes in postorder

Every node comes after its children, so the root is last

Fun factoid

Newick format is named after Newick’s seafood restaurant in Dover, NH, where a committee of phylogeneticists met in 1986 to agree on the format

Picture of Newick's restaurant

Newick’s restaurant

How many trees?

Tree space is enormous

Unrooted trees with \(n\) tips:1 \(\;1 \times 3 \times 5 \times \cdots \times (2n - 5)\)

Tips Unrooted trees At 1 million trees per second
4 3
10 2,027,025 2 seconds
15 \(7.9 \times 10^{12}\) 92 days
20 \(2.2 \times 10^{20}\) 7 million years
24 \(5.6 \times 10^{26}\) \(1.8 \times 10^{13}\) years

We can’t try them all! Programs like IQ-TREE and RAxML use heuristic search: Make small rearrangements, and keep the ones that improve the likelihood

The likelihood of a tree

We don’t have data from the ancestors

For two sequences, each site needed one lookup. E.g., \(P_{A \rightarrow G}(t)\)

On a tree, we have to add up over every possible combination of bases at the internal nodes

An unrooted tree of 24 genomes has 22 internal nodes:

\(4^{22} \approx 1.8 \times 10^{13}\) combinations, for every site

Question: Is there a way to avoid repeating work?

Felsenstein’s pruning algorithm1

For each node, from the tips to the root, calculate its conditional likelihoods: \(L_{\text{node}}(i)\) = probability of the data at the tips that descend from the node, if the node has base \(i\)

In doing this, the pruning algorithm follows 3 rules

Pruning algorithm Rule 1

At each tip node: the conditional likelihood is 1 for the base we observed and 0 for the other three bases

For example:

  • The probability we observe an “A” in our data, given the true base is A = 1
  • The probability we observe an “A” in our data, given the true base is C = 0

Pruning algorithm Rule 2

At each internal node: For each base \(i\), multiply together, over each child, the probability of getting the child’s data if the node has base \(i\):

\[ L_{\text{parent}}(i) = \prod_{\text{children}} \Big[ \sum_{j} P_{ij}(t_{\text{child}}) \, L_{\text{child}}(j) \Big] \]

For example, if the parent (\(i\)) has A and the child (\(j\)) has G, the term in the sum is the probability that A evolved to G along the branch connecting them (\(P_{A\rightarrow G}(t_{\text{child}})\)), times the probability of the data at the tips that descend from the child, if the child has G (\(L_{\text{child}}(G)\))

This is dynamic programming (like Needleman-Wunsch alignment!): Each \(L_{\text{child}}\) is calculated once, and re-used for every base of the parent

Pruning algorithm Rule 3

At the root, the likelihood of the site is the sum of the conditional likelihoods of the four bases, weighted by the base frequencies:

\(L_{\text{site}} = \sum_{i} \pi_i \, L_{\text{root}}(i)\)

This is a weighted average!

For example:

  • For \(i =\) A: The probability of the root having A (\(\pi_A\)), times the probability of the data at all the tips, if the root has A (\(L_{\text{root}}(A)\))
  • Add these products for all four bases

Pruning, step by step

Pruning one site

Node \(L(\texttt{A})\) \(L(\texttt{C})\) \(L(\texttt{G})\) \(L(\texttt{T})\)
Human, Chimp 1 0 0 0
Mouse, Rat 0 0 1 0
Primates 0.8215 0.0010 0.0010 0.0010
Rodents 0.0034 0.0034 0.6797 0.0034
Root 0.0266 0.0029 0.0385 0.0029

\(L_{\text{site}} = \frac{1}{4}(0.0266 + 0.0029 + 0.0385 + 0.0029) = 0.0177\)

Two shortcuts

  • Missing data: A gap at a tip gives conditional likelihoods of \((1, 1, 1, 1)\)
    • A site with no data has likelihood 1 (\(\log 1 = 0\))
  • Site patterns: Identical columns of the alignment have identical likelihoods
    • Four SARS-CoV-2 genomes: 29,903 sites, but only 56 patterns

Which tree is best?

SARS-CoV-2 lineages

Name Pango lineage1
Wuhan B (the reference genome)
B.1 B.1
Delta B.1.617.2
Alpha B.1.1.7
Gamma P.1 = B.1.1.28.1
Omicron BA.1.1 = B.1.1.529.1.1

Lineage names work like addresses: B.1.1 descends from B.1, which descends from B

Three trees (unrooted) for four genomes

Tree Log-likelihood Worse by
1: (Wuhan, Delta), (Alpha, Omicron) −42,433.11
2: (Wuhan, Alpha), (Delta, Omicron) −42,444.08 10.96
3: (Wuhan, Omicron), (Delta, Alpha) −42,447.15 14.04

How much does each kind of site favor Tree 1 over Tree 2?

Kind of site Sites Tree 1 − Tree 2
Two pairs (e.g., AATT) 4 +12.49
All four the same 29,750 +2.01
One genome differs 98 −4.44
Three or more bases 2 +0.93
Missing data 49 −0.02

Other methods

Parsimony

  • Maximum parsimony: Choose the tree that needs the fewest changes
  • Fast, and intuitive: For the four genomes in the exercise, it would use only the “two pairs” sites
    • 3 favor Tree 1 and 1 favors Tree 2
    • All other sites are ignored!
  • But it ignores branch lengths and multiple hits, so it can be positively misleading: Long branches attract each other1
  • Still useful when branches are very short: UShER2 uses parsimony to place millions of SARS-CoV-2 genomes on a tree

Distance-based algorithms

  • We will skip these today
  • But you used two of these methods in Capstone 1!
    • Neighbor-joining and UPGMA

Reconstructing ancestors

What base did the root have?

The conditional likelihoods at the root also give us the probability of each base:

\[ \Pr(\text{root has base } i \mid \text{data}) = \frac{\pi_i \, L_{\text{root}}(i)}{\sum_{k} \pi_k \, L_{\text{root}}(k)} \]

For our example site:

Base (\(i\)) \(L_\text{root}(i)\) \(\Pr(\text{root has base } i \mid \text{data})\)
A 0.0266 0.375
C 0.0029 0.042
G 0.0385 0.542
T 0.0029 0.042

Question: But what about the other ancestors?

The pulley principle1

  • For models like JC, K80, and HKY, the likelihood doesn’t depend on where the root is
  • We can slide the root anywhere, like a rope over a pulley
  • So, to reconstruct any ancestor: Reroot the tree at that ancestor, and use the probabilities at the root

The ancestors of SARS-CoV-2 lineages

The maximum likelihood tree of six genomes (JC model)

Reconstructed bases

Position Wuhan B.1 Delta Alpha Gamma Omicron B.1 ancestor B.1.1 ancestor
3,037 C T T T T T T T
14,408 C T T T T T T T
23,403 (D614G) A G G G G G G G
28,881–3 GGG GGG TGG AAC AAC AAC GGG AAC

Every reconstruction has a probability of about 0.9999 or higher

Confident is not the same as correct

  • Alpha, Gamma, and Omicron all have N501Y (position 23,063, in the spike protein)
  • So our reconstruction says the B.1.1 ancestor had it too (probability 0.9999)
  • But N501Y arose independently in each of these lineages, favored by natural selection (convergent evolution)
  • Our model knows nothing about selection, and our sample has no early B.1.1 genomes

Summary

  • A tree can be stored as a list of nodes in postorder
  • The number of trees grows astronomically, so programs use a heuristic search
  • Felsenstein’s pruning algorithm uses dynamic programming to calculate the likelihood of a tree from the tips to the root
  • A few sites strongly favor one tree, but every site informs the branch lengths; parsimony uses only those few sites, and can be misled by long branches
  • The root’s conditional likelihoods give ancestral base probabilities, and the pulley principle lets us reroot to reconstruct any ancestor

In the exercise

  1. get_tip_names: Loop over a tree’s list of nodes
  2. number_of_unrooted_trees: Count tree space
  3. tip_likelihoods, combine_children, get_root_likelihoods, site_log_likelihood: Build the pruning algorithm, one rule at a time
  4. Compare three trees for four SARS-CoV-2 genomes
  5. root_state_probabilities and reconstruct_ancestor: Reconstruct the ancestors of SARS-CoV-2 lineages

Watch for

Does rerooting the tree change its likelihood?

Exercise

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

References

Felsenstein, Joseph. 1978a. “Cases in Which Parsimony or Compatibility Methods Will Be Positively Misleading.” Systematic Zoology 27 (4): 401–10. https://doi.org/10.2307/2412923.
Felsenstein, Joseph. 1978b. “The Number of Evolutionary Trees.” Systematic Biology 27 (1): 27–33.
Felsenstein, Joseph. 1981. “Evolutionary Trees from DNA Sequences: A Maximum Likelihood Approach.” Journal of Molecular Evolution 17: 368–76.
Rambaut, Andrew, Edward C. Holmes, Áine O’Toole, et al. 2020. “A Dynamic Nomenclature Proposal for SARS-CoV-2 Lineages to Assist Genomic Epidemiology.” Nature Microbiology 5 (11): 1403–7. https://doi.org/10.1038/s41564-020-0770-5.
Turakhia, Yatish, Bryan Thornlow, Angie S. Hinrichs, et al. 2021. “Ultrafast Sample Placement on Existing tRees (UShER) Enables Real-Time Phylogenetics for the SARS-CoV-2 Pandemic.” Nature Genetics 53 (6): 809–16. https://doi.org/10.1038/s41588-021-00862-7.
Yang, Ziheng, Sudhir Kumar, and Masatoshi Nei. 1995. “A New Method of Inference of Ancestral Nucleotide and Amino Acid Sequences.” Genetics 141 (4): 1641–50. https://doi.org/10.1093/genetics/141.4.1641.

Acknowledgements

I drafted and iterated on these slides in collaboration with Anthropic’s Claude (Opus 5.5 model).