Introduction

These notes should enable the user to estimate phylogenetic trees from alignment data with different methods using the phangorn package (Schliep 2011) . Several functions of package are also described in more detail in (Paradis 2012). For more theoretical background on all the methods see e.g. [Felsenstein (2004)](Yang 2006). This document illustrates some of the package features to estimate phylogenetic trees using different reconstruction methods. Small adaptations to the scripts in section @ref(appendix) should enable the user to perform phylogenetic analyses.

Getting started

The first thing we have to do is to read in an alignment. Unfortunately there exists many different file formats that alignments can be stored in. The function read.phyDat is used to read in an alignment. There are several functions to read in alignments depending on the format of the data set (“nexus”, “phylip”, “fasta”) and the kind of data (amino acid or nucleotides) in the ape package (Paradis and Schliep 2019) and phangorn. The function read.phyDat calls these other functions and transform them into a phyDat object. For the specific parameter settings available look in the help files of the function read.dna (for phylip, fasta, clustal format), read.nexus.data for nexus files. For amino acid data additional read.aa is called.

We start our analysis loading the phangorn package and then reading in an alignment.

library(ape)
library(phangorn)
fdir <- system.file("extdata/trees", package = "phangorn")
primates <- read.phyDat(file.path(fdir, "primates.dna"),
                        format = "interleaved")

Distance based methods

After reading in the alignment we can build a first tree with distance based methods. The function dist.dna from the ape package computes distances for many DNA substitution models. To use the function dist.dna we have to transform the data to class DNAbin. For amino acids the function dist.ml offers common substitution models (for example “WAG”, “JTT”, “LG”, “Dayhoff”, “cpREV”, “mtmam”, “mtArt”, “MtZoa” or “mtREV24”).

After constructing a distance matrix we reconstruct a rooted tree with UPGMA and alternatively an unrooted tree using Neighbor Joining [Saitou and Nei (1987)](Studier and Keppler 1988). More distance methods like fastme are available in the ape package.

dm  <- dist.ml(primates)
treeUPGMA  <- upgma(dm)
treeNJ  <- NJ(dm)

We can plot the trees treeUPGMA and treeNJ with the commands:

plot(treeUPGMA, main="UPGMA")

Rooted UPGMA tree.

plot(treeNJ, "unrooted", main="NJ")

Unrooted NJ tree.

Bootstrap

To run the bootstrap we need to first write a function which computes a tree from an alignment. So we first need to compute a distance matrix and afterwards compute the tree. This function we can than give to the bootstrap.phyDat function.

fun <- function(x) upgma(dist.ml(x))
bs_upgma <- bootstrap.phyDat(primates,  fun)

With the new syntax of R 4.1 this can be written a bit shorter:

bs_upgma <- bootstrap.phyDat(primates,  \(x){dist.ml(x) |> upgma})

Finally we can plot the tree with bootstrap values added:

plotBS(treeUPGMA, bs_upgma, main="UPGMA")

Rooted UPGMA tree.

Distance based methods are very fast and we will use the UPGMA and NJ tree as starting trees for the maximum parsimony and maximum likelihood analyses.

Parsimony

The function parsimony returns the parsimony score, that is the number of changes which are at least necessary to describe the data for a given tree. We can compare the parsimony score or the two trees we computed so far:

parsimony(treeUPGMA, primates)
## [1] 751
parsimony(treeNJ, primates)
## [1] 746

The function most users want to use infer phylogenies with MP is pratchet, an implementation of the parsimony ratchet (Nixon 1999). This allows to escape local optima and find better trees than performing only NNI / SPR rearrangements.

The current implementation is

  1. Create a bootstrap data set \(D_b\) from the original data set.
  2. Take the current best tree and perform tree rearrangements on \(D_b\) and save bootstrap tree as \(T_b\).
  3. Use \(T_b\) and perform tree rearrangements on the original data set. If this trees has a lower parsimony score than the currently best tree replace it.
  4. Iterate 1:3 until either a given number of iteration is reached or no improvements have been recorded for a number of iterations.
treeRatchet  <- pratchet(primates, trace = 0, minit=100)
parsimony(treeRatchet, primates)
## [1] 746

