Slightly faster multiple pairwise RF distance

I just posted a slightly faster version of multiRF that takes advantage of the fact that the first element in the list of tree bipartitions from prop.part of the ape package is the "partition" with all the taxa in the tree - which can be ignored if our two input trees have the same taxa (they should). The result is about a 20% speed-up in the 100 × 99 ÷ 2 comparisons for 100, 100-taxon trees.

Here's a demo using my i3 Asus Vivobook (obviously - it would be much faster on a speedy desktop):

> library(phytools)
> library(phangorn)
> packageVersion("phytools")
[1] ‘0.3.53’
> ## circumvents rNNI from phangorn which is buggy
> ## in the current version
> randomNNI<-function(tree,moves=1){
  for(i in 1:moves) tree<-sample(nni(tree),1)[[1]]
  tree
 }
> ## random tree
> tree<-rtree(n=100,rooted=FALSE,br=NULL)
> ## 100 trees random NNIs
> trees<-lapply(round(runif(n=100,min=1,max=100)), function(x,y) randomNNI(y,x),y=tree)
> class(trees)<-"multiPhylo"
> ## this is a hack
> trees<-read.tree(text=write.tree(trees))
> ## compute all pairwise RF
> system.time(D1<-multiRF(trees))
   user  system elapsed
  16.41    0.02   16.58
> ## now compare to new version
> detach("package:phytools",unload=TRUE)
> install.packages("phytools_0.3-60.zip",repos=NULL)
package ‘phytools’ successfully unpacked and MD5 sums checked
> library(phytools)
> packageVersion("phytools")
[1] ‘0.3.60’
> system.time(D2<-multiRF(trees))
   user  system elapsed
  13.28    0.01   13.29
> plot(D1,D2)

(Just to verify that it still works!)

BTW, I just submitted phytools 0.3-60 to CRAN - so hopeful it is accepted & percolates through the mirror repositories this week.

Getting Robinson-Foulds distances for a set of trees

The phangorn function RF.dist - for computing the Robinson-Foulds distance between two trees - is very fast. However if we want to compare a set of trees stored in an object of class "multiPhylo" there is some opportunity to recycle internal calculations & thus improve on just calling RF.dist n × (n - 1) / 2 times for n trees. Specifically, here is some code that scales to multiple trees a little bit better than n × (n - 1) / 2 RF.dist by calling the ape function prop.part (which finds all bipartitions in the tree) only once for each tree:

multiRF<-function(trees){
 if(class(trees)!="multiPhylo")
  stop("trees should be an object of class \"multiPhylo\"")
 N<-length(trees)
 RF<-matrix(0,N,N)
 if(any(sapply(unclass(trees),is.rooted))){
  cat("Some trees are rooted. Unrooting all trees.\n")
  trees<-lapply(unclass(trees),unroot)
 }
 foo<-function(pp) lapply(pp,function(x,pp)
  sort(attr(pp,"labels")[x]),pp=pp)
 xx<-lapply(unclass(trees),function(x) foo(prop.part(x)))
 for(i in 1:(N-1)) for(j in (i+1):N)
  RF[i,j]<-RF[j,i]<-2*sum(!xx[[i]]%in%xx[[j]])
 RF
}

It works OK compared to calling RF.dist a whole bunch of times independently. First, here is an alternate version of the function which uses RF.dist:

multiRF2<-function(trees){
 N<-length(trees)
 RF<-matrix(0,N,N)
 for(i in 1:(N-1)) for(j in (i+1):N)
  RF[i,j]<-RF[j,i]<-RF.dist(trees[[i]],trees[[j]])
 RF
}

Now, a very simple comparison with random trees:

> N<-100
> n<-100
> trees<-rmtree(N,n,rooted=FALSE)
> system.time(X1<-multiRF(trees))
   user  system elapsed
   9.16    0.03    9.19
> system.time(X2<-multiRF2(trees))
   user  system elapsed
  19.98    0.00   19.99
> all(X1==X2)
[1] TRUE

Even more on plotting subtrees as triangles on a phylogenetic backbone tree

This is pretty much guaranteed to be the last post on this topic today, but I made some additional refinements of plot.backbonePhylo, my new plotting method to show subtrees as triangles.

The three main changes are as follows: (1) following Luke's suggestion, subtrees containing only 1 taxon are now plotted as a line (actually - a triangle with zero width, but I guarantee you won't be able to tell the difference), instead of as a triangle with unit width; (2) following Graham Slater's suggestion, I made a small change so that there is a space between adjacent triangles - the consequence is that clades with 2 taxa are of unit length, clades with 10 taxa are 9 units in width, etc., which actually lines up with tips containing 1 taxon, which now have unit width; and (3) I added the optional scaling factor vscale, which scales the vertical dimension so that large subtrees can be portrayed while still allowing singletons and spacing between triangles.

All of these updates are available here, but are also part of a new phytools build (phytools 0.3-51).

Here's a demo:

