Showing multiple conflicting phylogenetic trees in the same plot

Yesterday, an R-sig-phylo subscriber posted the following question:

“For comparison of two trees I can use very nice function cophyloplot plotting two trees (left and right) and connecting respective nodes by lines. Very nice and convenient to read. But only for two trees. Displaying multiple trees in multiple comparisons is not very convenient. To display dozens to hundreds of trees there is densitree. Also nice, but for this purpose I don't like its display. I have several trees (~5) and I wish to compare their topologies, show supports (at least for differing nodes) and highlight differences. I thought about some overlay/parallel plotting (similar to the attached image) where there would be complete topologies displayed and incongruences would be easily visible. It would be probably doable by plotting all separate trees by plot.phylo and then combining and tuning the figure in some vector editor (like Inkscape). But I hope there is some more automated way to do it.”

The poster also supplied the following image as example:

Here are a couple of possible solutions using phytools.

Solution 1: Overlay the trees with “inner” node position.

library(phytools)

## Loading required package: ape

## Loading required package: maps

trees

## 9 phylogenetic trees

## get depths of the tree
h<-sapply(trees,function(x) max(nodeHeights(x)))
## test plot the longest tree to get the maximum x dimension
plotTree(trees[[which(h==max(h))]],direction="leftwards")

xlim<-get("last_plot.phylo",envir=.PlotPhyloEnv)$x.lim
## set tip order based on a majority rule consensus tree
tips<-setNames(1:Ntip(trees[[1]]), 
    untangle(consensus(trees,p=0.5),"read.tree")$tip.label)
## plot our trees
for(i in 1:length(trees)){
    plotTree(trees[[i]],color=make.transparent("blue",0.2),
        tips=tips,add=i>1,direction="leftwards",xlim=xlim[2:1],
        mar=c(4.1,1.1,1.1,1.1),nodes="inner",
        ftype=if(i==1) "i" else "off")
    if(i==1) axis(1)
}

Solution 2: same, but fix to constant depth.

rescale.tree<-function(tree,scale){
    tree$edge.length<-tree$edge.length/max(nodeHeights(tree))*scale
    tree
}
rescaled<-lapply(trees,rescale.tree,scale=1)
plotTree(rescaled[[1]],direction="leftwards")

xlim<-get("last_plot.phylo",envir=.PlotPhyloEnv)$x.lim
## set tip order based on a majority rule consensus tree
tips<-setNames(1:Ntip(rescaled[[1]]), 
    untangle(consensus(rescaled,p=0.5),"read.tree")$tip.label)
## plot our trees
for(i in 1:length(rescaled)){
    plotTree(rescaled[[i]],color=make.transparent("blue",0.2),
        tips=tips,add=i>1,direction="leftwards",xlim=xlim[2:1],
        mar=c(4.1,1.1,1.1,1.1),nodes="inner",
        ftype=if(i==1) "i" else "off")
    if(i==1) axis(1)
}

Solution 3: same as #1, but using different colors.

## get depths of the tree
h<-sapply(trees,function(x) max(nodeHeights(x)))
## test plot the longest tree to get the maximum x dimension
plotTree(trees[[which(h==max(h))]],direction="leftwards")

xlim<-get("last_plot.phylo",envir=.PlotPhyloEnv)$x.lim
## set tip order based on a majority rule consensus tree
tips<-setNames(1:Ntip(trees[[1]]), 
    untangle(consensus(trees,p=0.5),"read.tree")$tip.label)
## set colors
colors<-rainbow(n=length(trees))
## plot our trees
for(i in 1:length(trees)){
    plotTree(trees[[i]],color=make.transparent(colors[i],0.4),
        tips=tips,add=i>1,direction="leftwards",xlim=xlim[2:1],
        mar=c(4.1,1.1,1.1,1.1),nodes="inner",
        ftype=if(i==1) "i" else "off")
    if(i==1) axis(1)
}

Solution 4: same as #3, but with a slight vertical offset of each tree.

ylim<-c(0,Ntip(trees[[1]])+1)
## get depths of the tree
h<-sapply(trees,function(x) max(nodeHeights(x)))
## test plot the longest tree to get the maximum x dimension
plotTree(trees[[which(h==max(h))]],direction="leftwards")

xlim<-get("last_plot.phylo",envir=.PlotPhyloEnv)$x.lim
## set tip order based on a majority rule consensus tree
tips<-setNames(1:Ntip(trees[[1]]), 
    untangle(consensus(trees,p=0.5),"read.tree")$tip.label)