Here we set the minimum iteration of the parsimony ratchet to 100 iterations. As the ratchet performed implicitly bootstrap resampling we already computed some branch support, in out case with at least 100 bootstrap iterations. The parameter trace=0 tells the function not write the current status to the console. The function may returns several best trees, but these trees have no branch length assigned to them yet. Now let’s do this:

treeRatchet  <- acctran(treeRatchet, primates)

After assigning edge weights we prune away internal edges of length, so our trees may contain multifurcations.

treeRatchet  <- di2multi(treeRatchet)

Some trees might have differed only between edges of length 0.

if(inherits(treeRatchet, "multiPhylo")){
  treeRatchet <- unique(treeRatchet)
}

The parsimony ratchet implicitly performs a bootstrap analysis (step 1). We make use of this and store the trees which where visited. This allows us to add bootstrap support values to the tree.

plotBS(midpoint(treeRatchet), type="phylogram")
add.scale.bar()

If treeRatchet is a list of trees, i.e. an object of class multiPhylo, we can subset the i-th trees with treeRatchet[[i]].

While in most cases pratchet will be enough to use, phangorn exports some function which might be useful. random.addition computes random addition and can be used to generate starting trees. The function optim.parsimony performs tree rearrangements to find trees with a lower parsimony score. The tree rearrangement implemented are nearest-neighbor interchanges (NNI) and subtree pruning and regrafting (SPR). The later one only works so far with the fitch algorithm.

treeRA <- random.addition(primates)
treeSPR  <- optim.parsimony(treeRA, primates)
## Final p-score 746 after  0 nni operations
parsimony(c(treeRA, treeSPR), primates)
## [1] 746 746

Branch and bound

For data sets with few species it is also possible to find all most parsimonious trees using a branch and bound algorithm (Hendy and Penny 1982). For data sets with more than 10 taxa this can take a long time and depends strongly on how tree like the data are. And for more than 20-30 taxa this will take almost forever.

(trees <- bab(primates[1:10,]))
## 1 phylogenetic tree

Maximum likelihood

The last method we will describe in this vignette is Maximum Likelihood (ML) as introduced by Felsenstein (Felsenstein 1981). We can compute the likelihood for a tree given the data

fit <- pml(treeNJ, data=primates)
fit
## model: JC 
## loglikelihood: -3075 
## unconstrained loglikelihood: -1230 
## 
## Rate matrix:
##   a c g t
## a 0 1 1 1
## c 1 0 1 1
## g 1 1 0 1
## t 1 1 1 0
## 
## Base frequencies:  
##    a    c    g    t 
## 0.25 0.25 0.25 0.25

The function pml returns an object of class pml. This object contains the data, the tree and many different parameters of the model like the likelihood. There are many generic functions for the class pml available, which allow the handling of these objects.

methods(class="pml")
## [1] AICc   BIC    anova  logLik plot   print  simSeq update vcov  
## see '?methods' for accessing help and source code

The object fit just estimated the likelihood for the tree it got supplied, but the branch length are not optimized for the Jukes-Cantor (Jukes and Cantor 1969) model yet, which can be done with the function optim.pml.

fitJC  <- optim.pml(fit, rearrangement="NNI")
## optimize edge weights:  -3075 --> -3068 
## optimize edge weights:  -3068 --> -3068 
##  optimize topology:  -3068 --> -3068 
##  NNI moves:  1 
##  optimize topology:  -3068 --> -3068 
##  NNI moves:  0 
## NNI moves:  1 
## optimize edge weights:  -3068 --> -3068 
##  optimize topology:  -3068 --> -3068 
##  NNI moves:  0 
## NNI moves:  0
logLik(fitJC)
## 'log Lik.' -3068 (df=25)

With the default values pml will estimate a Jukes-Cantor model. That means equal base frequencies and all transition rates are equal. The generic function update allows to change parameters manually. This is not what we usually want to do. However we might want to supply a different tree or change the number of rate categories.

