Skip to contents

Introduction

These notes describe the ancestral sequence reconstruction using the phangorn package (Schliep 2011). phangorn provides several methods to estimate ancestral character states with either Maximum Parsimony (MP) or Maximum Likelihood (ML). For more background on all the methods see e.g. (Felsenstein 2004) or (Yang 2006).

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. “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). The function acctran (accelerated transformation) assigns edge length and internal nodes to the tree (Swofford and Maddison 1987).

anc.acctran <- ancestral.pars(tree, primates, "ACCTRAN")
anc.mpr <- ancestral.pars(tree, primates, "MPR")

All the ancestral reconstructions for parsimony are based on the fitch algorithm and so far only bifurcating trees are allowed. However trees can get pruned afterwards using the function multi2di from ape.

The seqLogo function from the seqLogo package from Bioconductor provides a neat way to show proportions of a nucleotides of ancestral states (see figure 1).

library(seqLogo)
seqLogo( t(subset(anc.mpr, getRoot(tree), 1:20)[[1]]), ic.scale=FALSE)

You may need to install seqLogo before

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("seqLogo")
plotAnc(tree, anc.mpr, 17)
title("MPR")
Fig 2. Ancestral reconstruction using MPR.

Fig 2. Ancestral reconstruction using MPR.

plotAnc(tree, anc.acctran, 17)
title("ACCTRAN")
Fig 3. Ancestral reconstruction using ACCTRAN.

Fig 3. Ancestral reconstruction using ACCTRAN.

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. So far only a marginal construction is implemented (see [Yang (2006)](Koshi and Goldstein 1996)) and no joint reconstruction (Pupko et al. 2000).

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.

anc.ml <- ancestral.pml(fit, "ml")
anc.bayes <- ancestral.pml(fit, "bayes")

The differences of the two approaches for a specific site (17) are represented in the following figures.

plotAnc(tree, anc.ml, 17)
title("ML")
Fig 4. Ancestral reconstruction the using the maximum likelihood.

Fig 4. Ancestral reconstruction the using the maximum likelihood.

plotAnc(tree, anc.bayes, 17)
title("Bayes")
Fig 5. Ancestral reconstruction using (empirical) Bayes.

Fig 5. Ancestral reconstruction using (empirical) 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

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

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 produced

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.

fit_SYM$logLik
## [1] -15.31
fit_ace$loglik+log(1/3)
## [1] -15.31
all.equal(fit_SYM$logLik, fit_ace$loglik+log(1/3))
## [1] "Mean relative difference: 1.229e-07"
anc_SYM <- ancestral.pml(fit_SYM, "ml")
plotAnc(bird.orders, 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.

Session info

## R version 4.3.3 (2024-02-29)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.4 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so;  LAPACK version 3.10.0
## 
## locale:
##  [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8       
##  [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
##  [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
## [10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   
## 
## time zone: UTC
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] phangorn_3.0.0.0 ape_5.7-1.9     
## 
## loaded via a namespace (and not attached):
##  [1] vctrs_0.6.5       nlme_3.1-164      cli_3.6.2         knitr_1.45       
##  [5] rlang_1.1.3       xfun_0.43         highr_0.10        purrr_1.0.2      
##  [9] generics_0.1.3    textshaping_0.3.7 jsonlite_1.8.8    htmltools_0.5.8  
## [13] ragg_1.3.0        sass_0.4.9        rmarkdown_2.26    quadprog_1.5-8   
## [17] grid_4.3.3        evaluate_0.23     jquerylib_0.1.4   fastmap_1.1.1    
## [21] yaml_2.3.8        lifecycle_1.0.4   memoise_2.0.1     compiler_4.3.3   
## [25] codetools_0.2-19  igraph_2.0.3      fs_1.6.3          pkgconfig_2.0.3  
## [29] fastmatch_1.1-4   Rcpp_1.0.12       lattice_0.22-5    systemfonts_1.0.6
## [33] digest_0.6.35     R6_2.5.1          parallel_4.3.3    magrittr_2.0.3   
## [37] Matrix_1.6-5      bslib_0.7.0       tools_4.3.3       pkgdown_2.0.7    
## [41] cachem_1.0.8      desc_1.4.3

References

Felsenstein, Joseph. 2004. Inferring Phylogenies. Sunderland: Sinauer Associates.
Koshi, Jeffrey M., and Richard A. Goldstein. 1996. Probabilistic reconstruction of ancestral protein sequences.” Journal of Molecular Evolution 42 (2): 313–20. https://doi.org/10.1007/BF02198858.
Paradis, Emmanuel, and Klaus Schliep. 2019. “Ape 5.0: An Environment for Modern Phylogenetics and Evolutionary Analyses in r.” Bioinformatics 35 (3): 526–28. https://doi.org/10.1093/bioinformatics/bty633.
Pennell, M. W., J. M. Eastman, G. J. Slater, J. W. Brown, J. C. Uyeda, R. G. Fitzjohn, M. E. Alfaro, and L. J. Harmon. 2014. “Geiger V2.0: An Expanded Suite of Methods for Fitting Macroevolutionary Models to Phylogenetic Trees.” Bioinformatics 30: 2216–18. https://doi.org/10.1093/bioinformatics/btu181.
Pupko, Tal, Itsik Pe, Ron Shamir, and Dan Graur. 2000. A Fast Algorithm for Joint Reconstruction of Ancestral Amino Acid Sequences.” Molecular Biology and Evolution 17 (6): 890–96. https://doi.org/10.1093/oxfordjournals.molbev.a026369.
Revell, Liam J. 2012. “Phytools: An r Package for Phylogenetic Comparative Biology (and Other Things).” Methods in Ecology and Evolution 3: 217–23.
Schliep, Klaus Peter. 2011. “Phangorn: Phylogenetic Analysis in R.” Bioinformatics 27 (4): 592–93. https://doi.org/10.1093/bioinformatics/btq706.
Swofford, D. L., and W. P. Maddison. 1987. “Reconstructing Ancestral Character States Under Wagner Parsimony.” Math. Biosci. 87: 199–229.
Yang, Ziheng. 2006. Computational Molecular Evolution. Oxford: Oxford University Press.