Wednesday, July 10, 2013

Rotating all the nodes in a tree (or a set of nodes)

I just wrote a utility function that wraps around rotate in the ape package to rotate a set of nodes in a tree (including nodes="all"). It also addresses this weird problem that I discovered that results when a bunch of rotations are applied to a tree the "phylo" can become non-compliant with certain ape & phytools functions because the order of the tip numbers in tree$edge is not 1:n for n tips.

Here's the function:

rotateNodes<-function(tree,nodes,polytom=c(1,2),...){
  n<-length(tree$tip.label)
  if(nodes[1]=="all") nodes<-1:tree$Nnode+n
  for(i in 1:length(nodes))
    tree<-rotate(tree,nodes[i],polytom)
  if(hasArg(reversible)) reversible<-list(...)$reversible
  else reversible<-TRUE
  if(reversible){
    ii<-which(tree$edge[,2]<=n)
    jj<-tree$edge[ii,2]
    tree$edge[ii,2]<-1:n
    tree$tip.label<-tree$tip.label[jj]
  }
  return(tree)
}

Let's try it:

> tree<-pbtree(n=26)
> tree$tip.label<-LETTERS
> plotTree(tree)
> tree<-rotateNodes(tree,"all")
> plotTree(tree)

Adding a scale bar to a partial circular tree

Here's a neat trick that I thought of yesterday. We can pretty easily add a scale bar to a partial circular tree created using plotTree(...,type="fan") in the phytools package, and here's how:

> library(phytools)
> packageVersion("phytools")
[1] ‘0.2.93’
> # for this demo - simulate a tree of length 100
> tree<-pbtree(n=100,scale=100)
> # plot it
> # (you'll probably need to tweak 'part' to fit the axis)
> plotTree(tree,type="fan",part=0.93,fsize=0.9)

Note: type="fan" is in development.
Many options of type="phylogram" are not yet available.

> # ok, now plot the axis
> h<-max(nodeHeights(tree))
> axis(1,pos=-0.05*h,at=seq(0,h,by=h/2),lwd=2)
> text(x=0.5*h,y=-0.25*h,"time")

That's pretty cool.

Phylomorphospace with time since the root projected using a color gradient

Two different people (Marcio Pie and a student whose name has somehow slipped my mind - please identify yourself!) independently suggested to me recently that I could use a similar approach to that employed when overlaying a posterior density map from stochastic mapping on a phylomorphospace to project time since the root of the tree onto a phylomorphospace. In theory, this would help overcome the difficulty that all temporal information about the phylogeny is lost when it is projected into morphospace.

Well, this is quite easy** (**that is, easy only because I've already programmed all the pieces in phytools) to do. Here's a hack using simulated data.

First let's load phytools & simulate tree & data:

> require(phytools)
Loading required package: phytools
> packageVersion("phytools")
[1] ‘0.2.93’
> # simulate tree & data
> tree<-pbtree(n=30,scale=100)
> X<-fastBM(tree,nsim=2)
Next, we're going to create an object of class "contMap". This is just a placeholder into which we'll slot our temporal data:
# unfortunately, there's no way to prevent this from plotting
AA<-contMap(tree,X[,1])
Now, we'll swap our color palette (although we don't need to - this is just an arbitrary decision); and then we'll substitute information about the height from the root for our trait data mapped onto the tree:
> AA$cols[]<-rainbow(1001,start=0.7,end=0)
> H<-nodeHeights(tree)
> h<-max(H)
> for(i in 1:nrow(H)) names(AA$tree$maps[[i]])<- round((H[i,1]+cumsum(AA$tree$maps[[i]]))/h*1000)
> # check to verify that temporal information is correct
> plot(AA,legend=FALSE)
Finally, let's plot our phylomorphospace with color as a temporal axis:
> phylomorphospace(AA$tree,X,colors=AA$cols,lwd=3, node.by.map=TRUE,xlab="trait 1",ylab="trait 2")
> add.color.bar(2,AA$cols,title="time from the root", lims=c(0,h),digits=1)
Click where you want to draw the bar
(Click here for full size.)

The final step (add.color.bar(...)) is the biggest hack of all as it only works because I artfully positioned it close enough to the x axis so that the bottom legend text is hidden!

Try it out!

Monday, July 8, 2013

New version of plotBranchbyTrait with full user control of plotting options

Today I received the following request from a colleague regarding the phytools function plotBranchbyTrait:

A PhD student of mine is trying to use your function plotBranchbyTrait to plot a phylogeny where the branch colors reflect the rate of molecular evolution .... He has managed to plot the tree with the branch colors and that part works very nicely, but when trying to offset the tip labels (species names) or modifying the edge.width he gets a message saying that the argument "label.offset" (or "edge.width") is matched by multiple actual arguments.... Can you assist us as why you think we're getting these messages? Any tricks to offset the species names?

plotBranchbyTrait differs from other similar functions in phytools (such as densityMap and contMap) in that it requires a single input value for the mapped trait per edge, and then plots that trait value on a color gradient scale - but uniformly for each branch.