fitF81 <- update(fitJC, k=4, inv=0.2, bf=baseFreq(primates))
fitF81
## model: F81+G(4)+I 
## loglikelihood: -3037 
## unconstrained loglikelihood: -1230 
## Proportion of invariant sites: 0.2 
## Discrete gamma model
## Number of rate categories: 4 
## Shape parameter: 1 
## 
## Rate matrix:
##   a c g t
## a 0 1 1 1
## c 1 0 1 1
## g 1 1 0 1
## t 1 1 1 0
## 
## Base frequencies:  
##       a       c       g       t 
## 0.37481 0.40160 0.03911 0.18448

In the line above we changed the th model a (discrete) rate across site model with 4 rate categories (using the default shape parameter of 1), to 0.2 invariant sites and and supply empirical base frequencies.

fitGTR <- optim.pml(fitF81, model="GTR", optInv=TRUE, optGamma=TRUE,
    rearrangement = "NNI", control = pml.control(trace = 0))
fitGTR
## model: GTR+G(4)+I 
## loglikelihood: -2610 
## unconstrained loglikelihood: -1230 
## Proportion of invariant sites: 0.005846 
## Discrete gamma model
## Number of rate categories: 4 
## Shape parameter: 2.89 
## 
## Rate matrix:
##        a         c         g      t
## a  0.000  1.211438 74.564568  1.025
## c  1.211  0.000000  0.006716 31.278
## g 74.565  0.006716  0.000000  1.000
## t  1.025 31.278075  1.000000  0.000
## 
## Base frequencies:  
##       a       c       g       t 
## 0.39401 0.37870 0.04036 0.18692

We will change the model to the GTR + \(\Gamma(4)\) + I model and then optimize all the parameters.

With the control parameters the thresholds for the fitting process can be changed. Here we want just to suppress output during the fitting process. For larger trees the NNI rearrangements often get stuck in a local maximum. We added two stochastic algorithms to improve topology search. The first (set rearrangement="stochastic") performs stochastic rearrangements similar as in (Nguyen et al. 2015), which makes random NNI permutation to the tree, which than gets optimized to escape local optima. The second option (rearrangement=“ratchet”) perform the likelihood ratchet (Vos 2003).

While these algorithms may find better trees they will also take more time.

fitGTR <- optim.pml(fitGTR, model="GTR", optInv=TRUE, optGamma=TRUE,
    rearrangement = "stochastic", control = pml.control(trace = 0))
fitGTR
## model: GTR+G(4)+I 
## loglikelihood: -2608 
## unconstrained loglikelihood: -1230 
## Proportion of invariant sites: 0.005529 
## Discrete gamma model
## Number of rate categories: 4 
## Shape parameter: 2.628 
## 
## Rate matrix:
##         a         c         g      t
## a   0.000  1.782383 1.748e+02  1.873
## c   1.782  0.000000 5.291e-03 67.293
## g 174.791  0.005291 0.000e+00  1.000
## t   1.873 67.293164 1.000e+00  0.000
## 
## Base frequencies:  
##       a       c       g       t 
## 0.39731 0.37827 0.04023 0.18420

Recently we added a new function pml_bb which hides a few steps from the user as long as only some standard phylogentic analysis are performed.

fitGTR <- pml_bb(primates, model="GTR+G(4)+I", control = pml.control(trace = 0))
fitGTR
## model: GTR+G(4)+I 
## loglikelihood: -2608 
## unconstrained loglikelihood: -1230 
## Proportion of invariant sites: 0.005622 
## Discrete gamma model
## Number of rate categories: 4 
## Shape parameter: 2.937 
## 
## Rate matrix:
##         a         c         g       t
## a  0.0000  0.706497 43.157934  0.4464
## c  0.7065  0.000000  0.002709 17.9853
## g 43.1579  0.002709  0.000000  1.0000
## t  0.4464 17.985260  1.000000  0.0000
## 
## Base frequencies:  
##       a       c       g       t 
## 0.39388 0.37904 0.04027 0.18681

In this case it generates a starting tree, an pml object and optimizes the parameters given the model term using optim.pml.

Model selection

Usually we will try to find the best fitting model beforehand. For this we use the function modelTest to compare different nucleotide or protein models the AIC, AICc or BIC, similar to popular programs ModelTest and ProtTest (D. Posada and Crandall 1998), (David Posada 2008), (Abascal, Zardoya, and Posada 2005). By default available nucleotide or amino acid models available are compared.

