Sunday, January 31, 2016

Plotting continuous traits on the tree using varied dot sizes at the tips

There may be code in another package to automate this more easily, but in response to a user question I thought I'd give a quick demo on one way that we can plot circles to represent the phenotypic trait values of species at the tips of a tree.

This example uses draw.circle in the plotrix library, and it might be necessary to modify it in various ways for your own purposes.

I have written a custom function for this - called dotTree. The kinks probably still need to be ironed out. Here is a demo.

First load libraries & simulate some data

## load libraries
library(plotrix)
library(phytools)
## to simulate data
tip.labels<-paste(LETTERS,replicate(26,
    paste(sample(letters,round(runif(1,4,8))),collapse="")),
    sep="._") ## realistic looking tip labels
tree<-pbtree(n=26,tip.label=tip.labels,scale=1)
x<-fastBM(tree,sig2=1)

Now the functions themselves:

dotTree<-function(tree,x,...){
    ## plot tree
    tree<-reorder(tree,"cladewise")
    plotTree(tree,offset=1.7,ylim=c(-1/25*Ntip(tree),
        Ntip(tree)),...)
    ## get last phylo plot parameters
    obj<-get("last_plot.phylo",envir=.PlotPhyloEnv)
    x.tip<-obj$xx[1:obj$Ntip]
    y.tip<-obj$yy[1:obj$Ntip]
    ## in case any x<0
    min.x<-min(x)
    max.x<-max(x)
    if(any(x<0)) x<-x-min(x)
    ## plot points
    x<-x[tree$tip.label]
    draw.circle(x.tip+1.2*strwidth("W"),y.tip,nv=200,
        radius=(0.8*x/max(x)+0.1)/2*diff(par()$usr[1:2])/
        diff(par()$usr[3:4]),col="blue")
    ## add legend
    circle.legend(x=par()$usr[1]+0.1*max(nodeHeights(tree)),
        y=0.1*(1+par()$usr[3]),min.x,max.x,length=5,...)
}

circle.legend<-function(x,y,min,max,length=4,prompt=FALSE,...){
    if(hasArg(cex)) cex<-list(...)$cex
    else cex<-1
    if(prompt){
        obj<-locator(1)
        x<-obj$x
        y<-obj$y
    }
    text(x,y-0.5,round(min,2),pos=1,cex=cex)
    s<-(0.8*max(min,0)/min(max,max+min)+0.1)/2*diff(par()$usr[1:2])/
        diff(par()$usr[3:4])
    e<-(0.8/2*diff(par()$usr[1:2])+0.1)/diff(par()$usr[3:4])
    rr<-seq(s,e,length.out=length)
    temp<-c(0,cumsum(1.1*rep(2*max(rr),length-1)))
    draw.circle(x+temp,rep(y,length),nv=200,radius=rr,col="blue")
    text(max(x+temp),y-0.5,round(max,2),pos=1,cex=cex)
    y1<-0.1/25*Ntip(tree)
    lines(c(x,max(x+temp)),rep(y-0.5-y1,2))
    lines(c(x,x),y-c(y1+0.5,2*y1+0.5))
    lines(c(max(x+temp),max(x+temp)),y-c(y1+0.5,2*y1+0.5))  
}

Now let's try it:

dotTree(tree,x,ftype="i")

plot of chunk unnamed-chunk-3

The visualization seems to work OK in this sized plotting device up to about 60 or 70 tips - after that a different method (or a bigger plotting device) is needed. For instance:

tree<-pbtree(n=60,scale=1)
x<-fastBM(tree,a=10,sig2=2)
dotTree(tree,x,lwd=1,fsize=0.6,cex=0.6)

plot of chunk unnamed-chunk-4

That's it.

Friday, January 29, 2016

phylo.heatmap standardizing by column

In tinkering with the visualization function I just added to phytools to plot a phylogenetic heat map, I realized that an important limitation of the function as written was that all columns of X are plotted on the same scale. That means that for traits with low among-species variability, it might be difficult to detect otherwise strong phylogenetic patterns in the data.

For example, here I will simulate a tree & then trait data for 10 traits with various values of σ.

library(phytools)
tree<-rcoal(n=30)
X<-sapply(setNames((1:10)^2,paste("sig=",1:10,sep="")),
    function(sig2,tree) fastBM(tree,sig2=sig2),tree=tree)
phylo.heatmap(tree,X)

plot of chunk unnamed-chunk-1

The problem with this is fairly obvious. Traits with low σ show little evidence for phylogenetic structure because it is just not captured by the color ramp that is specified.

An easy solution to this is to first standardize the columns to have the same variance and then re-run the analysis. Here is what that looks like:

phylo.heatmap(tree,X,standardize=TRUE)

plot of chunk unnamed-chunk-2

I also added some other updates to the legend placement & the labels. The update can be seen here or installed directly from GitHub using detools.

New function for phylogenetic heat map

I just added a preliminary version of a new function to phytools for plotting a “phylogenetic heat map.” This visualization is roughly akin to the R stats function heatmap: although instead of plotting a dendrogram next to the plotted data, it plots an input tree. This was something requested of me by a colleague Nate Swenson some time ago, but unfortunately I didn't get around to trying it until today.

Here's a quick demo of how it works. To follow, one would have to install the latest version of phytools from GitHub:

library(devtools)
install_github("liamrevell/phytools")

Now, let's try it:

library(phytools)
packageVersion("phytools")
## [1] '0.5.13'
## simulate a tree & some data
tree<-pbtree(n=26,tip.label=LETTERS[26:1])
X<-fastBM(tree,nsim=12)
phylo.heatmap(tree,X)

plot of chunk unnamed-chunk-2

A few different attributes can be customized about the plot, and I'm open to making it more flexible if there is interest. For instance, it is quite straitforward to change the palette:

## terrain colors
phylo.heatmap(tree,X,colors=terrain.colors(20))

plot of chunk unnamed-chunk-3

## interpolated blue->red
colors<-colorRampPalette(colors=c("blue","red"))(100)
phylo.heatmap(tree,X,colors=colors)

plot of chunk unnamed-chunk-3

It also works & looks pretty good for non-ultrametric trees:

tree<-rtree(n=40)
X<-fastBM(tree,nsim=20)
colnames(X)<-paste("trait",1:20)
phylo.heatmap(tree,X,fsize=c(0.8,0.8,1))

plot of chunk unnamed-chunk-4

That's it for now.