Wednesday, June 12, 2019

Aesthetic & functional updates to the function ltt95 for visualizing a lineage-through-time plots for a set of trees

I just pushed some updates & additions (1, 2, 3, 4) to the phytools function ltt95.

ltt95 (e.g., here) is a relatively simple function for drawing (1-α) × 100% CIs around a lineage through time plot. It has two different modes - one in which the mean or median & confidence intervals are computed in terms of the number of lineages at a given time (method="lineages") and the other other in which the mean or median & confidence intervals are calculated for the amount of time that has elapsed, given a certain number of lineages (method="times").

I mainly added some improved visualization - for instance, using a shaded polygon instead of simple dotted lines to demarcate the CIs. I'll be the first to admit that the style is unoriginal. I believe that the paleotree package produces a plot that is stylistically highly similar, if not the same. Nonetheless, the function has some potentially useful features so I thought it was worth doing. I also corrected the lower confidence bounds - previously it was one element too low making the CIs slightly too large.

Here's a demo using some phylogenies of primates:

library(phytools)
primates<-read.nexus(file=
    "http://www.phytools.org/TS2018/data/10kTrees_Primates.nex")
ltt95(primates,log=TRUE)

plot of chunk unnamed-chunk-1

The weird effect where the number of accumulated lineages drops off to zero at the end of the plot is a characteristic sign of some of the trees not being ultrametric to the precision of the machine. This is because phytools::ltt (used internally by ltt95) intentionally does not assume that the trees are ultrametric. In our case, they should be - and the small deviations that we see are merely due to rounding of the edge lengths when our trees were written to file.

To fix that, let's use force.ultrametric as follows:

primates<-lapply(primates,force.ultrametric)
class(primates)<-"multiPhylo"

Now, let's proceed to re-compute our "ltt95" object, this time without plotting:

object<-ltt95(primates,log=TRUE,plot=FALSE)
object
## Object of class "ltt95".
##   alpha: 0.05
##   method:    lineages
##   mode:      median
##   N(lineages):   273
par(bty="l")
plot(object,shaded=TRUE,xaxis="flipped")

plot of chunk unnamed-chunk-3

That's cool. Now let's switch to mode="times", just to see how it looks:

object<-ltt95(primates,log=TRUE,plot=FALSE,method="times")
object
## Object of class "ltt95".
##   alpha: 0.05
##   method:    times
##   mode:      median
##   N(lineages):   273
par(bty="l")
plot(object,shaded=TRUE,xaxis="flipped",bg=rgb(1,0,0,0.25))

plot of chunk unnamed-chunk-4

Finally, here is what it looks like if we compute the average number of lineages at each time point, rather than the median. I'm also going to set α to 0.0 for the confidence intervals - essentially, I want my CI to include all trees in the sample:

object<-ltt95(primates,log=TRUE,plot=FALSE,mode="mean",alpha=0)
par(bty="l")
plot(object,shaded=TRUE,xaxis="flipped")
legend(x="topleft",c("mean number of lineages",
    "total range of lineages across trees"),
    lty=c(1,0),lwd=c(2,0),pch=c(NA,22),pt.cex=2.5,
    pt.bg=c("black",rgb(0,0,1,0.25)),col=c("black","transparent"),
    bty="n")

plot of chunk unnamed-chunk-5

Finally, here are all the LTTs plotted for comparison:

object<-ltt(primates,plot=FALSE)
max.t<-max(sapply(object,function(x) max(x$times)))
for(i in 1:length(object)) 
    object[[i]]$times<-object[[i]]$times+
        diff(c(max(object[[i]]$times),max.t))
plot(object,col=rgb(0,0,1,0.5))

plot of chunk unnamed-chunk-6

Neat.

Monday, June 10, 2019

Using geo.palette( ) to map geological periods on a phylomorphospace plot

Yesterday I added a function to the namespace of phytools called geo.palette that returns the standard geological time scale color palette.

The following gives a demo on how this can be mapped onto a phylomorphospace:

## first compute total tree depth
h<-max(nodeHeights(tree))
## extract the temporal limits of your geological periods, but from
## the root of the tree forward
limits<-h-geo.palette()$period[,"start"]
limits<-limits[limits>0]
limits<-c(0,sort(limits))
names(limits)[1]<-rownames(geo.palette()$period)[length(limits)]
limits
## Carboniferous       Permian      Triassic      Jurassic    Cretaceous 
##         0.000        21.100        67.830       118.700       175.000 
##     Paleogene       Neogene    Quaternary 
##       254.000       296.970       317.412
## pull the colors from geo.palette( )
cols<-setNames(geo.palette()$cols[length(limits):1],
    names(limits))
## create our period-mapped tree using make.era.map
period.tree<-make.era.map(tree,limits)
## now plot our tree
## here I plot the phylomorphospace twice to get a black outline,
## but this is not necessary.
par(lend=1)
phylomorphospace(tree,X,lwd=7,node.size=c(0,0),
    ftype="off",xlab="trait 1",ylab="trait 2")
phylomorphospace(period.tree,X,colors=cols,
    node.size=c(0,0),lwd=5,add=TRUE,ftype="off")
points(X,pch=21,cex=1.5,bg="grey")
legend(x="topright",legend=names(cols)[length(cols):1],
    pch=22,pt.bg=cols[length(cols):1],pt.cex=2)

plot of chunk unnamed-chunk-1

My data for this exercise were simulated as follows:

set.seed(999)
tree<-pbtree(n=40,scale=320)
X<-fastBM(tree,nsim=2)