On a strange algorithm for plotting trees in R

Here’s a weird way to plot a phylogenetic tree in R that that I haven’t found a use for yet (but still might).

First, get the coordinates of the nodes of the tree we’re going to plot.

Then, iterating over all the tips of the tree, compute the set of ancestors from the tip to the global root. For this, we’ll use the phangorn package, but we could’ve also done it via a post-order traversal.

Finally, graph a line (with type="s" for a square phylogram) from through the coordinates on each path from the root to each tip.

This is what this looks like in R.

## load libraries
library(phytools)
library(phangorn)
## load an example tree
data(mammal.tree)
## compute coordinates of tree to plot
plotTree(mammal.tree,color="transparent",fsize=0.7,
	ftype="i")
## get coordinates from "last_plot.phylo"
obj<-get("last_plot.phylo",envir=.PlotPhyloEnv)
for(i in 1:Ntip(mammal.tree)){
	## get the ancestors of every tip
	aa<-c(i,Ancestors(mammal.tree,i))
	## plot the line segments, using type="s", from
	## the root
	lines(obj$xx[aa],obj$yy[aa],type="s",lwd=5,
		col=palette()[2])
}

To see that the tree has been graphed this way, let’s assign an α transparency value to our color and see that edges near the root are plotted many times, while edges near the tips are plotted few times or only once.

plotTree(mammal.tree,color="transparent",fsize=0.7,
	ftype="i")
obj<-get("last_plot.phylo",envir=.PlotPhyloEnv)
for(i in 1:Ntip(mammal.tree)){
	aa<-c(i,Ancestors(mammal.tree,i))
	lines(obj$xx[aa],obj$yy[aa],type="s",lwd=5,
		col=make.transparent(palette()[2],0.1))
}

The only obvious advantage of using this algorithm, compared to plotting the edges and vertical connecting lines, one by one is that it allows us to more easily use “beveled” line join style. (This has no meaning if each line segment is plotted separately.)

For example:

par(ljoin=2,lend=1)
plotTree(mammal.tree,color="transparent",fsize=0.7,
	ftype="i")
obj<-get("last_plot.phylo",envir=.PlotPhyloEnv)
for(i in 1:Ntip(mammal.tree)){
	aa<-c(i,Ancestors(mammal.tree,i))
	lines(obj$xx[aa],obj$yy[aa],type="s",lwd=7,
		col=palette()[4])
}

This compared to with phytools::plotTree.

par(ljoin=2,lend=1)
plotTree(mammal.tree,fsize=0.7,ftype="i",lwd=7,
	color=palette()[4])

Lastly, just to see what accidental creation results, let’s try type="fan".

data(salamanders)
par(ljoin=2,lend=1)
plotTree(salamanders,color="transparent",
	ftype="i",type="fan")
obj<-get("last_plot.phylo",envir=.PlotPhyloEnv)
for(i in 1:Ntip(salamanders)){
	aa<-c(i,Ancestors(salamanders,i))
	lines(obj$xx[aa],obj$yy[aa],type="s",lwd=3,
		col=palette()[4])
}

Yikes.

Maybe a “radial” style will come out better.

par(ljoin=2,lend=1)
plotTree(salamanders,color="transparent",
	ftype="i",type="fan")
obj<-get("last_plot.phylo",envir=.PlotPhyloEnv)
for(i in 1:Ntip(salamanders)){
	aa<-c(i,Ancestors(salamanders,i))
	lines(obj$xx[aa],obj$yy[aa],lwd=3,
		col=palette()[4])
}

Fix to bug in phytools force.ultrametric function to coerce a tree to be ultrametric

R phylogenetics power-user & comparative methods star Jacob Berv contacted me this afternoon with a bug in the phytools method force.ultrametric in which the function appeared to be creating a basal polytomy where none had existed in the input tree.

