Friday, September 20, 2013

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.

6 comments:

  1. That's great! Good call on the vertical scaling parameter.

    ReplyDelete
    Replies
    1. There are circumstances in which you might what to scale non-linearly - for instance if you have 10,000 of one thing & 10 of another. One solution in that case is to log-tranform the translation matrix for phylo.toBackbone.

      Delete
  2. Hi Liam,

    I've just been using phylo.toBackbone with ladderized tree and noticed some strange behavior. It seems that with a ladderized tree, tips are not appopriately spaced in relation to their diversity, with the result that triangles overlap. Do you have any ideas why before I dig into the code?

    You can replicate the issue following the example you give above by adding

    lad.tree <- ladderize(tree)
    tt2 <- phylo.toBackbone(lad.tree, trans)
    plot(tt2)

    ReplyDelete
    Replies
    1. Hi Graham.

      Weird. This seems to have something to do with the object that is produced by ladderize. May be badly conformed - but I'm not sure in what way.

      You can circumvent by doing:

      lad.tree <- ladderize(tree)
      lad.tree <- read.tree(text=write.tree(lad.tree))
      tt2 <- phylo.toBackbone(lad.tree, trans)
      plot(tt2)

      - Liam

      Delete
    2. Cool - thanks Liam. Ladderized trees are a little problematic in other ways in my experience, so i suspect you're right (although I also don't know why)

      Delete
  3. The blog is absolutely fantastic.

    ReplyDelete