Requirements

Before diving into this slide deck, we recommend you to have a look at:

• Introduction to Galaxy Analyses

Likelihood

A way to select models, given data

Why Likelihood?

• Statistical rigour! Likelihood requires a statistical model and gains scientific rigour to defend hypotheses made.
• Comparing hypotheses! We can compare models of evolution and determine which is better supported by the data.
• Accuracy! Likelihood is generally the most accurate method to estimate phylogenies.

Likelihood Definition

The likelihood of a model is a quantity that is proportional to the probability that the model gave you the data.

It is easy(ish) to calculate the probability that a given model yields a particular outcome -- the data set.

Maximum Likelihood (ML) seeks to find the model for which this probability is the highest.

Dice example

Imagine I have a set of very strange 10-sided dice, where the sides of my blue die are labelled (1,1,2,3,4,5,5,6,6,6), and the sides of the red die are labelled (1,1,2,2,3,3,4,4,5,6).

Both dice have the same range of possible outcomes, 1 through 6, but the relative probabilities are different.

score	Pred	Pblue
1	0.2	0.2
2	0.1	0.2
3	0.1	0.2
4	0.1	0.2
5	0.2	0.1
6	0.3	0.1
Total	1.0	1.0

Dice example

Suppose I choose a die and roll it three times, getting values 3, 6, and 2.

Which die is most likely to have been used?

If it were the blue die then the probability of this outcome is $0.1\times 0.3 \times 0.1 = 0.003$.

If it were the red die then this probability is $0.2\times 0.1 \times 0.2 = 0.004$.

Neither of these probabilities is high! But the larger of the two is for the red die: thus, given the data, we would say the maximum likelihood estimate of the "model" (here, which die was used) is the red die.

Slightly formal, with Bayes

We could write this example more formally as

• Data $D$ is the outcome $(3,6,2)$;

• $P(\color{blue}{B})$ and $P(\color{red}{R})$ are the probabilities of choosing the blue and red dice (say 0.5 each);

• $P(D|\color{blue}{B})$ and $P(D|\color{red}{R})$ are the probabilities of the data given the model "blue die" and "red die";

• $P(\color{blue}{B}|D)$ is the probability of the model given the data, and it can be written

$P(\color{blue}{B}|D) = \frac{P(D|\color{blue}{B})P(\color{blue}{B})}{P(D)} = \frac{P(D|\color{blue}{B})P(\color{blue}{B})}{P(D|\color{blue}{B})P(\color{blue}{B})+P(D|\color{red}{R})P(\color{red}{R})}$

-- this is Bayes' formula.

• Similarly we can write $P(\color{red}{R}|D)$ as

$P(\color{red}{R}|D) = \frac{P(D|\color{red}{R})P(\color{red}{R})}{P(D)} = \frac{P(D|\color{red}{R})P(\color{red}{R})}{P(D|\color{blue}{B})P(\color{blue}{B})+P(D|\color{red}{R})P(\color{red}{R})}$

• These conditional probabilities can be thought of as the relative likelihoods of the model, given the data.

Important Concepts learned

The dice example was very simple but it showed us important concepts:

1. We could turn around the probability of the data given a model to the probability of a model, given the data;
   
2. Of the models considered, it was possible to choose one that was more likely;
   
3. These likelihoods can be really small, and don't have to sum to 1 (so they're not probabilities really).

Likelihood values in practice

In phylogenetics the probability of an alignment given a tree and its probabilities become really tiny so in order to avoid underflow numerical errors we use the natural log of the probabilities, and call it the log-likelihood.

Probability of a substitution

Now suppose we have a probability matrix whose entries are the probabilities of nucleotide substitution in some fixed time, say 1 million years:



We could then calculate the probability of going from one nucleotide to another, say A to G, in 1M years by reading it from the graph:

$Pr(A\to G) = P_{A,G} = 0.01.$

Probability of sequences changing