> library(phytools)
> packageVersion("phytools")
[1] ‘0.3.51’
> library(geiger) # to help us build our tree
> ## now let's create our backbone tree with
> ## random subtree diversities
> tree<-transform(pbtree(n=10),model="lambda",lambda=0.5)
> ## for old versions of geiger, use lambdaTree
>
> ## create a translation table
> ## leaving a couple of single-taxon clades for fun
> tip.label<-sample(tree$tip.label,8)
> clade.label<-LETTERS[1:8]
> N<-ceiling(runif(n=8,min=1,max=20))
> ## set crown node depth to 1/2 the maximum depth
> depth<-sapply(tip.label,function(x,y) 0.5*y$edge.length[which(tree$edge[,2]==which(y$tip.label== x))],y=tree)
> trans<-data.frame(tip.label,clade.label,N,depth)
> rownames(trans)<-NULL
> rm(tip.label,clade.label,N,depth)
>
> ## here's what trans looks like
> trans
  tip.label clade.label  N     depth
1        t4           A  3 1.0079214
2        t8           B 13 0.8367549
3        t1           C 15 1.1433276
4        t3           D  8 1.0129504
5        t6           E 12 0.9687320
6        t5           F  2 0.9687320
7        t9           G 14 0.7253729
8       t10           H 17 0.7253729
> 
> ## convert
> tt<-phylo.toBackbone(tree,trans)
>
> ## plot
> plot(tt)

> ## now let's make our clades much bigger
> trans$N<-trans$N*10
> trans
  tip.label clade.label   N     depth
1        t4           A  30 1.0079214
2        t8           B 130 0.8367549
3        t1           C 150 1.1433276
4        t3           D  80 1.0129504
5        t6           E 120 0.9687320
6        t5           F  20 0.9687320
7        t9           G 140 0.7253729
8       t10           H 170 0.7253729
> tt<-phylo.toBackbone(tree,trans)
> plot(tt)

> ## now let's use the new parameter vscale
> plot(tt,vscale=0.1)

Pretty cool.

We have to be careful with this, though, because if we adjust our clade widths such that any are less than 1 unit then we are going to get some really weird plots.

More on plotting backbone trees with triangles as subtrees

Yesterday (more accurately, early this morning), I posted about a new plotting method for plotting a backbone phylogeny with diversities for each subtree represented as a triangle. This type of visualization is quite common in the literature; however, I'm told (by Luke Mahler, so blame him if this is incorrect) that there is no such plotting method so far in R.

The way I did this was by first creating a new type of phylogeny object, of class "backbonePhylo". This object is stored in a very similar way to the typical phylogeny object in R (class "phylo"), which the exception that we have replaced the vector tip.label with the list tip.clade. tip.clade is a list in which each element contains a label, a diversity (i.e., the number of tips hidden in the subtree), and a crown depth.

I have written some of the methods around this object as S3 methods, which means they can be dispatched using (for instance) calls to the generic plot, reorder, print, etc.

This morning, I did a little more work on this suite of functions. Most importantly, I created a new function (phylo.toBackbone) that converts an object of class "phylo" to an object of class "backbonePhylo" with the help of a translation table. The translation table is a data frame with four variables: tip.label (containing the names of the tips in the input tree - these might be exemplar taxa); label (containing the names of the subtrees - these might be the names of higher taxa); N (obviously, the diversities of each subtree); and depth (the age of each crown group, in terms of the branch lengths of the input tree - note that this cannot be longer then the corresponding tip edge). Note that not all tips in the original tree need to be in the translation table. Those that are not will be treated as clades containing one taxon in the backbone tree.

Here's a demo using random subtree sizes on a stochastic backbone tree. New code for this is here.

> library(phytools)
> library(geiger) # to help us build our tree
> source("backbonePhylo.R")
> ## now let's create our backbone tree
> tree<-transform(pbtree(n=10),model="lambda",lambda=0.5)
> ## for old versions of geiger, use lambdaTree
>
> ## create a translation table
> ## leaving a couple of single-taxon clades for fun
> tip.label<-sample(tree$tip.label,8)
> clade.label<-LETTERS[1:8]
> N<-ceiling(runif(n=8,min=1,max=20))
> ## set crown node depth to 1/2 the maximum depth
> depth<-sapply(tip.label,function(x,y) 0.5*y$edge.length[which(tree$edge[,2]==which(y$tip.label== x))],y=tree)
> trans<-data.frame(tip.label,clade.label,N,depth)
> rownames(trans)<-NULL
> rm(tip.label,clade.label,N,depth)
>
> ## here's what trans looks like
> trans
  tip.label clade.label  N    depth
1        t2           A  6 0.8992932
2        t4           B  5 0.7617329
3        t3           C 16 0.8674696
4        t5           D  9 0.5984337
5        t1           E  5 0.8992932
6       t10           F  7 0.5554955
7        t7           G  8 0.5664347
8        t9           H  7 0.5554955
>
> ## convert
> tt<-phylo.toBackbone(tree,trans)
>
> ## plot
> plot(tt)

That's pretty good.

Plotting phylogenies with triangles for subtrees

Here's something new. Luke Mahler asked me if I could create a plotting method that could take a backbone phylogeny and plot subtrees using triangles that are scaled proportionally to the number of species in the subtree. I'm still working on this, but I posted some preliminary code here.