For context, phytools::force.ultrametric is a not a formal method for ultrametricizing a phylogenetic tree (like non-parametric rate smoothing, for example). Instead, it is designed to implement a couple of different simple procedures to coerce a tree to pass ape::is.ultrametric even if it should fail it due to numerical precision (e.g., rounding) issues with the edge lengths.

This is done using two main algorithms. The default (denominated method="nnls") uses non-negative least squares to find the ultrametric set of edge lengths that results in the set of patristic distances with the minimum sum of squared difference to the patristic distances implied by the original tree. (The other just extends all the tips – but this method was unaffected by the bug discovered by Jacob.)

Here is a demonstration of the bug, using the example tree supplied by Jacob.

library(phytools)
library(phangorn)
tree<-read.tree(file="test.tre")
par(mfrow=c(1,2))
plotTree(tree,ftype="off")
plotTree(phytools::force.ultrametric(tree),ftype="off")

## ***************************************************************
## *                          Note:                              *
## *    force.ultrametric does not include a formal method to    *
## *    ultrametricize a tree & should only be used to coerce    *
## *   a phylogeny that fails is.ultramtric due to rounding --   *
## *    not as a substitute for formal rate-smoothing methods.   *
## ***************************************************************

OK. Something is obviously wrong here.

Since method="nnls" actually uses phangorn::nnls.tree internally, I was able to quickly trace the issue to this phangorn update by Klaus Schliep.

Here is the new, fixed version of force.ultrametric which I have already pushed to GitHub.

force.ultrametric<-function(tree,method=c("nnls","extend"),...){
	if(hasArg(message)) message<-list(...)$message
	else message<-TRUE
	if(message){
		cat("***************************************************************\n")
		cat("*                          Note:                              *\n")
		cat("*    force.ultrametric does not include a formal method to    *\n")
		cat("*    ultrametricize a tree & should only be used to coerce    *\n")
		cat("*   a phylogeny that fails is.ultramtric due to rounding --   *\n")
		cat("*    not as a substitute for formal rate-smoothing methods.   *\n")
		cat("***************************************************************\n")
	}
	method<-method[1]
	if(method=="nnls") tree<-nnls.tree(cophenetic(tree),tree,
		method="ultrametric",rooted=is.rooted(tree),trace=0)
	else if(method=="extend"){
		h<-diag(vcv(tree))
		d<-max(h)-h
		ii<-sapply(1:Ntip(tree),function(x,y) which(y==x),
			y=tree$edge[,2])
		tree$edge.length[ii]<-tree$edge.length[ii]+d
	} else 
		cat("method not recognized: returning input tree\n\n")
	tree
}

Let’s make sure it works by repeating our plot from above.

par(mfrow=c(1,3))
plotTree(tree,ftype="off",mar=c(0.1,1.1,2.1,1.1))
mtext("a) input tree",line=0,adj=0.1)
plotTree(phytools::force.ultrametric(tree),ftype="off",
	mar=c(0.1,1.1,2.1,1.1))

## ***************************************************************
## *                          Note:                              *
## *    force.ultrametric does not include a formal method to    *
## *    ultrametricize a tree & should only be used to coerce    *
## *   a phylogeny that fails is.ultramtric due to rounding --   *
## *    not as a substitute for formal rate-smoothing methods.   *
## ***************************************************************

mtext("b) CRAN force.ultrametric",line=0,adj=0.1)
plotTree(force.ultrametric(tree),ftype="off",
	mar=c(0.1,1.1,2.1,1.1))

## ***************************************************************
## *                          Note:                              *
## *    force.ultrametric does not include a formal method to    *
## *    ultrametricize a tree & should only be used to coerce    *
## *   a phylogeny that fails is.ultramtric due to rounding --   *
## *    not as a substitute for formal rate-smoothing methods.   *
## ***************************************************************

mtext("c) GitHub force.ultrametric",line=0,adj=0.1)

Obviously, this bug will be fixed in the next CRAN release of phytools but can already be obtained by updating phytools from GitHub.

Happy holidays.