To calculate the probability of multiple nucleotides changing we multiply the individual probabilities together -- which is making direct use of the assumption that each nucleotide site evolves independently of the others.

We would then calculate the probability of going from sequence AAGT to AAGA as

$$\begin{align} prob & = P(A\to A)\times P(A\to G) \times P(G\to G) \times P(T\to A) \\ & = 0.97 \times 0.01 \times 0.97 \times 0.01 \\ & \approx 9.4\times 10^{-6} \end{align}$$

Changing the time

To calculate the probability of sequences changing in our toy example for say two million years, we would square the probability matrix $P$:

Notice that the probability of going from A to G is not quite twice what it was before, since now there is a small chance that A will change to G and then back again.

So if we have a probability matrix for a fixed time period we can convert it to a probability matrix for multiples of that time.

Probabilities to Rates

This doesn't generalise to arbitrary time periods, so instead we use a rate matrix, like this:



Converting a rate matrix to a probability requires a matrix exponentiation, giving formulae like this:

$P(t) = e^{Qt}.$

We won't worry here about technical details, but note that some classes of rate matrix are easier to work with than others, the easiest being the JC69 (one parameter) matrix.

Likelihood uses rate matrices

• Many models of sequence evolution available, including different base (nucleotide) frequencies, and different sites being allowed different rates.

• Likelihoods will be expressed as a (negative) log-likelihood, and these are used to compare models.

• Typically these will be quite large (e.g., -37000, -80000), but it is the relative likelihoods that matter, so the difference between log likelihoods.

A model is selected to give the best chance of getting the phylogeny right.

A huge range of models

• More parameter-rich models fit better; we need to avoid over-fitting.

• The (log)-likelihood value can be adjusted to penalise the model, based on the number of parameters.

• Common penalties are the Akaike Information Criterion (AIC).

• There is also a "corrected" AIC (AICc), and a "Bayesian Information Criterion" penalty (BIC).

Penalised likelihood

For a model with $k$ free parameters and a log likelihood $L$, here are the penalty functions:

AIC is the original information-theoretic correction: $AIC = 2k-2\ln(L)$

AICc is a correction for small sample sizes and is equal to $AICc = AIC + \frac{2k^{2}+2k}{n-k-1}$

BIC is defined by $BIC = k\ln(n) - 2\ln(L)$ (where $n$ is a measure of sample size).

All of these are in common use but the commonest is AIC, and that is the method we use in the tutorial.

Phylogenetic likelihood is mostly directionless

The likelihood calculations in this tutorial are based on reversible models: that means that there is no difference between the likelihood of going from sequence A to sequence B of from B to A.

This means that the output trees are unrooted.

(There are some asymmetric models included in IQTree but they are well beyond the scope of this tutorial and are not in common use.)

Likelihood of a tree and model

Now we have the huge task:

 Copy
ML Search:
    for every tree $T$ in tree space {
        for each parameter in the substitution rate matrix {
            for each branch length {
                Calculate the (log-)likelihood of the alignment
            }
        }
    }

... and find the best combination!

All this means that Maximum Likelihood is slow.

Likelihood under JC69

where the asterisk is a short-hand to make the row-sums equal 0.

• This model has only one parameter.
• Note: we cannot separate this from the overall rate (just as you can't tell people how "far" somewhere is by assuming how fast they'll travel).
• Thus tree search is "only" over all trees and all branch lengths.

Likelihood under HKY85

• The HKY85 model has four parameters: $\pi_{A}$, $\pi_{G}$, $\pi_{C}$ (and this fixes $\pi_{T}$), and a transition/ transversion rate ratio $\kappa$.

• Can fix $\pi$ or estimate it from the observed frequencies, or even estimate it by searching for the best values (which means it takes longer of course).

GTR

The General Time Reversible model has nine parameters: three from the base frequencies and six more from the parameters in the substitution rate matrix:



(The Greek letters are not standardised: they just show the structure.)

Rates across sites

• Some sites are non-coding: they evolve faster than others.

