Monday, February 13, 2017

Plotting up- and down-facing "contMap" & "densityMap" trees

Yesterday, I posted on plotting upward & downward facing trees with plotSimmap (and the wrapping function plotTree).

The reason I programmed this is to eventually extend this functionality to the methods that use plotSimmap internally - a long list that includes such functions as plot.contMap (& contMap) and plot.densityMap (& densityMap). These functions seem to be quite popular so I thought I'd do them first.

The biggest pain in the butt for these methods was not getting the phylogeny oriented correctly - after yesterday's phytools update, that part was easy. The hard part was getting the color-bar legend correctly positioned & sized. The reason this turned out to be challenging was because our plotting region is not square - meaning that x & y dimensions can have different scales. When our legend & branch lengths point in the same direction (right-to-left or left-to-right) we have nothing to worry about; however if we want to plot our tree upward or downward, but our legend (for the vertical branches) horizontally, then we start to have problems. To make this happen, we have to take advantage of par("pin") which gives us the relative dimensions of our plotting area in x & y - but in inches, not in our plotting units (these are given by par("usr"), by the way).

The updates can be seen on GitHub here.

Let's try a demo. Here I have simulated a simple tree, but added one longish tip label to make sure that the sizing of our plotting area to accomodate tip labels is working properly. (As mentioned in an earlier post, this is one of the other challenging thing about plotting phylogenies.)

library(phytools)
packageVersion("phytools")
## [1] '0.5.76'
tree
## 
## Phylogenetic tree with 26 tips and 25 internal nodes.
## 
## Tip labels:
##  A longish species name, B, C, D, E, F, ...
## 
## Rooted; includes branch lengths.
x
## A longish species name                      B                      C 
##             -1.3444974             -2.7182556             -3.6415576 
##                      D                      E                      F 
##             -2.7367692             -2.9094996              1.0659075 
##                      G                      H                      I 
##              0.4134942              0.6560451              0.9139865 
##                      J                      K                      L 
##              1.3646047              1.3643494              0.4296533 
##                      M                      N                      O 
##              0.9882680              0.8843809              0.7174877 
##                      P                      Q                      R 
##              1.0307320              1.3304640              2.8291197 
##                      S                      T                      U 
##              0.6508031              0.1046013              0.7843597 
##                      V                      W                      X 
##              0.8382973              0.4644528              0.5338497 
##                      Y                      Z 
##              2.5244863              2.2260855
obj<-contMap(tree,x,plot=FALSE)
obj
## Object of class "contMap" containing:
## 
## (1) A phylogenetic tree with 26 tips and 25 internal nodes.
## 
## (2) A mapped continuous trait on the range (-3.641558, 2.82912).
plot(obj,direction="upwards")

plot of chunk unnamed-chunk-1

Next, using densityMap:

trees<-make.simmap(tree,y,nsim=100)
## make.simmap is sampling character histories conditioned on the transition matrix
## 
## Q =
##            a          b
## a -0.9922819  0.9922819
## b  0.9922819 -0.9922819
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##   a   b 
## 0.5 0.5
## Done.
obj<-densityMap(trees,plot=FALSE)
## sorry - this might take a while; please be patient
obj
## Object of class "densityMap" containing:
## 
## (1) A phylogenetic tree with 26 tips and 25 internal nodes.
## 
## (2) The mapped posterior density of a discrete binary character with states (a, b).
plot(obj,direction="downwards",lwd=5,outline=TRUE)

plot of chunk unnamed-chunk-2

obj<-setMap(obj,c("white","black"))
plot(obj,direction="upwards",lwd=5,outline=TRUE)

plot of chunk unnamed-chunk-2

Neat.

The data for this exercise were simulated as follows:

tree<-pbtree(n=26,tip.label=LETTERS)
tree$tip.label[1]<-"A longish species name"
x<-fastBM(tree)
Q<-matrix(c(-1,1,1,-1),2,2)
rownames(Q)<-colnames(Q)<-letters[1:2]
y<-sim.history(tree,Q)$states

Sunday, February 12, 2017

Plotting up- & down-facing phylograms in phytools::plotSimmap

I just pushed a significant update to the workhorse phytools plotting function plotSimmap to permit plotting in both up- and down-facing square phylograms, in addition to the right & left phylograms and the circular fan trees already available.

When I call plotSimmap a “workhorse,” what I mean is that it is used internally by a whole bunch of other functions & methods in the package. For instance, in addition to the S3 plot method for "simmap" and "multiSimmap" object classes, it is also used by bind.tip (in interactive mode), contMap (& plot.contMap), densityMap (& plot.densityMap), densityTree, dotTree (for only one character), fancyTree (for various methods), ltt, plot.describe.simmap, plotTree, plotTree.wBars (along with both plotTree.barplot, and plotTree.boxplot), reroot (for interactive mode), and treeSlice (in prompt=TRUE mode), among other things. Adding the up- & down-facing directions to this method will eventually allow me to add this capability to many of these functions as well.

