Summary

• We’ve extensively discussed population structure, but now we want to capture the history of shared ancestry among populations and species.
• A phylogenetic tree is a hypothesis of shared ancestry, estimated from characters (usually DNA) under an explicit model.
• Today: the likelihood of a tree from first principles, models of sequence evolution, model selection and ML.
• Thursday: the multispecies coalescent, species-tree estimation, coalescent species delimitation, and introgression.

Anatomy of a phylogeny

Rooting gives time a direction

Placing a root r on the internal branch gives one rooted tree: \(((A,B),(C,D))\).

Same unrooted tree, different root

Root on the branch to A instead and the same unrooted tree reads \((A,(B,(C,D)))\). An unrooted four-taxon tree has five branches, so five possible rootings — the outgroup picks one.

The inference pipeline

From sequences to a supported tree — model choices enter at every step:

1. Homology / alignment — orthology assessment, multiple sequence alignment.
2. Model — a substitution model (partitions, rate variation).
3. Search / inference — optimize or sample trees (ML, Bayesian).
4. Support — bootstrap or posterior probabilities.

Distance & parsimony

• Distance methods (e.g. neighbor-joining): reduce sequences to pairwise distances, build a tree by clustering. Fast, but limitations beyond populations or closely related species.

• Parsimony: choose the tree minimizing the number of changes. Intuitive, but can fail under various scenarios such as high rates of evolution (“long-branch attraction”).

• Parsimony built the first broad molecular trees of plants: 475 and 499 seed-plant sequences of the plastid gene rbcL (Chase et al. 1993).

Two rules of probability

Roll a fair 10-sided die twice.

Every likelihood on a tree is built from these two rules: multiply along a history, add over the histories we can’t observe.

Probability vs likelihood

I rolled 2, 1, 1, 4. Is that surprising?

\[\Pr\{2,1,1,4 \mid \text{10-sided die}\} = \left(\tfrac{1}{10}\right)^4 = \tfrac{1}{10\,000}\]

The likelihood function

The likelihood of a model is the probability of the observed data under it:

\[ \begin{aligned} L(\text{4-sided die}) &= f(X \mid \text{4-sided die}) = \prod_{h=1}^{N} f(x_h \mid \text{4-sided die}) \\[4pt] L(\theta) &= f(X \mid \theta) = \prod_{h=1}^{N} f(x_h \mid \theta) \end{aligned} \]

• \(X = (x_1,\dots,x_N)\) are the \(N\) observations (for an alignment, its sites); \(\theta\) is the model parameters. Do not confuse with our popgen parameter \(\theta\)
• The product is the AND rule — observations are assumed independent.
• Maximum likelihood: pick the \(\theta\) (a combination of many parameters) that makes the observed data least surprising.

Why programs always report \(\log L\)

Logs turn products into sums

\[\log(x\times y) = \log x + \log y \quad\Longrightarrow\quad \ell(\theta) = \log L(\theta) = \sum_{h=1}^{N}\log f(x_h\mid\theta)\]

The log-likelihood is maximized by the model that surprises us the least.

Pretend we know the ancestors

The root: \(\pi_A\)

Branch \(v_1\): root → tip A

Branch \(v_2\): root → internal node

Branches \(v_3\) and \(v_4\): A → T, twice

Two probability components to keep track of

But we never know the ancestors

One term at a time

Each choice of ancestral states is one history; its probability is a product. The site likelihood is the sum over all of them.

All 16 histories

With \(n\) tips there are \(n-1\) internal nodes, so \(4^{\,n-1}\) histories per site: 20 taxa → \(4^{19} \approx 2.7\times10^{11}\).

A lot of repetition

The lower-subtree term \(P_{yT}(v_3)P_{yT}(v_4)\) depends only on \(y\) — yet it’s recomputed for every root state.

Pruning: do the repeated part once

Same sum, grouped so the repeated part appears only once (\(ab+ac=a(b+c)\)):

\[ L = \sum_{x}\pi_x\ \underbrace{P_{xA}(v_1)}_{\text{tip 1}}\ \sum_{y}P_{xy}(v_2)\ \underbrace{\big[P_{yT}(v_3)\,P_{yT}(v_4)\big]}_{\text{the two T tips, if the internal node is } y} \]

The bracket depends only on \(y\), not on the root state \(x\). So compute it once for each of the 4 values of \(y\), write the 4 numbers down, and reuse them for every root state.

Pruning, step 1: start at the tips

A tip’s base is known, so its list is simple: the observed base explains the data perfectly (1), and any other base cannot (0).

Our column, tips A, T, T:

Pruning, step 2: move up one node

For each possible base at a node, handle the two children separately, then combine:

1. One child: the chance the node’s base turns into each possible child base, times that child’s list, added up over the child’s bases.
2. Two children: multiply their answers (given the node’s base, the two sides evolve independently).

Example: node \(S\) (children: tips 3 and 4, both T). If \(S\) has base \(y\):

