The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.
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 this 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’s features to estimate
phylogenetic trees using different reconstruction methods.
The first thing we have to do is to read in an alignment.
Unfortunately there exist many different file formats that alignments
can be stored in. In most cases, the function read.phyDat
is used to read in an alignment. In the ape package (Paradis and Schliep 2019) and
phangorn, 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, nucleotides, morphological data). The
function read.phyDat
calls these other functions and
transforms 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. Morphological data will be
shown later in the vignette Phylogenetic trees from morphological
data.
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")
After reading in the nucleotide 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, but to use the function dist.dna
, we have to
transform the data to class DNAbin. The function dist.ml
from phangorn offers the substitution models “JC69” and “F81”
for DNA, and also common substitution models for amino acids
(e.g. “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")
Rooted UPGMA tree.
plot(treeNJ, "unrooted", main="NJ")
Unrooted NJ tree.
To run the bootstrap we first need to write a function which computes
a tree from an alignment. So we first need to compute a distance matrix
and afterwards compute the tree. We can then give this function 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")
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.
The function parsimony returns the parsimony score, that is the minimum number of changes necessary to describe the data for a given tree. We can compare the parsimony score for the two trees we computed so far:
parsimony(treeUPGMA, primates)
## [1] 751
parsimony(treeNJ, primates)
## [1] 746
The function most users want to use to infer phylogenies with MP
(maximum parsimony) is pratchet
, an implementation of the
parsimony ratchet (Nixon 1999). This
allows to escape local optima and find better trees than only performing
NNI / SPR rearrangements.
The current implementation is
minit
) or no improvements have been recorded for a number
of iterations (k
).<- pratchet(primates, trace = 0, minit=100)
treeRatchet parsimony(treeRatchet, primates)
## [1] 746
Here we set the minimum iteration of the parsimony ratchet
(minit
) to 100 iterations, the default number for
k
is 10. As the ratchet implicitly performs bootstrap
resampling, we already computed some branch support, in our 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 return 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
tol
(default = 1e-08), 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 }
As mentioned above, 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 rearrangements implemented are nearest-neighbor interchanges
(NNI) and subtree pruning and regrafting (SPR). The latter so far only
works 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
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 is. And for more than 20-30 taxa this will take almost forever.
<- bab(primates[1:10,])) (trees
## 1 phylogenetic tree
The last method we will describe in this vignette is Maximum Likelihood (ML) as introduced by Felsenstein (Felsenstein 1981).
Usually, as a first step, we will try to find the best fitting model.
For this we use the function modelTest
to compare different
nucleotide or protein models with 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 are compared.
The Vignette Markov models and transition rate matrices gives further background on those models, how they are estimated and how you can work with them.
<- modelTest(primates) mt
It’s also possible to only select 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 | 5316 | 0.00 | 5408 |
HKY+I | 30 | -2621 | 5302 | 0.00 | 5311 | 0.00 | 5406 |
HKY+G(4) | 30 | -2613 | 5285 | 0.18 | 5294 | 0.46 | 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 | -2618 | 5303 | 0.00 | 5314 | 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 | 5400 |
GTR+G(4)+I | 35 | -2607 | 5284 | 0.27 | 5297 | 0.11 | 5405 |
To speed computations up the thresholds for the optimizations in
modelTest
are not as strict as for optim.pml
(shown in the coming vignettes) 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. This can either be done for a specific model, or for a
specific criterion.
<- as.pml(mt, "HKY+G(4)+I")
fit <- as.pml(mt, "BIC") fit
To simplify the workflow, we can give the result of
modelTest
to the function pml_bb
and optimize
the parameters taking the best model according to BIC. Ultrafast
bootstrapping (Minh, Nguyen, and Haeseler
2013) is conducted automatically if the default
rearrangements="stochastic"
is used. If
rearrangements="NNI"
is used, no bootstrapping is
conducted.
<- 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.1
##
## Rate matrix:
## a c g t
## a 0 1 57 1
## c 1 0 1 57
## g 57 1 0 1
## t 1 57 1 0
##
## Base frequencies:
## a c g t
## 0.421 0.362 0.044 0.173
We can also use pml_bb
with a defined model to infer a
phylogenetic tree.
<- pml_bb(primates, model="GTR+G(4)+I") fitGTR
If we instead want to conduct standard bootstrapping (Felsenstein 1985; Penny and Hendy 1985), we can
do so with the function bootstrap.pml
:
<- bootstrap.pml(fit_mt, 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 compare the standard bootstrap values to the ultrafast
bootstrap values. With the function plotBS
it is not only
possible to plot these two, but also the transfer bootstraps (Lemoine et al. 2018) which are especially
useful for large data sets.
plotBS(midpoint(fit_mt$tree), p = .5, type="p", digits=2, main="Ultrafast bootstrap")
plotBS(midpoint(fit_mt$tree), bs, p = 50, type="p", main="Standard bootstrap")
plotBS(midpoint(fit_mt$tree), bs, p = 50, type="p", digits=0, method = "TBE", main="Transfer bootstrap")
Unrooted tree (midpoint rooted) with ultrafast, standard and transfer bootstrap support values.
If we want to assign the standard or transfer bootstrap values to the
node labels in our tree instead of plotting it (e.g. to export the tree
somewhere else), plotBS
gives that option with
type = "n"
:
# assigning standard bootstrap values to our tree; this is the default method
<- plotBS(fit_mt$tree, bs, type = "n")
tree_stdbs
# assigning transfer bootstrap values to our tree
<- plotBS(fit_mt$tree, bs, type = "n", method = "TBE") tree_tfbs
It is also possible to look at consensusNet
to identify
potential conflict.
<- consensusNet(bs, p=0.2)
cnet plot(cnet, show.edge.label=TRUE)
ConsensusNet from the standard 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.
Now that we have our tree with bootstrap values, we can easily write it to a file in Newick-format:
# tree with ultrafast bootstrap values
write.tree(fit_mt$tree, "primates.tree")
# tree with standard bootstrap values
write.tree(tree_stdbs, "primates.tree")
# tree with transfer bootstrap values
write.tree(tree_tfbs, "primates.tree")
When we assume a “molecular clock” phylogenies can be used to infer
divergence times (Zuckerkandl and Pauling
1965). 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 the tree has to be ultrametric, or
the constraints given by the tip dates.
For an ultrametric starting tree we can use an UPGMA or WPGMA tree.
<- pml_bb(primates, model="HKY+G(4)", method="ultrametric",
fit_strict rearrangement="NNI", control = pml.control(trace = 0))
plot(fit_strict)
With phangorn we also can estimate tipdated phylogenies. Here we use a H3N2 virus data set from treetime (Sagulenko, Puller, and Neher 2018) as an example. 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 in 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 the 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 at 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 only estimate a parameter for each internal node and for tipdated trees we have one additional parameter for the rate.
## R version 4.2.2 Patched (2022-11-10 r83330)
## 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.41 phangorn_2.11.1 ape_5.6-2
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.9 bslib_0.4.2 compiler_4.2.2 jquerylib_0.1.4
## [5] highr_0.10 tools_4.2.2 digest_0.6.31 jsonlite_1.8.4
## [9] evaluate_0.20 lifecycle_1.0.3 nlme_3.1-161 lattice_0.20-45
## [13] pkgconfig_2.0.3 rlang_1.0.6 Matrix_1.5-3 fastmatch_1.1-3
## [17] igraph_1.3.5 cli_3.6.0 rstudioapi_0.14 yaml_2.3.6
## [21] parallel_4.2.2 xfun_0.36 fastmap_1.1.0 stringr_1.5.0
## [25] generics_0.1.3 vctrs_0.5.1 sass_0.4.4 grid_4.2.2
## [29] glue_1.6.2 R6_2.5.1 rmarkdown_2.20 magrittr_2.0.3
## [33] codetools_0.2-18 htmltools_0.5.4 quadprog_1.5-8 stringi_1.7.12
## [37] cachem_1.0.6
These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.