Parsimony reconstructions
To reconstruct ancestral sequences we first load some data and
reconstruct a tree:
library(phangorn)
fdir <- system.file("extdata/trees", package = "phangorn")
primates <- read.phyDat(file.path(fdir, "primates.dna"),
format = "interleaved")
tree <- pratchet(primates, trace=0) |> acctran(primates) |> makeNodeLabel()
parsimony(tree, primates)
## [1] 746
For parsimony analysis of the edge length represent the observed
number of changes. Reconstructing ancestral states therefore defines
also the edge lengths of a tree. However there can exist several equally
parsimonious reconstructions or states can be ambiguous and therefore
edge length can differ. In phangorn all the ancestral
reconstructions for parsimony used to be based on the fitch algorithm
(below version 3.0) and needed bifurcating trees. However trees can get
pruned afterwards using the function multi2di from
ape. Recently we replaced the acctran routine with a method
based on the sankoff algorithm adopting the algorithm for joint
reconstruction (Pupko et al. 2000) and
breaking ties at random. This has the additional benefit that to infer
phylogenies with multifurcations.
“MPR” reconstructs the ancestral states for each (internal) node as
if the tree would be rooted in that node. However the nodes are not
independent of each other. If one chooses one state for a specific node,
this can restrict the choice of neighboring nodes (figures 2 and 3).
There is also an option “POSTORDER” which is only a one pass algotrithm,
which we use for teaching purposes. The function acctran (accelerated
transformation) assigns edge length and internal nodes to the tree (Swofford and Maddison 1987).
anc.pars <- anc_pars(tree, primates)
The plotSeqLogo function is a wrapper around the from
the ggseqlogo function in the ggseqlogo package (Wagih 2024) and provides a simple way to show
proportions of a nucleotides of ancestral states (see figure 1).
plotSeqLogo(anc.pars, node=getRoot(tree), 1, 20)
Fig 1. Ancestral reconstruction for a
node.
plotAnc(anc.pars, 17)
title("MPR")
Fig 2. Ancestral reconstruction using MPR.
Likelihood reconstructions
phangorn also offers the possibility to estimate ancestral
states using ML. The advantages of ML over parsimony is that the
reconstruction accounts for different edge lengths. Currently a marginal
construction is implemented (see [Yang
(2006)](Koshi and Goldstein 1996)),
but the joint reconstruction (Pupko et al.
2000) only for models without rate variation (e.g. gamma models)
or invariant sites.
fit <- pml(tree, primates)
fit <- optim.pml(fit, model="F81", control = pml.control(trace=0))
We can assign the ancestral states according to the highest
likelihood (“ml”): \[
P(x_r = A) = \frac{L(x_r=A)}{\sum_{k \in \{A,C,G,T\}}L(x_r=k)}
\] and the highest posterior probability (“bayes”) criterion:
\[
P(x_r=A) = \frac{\pi_A L(x_r=A)}{\sum_{k \in \{A,C,G,T\}}\pi_k
L(x_r=k)},
\] where \(L(x_r)\) is the joint
probability of states at the tips and the state at the root \(x_r\) and \(\pi_i\) are the estimated base frequencies
of state \(i\). Both methods agree if
all states (base frequencies) have equal probabilities.
The differences of the two approaches for a specific site (17) are
represented in the following figures.
plotAnc(anc.ml, 17)
title("ML")
Fig 4. Ancestral reconstruction the using the
maximum likelihood.
#plotAnc(anc.bayes, 17)
#title("Bayes")
Fitting for discrete comparative data
Often have already a phylogeny and only want estimate the ancestral
reconstruction for this tree. This is a common problem in phylogentic
comparative methods and we can use the function ace in the ape
(Paradis and Schliep 2019),
fitDiscrete in the geiger (Pennell et
al. 2014) or fitMK in the phytools (Revell 2012) package. Here we want to show how
to fit these models using optim.pml.
First we load a tree and create some data.
data("bird.orders")
x <- c(rep(0, 5), rep(1, 18))
x[c(20,22,23)] <- 2
x <- factor(x)
names(x) <- bird.orders$tip.label
dat <- phyDat(x, "USER", levels=c(0,1,2))
We than set up the pml object and optimize the model.
Instead of optimizing the edge length we only optimize the rate.
fit <- pml(bird.orders, dat)
fit_ER <- optim.pml(fit, optEdge = FALSE, optRate=TRUE,
control = pml.control(trace=0))
fit_ER
## model: Mk
## loglikelihood: -16.47
## unconstrained loglikelihood: 0
##
## Rate matrix:
## 0 1 2
## 0 0 1 1
## 1 1 0 1
## 2 1 1 0
##
## Base frequencies:
## 0 1 2
## 0.3333 0.3333 0.3333
##
## Rate: 0.007846
We can also fit the symmetric (model=“SYM”) or ordered
metristic model (model=“ORDERED”).
fit_SYM <- optim.pml(fit, optEdge = FALSE, optRate=TRUE, model="SYM",
control = pml.control(trace=0))
fit_SYM
## model: SYM
## loglikelihood: -15.31
## unconstrained loglikelihood: 0
##
## Rate matrix:
## 0 1 2
## 0 0.000e+00 0.2747 1.604e-06
## 1 2.747e-01 0.0000 1.000e+00
## 2 1.604e-06 1.0000 0.000e+00
##
## Base frequencies:
## 0 1 2
## 0.3333 0.3333 0.3333
##
## Rate: 0.00678
We can compare the estimate with the one from ace from
ape.
fit_ace <- ace(x, bird.orders, model="SYM", type = "d")
## Warning in sqrt(diag(solve(h))): NaNs wurden erzeugt
The log-likelihood values differ slightly as in phangorn the values
get multiplied by the state frequencies. Thus if we add
log(1/3) as we have three states to ace estimate the two
estimates are almost identical.
## [1] -15.31
## [1] -15.31
all.equal(fit_SYM$logLik, fit_ace$loglik+log(1/3))
## [1] "Mean relative difference: 1.229e-07"
anc_SYM <- anc_pml(fit_SYM)
plotAnc(anc_SYM)
More complicated models can be applied using defining the rate matrix
as shown in the vignette Markov models and transition rate
matrices. The “ARD” model is currently not available as
phangorn only fits reversible models.