Loading treescape and data:
Load the required packages:
library("treescape")
library("phangorn")
Load BEAST trees:
We load a random sample of 500 of the trees (from the second half of the posterior) produced using BEAST v1.8 with xml file 4 from Drummond and Rambaut (2007). It uses the standard GTR + Gamma + I substitution model with uncorrelated lognormal-distributed relaxed molecular clock. Each tree has 17 tips.
For convenience in our initial analysis we will take a random sample of 200 of these trees; sample sizes can be increased later.
set.seed(123)
BEASTtrees <- DengueTrees[sample(1:length(DengueTrees),200)]
Load nucleotide sequences:
Creating neighbour-joining and maximum likelihood trees:
Create a neighbour-joining (NJ) tree using the Tamura and Nei (1993) model (see ?dist.dna for more information) and root it on the outgroup "D4Thai63":
Dnj <- nj(dist.dna(DengueSeqs, model = "TN93"))
DnjRooted <- root(Dnj, resolve.root=TRUE, outgroup="D4Thai63")
plot(DnjRooted)
We note that bootstrapping gives 100% support:
Dnjboots <- boot.phylo(Dnj, as.matrix(DengueSeqs, model = "TN93"), B=100, function(e)
nj(dist.dna(DengueSeqs, model = "TN93")), trees=TRUE)
## $BP
## [1] 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100
##
## $trees
## 100 phylogenetic trees
Create a maximum-likelihood (ML) tree
Dfit.ini <- pml(DnjRooted, as.phyDat(DengueSeqs), k=4)
Dfit <- optim.pml(Dfit.ini, optNni=TRUE, optBf=TRUE,
optQ=TRUE, optGamma=TRUE, model="GTR")
## Warning: I unrooted the tree
# root:
DfitTreeRooted <- root(Dfit$tree, resolve.root=TRUE, outgroup="D4Thai63")
# bootstrap supports:
DMLboots <- bootstrap.pml(Dfit, optNni=TRUE)
# root:
DMLbootsrooted <- lapply(DMLboots, function(x) root(x, resolve.root=TRUE, outgroup="D4Thai63"))
Using treescape to compare trees
We now use the function treescape to find and plot distances between all these trees:
# collect the trees into a single object of class multiPhylo:
DengueTrees <- c(BEASTtrees,DMLbootsrooted,list(DfitTreeRooted),list(DnjRooted))
class(DengueTrees) <- "multiPhylo"
# add tree names:
names(DengueTrees)[1:200] <- paste0("BEAST",1:200)
names(DengueTrees)[201:300] <- paste0("ML_boot",1:100)
names(DengueTrees)[[301]] <- "ML"
names(DengueTrees)[[302]] <- "NJ"
# create vector corresponding to tree inference method:
Dtype <- c(rep("BEAST",200),rep("MLboots",100),"ML","NJ")
# use treescape to find and project the distances:
Dscape <- treescape(DengueTrees, nf=5)
# simple plot:
plotGrovesD3(Dscape$pco, groups=Dtype)
The function plotGrovesD3 produces interactive d3 plots which enable zooming, moving, tooltip text and legend hovering. We now refine the plot with colour-blind friendly colours (selected using ColorBrewer2), bigger points, varying point opacity, informative legend title and smaller legend width:
Dcols <- c("#abd9e9","#fdae61","#d7191c","#2c7bb6")
plotGrovesD3(Dscape$pco, groups=Dtype, colors=Dcols, point_size=200, point_opacity=c(rep(0.5,300),1,1), col_lab="Tree type", legend_width=80)
We can also add tree labels to the plot. Where these overlap, the user can use “drag and drop” to move them around for better visibility.
Dcols <- c("#abd9e9","#fdae61","#d7191c","#2c7bb6")
plotGrovesD3(Dscape$pco, groups=Dtype, treeNames = names(DengueTrees), colors=Dcols, point_size=200, point_opacity=c(rep(0.5,300),1,1), col_lab="Tree type", legend_width=80)
Alternatively, where labels are too cluttered, it may be preferable not to plot them but to make the tree names available as tooltip text instead:
Dcols <- c("#abd9e9","#fdae61","#d7191c","#2c7bb6")
plotGrovesD3(Dscape$pco, groups=Dtype, tooltip_text = names(DengueTrees), colors=Dcols, point_size=200, point_opacity=c(rep(0.5,300),1,1), col_lab="Tree type", legend_width=80)
The scree plot is available as part of the treescape output:
barplot(Dscape$pco$eig, col="navy")
We can also view the plot in 3D:
Dcols3D <- c(rep(Dcols[[1]],200),rep(Dcols[[2]],100),Dcols[[3]],Dcols[[4]])
rgl::plot3d(Dscape$pco$li[,1],Dscape$pco$li[,2],Dscape$pco$li[,3],
type="s", size=c(rep(1.5,300),2,2), col=Dcols3D,
xlab="", ylab="", zlab="")
You must enable Javascript to view this page properly.
treescape analysis
From these plots we can see that treescape has identified variation in the trees according to the Kendall Colijn metric (\(\lambda=0\), ignoring branch lengths). We can compare pairs of trees using the plotTreeDiff function to see exactly where these differences arise. Tips with identical ancestry in the two trees are coloured grey, whereas tips with differing ancestry are coloured peach-red, with the colour darkening according to the number of ancestral differences found at each tip. Since we are comparing the trees topologically (ignoring branch lengths, for the moment), we plot with constant branch lengths for clarity.
# comparing NJ and ML:
plotTreeDiff(DnjRooted,DfitTreeRooted, use.edge.length=FALSE)
treeDist(DnjRooted,DfitTreeRooted)
## [1] 13.93
For pairwise comparisons it is helpful to find a small number of representative trees. We can find a geometric median tree from the BEAST trees using the medTree function:
BEASTmed <- medTree(BEASTtrees)
There are two median trees, with identical topology:
## 2 phylogenetic trees
treeDist(BEASTmed$trees[[1]],BEASTmed$trees[[2]])
## [1] 0
so we may select one of them as a BEAST representative tree. Note that for a more thorough analysis it may be appropriate to identify clusters among the BEAST trees and select a summary tree from each cluster: we demonstrate this approach later in the vignette.
BEASTrep <- BEASTmed$trees[[1]]
# comparing BEAST median and NJ:
plotTreeDiff(BEASTrep,DnjRooted, use.edge.length=FALSE)
treeDist(BEASTrep,DnjRooted)
## [1] 13.27
# comparing BEAST median and ML:
plotTreeDiff(BEASTrep,DfitTreeRooted, use.edge.length=FALSE)
treeDist(BEASTrep,DfitTreeRooted)
## [1] 9.487
# comparing BEAST median to a random BEAST tree:
num <- runif(1,1,200)
randomBEASTtree <- BEASTtrees[[num]]
plotTreeDiff(BEASTrep, randomBEASTtree, use.edge.length=FALSE)
treeDist(BEASTrep,randomBEASTtree)
## [1] 5.568
Using treescape to analyse the BEAST trees in more detail:
We used TreeAnnotator (Drummond and Rambaut, 2007) to create a Maximum Clade Credibility (MCC) tree from amongst the BEAST trees.
# load the MCC tree
data(DengueBEASTMCC)
# concatenate with other BEAST trees
BEAST201 <- c(BEASTtrees,list(DengueBEASTMCC))
# compare using treescape:
BEASTscape <- treescape(BEAST201, nf=5)
# simple plot:
plotGrovesD3(BEASTscape$pco)
There appear to be clusters of tree topologies within the BEAST trees. We can use the function findGroves to identify clusters:
# find clusters or 'groves':
BEASTGroves <- findGroves(BEASTscape,nclust=4,clustering="single")
and to find a median tree per cluster:
# find median tree(s) per cluster:
BEASTMeds <- medTree(BEAST201, groups=BEASTGroves$groups)
# for each cluster, select a single median tree to represent it:
BEASTMedTrees <- c(BEASTMeds$`1`$trees[[1]],
BEASTMeds$`2`$trees[[1]],
BEASTMeds$`3`$trees[[1]],
BEASTMeds$`4`$trees[[1]])
We can now make the plot again, highlighting the MCC tree and the four median trees:
# extract the numbers from the tree list 'BEASTtrees' which correspond to the median trees:
BEASTMedTreeNums <-c(which(BEASTGroves$groups==1)[[BEASTMeds$`1`$treenumbers[[1]]]],
which(BEASTGroves$groups==2)[[BEASTMeds$`2`$treenumbers[[1]]]],
which(BEASTGroves$groups==3)[[BEASTMeds$`3`$treenumbers[[1]]]],
which(BEASTGroves$groups==4)[[BEASTMeds$`4`$treenumbers[[1]]]])
# prepare a vector to highlight median and MCC trees
highlightTrees <- rep(1,201)
highlightTrees[[201]] <- 2
highlightTrees[BEASTMedTreeNums] <- 2
# plot:
plotGrovesD3(BEASTscape$pco,
groups=as.vector(BEASTGroves$groups),
colors=Dcols,
col_lab="Cluster",
symbol_var = highlightTrees,
size_range = c(60,600),
size_var = highlightTrees,
legend_width=0)
To understand the differences between the representative trees we can use plotTreeDiff again, for example:
# differences between the MCC tree and the median from the largest cluster:
treeDist(DengueBEASTMCC,BEASTMedTrees[[1]])
## [1] 2
plotTreeDiff(DengueBEASTMCC,BEASTMedTrees[[1]], use.edge.length=FALSE)
# differences between the median trees from clusters 1 and 2:
treeDist(BEASTMedTrees[[1]],BEASTMedTrees[[2]])
## [1] 10.63
plotTreeDiff(BEASTMedTrees[[1]],BEASTMedTrees[[2]], use.edge.length=FALSE)
References
[1] Drummond, A. J., and Rambaut, A. (2007) BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evolutionary Biology, 7(1), 214.
[2] Lanciotti, R. S., Gubler, D. J., and Trent, D. W. (1997) Molecular evolution and phylogeny of dengue-4 viruses. Journal of General Virology, 78(9), 2279-2286.