Load packages

We start by loading the phangorn package and setting a random seed:

library(phangorn)
set.seed(9)

Load data

The dataset we’re using contains morphological data for 12 mite species, with 79 encoded characters (Schäffer et al. 2010). When reading in the .csv file, row.names = 1 uses the first column (species) as row names. To get a phyDat object, we have to convert the dataframe into a matrix with as.matrix.

fdir <- system.file("extdata", package = "phangorn")
mm <- read.csv(file.path(fdir, "mites.csv"), row.names = 1)
mm_pd <- phyDat(as.matrix(mm), type = "USER", levels = 0:7)

The data can then be written into a nexus file:

write.phyDat(mm_pd, file.path(fdir, "mites.nex"), format = "nexus")

Reading in a nexus file is even easier than reading in a csv file:

mm_pd <- read.phyDat(file.path(fdir, "mites.nex"), format = "nexus", type = "STANDARD")

After reading in the nexus file, we have the states 0:9, but the data only has the states 0:7. Here is one possibility to change the contrast matrix:

contrast <- matrix(data = c(1,0,0,0,0,0,0,0,0,
    0,1,0,0,0,0,0,0,0,
    0,0,1,0,0,0,0,0,0,
    0,0,0,1,0,0,0,0,0,
    0,0,0,0,1,0,0,0,0,
    0,0,0,0,0,1,0,0,0,
    0,0,0,0,0,0,1,0,0,
    0,0,0,0,0,0,0,1,0,
    0,0,0,0,0,0,0,0,1,
    1,1,1,1,1,1,1,1,1),
    ncol = 9, byrow = TRUE)
dimnames(contrast) <- list(c(0:7,"-","?"),
    c(0:7, "-"))
contrast

##   0 1 2 3 4 5 6 7 -
## 0 1 0 0 0 0 0 0 0 0
## 1 0 1 0 0 0 0 0 0 0
## 2 0 0 1 0 0 0 0 0 0
## 3 0 0 0 1 0 0 0 0 0
## 4 0 0 0 0 1 0 0 0 0
## 5 0 0 0 0 0 1 0 0 0
## 6 0 0 0 0 0 0 1 0 0
## 7 0 0 0 0 0 0 0 1 0
## - 0 0 0 0 0 0 0 0 1
## ? 1 1 1 1 1 1 1 1 1

mm_pd <- phyDat(mm_pd, type="USER", contrast=contrast)

Now that we have our data, we can start the analyses.

Parsimony

For morphological data, one of the most frequently used approaches to conduct phylogenetic trees is maximum parsimony (MP). pratchet (as already described in Estimating phylogenetic trees with phangorn) implements the parsimony ratchet (Nixon 1999). To create a starting tree, we can use the function random.addition:

mm_start <- random.addition(mm_pd)

This tree can then be given to pratchet:

mm_tree <- pratchet(mm_pd, start = mm_start, minit = 1000, maxit = 10000,
                    all = TRUE, trace = 0)
mm_tree

## 19 phylogenetic trees

With all=TRUE we get all (in this case 19) trees with lowest parsimony score in a multiPhylo object. Since we we did a minimum of 1000 iterations, we already have some edge support. Now we can assign the edge lengths.

mm_tree <- acctran(mm_tree, mm_pd)

mm_bab <- bab(mm_pd, trace = 0)
mm_bab

## 37 phylogenetic trees

mm_tree_rooted <- root(mm_tree, outgroup = "C._cymba", resolve.root = TRUE,
                       edgelabel = TRUE)

# subsetting for tree nr. 9
plotBS(mm_tree_rooted[[9]], digits = 2)

# save plot as pdf
pdf(file = "mm_rooted.pdf")
plotBS(mm_tree_rooted[[9]], digits = 2)
dev.off()

Consensus tree

To look at the consensus tree of our 19 trees from pratchet, or of our 37 most parsimonious trees from bab, we can use the consensus function from ape.

# unrooted pratchet tree
mm_cons <- consensus(mm_tree)

# rooted pratchet tree
mm_cons_root <- consensus(mm_tree_rooted, rooted = TRUE)

# branch and bound, we root the consensus tree in the same step
mm_bab_cons <- root(consensus(mm_bab), outgroup = "C._cymba",
                    resolve.root = TRUE, edgelabel = TRUE)

plot(mm_cons, main="Unrooted pratchet consensus tree")
plot(mm_cons_root, main="Rooted pratchet consensus tree")
plot(mm_bab_cons, main="Rooted bab consensus tree")

Unrooted and rooted consensus trees of the mites dataset with MP.

We can clearly see that, as expected, the two rooted trees have the same topology.

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] 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      knitr_1.48       
## [21] tibble_3.2.1      munsell_0.5.1     bslib_0.8.0       pillar_1.9.0     
## [25] rlang_1.1.4       utf8_1.2.4        fastmatch_1.1-4   cachem_1.1.0     
## [29] xfun_0.47         quadprog_1.5-8    sass_0.4.9        cli_3.6.3        
## [33] withr_3.0.1       magrittr_2.0.3    digest_0.6.37     grid_4.4.1       
## [37] rstudioapi_0.16.0 lifecycle_1.0.4   nlme_3.1-165      vctrs_0.6.5      
## [41] evaluate_0.24.0   glue_1.7.0        farver_2.1.2      codetools_0.2-20 
## [45] ggseqlogo_0.2     fansi_1.0.6       colorspace_2.1-1  rmarkdown_2.28   
## [49] tools_4.4.1       pkgconfig_2.0.3   htmltools_0.5.8.1