Felsenstein's tree-pruning algorithm

In statistical genetics, Felsenstein's tree-pruning algorithm (or Felsenstein's tree-peeling algorithm), attributed to Joseph Felsenstein, is an algorithm for efficiently computing the likelihood of an evolutionary tree from nucleic acid sequence data. ^{[1] [2]}

The algorithm is often used as a subroutine in a search for a maximum likelihood estimate for an evolutionary tree. Further, it can be used in a hypothesis test for whether evolutionary rates are constant (by using likelihood ratio tests). It can also be used to provide error estimates for the parameters describing an evolutionary tree.

Details

A simple phylogenetic tree example made from arbitrary data D

The likelihood of a tree ${\displaystyle T}$ is, by definition, the probability of observing certain data ${\displaystyle D}$ (${\displaystyle D}$ being a nucleotide sequence alignment for example i.e. a succession of ${\displaystyle n}$ DNA site ${\displaystyle s}$) given the tree. It is often written : ${\displaystyle P(D|T)}$.

Here is an example of an evolutionary tree on arbitrary sequence data ${\displaystyle D}$:

This is a key value and is often quite complicated to compute. To ease the computations, Felsenstein and his colleagues used several assumptions that are still widely used today. The main assumption is that mutations between DNA sites are independent of each other. This permits to compute the likelihood as a simple product of probabilities. Now you can divide the data ${\displaystyle D}$ between several ${\displaystyle D_{s}}$ for each nucleotide site ${\displaystyle s}$ inside of ${\displaystyle D}$. The global likelihood of the tree will be the product of the likelihoods of each site:

Same tree but made from D1, which consists in the first DNA sites from D

${\displaystyle P(D|T)=\prod _{s=1}^{n}{P(D_{s}|T)}}$

If I reuse the example above, ${\displaystyle D_{1}}$ tree would be:

The second assumption concerns the models of DNA sequence evolution. The computation of the likelihood needs the nucleotide frequencies as well as the transition probabilities (when a mutation occurs, probability of going from a nucleotide to another). The simplest model is the Jukes-Cantor model, assuming equal nucleotide frequencies ${\displaystyle \left(\pi _{A}=\pi _{G}=\pi _{C}=\pi _{T}={1 \over 4}\right)}$ and equal transition probabilities from ${\displaystyle X}$ to ${\displaystyle Y}$ (${\displaystyle p_{A\rightarrow T}=p_{A\rightarrow C}=p_{A\rightarrow G}={\frac {1}{4}}(1-e^{-\mu l})}$ and ${\displaystyle p_{A\rightarrow A}=e^{-\mu l}+{\frac {1}{4}}(1-e^{-\mu l})}$) and idem for the other bases. Here ${\displaystyle \mu }$ is the global mutation rate of the model.

Felsenstein proposed to decomposed computations even more by using "partial likelihoods" in the computation of ${\displaystyle P(D_{s}|T)}$. Here is how it works.

Assume we are on a node ${\displaystyle k}$ on the tree ${\displaystyle T}$. ${\displaystyle k}$ has two daughter nodes ${\displaystyle i}$ and ${\displaystyle j}$ and, for each DNA base ${\displaystyle X=\{A,T,C,G\}}$ present on ${\displaystyle k}$, we can define a partial likelihood ${\displaystyle w_{k}(X)}$ such as:

${\displaystyle w_{k}(X)=(\sum _{Y}p_{X\rightarrow Y}\centerdot w_{i}(Y))\centerdot (\sum _{Z}p_{X\rightarrow Z}\centerdot w_{j}(Z))}$

where ${\displaystyle Y}$ and ${\displaystyle Z}$ are also DNA bases. ${\displaystyle p_{X\rightarrow Y}}$ is the transition probability from nucleotide ${\displaystyle X}$ to nucleotide ${\displaystyle Y}$ (idem for ${\displaystyle p_{X\rightarrow Z}}$). ${\displaystyle w_{i}(Y)}$ is the partial likelihood of the daughter node ${\displaystyle i}$, evaluated on nucleotide ${\displaystyle Y}$ (idem for ${\displaystyle w_{j}(Z)}$).

Using this formula, one has to start from the tips of the tree ${\displaystyle T}$, then move towards the root and compute the partial likelihoods of each necessary node on the way (4 partial likelihoods per node). Having finished at the root of the tree, the likelihood of the tree (for this particular site) is then the sum of the partial likelihoods of the root times the appropriated nucleotide frequency.

${\displaystyle P(D_{s}|T)=\sum _{X}p_{X}\centerdot w_{r}(X)}$

After doing so for every site ${\displaystyle s}$, one can finally obtain the likelihood of the global evolutionary tree by multiplying each "sublikelihood".

References

1. Felsenstein, J. (1973). "Maximum Likelihood and Minimum-Steps Methods for Estimating Evolutionary Trees from Data on Discrete Characters". Systematic Biology. 22 (3): 240–249. doi:10.1093/sysbio/22.3.240.
2. Felsenstein, J. (1981). "Evolutionary trees from DNA sequences: A maximum likelihood approach". Journal of Molecular Evolution. 17 (6): 368–376. Bibcode:1981JMolE..17..368F. doi:10.1007/BF01734359. PMID   7288891. S2CID   8024924.