• Tip 3: \(\sum_x P_{yx}(v_3)\,L_3(x) = P_{yT}(v_3)\), because tip 3’s list is 0 except at T.
• Tip 4: likewise \(P_{yT}(v_4)\).
• So \(L_S(y) = P_{yT}(v_3)\times P_{yT}(v_4)\): the conditional likelihood we reused above.

In symbols, for any node \(k\) with children \(i\) and \(j\):

\[ L_k(s) = \Big[\sum_{x}P_{sx}(v_i)\,L_i(x)\Big]\times\Big[\sum_{x}P_{sx}(v_j)\,L_j(x)\Big] \]

Sum = OR, multiply = AND

Pruning, step 3: the root, then done

Treat the root \(R\) like any other node (children: tip 1 and \(S\)) to get its list \(L_R\). Nothing lies above it, so weight each base by its equilibrium frequency \(\pi_s\) and add:

\[ L_{\text{site}} = \sum_{s}\pi_s\,L_R(s) \]

The whole algorithm for one site:

1. Write the 0/1 lists for the tips.
2. Work upward, node by node, using the step 2 recipe.
3. At the root, weight by \(\pi\) and add.

Cost: one small calculation per internal node, so it grows with the number of taxa, versus \(4^{\,n-1}\) histories by brute force (Felsenstein 1981). Repeat for every site.

A worked number: the set-up

Along a branch, a base either stays the same or changes to one particular other base:

Check: 0.90638 + 3 × 0.03121 = 1 (stay, or change to one of the 3 other bases).

A worked example: node \(S\)

\(S\) has two T tips on branches of length 0.1. For each possible state of \(S\), how likely are two T tips?

If \(S\) were anything other than T, two separate changes would be needed: about 850 times less likely.

We keep these as a list with one entry per state of \(S\):

\[ L_S = (A,\ C,\ G,\ T) = (0.00097,\ 0.00097,\ 0.00097,\ 0.82152) \]

A worked example: the probability of the whole graph

Each root state is equally likely a priori (\(\tfrac14\)). Weight the products and add:

\[ L = \tfrac14\,(0.02191 + 0.00156 + 0.00156 + 0.04358) = \mathbf{0.01715}, \qquad \log L = -4.066 \]

• Check: adding up all 16 root/internal-state histories by brute force gives the same number. Pruning is the same sum, organized.
• Which root state is most plausible? Each one’s share of the total: T 64%, A 32%, C 2%, G 2%. The two T tips outvote the single A, because tip 1’s branch (0.2) is the longest, so it says the least about the root.
• That was one site. An alignment has many sites, and we combine their likelihoods next.

From one site to an alignment

Sites are assumed independent — the AND rule again. \(L_h\) is the likelihood of site \(h\) (the pruning result):

\[ L = \Pr(\text{data}\mid \tau, \text{model}) = \prod_{h=1}^{N} L_h \qquad\Longrightarrow\qquad \ell = \sum_{h=1}^{N}\log L_h \]

• Maximum likelihood (ML): find the topology \(\tau\), branch lengths and model parameters that maximize \(\ell\) (RAxML, IQ-TREE).
• Identical columns (site patterns) are computed once and weighted by their counts.

Where do the \(P_{ij}(v)\) come from?

Branch lengths confound rate and time

The substitution process is a Markov chain

Model each site as a continuous-time Markov chain over states {A, C, G, T}. The instantaneous rate matrix \(Q\) gives the rate of change from state \(i\) to \(j\):

\[ Q = \{q_{ij}\},\qquad q_{ii} = -\sum_{j\neq i} q_{ij}\ \text{(rows sum to 0)}. \]

The probability of ending in state \(j\) after time \(t\), starting in \(i\):

\[ P(t) = e^{Qt}, \qquad P_{ij}(t) = \Pr(j \text{ at } t \mid i \text{ at } 0). \]

• \(P(t)\) averages over all unobserved histories on a branch — including multiple hits.
• Markov process: only the current state matters, not how it got there. Named after Andrey Markov (1856–1922).
• Time-reversible: \(\pi_i q_{ij} = \pi_j q_{ji}\) — the process looks the same run backwards.

Jukes & Cantor (1969): the simplest \(Q\)

Deriving \(P(t)\) for JC69: the set-up

• Let “events” happen at rate \(\mu\). Each event draws a new base uniformly from all four — so an event can redraw A as A.
• A specific change, say A→C, happens at rate \(\mu/4\). That’s \(\lambda\), so \(\mu = 4\lambda\).

Deriving \(P(t)\) for JC69: the Poisson step

The number of events in time \(t\) is Poisson with mean \(\mu t\):

\[ \begin{aligned} \Pr(\text{no event}) &= e^{-\mu t} \qquad\qquad \Pr(\text{at least one event}) = 1 - e^{-\mu t}\\[4pt] \Pr(\text{last event drew a T}) &= \tfrac14 \end{aligned} \]

In branch-length units

Measure a branch by the expected number of real substitutions per site, \(v\): the rate of leaving a base times time, \(v = 3\lambda t\). So \(4\lambda t = \tfrac43 v\) and

\[ P_{ii}(v) = \tfrac14 + \tfrac34 e^{-4v/3}, \qquad P_{ij}(v) = \tfrac14 - \tfrac14 e^{-4v/3}. \]