mt <- modelTest(primates)

Here we can select only some common models:

mt <- modelTest(primates, model=c("JC", "F81", "K80", "HKY", "SYM", "GTR"), 
                control = pml.control(trace = 0))

The results of modelTest is illustrated in following table:

Model df logLik AIC AICw AICc AICcw BIC
JC 25 -3068 6187 0.00 6193 0.00 6273
JC+I 26 -3063 6177 0.00 6184 0.00 6267
JC+G(4) 26 -3067 6186 0.00 6193 0.00 6275
JC+G(4)+I 27 -3063 6179 0.00 6187 0.00 6272
F81 28 -2918 5892 0.00 5900 0.00 5989
F81+I 29 -2909 5876 0.00 5885 0.00 5976
F81+G(4) 29 -2913 5883 0.00 5892 0.00 5983
F81+G(4)+I 30 -2909 5877 0.00 5886 0.00 5980
K80 26 -2953 5958 0.00 5965 0.00 6048
K80+I 27 -2945 5943 0.00 5950 0.00 6036
K80+G(4) 27 -2945 5944 0.00 5951 0.00 6037
K80+G(4)+I 28 -2942 5941 0.00 5949 0.00 6037
HKY 29 -2625 5308 0.00 5317 0.00 5408
HKY+I 30 -2621 5303 0.00 5312 0.00 5406
HKY+G(4) 30 -2613 5285 0.18 5295 0.45 5389
HKY+G(4)+I 31 -2612 5287 0.08 5297 0.14 5394
SYM 30 -2814 5688 0.00 5697 0.00 5791
SYM+I 31 -2812 5685 0.00 5695 0.00 5792
SYM+G(4) 31 -2805 5671 0.00 5681 0.00 5778
SYM+G(4)+I 32 -2805 5673 0.00 5684 0.00 5784
GTR 33 -2619 5303 0.00 5315 0.00 5417
GTR+I 34 -2614 5295 0.00 5307 0.00 5412
GTR+G(4) 34 -2608 5283 0.47 5295 0.29 5401
GTR+G(4)+I 35 -2607 5284 0.27 5297 0.11 5405

To speed computations up the thresholds for the optimization in modelTest are not as strict as for optim.pml and no tree rearrangements are performed, which is the most time consuming part of the optimizing process. As modelTest computes and optimizes a lot of models it would be a waste of computer time not to save these results. The results are saved as call together with the optimized trees in an environment and the function as.pml evaluates this call to get a pml object back to use for further optimization or analysis.

(fit <- as.pml(mt, "HKY+G(4)+I"))
## model: HKY+G(4)+I 
## loglikelihood: -2612 
## unconstrained loglikelihood: -1230 
## Proportion of invariant sites: 0.002716 
## Discrete gamma model
## Number of rate categories: 4 
## Shape parameter: 2.127 
## 
## Rate matrix:
##       a     c     g     t
## a  0.00  1.00 56.01  1.00
## c  1.00  0.00  1.00 56.01
## g 56.01  1.00  0.00  1.00
## t  1.00 56.01  1.00  0.00
## 
## Base frequencies:  
##       a       c       g       t 
## 0.42042 0.36229 0.04394 0.17335
(fit <- as.pml(mt, "BIC"))
## model: HKY+G(4) 
## loglikelihood: -2613 
## unconstrained loglikelihood: -1230 
## Discrete gamma model
## Number of rate categories: 4 
## Shape parameter: 2.145 
## 
## Rate matrix:
##       a     c     g     t
## a  0.00  1.00 51.44  1.00
## c  1.00  0.00  1.00 51.44
## g 51.44  1.00  0.00  1.00
## t  1.00 51.44  1.00  0.00
## 
## Base frequencies:  
##       a       c       g       t 
## 0.41943 0.36287 0.04387 0.17383

To simplify the workflow we can give the result of modelTest to pml_bb function and optimize the parameters taking the best model according to BIC.