The function takes a modified tree as an object of class "backbonePhylo". I'm going to make some functions to convert from "phylo" objects, but for now this has to be built manually. The "backbonePhylo" object is like a "phylo" object, but with tip.clade substituted for tip.label. tip.clade is a list in which each subtree is represented by a label, a number of descendant taxa, and crown group depth.

Here's a quick demo of what can be done so far:

> library(phytools)
> library(geiger) # to help us build our tree
> source("backbonePhylo.R")
> ## now let's create our backbone tree with
> ## random subtree diversities
> tree<-transform(pbtree(n=10),model="lambda",lambda=0.5)
> tree$tip.clade<-list()
> for(i in 1:10){
  tree$tip.clade[[i]]<-
    list(label=paste("Group",i),
    N=ceiling(runif(n=1,min=0,max=20)),
    depth=tree$edge.length[which(tree$edge[,2]==i)])
  }
> class(tree)<-"backbonePhylo"
> tree

Backbone phylogenetic tree with 10 subtrees and 9 resolved internal nodes.

Labels: Group 1, Group 2, Group 3, Group 4, Group 5, ...
Diversities: 2, 18, 13, 7, 9, ...

> ## ok, let's plot it
> plot(tree)

Cool. That's more or less what we were going for.

Coloring tip labels by character state

Yesterday I received the query:

"I want to color my tree tip labels according to a particular character and provide a legend for the same."

Here's a demo of one way to do this, using the phytools data object anoletree:

require(phytools)
data(anoletree)
tree<-anoletree
x<-getStates(tree,"tips")
par(fg="transparent")
plotTree(tree,fsize=0.6,ylim=c(-1,length(tree$tip.label)))
lastPP<-get("last_plot.phylo",env=.PlotPhyloEnv)
ss<-sort(unique(x))
colors<-setNames(palette()[1:length(ss)],ss)
par(fg="black")
add.simmap.legend(colors=colors,vertical=FALSE,x=0.25,
  y=-1,prompt=FALSE)
colors<-sapply(x,function(x,y) y[which(names(y)==x)],
  y=colors)
tt<-gsub("_"," ",tree$tip.label)
text(lastPP$xx[1:length(tt)],lastPP$yy[1:length(tt)],
  tt,cex=0.6,col=colors,pos=4,offset=0.1,font=3)

Here's the result:

Obviously, some of this is idiosyncratic to this specific dataset. Hopefully, users can figure out how to modify this to their own data and problem.

Simulating DNA sequences with Γ rate heterogeneity among sites

I just added this to the (fairly maligned) function genSeq, for simulating DNA sequences by (inefficiently) wrapping the phytools function sim.history. Here's what the function looks like now:

genSeq<-function(tree,l=1000,Q=NULL,rate=1, format="DNAbin",...){
 if(is.null(Q)){
  Q<-matrix(1,4,4)
  rownames(Q)<-colnames(Q)<-c("a","c","g","t")
  diag(Q)<-0
  diag(Q)<--colSums(Q)
 }
 if(length(rate)!=l){
  if(length(rate)==1) rate<-rep(rate,l)
  else {
   cat("warning: length(rate) & l should match for      length(rate)>1\n")
   cat(" rate will be recycled.\n")
   rate<-rep(rate,ceiling(l/length(rate)))[1:l]
  }
 }
 cat("simulating sequences....\n")
 X<-sapply(rate,function(a,b,c)   sim.history(b,a*c)$states,b=tree,c=Q)
 if(format=="DNAbin") return(as.DNAbin(X))
 else if(format=="phyDat") return(as.phyDat(X))
 else if(format=="matrix") return(X)
}

One can easily see that the way rate heterogeneity is handled is by allowing the user to supply a site-specific vector of rates. Here's what this would look like if we wanted to simulate sequences under a continuous Γ rate heterogeneity model with a shape parameter, α = 0.25 and a mean rate of 1.0 substitution per site per unit of branch length:

> library(phytools)
> tree<-pbtree(n=26,tip.label=LETTERS,scale=0.1)
> gg<-rgamma(n=1000,shape=0.25,rate=0.25)
> X<-genSeq(tree,rate=gg)
simulating sequences....
> X
26 DNA sequences in binary format stored in a matrix.

All sequences of same length: 1000

Labels: A B C D E F ...

Base composition:
    a    c    g    t
0.264 0.245 0.242 0.249

We can do invariant sites the same way - although in this case it will break the function if we set any of the rates to zero, so we need to set them very slightly > 0:

> ii<-rep(1,1000)
> tol<-1e-12
> ii[sample(1:1000,100)]<-tol
> X<-genSeq(tree,rate=ii)
simulating sequences....
> X
26 DNA sequences in binary format stored in a matrix.

All sequences of same length: 1000

Labels: A B C D E F ...

Base composition:
    a    c    g    t
0.244 0.261 0.238 0.257

As Klaus pointed out, this method is very slow - so you probably have better options. But it's still a fun exercise to see how it works - and offers the interesting possibility that we could store and return all character histories along with the states for tips.