Adding a tip at random, but at a fixed height or depth

A phytools user/reader recently posted the following question:

“In the add.random function, using an ultrametric tree, is there any way to control the edge lengths? Ideally I wanna constrain where (timewise) the new tips get added. I don't care where they end up relationship- wise, just want to be able to add them all say <10 million years, or <5 million years, what have you."

There is not a way to do this automatically within add.random, but it is fairly straightforward with bind.tip.

Here I deal with the (fairly) straightforward problem of adding a tip randomly to a tree at a fixed depth. Perhaps in another post I will address the issue of adding a tip on a time interval.

library(phytools)
## for our purposes, simulate a tree of depth 100
tree<-pbtree(n=26,tip.label=LETTERS,scale=100)
## add a tip randomly with height above the root of 60:
h<-60
H<-nodeHeights(tree)
ii<-intersect(which(H[,1]<h),which(H[,2]>h))
edges<-tree$edge[ii,2]
depths<-H[ii,2]-h
## pick one at random:
jj<-sample(1:length(edges),1)
obj<-bind.tip(tree,"Randomly added tip",where=edges[jj],position=depths[jj])
plotTree(obj,mar=c(3.1,0.1,0.1,0.1))
axis(1,at=seq(0,100,by=20))
lines(rep(h,2),par()$usr[3:4],lty="dashed")

We can write a silly little function to repeat this, though here I will give it as "depth” (from the present), rather than height above the root:

bind.at.depth<-function(tree,tip.label,depth){
    H<-nodeHeights(tree)
    h<-max(H)-depth
    ii<-intersect(which(H[,1]<h),which(H[,2]>h))
    edges<-tree$edge[ii,2]
    depths<-H[ii,2]-h
    jj<-sample(1:length(edges),1)
    bind.tip(tree,tip.label,where=edges[jj],position=depths[jj])
}

and then replicate it a bunch of times. Here, when I plot them, for fun I'll paint the added edge red:

trees<-replicate(9,bind.at.depth(tree,"Taxon",40),simplify=FALSE)
class(trees)<-"multiPhylo"
par(mfrow=c(3,3))
for(i in 1:9){
    plot(paintBranches(trees[[i]],which(trees[[i]]$tip.label=="Taxon"),
        state="2"),mar=c(3.1,0.1,0.1,0.1),fsize=0.6,lwd=3,
        split.vertical=TRUE)
    axis(1,at=seq(0,100,by=20))
    lines(rep(60,2),par()$usr[3:4],lty="dashed")
}

## no colors provided. using the following legend:
##       1       2 
## "black"   "red"

## no colors provided. using the following legend:
##       1       2 
## "black"   "red"

## no colors provided. using the following legend:
##       1       2 
## "black"   "red"

## no colors provided. using the following legend:
##       1       2 
## "black"   "red"

## no colors provided. using the following legend:
##       1       2 
## "black"   "red"

## no colors provided. using the following legend:
##       1       2 
## "black"   "red"

## no colors provided. using the following legend:
##       1       2 
## "black"   "red"

## no colors provided. using the following legend:
##       1       2 
## "black"   "red"

## no colors provided. using the following legend:
##       1       2 
## "black"   "red"

Or, with even more replicates!

trees<-replicate(36,bind.at.depth(tree,"Taxon",40),simplify=FALSE)
class(trees)<-"multiPhylo"
par(mfrow=c(6,6))
for(i in 1:36){
    plot(paintBranches(trees[[i]],which(trees[[i]]$tip.label=="Taxon"),
        state="2"),ftype="off",lwd=2,split.vertical=TRUE,
        colors=setNames(c("blue","red"),1:2))
    lines(rep(60,2),par()$usr[3:4],lty="dashed")
}

That's all there is to it!

Average trees and maximum clade credibility trees

I have posted a couple of times on so-called average trees. This is the practice of identifying the phylogeny with the minimum sum of squares distances to the trees in a set. This is something I got into primarily because I was interested in summarizing distributions of trees from comparative methods such as stochastic character mapping.

Though I stumbled into it by accident, I'm discovering that there is a literature on average trees - although in phylogenetic biology, at least, it doesn't seem to be particularly large (compared to that on the different task of computing consensus trees). One significant contribution is an article dating to 1997 by Lapointe & Cucumel entitled “The Average Consensus Procedure: Combination of Weighted Trees Containing Identical or Overlapping Sets of Taxa.” This paper describes a procedure that I independently discovered (20 years later) and have implemented in the phytools function ls.consensus.

According to this consensus method, the 1 through n for n trees patristic distance matrices, D_i are computed and averaged to produce D, then a tree & branch lengths are identified in which the sum of squares difference between D and the patristic distance matrix of the consensus tree, D_consensus is minimized.

Here is a quick demo using the sample of really bad primate trees from before:

library(phytools)
library(phangorn)
ls.con<-ls.consensus(trees)

## Best Q = 0.0052509716725781

## Solution found after 1 set of nearest neighbor interchanges.

plotTree(ls.con<-midpoint(ls.con))

Now, average trees have been criticized - not least because (evidently) they are not guaranteed to have the Pareto property, meaning that it is not guaranteed that every bipartition present in the consensus can be found in at least one of the input trees!