fit_mt <- pml_bb(mt, control = pml.control(trace = 0))
fit_mt
## model: HKY+G(4) 
## loglikelihood: -2612 
## unconstrained loglikelihood: -1230 
## Discrete gamma model
## Number of rate categories: 4 
## Shape parameter: 2.07 
## 
## Rate matrix:
##       a     c     g     t
## a  0.00  1.00 56.55  1.00
## c  1.00  0.00  1.00 56.55
## g 56.55  1.00  0.00  1.00
## t  1.00 56.55  1.00  0.00
## 
## Base frequencies:  
##      a      c      g      t 
## 0.4209 0.3623 0.0438 0.1731

We can compare nested models for the JC and GTR + \(\Gamma(4)\) + I model using likelihood ratio statistic

anova(fitJC, fitGTR)
## Likelihood Ratio Test Table
##   Log lik. Df Df change Diff log lik. Pr(>|Chi|)    
## 1    -3068 25                                       
## 2    -2608 35        10           920     <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

with the Shimodaira-Hasegawa test

SH.test(fitGTR, fitJC)
##      Trees  ln L Diff ln L p-value
## [1,]     1 -2608       0.0  0.4993
## [2,]     2 -3068     459.8  0.0000

or with the AIC

AIC(fitJC)
## [1] 6187
AIC(fitGTR)
## [1] 5287
AICc(fitGTR)
## [1] 5300
BIC(fitGTR)
## [1] 5408

Bootstrap

At last we may want to apply bootstrap to test how well the edges of the tree are supported:

bs <- bootstrap.pml(fitJC, bs=100, optNni=TRUE,
    control = pml.control(trace = 0))

Now we can plot the tree with the bootstrap support values on the edges and also look at consensusNet to identify potential conflict.

plotBS(midpoint(fitJC$tree), bs, p = 50, type="p")

Tree with bootstrap support. Unrooted tree (midpoint rooted) with bootstrap support values.

cnet <- consensusNet(bs, p=0.2)
plot(cnet, show.edge.label=TRUE)

ConsensusNet from the bootstrap sample. Several analyses, e.g.bootstrap and modelTest, can be computationally demanding, but as nowadays most computers have several cores one can distribute the computations using the parallel package. However it is only possible to use this approach if R is running from command line (“X11”), but not using a GUI (for example “Aqua” on Macs) and unfortunately the parallel package does not work at all under Windows.

Molecular dating with a strict clock for ultrametric and tipdated phylogenies

We implemented a strict clock as described in (Felsenstein 2004), p. 266 allowing to infer ultrametric and tip-dated phylogenies. We need a starting tree that fulfills the assumptions, so either tree has to be ultrametric or the constraints given by the tip dates.
For a ultrametric starting tree starting tree we can use an UPGMA or WPGMA tree. When we optimize the tree with optim.pml we have to make sure set we set optRooted = TRUE and for tip-dated we need additionally set optRate=TRUE!

fit_strict <- pml_bb(primates, model="HKY+G(4)", method="ultrametric",
                     control = pml.control(trace = 0))
plot(fit_strict)

With phangorn we also can estimate tipdated phylogenies. Here we use as an example a H3N2 virus data set from treetime (Sagulenko, Puller, and Neher 2018). Additionally to the alignment we also need to read in data containing the dates of the tips.

fdir <- system.file("extdata/trees", package = "phangorn")
tmp <- read.csv(file.path(fdir,"H3N2_NA_20.csv"))
H3N2 <- read.phyDat(file.path(fdir,"H3N2_NA_20.fasta"), format="fasta")

We first process the sampling dates and create a named vector. The lubridate package (Grolemund and Wickham 2011) comes very handy dates in case one has to recode dates, e.g. days and months.

dates <- setNames(tmp$numdate_given, tmp$name)
head(dates)
##               A/Hawaii/02/2013|KF789866|05/28/2013|USA|12_13|H3N2/1-1409 
##                                                                     2013 
##         A/Boston/DOA2_107/2012|CY148382|11/01/2012|USA|12_13|H3N2/1-1409 
##                                                                     2013 
##               A/Oregon/15/2009|GQ895004|06/25/2009|USA|08_09|H3N2/1-1409 
##                                                                     2009 
## A/Hong_Kong/H090_695_V10/2009|CY115546|07/10/2009|Hong_Kong||H3N2/8-1416 
##                                                                     2010 
##            A/New_York/182/2000|CY001279|02/18/2000|USA|99_00|H3N2/1-1409 
##                                                                     2000 
##        A/Canterbury/58/2000|CY009150|09/05/2000|New_Zealand||H3N2/8-1416 
##                                                                     2001