## set colors
colors<-rainbow(n=length(trees))
## plot our trees
for(i in 1:length(trees)){
    y.shift<-(i-median(1:length(trees)))/length(trees)/2
    plotTree(trees[[i]],color=make.transparent(colors[i],0.4),
        tips=tips+y.shift,add=i>1,direction="leftwards",xlim=xlim[2:1],
        mar=c(4.1,1.1,1.1,1.1),nodes="inner",
        ftype=if(i==floor(median(1:length(trees)))) "i" else "off",
        ylim=ylim)
    if(i==1) axis(1)
}

Solution 5: same as #4, but with a fixed total depth.

rescale.tree<-function(tree,scale){
    tree$edge.length<-tree$edge.length/max(nodeHeights(tree))*scale
    tree
}
rescaled<-lapply(trees,rescale.tree,scale=1)
ylim<-c(0,Ntip(rescaled[[1]])+1)
## get depths of the tree
h<-sapply(rescaled,function(x) max(nodeHeights(x)))
## test plot the longest tree to get the maximum x dimension
plotTree(rescaled[[i]],direction="leftwards")

xlim<-get("last_plot.phylo",envir=.PlotPhyloEnv)$x.lim
## set tip order based on a majority rule consensus tree
tips<-setNames(1:Ntip(rescaled[[1]]), 
    untangle(consensus(rescaled,p=0.5),"read.tree")$tip.label)
## set colors
colors<-rainbow(n=length(rescaled))
## plot our trees
for(i in 1:length(rescaled)){
    y.shift<-(i-median(1:length(rescaled)))/length(rescaled)/2
    plotTree(rescaled[[i]],color=make.transparent(colors[i],0.4),
        tips=tips+y.shift,add=i>1,direction="leftwards",xlim=xlim[2:1],
        mar=c(4.1,1.1,1.1,1.1),nodes="inner",
        ftype=if(i==floor(median(1:length(rescaled)))) "i" else "off",
        ylim=ylim)
    if(i==1) axis(1)
}

We don't need to use “inner” node positions, I just prefer that because I feel that it makes it easier to visualize similarity of topology. Here is one final solution (#6) using “standard” node positions:

rescale.tree<-function(tree,scale){
    tree$edge.length<-tree$edge.length/max(nodeHeights(tree))*scale
    tree
}
rescaled<-lapply(trees,rescale.tree,scale=1)
ylim<-c(0,Ntip(trees[[1]])+1)
## get depths of the tree
h<-sapply(rescaled,function(x) max(nodeHeights(x)))
## test plot the longest tree to get the maximum x dimension
plotTree(rescaled[[i]],direction="leftwards")

xlim<-get("last_plot.phylo",envir=.PlotPhyloEnv)$x.lim
## set tip order based on a majority rule consensus tree
tips<-setNames(1:Ntip(rescaled[[1]]), 
    untangle(consensus(rescaled,p=0.5),"read.tree")$tip.label)
## set colors
colors<-rainbow(n=length(rescaled))
## plot our trees
for(i in 1:length(rescaled)){
    y.shift<-(i-median(1:length(rescaled)))/length(rescaled)/2
    plotTree(rescaled[[i]],color=make.transparent(colors[i],0.4),
        tips=tips+y.shift,add=i>1,direction="leftwards",xlim=xlim[2:1],
        mar=c(4.1,1.1,1.1,1.1),
        ftype=if(i==floor(median(1:length(rescaled)))) "i" else "off",
        ylim=ylim)
    if(i==1) axis(1)
}

That's it.

BTW, the trees used in this exercise were simulated to be similar but different using the following code with phangorn:

library(phangorn)
true.tree<-pbtree(n=10,tip.label=LETTERS[1:10])
trees<-untangle(rNNI(true.tree,n=9),"read.tree")
D<-cophenetic(true.tree)
trees<-lapply(trees,nnls.tree,dm=D,rooted=TRUE)
class(trees)<-"multiPhylo"
trees<-untangle(trees,"read.tree")
rescale.tree<-function(tree,scale){
    tree$edge.length<-tree$edge.length/max(nodeHeights(tree))*scale
    tree
}
h<-rnorm(n=length(trees),mean=1,sd=0.1)
trees<-mapply(rescale.tree,tree=trees,scale=h,SIMPLIFY=FALSE)
class(trees)<-"multiPhylo"

[Ed. - A previous version of this post had ylim=c(0,11) rather than ylim<-c(0,Ntip(trees[[1]])+1) & ylim=ylim as in the current edit. Sorry about that!]

[Ed. #2 - I have also changed i==median(1:length(trees)) to i==floor(median(1:length(trees))) to address cases of an even number of trees in which the median^th tree is not an integer!]