Description: This script provides examples of the new functions available in the phangorn library to ‘intertwine’ trees and networks, i.e. compare trees and networks and data transferrance. It also provides a step-by-step guide for users new to R.
Methodological advancement: The major advancement in this phangorn update is the introduction of a generic network object with a wide range of related transfer and analysis functions. These new functions provide the first means to directly transfer information amongst a wide range of phylogenetic trees and networks, as well as means to visualise and further analyse this information. This should provide a platform for individuals to easily conduct tree-network comparisons and stimulate further function development by the community.
What next?: By implementing full network handling compatibility in R, and providing exemplar scripts (such as this) and support, the scientific community now has an easy means to analyse and compare the results of phylogenetic trees vs. network approaches. We hope this will open a largely unexplored world to the general phylogenetic audience.
Open R and run the two lines of code below in the command line.
(You will need to select a region from which to download the library)
install.packages("phangorn", dependencies=TRUE)
# install latest development version needs devtools
install.packages("devtools", dependencies=TRUE)
library(devtools)
install_github("KlausVigo/phangorn")library(phangorn) # load the phangorn libraryThis is often a major stumbling block for new R users. You need to specify where in your folder structure you wish to work. i.e, where are the files stored that you wish to input?
This is done using the setwd() function, e.g. setwd(“C:/TreesNetworks/Example Files”).
We now set it to the folder of the phangorn package, which contains the files we want to load for this example.
This example files are based on the woodmouse dataset available in the ape library. Ultimately, this dataset is from this study: Michaux, J. R., Magnanou, E., Paradis, E., Nieberding, C. and Libois, R. (2003) Mitochondrial phylogeography of the Woodmouse (Apodemus sylvaticus) in the Western Palearctic region. Molecular Ecology, 12, 685-697.)
## automatically set the correct working directory for the examples below
# setwd(system.file("extdata/trees", package = "phangorn"))
# for this vignette we create a path to the files we want to load
fdir <- system.file("extdata/trees", package = "phangorn")
## in your case it may look something like this...
# setwd("C:/TreesNetworks/Example Files")
##DNA Matrix, maybe not needed
woodmouse <- read.phyDat(file.path(fdir, "woodmouse.fasta"),format="fasta")
## RAxML best-known tree with bipartition support (from previous analysis)
raxml.tree <- read.tree(file.path(fdir,"RAxML_bipartitions.woodmouse"))
## RAxML bootstrap trees (from previous analysis)
raxml.bootstrap <- read.tree(file.path(fdir,"RAxML_bootstrap.woodmouse"))
## MrBayes consensus tree (50% majority rule) (from previous analysis)
mrbayes.tree <- read.nexus(file.path(fdir,"woodmouse.mrbayes.nex.con"))
## MrBayes sample runs 1 and 2 (from previous analysis)
run1 <- read.nexus(file.path(fdir,"woodmouse.mrbayes.nex.run1.t"))
run2 <- read.nexus(file.path(fdir,"woodmouse.mrbayes.nex.run2.t"))
## How many trees are in the MrBayes tree sample?
run1## 1001 phylogenetic treesrun2## 1001 phylogenetic trees## Combining the two runs and removing 25% burn-in
mrbayes.trees <- c(run1[251:1001],run2[251:1001])
## NeigbourNet Nexus file generated by SplitsTree (from previous analysis)
Nnet <- read.nexus.networx(file.path(fdir,"woodmouse.nxs"))All example files read into R.
par(mfrow=c(1,2), mar=c(2,2,2,2)) # Setting plot parameters
### Plotting trees with support values:
## RAxML
plot(raxml.tree)
nodelabels(raxml.tree$node.label, adj = c(1, 0), frame = "none")
## MrBayes
plot(mrbayes.tree)
nodelabels(mrbayes.tree$node.label, adj = c(1, 0), frame = "none")
par(mfrow=c(1,1)) # Setting plot parameters
# NeighbourNet
plot.networx(Nnet,"2D")
## alternatively,
# plot(Nnet,"2D")Identification of edge bundles (in black) in a neighbour-net (NN) network that correspond to branches (labelled 1-12) in a tree (a maximum likelihood tree, in this case).
# create a vector of labels for the network corresponding to edges in the tree
edge.lab <- createLabel(Nnet, raxml.tree, raxml.tree$edge[,2], "edge")
# could be also 1:27 instead of raxml.tree$edge[,2]
# Show the correspondingly labelled tree and network in R
par(mfrow=c(1,2), mar=c(1,1,1,1))
#plotBS(raxml.tree, rotate.tree = 180)
plot(raxml.tree, "u", rotate.tree = 180)
edgelabels(raxml.tree$edge[,2],col="blue", frame="none")
# find edges that are in the network but not in the tree
edge.col <- rep("black", nrow(Nnet$edge))
edge.col[ is.na(edge.lab) ] <- "red"
# or a simpler alternative...
edge.col <- createLabel(Nnet, raxml.tree, "black", nomatch="red")
x <- plot.networx(Nnet, edge.label = edge.lab, show.edge.label = T, "2D", edge.color = edge.col,
col.edge.label = "blue")
# the above plot function returns an invisible networx object and this object also
# contains the colors for the edges.Bootstrap support for all branches (branch labels) in the ML tree mapped on the corresponding edge bundles of the NN network. The edges in the network which are not found as ML tree branches are highlighted in red.
x <- addConfidences(Nnet,raxml.tree)
# find splits that are in the network but not in the tree
split.col <- rep("black", length(x$splits))
split.col[ !matchSplits(as.splits(x), as.splits(raxml.tree)) ] <- "red"
# simpler alternative...
split.col2 <- createLabel(x, raxml.tree, label="black", "split", nomatch="red")
# Plotting in R
par(mfrow=c(1,1))
out.x <- plot.networx(x,"2D",show.edge.label=TRUE, split.color=split.col, col.edge.label = "blue")
Again we can write to SplitsTree for viewing…
# write.nexus.networx(out.x,"woodmouse.tree.support.nxs")
## or we can also export the splits alone (for usage in software other than SplitsTree)
# write.nexus.splits(as.splits(out.x),"woodmouse.splits.support.nxs")Frequencies of bipartitions found in the bootstrap pseudoreplicates mapped on the corresponding edge bundles of the NN network using a threshold of 10% (i.e. any edge is labelled that occurs in at least 100 of the 1000 ML-BS pseudoreplicates). Edge bundles not found in the ML tree are labelled using grey edges.
y <- addConfidences(Nnet, as.splits(raxml.bootstrap))
edge.col <- createLabel(y, raxml.tree, label="black", "edge", nomatch="grey")
y <- plot(y,"2D",show.edge.label=TRUE, edge.color=edge.col)
## Write to SplitsTree for viewing
# write.nexus.networx(y,"NN.with.bs.support.nxs")We can also compare the neighborNet with a consensusNet (Holland BR, Huber KT, Moulton V, Lockhart PJ,2004, Using consensus networks to visualize contradictory evidence for species phylogeny. Molecular Biology and Evolution, 21, 1459-1461). Furthermore, we can extract the support values from the consensus network, and place these on to the NeighbourNet (this is similar to the process explained in 1C above).
cnet <- consensusNet(raxml.bootstrap,prob=0.10)
edge.col <- createLabel(cnet, Nnet, label="black", "edge", nomatch="grey")
cnet <- plot.networx(cnet, "2D", show.edge.label = TRUE, edge.color=edge.col)
edge.col <- createLabel(Nnet, cnet, label="black", "edge", nomatch="grey")
z <- plot.networx(Nnet, "2D", show.edge.label = TRUE, edge.color=edge.col)
obj <- addConfidences(Nnet,cnet)
plot.networx(obj,"2D",show.edge.label=T, edge.color=edge.col, col.edge.label = "blue")
## Write to SplitsTree for viewing
# write.nexus.networx(obj,"Nnet.with.ML.Cnet.Bootstrap.nxs")