Again we use the pml_bb function, which optimizes the tree given th constraints of the tip.dates vector.

fit_td <- pml_bb(H3N2, model="GTR+G(4)", method="tipdated", tip.dates=dates, 
               rearrangement="NNI", control = pml.control(trace = 0))

And we last we plot the tree with a timescale.

tree_td <- fit_td$tree
root_time <- max(dates) - max(node.depth.edgelength(tree_td))
plot(tree_td, show.tip.label = FALSE)
axisPhylo(root.time = root_time, backward = FALSE)

While the loglikelihood is lower than for an unrooted tree, we have to keep in mind that rooted trees use less parameters. In unrooted trees we estimate one edge length parameter for each tree, for ultrametric trees we estimate only a parameter for each internal node and for tipdated trees we have one additional parameter for the rate.

Appendix

Standard scripts for nucleotide analysis

Here we provide two standard scripts which can be adapted for the most common tasks. Most likely the arguments for read.phyDat have to be adapted to accommodate your file format. Both scripts assume that the parallel package works on your platform, see comments above.

library(phangorn)
file <- "myfile"
dat <- read.phyDat(file, format="fasta")
dm <- dist.ml(dat, "F81")
tree <- NJ(dm)
# as alternative for a starting tree:
tree <- pratchet(dat)          # parsimony tree
tree <- nnls.phylo(tree, dm)   # need edge weights

# 1. alternative: quick and dirty: GTR + G
(fit <- pml_bb(dat, model="GTR+G"))

# 2. alternative: choose with modelTest
mt <- modelTest(dat, multicore=TRUE)
mt[order(mt$BIC),]
# chooses best model from the table according to BIC (default)
fit <- pml_bb(mt)
bs <- bootstrap.pml(fit, bs=100, optNni=TRUE, multicore=TRUE)

tree_ml <- plotBS(fit$tree, bs)
# export the tree 
write.tree(tree_ml, file="tree_ml.nwk")

Standard scripts for amino acid analysis

You can specify different several models build in which you can specify, e.g. “WAG”, “JTT”, “Dayhoff”, “LG”. Optimizing the rate matrix for amino acids is possible, but would take a long, a very long time and you will need to have a large alignment to estimate all the parameters. So make sure to set optBf=FALSE and optQ=FALSE in the function optim.pml, which is also the default.

library(phangorn)
file <- "myfile"
dat <- read.phyDat(file, type = "AA", format="fasta")

# compute a neighbor joining tree
dm <- dist.ml(dat, model="JTT")
tree <- NJ(dm)

# parallel will only work safely from command line and not at all windows
# Select specific models. 
(mt <- modelTest(dat, model=c("JTT", "LG", "WAG"), multicore=TRUE))
# run all available amino acid models
(mt <- modelTest(dat, model="all", multicore=TRUE))

(fit <- pml_bb(mt))
# non-parametric bootstrap
bs <- bootstrap.pml(fit, bs=100, optNni=TRUE, multicore=TRUE)

tree_ml <- plotBS(fit$tree, bs)
# export the tree 
write.tree(tree_ml, file="tree_ml.nwk")

Session info