The model corrects for multiple hits

HKY85: transitions ≠ transversions

GTR: let the data decide

\[ Q = \begin{array}{c|cccc} & A & C & G & T\\ \hline A & - & \pi_C a\mu & \pi_G b\mu & \pi_T c\mu\\ C & \pi_A a\mu & - & \pi_G d\mu & \pi_T e\mu\\ G & \pi_A b\mu & \pi_C d\mu & - & \pi_T f\mu\\ T & \pi_A c\mu & \pi_C e\mu & \pi_G f\mu & - \end{array} \]

• Six exchangeabilities \(a\ldots f\) (five free, one fixed for scale) plus three free base frequencies: 8 parameters (Tavaré 1986).
• Most models in use sit between JC69 and GTR. Classic RAxML uses GTR for DNA; IQ-TREE picks a model for you.
• Amino-acid and codon models exist too, with their own complications.

A hierarchy of substitution models

Models differ in how much of \(Q\) they let vary (base frequencies \(\pi\), exchangeabilities):

Each model is a special case of the one below it: they are nested.

Rate variation among sites

Every model above can add rate variation: some sites evolve quickly, others hardly at all.

Even RuBisCO, encoded by rbcL — the workhorse gene of early plant phylogenetics — combines a mass of strongly constrained sites with a few residues under positive selection in most land-plant lineages (Kapralov and Filatov 2007).

The distribution of sites under selection in rbcL and their relevance in the context of a 3-D protein structure. Figs. 1 and 2 from Kapralov and Filatov (2007)

Approach 1: partition the alignment

Each region gets its own rate multiplier on the same \(Q\). For JC69 in region \(k\):

\[ P_{ii}(v) = \tfrac14 + \tfrac34 e^{-4 r_k v/3}, \qquad P_{ij}(v) = \tfrac14 - \tfrac14 e^{-4 r_k v/3} \]

• Pro: each site’s likelihood is computed once.
• Con: the regions must be chosen a priori, and rates still vary within them.

Approach 2: a mixture at every site

Let every site have any of \(K\) rates, and average (the OR rule) over them. With \(K=4\) equally likely rates:

\[ L(\theta) = \prod_{h=1}^{N}\left[\tfrac14\Pr(x_h\mid r_1) + \tfrac14\Pr(x_h\mid r_2) + \tfrac14\Pr(x_h\mid r_3) + \tfrac14\Pr(x_h\mid r_4)\right] \]

Here \(x_h\) is alignment column \(h\). No partitions needed. But where do \(r_1,\dots,r_4\) come from?

Rates from a gamma distribution

The discrete gamma

Beyond the gamma

• Free-rate model: estimate the \(K\) rates and their weights directly, with no gamma shape assumed (Yang 1995). Usually written \(+R_K\) in IQ-TREE.
• Needs enough data to estimate the categories; \(+\Gamma\) is still widely used.
• \(+I\): a proportion of sites that never change.

The Akaike Information Criterion

\[ \mathrm{AIC} = 2k - 2\log \hat L \]

• \(k\) = number of free parameters; \(\hat L\) = the maximum likelihood.
• More parameters always raise \(\hat L\); the \(2k\) term charges for them (Akaike 1974).
• Goal: the best — or least-worst — model among those we tried.

A true model exists

Our models live in a smaller space

An approximating model

The KL divergence contains the AIC

Several candidate models

AIC estimates each model’s distance from the truth, up to a constant shared by all of them.

Alternatives to AIC

\[ \mathrm{AIC_c} = 2k - 2\log\hat L + \frac{2k^2 + 2k}{N - k - 1} \qquad\qquad \mathrm{BIC} = k\log N - 2\log\hat L \]

• AICc corrects for small sample size \(N\) (alignment length).
• BIC penalizes parameters more heavily as \(N\) grows (Schwarz 1978).
• There’s no consensus on which is best; programs report all three. ModelFinder uses BIC by default.

Many loci, one analysis

If you can analyze one locus well, you can analyze a hundred — usually by concatenating them into one supermatrix.

Should each gene have its own model?

Partitioned analyses give each block its own \(Q\), rate variation and/or relative rate. Concatenated ML with partitioned mixture models is in almost every phylogenomics paper.

Same data, different models, different trees

What is sister to all other flowering plants? Disagreement between the slowly evolving sites and fast sites (Xi et al. 2014).

Conflicting evidence for early angiosperm evolution based on partitioning strategies and model choice. Fig. 1 from Xi et al. (2014)

The Amborella genome paper (Amborella Genome Project 2013) has Amborella alone as sister to all other angiosperms, the more accepted hypothesis.

Target capture: conserved exons, variable flanks

Data-driven partitioning and model selection

Support for ML trees: the bootstrap

• Resample alignment columns with replacement, re-estimate the tree, repeat 100–1000×.
• Bootstrap support for a clade = the fraction of replicate trees that contain it (Felsenstein 1985).
• It measures how consistently the data point to a clade, given the model. It is not the probability that the clade is true.