# 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)
<- system.file("extdata/trees", package = "phangorn")
fdir <- read.phyDat(file.path(fdir, "primates.dna"),
primates 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.

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

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

`plot(treeUPGMA, main="UPGMA")`

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

## 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.

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

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

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

Finally we can plot the tree with bootstrap values added:

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

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

- Create a bootstrap data set \(D_b\) from the original data set.
- Take the current best tree and perform tree rearrangements on \(D_b\) and save bootstrap tree as \(T_b\).
- 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.
- Iterate 1:3 until either a given number of iteration is reached or no improvements have been recorded for a number of iterations.

```
<- pratchet(primates, trace = 0, minit=100)
treeRatchet 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:

`<- acctran(treeRatchet, primates) treeRatchet `

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

`<- di2multi(treeRatchet) treeRatchet `

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

```
if(inherits(treeRatchet, "multiPhylo")){
<- unique(treeRatchet)
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.

```
<- random.addition(primates)
treeRA <- optim.parsimony(treeRA, primates) treeSPR
```

`## 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.

`<- bab(primates[1:10,])) (trees `

`## 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

```
<- pml(treeNJ, data=primates)
fit 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`

.

`<- optim.pml(fit, rearrangement="NNI") fitJC `

```
## 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.

```
<- update(fitJC, k=4, inv=0.2, bf=baseFreq(primates))
fitF81 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.

```
<- optim.pml(fitF81, model="GTR", optInv=TRUE, optGamma=TRUE,
fitGTR 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.

```
<- optim.pml(fitGTR, model="GTR", optInv=TRUE, optGamma=TRUE,
fitGTR 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.

```
<- pml_bb(primates, model="GTR+G(4)+I", control = pml.control(trace = 0))
fitGTR 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.

`<- modelTest(primates) mt `

Here we can select only some common models:

```
<- modelTest(primates, model=c("JC", "F81", "K80", "HKY", "SYM", "GTR"),
mt 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.

`<- as.pml(mt, "HKY+G(4)+I")) (fit `

```
## 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
```

`<- as.pml(mt, "BIC")) (fit `

```
## 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.

```
<- pml_bb(mt, control = pml.control(trace = 0))
fit_mt 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:

```
<- bootstrap.pml(fitJC, bs=100, optNni=TRUE,
bs 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")`

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

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`

!

```
<- pml_bb(primates, model="HKY+G(4)", method="ultrametric",
fit_strict 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.

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

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.

```
<- setNames(tmp$numdate_given, tmp$name)
dates 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.

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

And we last we plot the tree with a timescale.

```
<- fit_td$tree
tree_td <- max(dates) - max(node.depth.edgelength(tree_td))
root_time 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)
<- "myfile"
file <- read.phyDat(file, format="fasta")
dat <- dist.ml(dat, "F81")
dm <- NJ(dm)
tree # as alternative for a starting tree:
<- pratchet(dat) # parsimony tree
tree <- nnls.phylo(tree, dm) # need edge weights
tree
# 1. alternative: quick and dirty: GTR + G
<- pml_bb(dat, model="GTR+G"))
(fit
# 2. alternative: choose with modelTest
<- modelTest(dat, multicore=TRUE)
mt order(mt$BIC),]
mt[# chooses best model from the table according to BIC (default)
<- pml_bb(mt)
fit <- bootstrap.pml(fit, bs=100, optNni=TRUE, multicore=TRUE)
bs
<- plotBS(fit$tree, bs)
tree_ml # 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)
<- "myfile"
file <- read.phyDat(file, type = "AA", format="fasta")
dat
# compute a neighbor joining tree
<- dist.ml(dat, model="JTT")
dm <- NJ(dm)
tree
# parallel will only work safely from command line and not at all windows
# Select specific models.
<- modelTest(dat, model=c("JTT", "LG", "WAG"), multicore=TRUE))
(mt # run all available amino acid models
<- modelTest(dat, model="all", multicore=TRUE))
(mt
<- pml_bb(mt))
(fit # non-parametric bootstrap
<- bootstrap.pml(fit, bs=100, optNni=TRUE, multicore=TRUE)
bs
<- plotBS(fit$tree, bs)
tree_ml # 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

*Bioinformatics*21 (9): 2104–5. https://doi.org/10.1093/bioinformatics/bti263.

*Journal of Molecular Evolution*17: 368–76.

*Inferring Phylogenies*. Sunderland: Sinauer Associates.

*Journal of Statistical Software*40 (3): 1–25. https://www.jstatsoft.org/v40/i03/.

*Math. Biosc.*59: 277–90.

*Mammalian Protein Metabolism*, edited by H. N. Munro, 21–132. Academic Press.

*Molecular Biology and Evolution*32 (1): 268–74. https://doi.org/10.1093/molbev/msu300.

*Cladistics*15: 407–14.

*Analysis of Phylogenetics and Evolution with r*. Second. New York: Springer.

*Bioinformatics*35 (3): 526–28. https://doi.org/10.1093/bioinformatics/bty633.

*Molecular Biology and Evolution*25 (7): 1253–56. https://doi.org/10.1093/molbev/msn083.

*Bioinformatics*14 (9): 817–18.

*Virus Evolution*4 (1): vex042.

*Molecular Biology and Evolution*4 (4): 406–25.

*Bioinformatics*27 (4): 592–93. https://doi.org/10.1093/bioinformatics/btq706.

*Molecular Biology and Evolution*5 (6): 729–31.

*Systematic Biology*52 (3): 368–73.

*Computational Molecular Evolution*. Oxford: Oxford University Press.