Estimating phylogenetic trees with phangorn

Klaus Schliep, Iris Bardel-Kahr

Graz University of Technology, University of Graz
klaus.schliep@gmail.com

2024-09-17

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

Getting started

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

Distance based methods

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.

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.
Rooted UPGMA tree.
plot(treeNJ, "unrooted", main="NJ")
Unrooted NJ tree.
Unrooted NJ tree.

Bootstrap

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.

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

  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 tree 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 (minit) or no improvements have been recorded for a number of iterations (k).
treeRatchet  <- pratchet(primates, trace = 0, minit=100)
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:

treeRatchet  <- acctran(treeRatchet, primates)

After assigning edge weights, we prune away internal edges of length tol (default = 1e-08), 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)
}

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.

treeRA <- random.addition(primates)
treeSPR  <- optim.parsimony(treeRA, primates)
## Final p-score 746 after  1 nni operations
parsimony(c(treeRA, treeSPR), primates)
## [1] 756 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 is. 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).

Model selection

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.

mt <- modelTest(primates)

It’s also possible to only select 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 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.

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

Conducting a ML tree

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.

fit_mt <- pml_bb(mt, control = pml.control(trace = 0))
fit_mt
## model: HKY+G(4) 
## loglikelihood: -2615 
## unconstrained loglikelihood: -1230 
## Model of rate heterogeneity: Discrete gamma model
## Number of rate categories: 4 
## Shape parameter: 2.3 
##   Rate Proportion
## 1 0.32       0.25
## 2 0.68       0.25
## 3 1.08       0.25
## 4 1.92       0.25
## 
## Rate matrix:
##    a  c  g  t
## a  0  1 55  1
## c  1  0  1 55
## g 55  1  0  1
## t  1 55  1  0
## 
## Base frequencies:  
##     a     c     g     t 
## 0.375 0.402 0.039 0.184

We can also use pml_bb with a defined model to infer a phylogenetic tree.

fitGTR <- pml_bb(primates, model="GTR+G(4)+I")

Bootstrap

If we instead want to conduct standard bootstrapping (Felsenstein 1985; Penny and Hendy 1985), we can do so with the function bootstrap.pml:

bs <- bootstrap.pml(fit_mt, bs=100, optNni=TRUE,
    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.Unrooted tree (midpoint rooted) with ultrafast, standard and transfer bootstrap support values.Unrooted tree (midpoint rooted) with ultrafast, standard and transfer bootstrap support values.

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
tree_stdbs <- plotBS(fit_mt$tree, bs, type = "n")

# assigning transfer bootstrap values to our tree
tree_tfbs <- plotBS(fit_mt$tree, bs, type = "n", method = "TBE")

It is also possible to look at consensusNet to identify potential conflict.

cnet <- consensusNet(bs, p=0.2)
plot(cnet, show.edge.label=TRUE)
ConsensusNet from the standard bootstrap sample.
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.

Exporting a tree

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")

Molecular dating with a strict clock for ultrametric and tipdated phylogenies

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. The function pml_bb ensures that the tree is ultrametric, or the constraints given by the tip dates are fulfilled. That differs from the function optim.pml where th tree supplied to the function has to fulfill the constraints. In this case for an ultrametric starting tree we can use an UPGMA or WPGMA tree.

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

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 in 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 the constraints of the tip.dates vector.

fit_td <- pml_bb(H3N2, model="HKY+I", method="tipdated", tip.dates=dates, 
               rearrangement="NNI", control = pml.control(trace = 0))
fit_td
## model: HKY+I 
## loglikelihood: -3118 
## unconstrained loglikelihood: -2884 
## Proportion of invariant sites: 0.69 
## 
## Rate: 0.0025 
## 
## Rate matrix:
##     a   c   g   t
## a 0.0 1.0 9.9 1.0
## c 1.0 0.0 1.0 9.9
## g 9.9 1.0 0.0 1.0
## t 1.0 9.9 1.0 0.0
## 
## Base frequencies:  
##    a    c    g    t 
## 0.31 0.19 0.24 0.26 
## 
## Rate: 0.0025

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. The rate is here comparable to the slope fo the tip-to-root regression in programs like TempEst (Rambaut et al. 2016).

And at last we plot the tree with a timescale.

plot(fit_td, align.tip.label=TRUE)

Session info

## R version 4.4.1 (2024-06-14)
## Platform: x86_64-pc-linux-gnu
## Running under: Debian GNU/Linux 12 (bookworm)
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.11.0 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.11.0
## 
## 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       
## 
## time zone: Europe/Vienna
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] knitr_1.48      phangorn_2.12.1 ape_5.8-0.1    
## 
## loaded via a namespace (and not attached):
##  [1] Matrix_1.7-0      gtable_0.3.5      jsonlite_1.8.8    dplyr_1.1.4      
##  [5] compiler_4.4.1    highr_0.11        tidyselect_1.2.1  Rcpp_1.0.13      
##  [9] parallel_4.4.1    jquerylib_0.1.4   scales_1.3.0      yaml_2.3.10      
## [13] fastmap_1.2.0     lattice_0.22-6    ggplot2_3.5.1     R6_2.5.1         
## [17] labeling_0.4.3    generics_0.1.3    igraph_2.0.3      tibble_3.2.1     
## [21] munsell_0.5.1     bslib_0.8.0       pillar_1.9.0      rlang_1.1.4      
## [25] utf8_1.2.4        fastmatch_1.1-4   cachem_1.1.0      xfun_0.47        
## [29] quadprog_1.5-8    sass_0.4.9        cli_3.6.3         withr_3.0.1      
## [33] magrittr_2.0.3    digest_0.6.37     grid_4.4.1        rstudioapi_0.16.0
## [37] lifecycle_1.0.4   nlme_3.1-165      vctrs_0.6.5       evaluate_0.24.0  
## [41] glue_1.7.0        farver_2.1.2      codetools_0.2-20  ggseqlogo_0.2    
## [45] fansi_1.0.6       colorspace_2.1-1  rmarkdown_2.28    tools_4.4.1      
## [49] pkgconfig_2.0.3   htmltools_0.5.8.1

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.
———. 1985. “Confidence Limits on Phylogenies. An Approach Using the Bootstrap.” Evolution 39: 783–91.
———. 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.
Lemoine, Fréderic, J-B Domelevo Entfellner, Eduan Wilkinson, Damien Correia, M Dávila Felipe, Tulio De Oliveira, and Olivier Gascuel. 2018. “Renewing Felsenstein’s Phylogenetic Bootstrap in the Era of Big Data.” Nature 556 (7702): 452–56.
Minh, Bui Quang, Minh Anh Thi Nguyen, and Arndt von Haeseler. 2013. “Ultrafast Approximation for Phylogenetic Bootstrap.” Molecular Biology and Evolution 30 (5): 1188–95.
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.
Penny, D., and M. D. Hendy. 1985. “Testing Methods Evolutionary Tree Construction.” Cladistics 1: 266–78.
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.
Rambaut, Andrew, Tommy T. Lam, Luiz Max Carvalho, and Oliver G. Pybus. 2016. “Exploring the temporal structure of heterochronous sequences using TempEst (formerly Path-O-Gen).” Virus Evolution 2 (1). https://doi.org/10.1093/ve/vew007.
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.
Yang, Ziheng. 2006. Computational Molecular Evolution. Oxford: Oxford University Press.
Zuckerkandl, Emile, and Linus Pauling. 1965. “Molecules as Documents of Evolutionary History.” Journal of Theoretical Biology 8 (2): 357–66.