The truth is, I developed this function primarily to help a colleague - so it is not endowed with the full range of plotting options. Unlike the vast majority of tree plotting functions in phytools, this one calls plot.phylo from the ape package internally - rather than plotSimmap in phytools.

This afternoon, I took a the twenty minutes (or so) required to pass control of all* (*with the exception of a couple, such as edge.color, which don't make sense) to the user. In so doing, I've kept the non-optional arguments of plotBranchbyTrait the same, although I've changed a few of the internal defaults - just to make for nicer plots. (For instance, I discovered that the default x & y limits meant that the legend was being plotting out of bounds in many instances.)

The updated code for plotBranchbyTrait is here, but I have also posted a new version of phytools (phytools 0.2-93), which can be downloaded and installed from source. (For more information on installing packages from source, search the web.)

Here's the usual demo:

> library(phytools)
> packageVersion("phytools")
[1] ‘0.2.93’
> tree<-pbtree(n=50)
> x<-fastBM(tree)
> plotBranchbyTrait(tree,x,mode="tips",edge.width=4, prompt=TRUE)
Click where you want to draw the bar

For fun, we can compare to contMap, which interpolates smoothly along each internal & terminal edge. (Of course, plotBranchbyTrait can also be used for any arbitrary trait in which we have a state on a continuous scale at each edge of the tree.) These two plots should look very similar - it is just a matter of whether or not our eyes can easily distinguish between continuously & abruptly changing colors.

> XX<-contMap(tree,x,lims=c(-3.5,1.5))
> XX$cols[]<-rainbow(1001,start=0.7,end=0)
> plot(XX,outline=F,lwd=4,legend=0.82)

Virtually the same to my eye!

Sunday, July 7, 2013

Plotting 'partial' circular trees

At this year's Evolution meeting in Snowbird, Utah, I saw a very nice poster by Patrick Fuller of the Wainwright Lab at UCDavis. One cool attribute of the poster was that he had a half circular phylogeny aligned to the bottom of this poster. My immediate thought, naturally, was that phytools should do this.

I have now added this capability to the functions plotSimmap (as well as plotTree, which uses plotSimmap internally). I have included these updates in a new non-CRAN release of phytools (phytools 0.2-91), which can be downloaded and installed from source.

Doing this was not hard. Basically, when we are deciding the circular positions of all the terminal edges of the tree, we start by ordering them cladewise, and then we evenly space the tips from 0 through 2π before wrapping them into our circular space. (I describe this a little more here.) To plot a partial fan tree, we just decide what fraction, θ, of the circle we want to use, then we space our terminal edges, instead of evenly from 0 to 2π, evenly from 0 through 2πθ. Aside from resizing our plot axes to match - that's pretty much all there is to it.

Here's a quick demo of the results:

> library(phytools)
> packageVersion("phytools")
[1] ‘0.2.91’
> data(anoletree)
> plotSimmap(anoletree,type="fan",part=0.5,fsize=0.9, ftype="i")
no colors provided. using the following legend:
    CG     GB      Non-     TC     TG        Tr        TW
"black"  "red"  "green3" "blue" "cyan" "magenta"  "yellow"

Note: type="fan" is in development.
Many options of type="phylogram" are not yet available.

> ss<-sort(unique(getStates(anoletree,"tips")))
> add.simmap.legend(colors=setNames(palette()[1:length(ss)],ss))
Click where you want to draw the legend
(Click here for full resolution.)

Pretty cool. That's about what we were going for. Here's a smaller, simulated 1/4 fan tree - again, with a mapped discrete character:

> Q<-matrix(c(-1,1,1,-1),2,2)
> tree<-sim.history(pbtree(n=40,scale=1),Q,anc=1)
> plotSimmap(tree,colors=setNames(c("blue","red"),1:2), type="fan",part=0.25,lwd=4)

Note: type="fan" is in development.
Many options of type="phylogram" are not yet available.

> add.simmap.legend(colors=setNames(c("blue","red"),1:2))
Click where you want to draw the legend

That's it.

Thursday, July 4, 2013

phytools page now hosted on 3rd party server

Happy 4th of July - the faculty pages web server at UMass-Boston has been suffering a lot of service interruptions over the past few days, so I have now moved the phytools page (http://www.phytools.org) to a new 3rd party provider. Hopefully that will take care of these service interruptions - and, of course, phytools is always available through CRAN.

Correction: the link above to http://www.phytools.org originally pointed to the wrong web address. Fixed.

Monday, July 1, 2013

Static help pages for phytools

The html help pages in R are created dynamically from the .Rd manual pages; however it is also possible to create static .html help pages to post on the web (or have available for review without R open). I have done this now & posted them to the phytools page. For example, the static help page for make.simmap can be seen here. You can also see an index of all functions available in the latest version of phytools. Check it out.

For more information about building static html help pages in R, see the following link.