Phylogenetic trees from morphological data
Iris Bardel-Kahr, Klaus Schliep
University of Graz, Graz University of Technologyklaus.schliep@gmail.com
2024-12-05
Source:vignettes/Morphological.Rmd
Morphological.Rmd
In this vignette, we will show how to work with morphological data in phangorn (Schliep 2011). In most cases the different morphological characters or character states are encoded with the numbers 0:9 (or less, if there are less differences). Morphological data can come in different formats. The most common ones are .csv and .nexus.
Load packages
We start by loading the phangorn package and setting a random seed:
## Loading required package: ape
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)
Branch and bound
In the case of our mites-dataset with 12 sequences, it’s also
possible to use the branch and bound algorithm (Hendy and Penny 1982) to find all most
parsimonious trees. With bigger datasets it is definitely recommended to
use pratchet
.
mm_bab <- bab(mm_pd, trace = 0)
mm_bab
## 37 phylogenetic trees
Root trees
If we want our unrooted trees to be rooted, we have the possibility
to use midpoint
to perform midpoint rooting. Rooting the
trees with a specific species (we chose C. cymba here) can be
done with the function root
from the ape package
(Paradis and Schliep 2019). To save the
correct node labels (edge support), it’s important to set
edgelabel=TRUE
.
mm_tree_rooted <- root(mm_tree, outgroup = "C._cymba", resolve.root = TRUE,
edgelabel = TRUE)
Plot trees
With plotBS
we plot a tree with their respective edge
support. It is also possible to save the plots as .pdf (or
various other formats, e.g. svg, png, tiff) file. digits
is
an argument to determine the number of digits shown for the bootstrap
values.
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")
We can clearly see that, as expected, the two rooted trees have the same topology.
Session info
## R version 4.4.2 (2024-10-31)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 22.04.5 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; LAPACK version 3.10.0
##
## locale:
## [1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8
## [4] LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8
## [7] LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C
## [10] LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C
##
## time zone: UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] phangorn_2.12.1.1 ape_5.8-0.5
##
## loaded via a namespace (and not attached):
## [1] nlme_3.1-166 cli_3.6.3 knitr_1.49 rlang_1.1.4
## [5] xfun_0.49 generics_0.1.3 textshaping_0.4.0 jsonlite_1.8.9
## [9] htmltools_0.5.8.1 ragg_1.3.3 sass_0.4.9 rmarkdown_2.29
## [13] quadprog_1.5-8 grid_4.4.2 evaluate_1.0.1 jquerylib_0.1.4
## [17] fastmap_1.2.0 yaml_2.3.10 lifecycle_1.0.4 compiler_4.4.2
## [21] codetools_0.2-20 igraph_2.1.1 fs_1.6.5 pkgconfig_2.0.3
## [25] Rcpp_1.0.13-1 htmlwidgets_1.6.4 systemfonts_1.1.0 lattice_0.22-6
## [29] digest_0.6.37 R6_2.5.1 parallel_4.4.2 magrittr_2.0.3
## [33] Matrix_1.7-1 bslib_0.8.0 tools_4.4.2 pkgdown_2.1.1
## [37] cachem_1.1.0 desc_1.4.3