However, in this simple case I thought it might be interesting to see how our average tree by this method compared to, say, a maximum clade credibility tree from the same input set. Note that this is not a maximum clade credibility from a Bayesian posterior distribution - however these trees do come from a probability distribution with expectation equal to the true tree and broad variance, so this could be equivalent to a sample from the posterior distribution under some circumstances, such as when there is a lot of uncertainty regarding phylogenetic relationships in the data.

We can compute the maximum clade credibility tree using phangorn - however to do this one will have to install the latest phangorn version from GitHub.

mcc<-maxCladeCred(trees,rooted=FALSE)
nn<-fastMRCA(mcc,"Lemur","Tarsier")
plotTree(mcc<-reroot(mcc,nn,mcc$edge.length[which(mcc$edge[,2]==nn)]/2))

Topologically this tree is a pretty good tree, but it is not the correct tree.

(To see how it is wrong, we can do this:

obj<-cophylo(ls.con,mcc)

## Rotating nodes to optimize matching...
## Done.

plot(obj)

Cool, right?)

Something that we usually cannot do is compute the clade credibility of a tree not in the set that we have drawn from our (posterior) probability distribution. However, here we have a tree who's credibility is of interest, and so we can compute its clade credibility* simply by tacking it on to the end of our "multiPhylo" object. (*Note that this is not precisely the clade credibility because we have included the clades from our reference tree in the calculation of the clade frequencies - but it is pretty close.)

obj<-list(unroot(ls.con))
obj<-c(trees,obj)
cc<-maxCladeCred(obj,tree=FALSE,rooted=FALSE)
obj<-hist(cc[1:100],breaks=seq(-34,-8,by=2),col="gray",main="",
    xlab="log(clade credibility)")
arrows(x0=cc[101],y0=0.5*max(obj$counts),x1=cc[101],y1=0,col="red",
    lwd=2,length=0.15,angle=20)
text(x=cc[101],y=0.5*max(obj$counts),
    "clade credibility of \"average\" tree",adj=c(0,0.3),srt=90)

So that means the average tree has a higher probability (credibility) than any tree in the set. Interesting.

One could reasonably point out that if we'd merely sampled more trees we would have surely sampled the tree with highest credibility, and then we would have found it to have the same credibility as our LS consensus tree. That's probably true in this case - but, remember, the number of trees for many taxa is vast and so it may not be the case - depending on the shape of our posterior density - that we sample the highest probability topology in our MCMC. (I'm not saying that we don't - just that it cannot be assumed to be guaranteed.)

Drawing a line on (& slicing) a tree to create k subtrees

A phytools user recently asked the following:

I know this is very old post, but I am struggling to make a vertical line in my phylo plot for a given cut. The idea is drawing the line to show a cut of the dendrogram in k groups, rather than a given height.

OK, so we want to show a vertical line on a rooted tree at a point that demarcates k descendant subgroups? This is pretty easy, in fact. We can take advantage of the phytools function ltt as follows:

library(phytools)
tree

## 
## Phylogenetic tree with 26 tips and 25 internal nodes.
## 
## Tip labels:
##  A, B, C, D, E, F, ...
## 
## Rooted; includes branch lengths.

obj<-ltt(tree,plot=FALSE)
obj

## Object of class "ltt" containing:
## 
## (1) A phylogenetic tree with 26 tips and 25 internal nodes.
## 
## (2) Vectors containing the number of lineages (ltt) and branching times (times) on the tree.
## 
## (3) A value for Pybus & Harvey's "gamma" statistic of 0.4726, p-value = 0.6365.

## set k
k<-6
## find the height when there are k lineages
h<-obj$times[which(obj$ltt==k)]
## plot our tree
plotTree(tree,lwd=1,mar=c(4.1,1.1,1.1,1.1))
axis(1)
lines(rep(h,2),par()$usr[3:4],col="red",lty="dashed")

The only trouble with this tactic is that because we have k subgroups at the event leading to the kth group, our line overlaps exactly the vertical line connecting the clades. Alternatively, we could split the difference between the 6th & 7th event and plot that line!

h<-mean(obj$times[c(which(obj$ltt==k),which(obj$ltt==(k+1)))])
plotTree(tree,lwd=1,mar=c(4.1,1.1,1.1,1.1))
axis(1)
lines(rep(h,2),par()$usr[3:4],col="red",lty="dashed")

Now we can clearly see that our vertical line cuts the tree at a point where there are k=6 descendant subtrees.

We can also physically slice the tree at this point if we are so inclined - using the phytools function treeSlice:

subtrees<-treeSlice(tree,h,trivial=TRUE)
print(subtrees,details=TRUE)

## 6 phylogenetic trees
## tree 1 : 2 tips
## tree 2 : 7 tips
## tree 3 : 4 tips
## tree 4 : 4 tips
## tree 5 : 8 tips
## tree 6 : 1 tips

par(mfrow=c(3,2))
lapply(subtrees,function(x) plotTree.singletons(rootedge.to.singleton(x)))

## Error in if (attr(tree, "order") == "cladewise") cw <- tree else stop("tree must be in \"cladewise\" order."): argument is of length zero

Unfortunately, none of the plotting functions seem to like plotting trees with only one tip so we only get to see five of our six subtrees here. Nonetheless, you get the idea.