plotSimmap is also used inside the new functions of the physketch package, such as in the phylo.tracer function to draw a "phylo" object on top of the image of a plotted tree. When I extend the direction argument to this function it will be possible to trace a phylogeny plotted in any direction.

Here's a quick demo of up- & down-facing trees:

library(phytools)
packageVersion("phytools")
## [1] '0.5.75'
plotTree(tree,direction="upwards",mar=c(0.1,4.1,0.1,0.1))
axis(side=2,at=seq(0,100,by=20))
title(ylab="time since the root")

plot of chunk unnamed-chunk-1

We can flip the direction of the axis - that is, plotting time from the present instead of time since the root - by plotting our tree downwards, but then switch ylim{1] & ylim{2].

For instance:

obj<-get("last_plot.phylo",envir=.PlotPhyloEnv)
ylim<-obj$y.lim[2:1]-(obj$y.lim[2]-max(nodeHeights(tree)))
plotTree(tree,direction="downwards",mar=c(0.1,4.1,0.1,0.1),ylim=ylim)
axis(side=2,at=seq(0,100,by=20))
title(ylab="time since the present")

plot of chunk unnamed-chunk-2

One (surprising?) challenge with tree plotting is ensuring that there is enough, but not too much, space left for tip labels. The reason this is such a pain is because we don't know how much space the labels will take up in the units of our tree until we've already plotted the tree! Let's make one crazy long label to see if I got it right this time:

tree$tip.label[1]<-"The quick brown fox jumped over the lazy dog"
plotTree(tree,direction="upwards")

plot of chunk unnamed-chunk-3

There are still some bugs to be worked it before I can add this option to all the other functions listed above - so please be patient.

That's it.

Friday, February 10, 2017

Additional function for interactive drawing of a cladogram or ultrametric phylogram

I just wrote another function designed to permit users to interactively draw phylogenies in the R environment.

This one is based closely on a demo I posted to my blog a couple of weeks ago showing how to draw a tree using phytools::bind.tip. This one uses bind.tip in a couple of different ways - either interactively (as in the case for method="cladogram" or non-interactively, but via the interactive phytools (mostly internal - although now exported to the namespace) function get.treepos.

The following is a quick demo:

library(phytools)

tips<-c("lemur","robin","whale","coelecanth","bat","cow",
    "goldfish","pig","iguana","human")
tips
##  [1] "lemur"      "robin"      "whale"      "coelecanth" "bat"       
##  [6] "cow"        "goldfish"   "pig"        "iguana"     "human"
outgroup<-"shark"
outgroup
## [1] "shark"
draw.ultrametric<-function(ingroup,outgroup=NULL,depth=1.0,
    method=c("phylogram","cladogram")){
    method<-method[1]
    if(is.null(outgroup)) out<-"OUTGROUP"
    else out<-outgroup
    tree<-pbtree(n=2,tip.label=c(ingroup[1],out),scale=depth)
    if(method=="cladogram"){
        for(i in 2:length(ingroup)) 
            tree<-bind.tip(tree,ingroup[i],interactive=TRUE)
        tree$edge.length<-NULL
    } else if(method=="phylogram"){
        dev.hold()
        plotTree(tree,mar=c(0.1,0.1,3.1,0.1))
        v<-seq(0,depth,by=depth/10)
        axis(3,at=v)
        abline(v=v,lty="dashed",col=make.transparent("grey",0.7))
        dev.flush()
        cat(paste("Click where you would like to bind the tip \"",
            ingroup[2],"\"\n",sep=""))
        flush.console()
        for(i in 2:length(ingroup)){
            obj<-get.treepos(message=FALSE)
            tree<-bind.tip(tree,ingroup[i],where=obj$where,
                pos=obj$pos)
            dev.hold()
            plotTree(tree,mar=c(0.1,0.1,3.1,0.1))
            axis(3,at=v)
            abline(v=v,lty="dashed",col=make.transparent("grey",0.7))
            dev.flush()
            if(i<length(ingroup))
                cat(paste("Click where you would like to bind the tip \"",
                    ingroup[i+1],"\"\n",sep=""))
            flush.console()
        }
    }
    tree
}
vertebrates<-draw.ultrametric(tips,outgroup)
plotTree(vertebrates)

plot of chunk unnamed-chunk-3

[The above example is based on one I use in undergrad classes in which I try to get them to come up with a 'phylogeny of vertebrate' from the following slide:

from which they are supposed to draw a tree.]

This function will most likely go in the physketch package.

That's all.