## R version 4.2.1 (2022-06-23)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.1 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
## 
## locale:
##  [1] LC_CTYPE=de_AT.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=de_AT.UTF-8        LC_COLLATE=C              
##  [5] LC_MONETARY=de_AT.UTF-8    LC_MESSAGES=de_AT.UTF-8   
##  [7] LC_PAPER=de_AT.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=de_AT.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] knitr_1.40      phangorn_2.10.0 ape_5.6-3      
## 
## loaded via a namespace (and not attached):
##  [1] igraph_1.3.4     Rcpp_1.0.9       rstudioapi_0.14  magrittr_2.0.3  
##  [5] lattice_0.20-45  R6_2.5.1         quadprog_1.5-8   rlang_1.0.5     
##  [9] fastmatch_1.1-3  fastmap_1.1.0    highr_0.9        stringr_1.4.1   
## [13] tools_4.2.1      parallel_4.2.1   grid_4.2.1       nlme_3.1-159    
## [17] xfun_0.33        cli_3.4.0        jquerylib_0.1.4  htmltools_0.5.3 
## [21] yaml_2.3.5       digest_0.6.29    Matrix_1.5-1     codetools_0.2-18
## [25] sass_0.4.2       prettydoc_0.4.1  cachem_1.0.6     evaluate_0.16   
## [29] rmarkdown_2.16   stringi_1.7.8    compiler_4.2.1   bslib_0.4.0     
## [33] generics_0.1.3   jsonlite_1.8.0   pkgconfig_2.0.3

References

Abascal, Federico, Rafael Zardoya, and David Posada. 2005. “ProtTest: Selection of Best-Fit Models of Protein Evolution.” Bioinformatics 21 (9): 2104–5. https://doi.org/10.1093/bioinformatics/bti263.
Felsenstein, Joseph. 1981. “Evolutionary Trees from DNA Sequences: A Maxumum Likelihood Approach.” Journal of Molecular Evolution 17: 368–76.
———. 2004. Inferring Phylogenies. Sunderland: Sinauer Associates.
Grolemund, Garrett, and Hadley Wickham. 2011. “Dates and Times Made Easy with lubridate.” Journal of Statistical Software 40 (3): 1–25. https://www.jstatsoft.org/v40/i03/.
Hendy, M. D., and D. Penny. 1982. “Branch and Bound Algorithms to Determine Minimal Evolutionary Trees.” Math. Biosc. 59: 277–90.
Jukes, Thomas H., and Charles R. Cantor. 1969. “{CHAPTER} 24 - Evolution of Protein Molecules.” In Mammalian Protein Metabolism, edited by H. N. Munro, 21–132. Academic Press.
Nguyen, Lam-Tung, Heiko A. Schmidt, Arndt von Haeseler, and Bui Quang Minh. 2015. “IQ-TREE: A Fast and Effective Stochastic Algorithm for Estimating Maximum-Likelihood Phylogenies.” Molecular Biology and Evolution 32 (1): 268–74. https://doi.org/10.1093/molbev/msu300.
Nixon, K. 1999. “The Parsimony Ratchet, a New Method for Rapid Rarsimony Analysis.” Cladistics 15: 407–14.
Paradis, Emmanuel. 2012. Analysis of Phylogenetics and Evolution with r. Second. New York: Springer.
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.
Posada, David. 2008. jModelTest: Phylogenetic Model Averaging.” Molecular Biology and Evolution 25 (7): 1253–56. https://doi.org/10.1093/molbev/msn083.
Posada, D., and K. A. Crandall. 1998. MODELTEST: Testing the Model of DNA Substitution.” Bioinformatics 14 (9): 817–18.
Sagulenko, Pavel, Vadim Puller, and Richard A Neher. 2018. “TreeTime: Maximum-Likelihood Phylodynamic Analysis.” Virus Evolution 4 (1): vex042.
Saitou, N., and M. Nei. 1987. “The Neighbor-Joining Method - a New Method for Reconstructing Phylogenetic Trees.” Molecular Biology and Evolution 4 (4): 406–25.
Schliep, Klaus Peter. 2011. “Phangorn: Phylogenetic Analysis in R.” Bioinformatics 27 (4): 592–93. https://doi.org/10.1093/bioinformatics/btq706.
Studier, J. A., and K. J. Keppler. 1988. “A Note on the Neighbor-Joining Algorithm of Saitou and Nei.” Molecular Biology and Evolution 5 (6): 729–31.
Vos, R. A. 2003. “Accelerated Likelihood Surface Exploration: The Likelihood Ratchet.” Systematic Biology 52 (3): 368–73.
Yang, Ziheng. 2006. Computational Molecular Evolution. Oxford: Oxford University Press.