• Third codon position sites are largely redundant: for most amino acids, the third nucleotide is free to vary.

• Some are under strong selection and evolve slowly.

• Some are so important that they must stay in the same state else a protein will not fold properly: these are fixed.
  
  This means that even if the same rate matrix might be applicable at different sites, the overall substitution rate might be hundreds of times faster at one site than another.
  
  Models that accommodate this are called "Rates Across Sites" or RAS models.
  
  By default, IQTree includes RAS models.

Rates across sites

While the proportions of substitution types might be the same, the overall rate at these two sites can differ.

IQTree

IQTree is at the forefront of maximum likelihood phylogenetic tools. It is by far the most advanced such tool and is under active development, adding new models and features.

It is available from IQtree.org.

To cite IQTree2, please use

B.Q. Minh, H.A. Schmidt, O. Chernomor, D. Schrempf, M.D. Woodhams, A. von Haeseler, R. Lanfear (2020) IQ-TREE 2: New models and efficient methods for phylogenetic inference in the genomic era. Mol. Biol. Evol., 37:1530-1534. https://doi.org/10.1093/molbev/msaa015

Tree Space

A way to navigate among trees

What is "tree space"?

A network:

• whose nodes are trees, and
• whose edges are adjacencies between trees, where
• each node has some score or weight: here, it will be likelihood.
  
  As we cannot check every tree we must wander through the set of trees in the hope of finding "good" ones, under the assumption that the best (most likely) trees will have a lot in common with each other.
  
  Fortunately this turns out to be the case! The best trees tend to clump together in tree space.

Tree perturbations 1: NNI

Nearest Neighbour Interchange is a way to transform one tree into another -- so we can search for the best tree.

An NNI move consists of effectively collapsing an internal branch of the tree and re-expanding it out in either of two possible ways.

For example, given the original tree

there are two internal edges: $x-y$ and $y-z$. We can do an NNI move on each of them.

Tree perturbations 1: NNI

Select edge $xy$ to collapse.

Tree perturbations 1: NNI

$xy$ is collapsed.

Tree perturbations 1: NNI

Expand $xy$ in either of two ways.

For an unrooted binary tree with $n$ leaves, there are $2(n-3)$ neighbouring trees under this move.

Simple example: five taxa

Tree perturbations 2: SPR

Subtree Pruning and Regrafting works by choosing an edge, cutting the tree at that edge, and then reconnecting the smaller part anywhere else in the larger part.

Tree perturbations 2: SPR

Trim the root of the subtree...

Tree perturbations 2: SPR

Regraft to $y-D$,

Tree perturbations 2: SPR

... or $y-z$,

Tree perturbations 2: SPR

... or $z-E$,

Tree perturbations 2: SPR

... or $z-F$.

Under this move there are more neighbours: $4(n-3)(n-2)$ possible moves.

Tree perturbations 3: TBR

Tree Bisection and Reconnection works by cutting the tree at any internal edge, then re-rooting each of the parts just formed, and reconnecting the roots:

Tree perturbations 3: TBR

Tree perturbations 3: TBR

Selecting edges as roots of the new subtrees

Tree perturbations 3: TBR

And reconnect!

Note: There is no neat form for the number of neighbours of a given tree under the SPR moves, but it's more than TBR.

Painted Hills near Mitchell, Oregon.

This file is licensed under the Creative Commons Attribution-Share Alike 3.0 Unported license.

Hill Climbing

Tree Search algorithms work essentially by

1. Starting with a reasonably good estimate, e.g., using Neighbour-Joining;

2. Checking all the neighbouring trees under whichever perturbation is being used for an acceptable new tree;

3. If an acceptable new tree is available, move to it! If not, then stop.

   • Most simply, "acceptable" means "better than this one".
     
     This kind of heuristic search is called hill climbing, from the obvious analogy.

This figure from Cao et al. (2019) illustrates moving around in tree space.

Source: https://doi.org/10.1101/746362

Anolis phylogeny example