New version of plotTree.wBars that permits positive & negative phenotypic trait values

A couple of weeks ago I posted a new function (1, 2) to plot bars at the tips of a circular or square phylogram. One limitation of this function is that because the bars are plotted 'growing' out of each leaf of the tree, the values of the phenotypic trait data underlying the bars cannot be negative. Negative values would result in bars growing (in the case of a fan tree) towards the root of the tree - which, of course, does not look right at all.

Here's what I mean:

> tree<-pbtree(n=100)
> x<-fastBM(tree)
> plotTree.wBars(tree,x,scale=0.5)

> ## or setting type="fan"
> plotTree.wBars(tree,x,scale=0.5,type="fan")

Now let's try the new version:

> require(phytools)
Loading required package: phytools
> packageVersion("phytools")
[1] ‘0.4.11’
> plotTree.wBars(tree,x,scale=0.5)

> plotTree.wBars(tree,x,scale=0.5,type="fan")

I'm not sure it makes a great visual in this case - but, nonetheless, it's better.

When some values of x are negative and the tree is non-ultrametric, then the bars are also centered a constant distance from the maximum tip height in the tree. For instance:

> tree<-rtree(n=50)
> x<-fastBM(tree)
> plotTree.wBars(tree,x,scale=0.3)

The code for this new function version is here; but you can also install a new build of phytools from source with this update.

Some useful updates to contMap & densityMap

The most common error reported by users of contMap & densityMap is the following:

Error in while (x > trans[i]) { : missing value where TRUE/FALSE needed

This is usually (but not always) due to the tree having internal edges of zero length, and is just the result of how I set up the internal machinery of these functions. So, for instance:

> tree<-rtree(10)
> tree$edge.length[which(tree$edge[,2]==14)]<-0
> plotTree(tree)

> x<-fastBM(tree)
> obj<-contMap(tree,x)
Error in while (x > trans[i]) { : missing value where TRUE/FALSE needed

Obviously, this could be circumvented by first collapsing zero-length branches into polytomies, and this is typically what I have advised users to do - but there is no really good reason why contMap shouldn't be able to work with the zero-length branches still intact. I have just posted a new version with the bug fixed, and it is in a new release of phytools, which can be downloaded & installed from source here.

Here's a demo using the same tree as before:

> packageVersion("phytools")
[1] ‘0.4.10’
> obj<-contMap(tree,x)

I have also added a new function to automate the process of changing the color ramp, described here. This function is called setMap. Here's a demo:

> obj<-setMap(obj,colors=c("blue","green","yellow"), space="Lab")
> plot(obj)

Obviously, there is no way that R can make you use a reasonable color scheme here!

In addition to fixing the same branches of zero-length problem in densityMap, it also now can plot the posterior density of any two-state character, not just a binary character with states of "0" & "1". This required a few small modifications, including accommodating the possibility that someone might use one & only one of states "0" & "1" (and then some other state not "0" or "1"). Here's a demo:

> tree<-pbtree(n=26,tip.label=letters[26:1],scale=1)
> Q<-matrix(c(-1,1,1,-1),2,2)
> rownames(Q)<-colnames(Q)<-LETTERS[1:2]
> x<-sim.history(tree,Q)$states
> trees<-make.simmap(tree,x,nsim=100)
make.simmap is sampling character histories conditioned on the transition matrix
Q =
           A          B
A -0.5933734  0.5933734
B  0.5933734 -0.5933734
(estimated using likelihood);
and (mean) root node prior probabilities
pi =
  A   B
0.5 0.5
Done.
> obj<-densityMap(trees,plot=FALSE,res=500)
sorry - this might take a while; please be patient
> plot(obj,outline=TRUE,lwd=6)

That's it for now.

Changing the color ramp in contMap or densityMap

I recently received the following question:

"Would it be possible in contMap to specify our own color ramp instead of stuck with the default red-to-blue?"

This cannot yet be done automatically, but it is not too hard. Here are a couple of examples. Note that they also apply equally to objects created using densityMap.

> ## first create simulated tree & data
> tree<-pbtree(n=26,tip.label=LETTERS[26:1])
> x<-fastBM(tree)
> ## build object of class "contMap"
> 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 (-4.458, 2.5).

> ## default colors
> plot(obj)

> ## what is the length of the current color ramp?
> n<-length(obj$cols)
> ## change to grey scale
> obj$cols[1:n]<-grey(0:(n-1)/(n-1))
> plot(obj)

> ## change to blue -> red
> obj$cols[1:n]<-colorRampPalette(c("blue","red"), space="Lab")(n)
> plot(obj)

That's it!

Plotting bars at the tips of a tree, part II

Earlier today, I responded to an R-sig-phylo request to be able to plot bars showing phenotypic trait values for species at the times of a circular or 'fan' tree. I have now added this function (plotTree.wBars) to the phytools package. It can be downloaded & installed from source here.

In addition to plotting a circular tree, this function version also:

(1) Plots a square phylogram in "rightward" or "leftward" orientation.

(2) Plots a stochastic character mapped tree (using plotSimmap instead of plotTree internally).

(3) Can accept most of the arguments of plotTree & plotSimmap (see the documentation page for plotSimmap for more information.)

Here's a demo:

> packageVersion("phytools")
[1] ‘0.4.8’
> tree<-pbtree(n=50)
> x<-fastBM(tree,bounds=c(0,Inf))
> plotTree.wBars(tree,x)

> Q<-matrix(c(-1,1,1,-1),2,2)
> tree<-sim.history(tree,Q)
> plotTree.wBars(tree,x,method="plotSimmap", colors=setNames(c("blue","red"),1:2),lwd=3)

> plotTree.wBars(tree,x,type="fan",method="plotSimmap", colors=setNames(c("blue","red"),1:2),lwd=3,scale=0.5, width=0.2)

That's it.

Plotting bars at the tips of a circular tree

Today Matt Helmus posted the following request to the R-sig-phylo mailing list:

"Does anyone know of R code (or perhaps another program) to plot bars across the tips of a radial/fan phylogeny? Specifically, I have a large phylogeny and a corresponding vector of continuous trait values for the tips, and while I could use those values to vary the color and size of the tip labels, or also to plot points of varying color or size at those tips, I think a better depiction of the data would be to plot bars that vary in height."

Naturally, this appealed to me. The biggest challenge is resuscitating high-school (or perhaps elementary school) trig (think SOH-CAH-TOA) to get the position of the vertices for each tip bar/rectangle correct. Here is my attempt. Note that it doesn't contain much in the ways of bells-and-whistles, but these could easily be added:

plotFan.wBars<-function(tree,x,scale=1,width=NULL){
  lims=c(-max(nodeHeights(tree))-scale*max(x),
    max(nodeHeights(tree))+scale*max(x))
  obj<-plotTree(tree,type="fan",ftype="off",
    xlim=lims,ylim=lims)
  if(is.null(width)) width<-(par()$usr[2]-par()$usr[1])/100
  w<-width
  x<-x[tree$tip.label]*scale
  for(i in 1:length(x)){
    dx<-obj$xx[i]
    dy<-obj$yy[i]
    theta<-atan(dy/dx)
    x1<-dx+(w/2)*cos(pi/2-theta)
    y1<-dy-(w/2)*sin(pi/2-theta)
    x2<-dx-(w/2)*cos(pi/2-theta)
    y2<-dy+(w/2)*sin(pi/2-theta)
    s<-if(dx>0) 1 else -1
    x3<-s*x[i]*cos(theta)+x2
    y3<-s*x[i]*sin(theta)+y2
    x4<-s*x[i]*cos(theta)+x1
    y4<-s*x[i]*sin(theta)+y1
    polygon(c(x1,x2,x3,x4),c(y1,y2,y3,y4),col="grey")
  }
  invisible(obj)
}

The two inobvious input arguments, scale & width, control the scale of the tip bars relative to the total height of the tree & the width of bars, respectively.

Let's try it:

> tree<-pbtree(n=50,scale=1)
> x<-fastBM(tree)
> x<-abs(x)
> plotFan.wBars(tree,x,scale=0.2)

That's all.

contMap with missing data

Today a phytools user emailed with the following question:

"I am still having difficulty running an analysis on the evolution of a continuous character with a partial data set (I have character states for about 50% of the taxa). I would like to display the results using contMap, but since fastAnc is built into the code, it won't map the character because the # of taxa in the data set doesn't match the # of taxa in the tree. I attempted to replace fastAnc with the modified version of anc.ML from one of your blog posts (link), but I couldn't get that to work either. I received the same error message as with fastAnc:
Error in ace(x, a, method = "pic") : length of phenotypic and of phylogenetic data do not match
Is there any way to use contMap with missing data? Any help would be greatly appreciated. Thank you for your time."

They are totally on the right track, but things get a little trickier than they might have anticipated because whereas fastAnc (when used in it's default mode) just returns a vector of estimated character values at nodes, anc.ML (when some taxa are present in the tree but missing from the input data) returns a list with the vector missing.x containing the MLE phenotypic trait values of taxa with missing data. It also happens to be the case that I never got around to adding this updated version of anc.ML to the phytools package.

Well, I have now made both updates so it is possible to (1) specify an inference method (method="fastAnc" or method="anc.ML"); and (2) phytools function anc.ML now permits missing data, in which case the states for the tip taxa in the tree missing from the data vector will be estimated using ML. The code for these updates are here: anc.ML, contMap, but since both use internal functions not in the phytools namespace it is probably wise to download the latest phytools build and install from source.

Here's a quick demo:

> packageVersion("phytools")
[1] ‘0.4.7’
> ## simulate tree & data
> tree<-pbtree(n=26,tip.label=LETTERS)
> x<-fastBM(tree)
> ## compute limits on x to use in both plots
> lims<-c(floor(min(x)*10)/10,ceiling(max(x)*10)/10)
> ## simulate missing tip values
> y<-sample(x,20)
> ## which taxa are missing
> setdiff(tree$tip.label,names(y))
[1] "H" "I" "O" "T" "W" "Z"
> ## contMap from full dataset
> contMap(tree,x,lims=lims,legend=1.5)

> ## contMap from reduced dataset
> contMap(tree,y,lims=lims,method="anc.ML",legend=1.5)

OK - that's it.

New version of phytools submitted to CRAN

I just submitted a new version of phytools (phytools 0.4-05) to CRAN. Some updates in the current version relative to the previous CRAN version include the following:

1. A new function, cladelabels, to label membership in a clade.

2. A new function, fastHeight, to get the height above the root of the common ancestor of a pair of species.

3. A new function, roundPhylogram (1, 2), to plot a round phylogram.

4. A new function, paintBranches, to paint an edge or set of edges on the tree.

5. Some new options and S3 methods for the function ltt95 (described here).

6. A simple new function, nodeheight, to quickly compute the height above the root of a single node. (For large trees this method will be much quicker than nodeHeights if the height of only one or a small number of nodes are of interest).

7. New plot & print methods for describe.simmap (described here).

8. Some new plotting options for contMap and densityMap (described here).

9. Finally, last & possibly least, a new function drop.tip.singleton (1, 2) to drop one or multiple tips from the tree while retaining all ancestral nodes of remaining tips from the original tree as singleton nodes.

Hopefully, this new version will be accepted & posted to CRAN; but in the meantime it can be downloaded & installed from source from the phytools page. Check it out.

Fixed drop.tip.singleton

Here is a fixed version of the function I just posted to drop leaves while retaining all ancestral nodes of remaining extant taxa as singletons:

drop.tip.singleton<-function(tree,tip){
  N<-length(tree$tip.label)
  m<-length(tip)
  ii<-sapply(tip,function(x,y) which(y==x),y=tree$tip.label)
  tree$tip.label<-tree$tip.label[-ii]
  ii<-sapply(ii,function(x,y) which(y==x),y=tree$edge[,2])
  tree$edge<-tree$edge[-ii,]
  tree$edge.length<-tree$edge.length[-ii]
  tree$edge[tree$edge<=N]<-
    as.integer(rank(tree$edge[tree$edge<=N]))
  tree$edge[tree$edge>N]<-tree$edge[tree$edge>N]-m
  N<-N-m
  if(any(sapply(tree$edge[tree$edge[,2]>N,2],"%in%",
    tree$edge[,1])==FALSE)) internal<-TRUE
  else internal<-FALSE
  while(internal){
    ii<-which(sapply(tree$edge[,2],"%in%",
      c(1:N,tree$edge[,1]))==FALSE)[1]
    nn<-tree$edge[ii,2]
    tree$edge<-tree$edge[-ii,]
    tree$edge.length<-tree$edge.length[-ii]
    tree$edge[tree$edge>nn]<-tree$edge[tree$edge>nn]-1
    tree$Nnode<-tree$Nnode-length(ii)
    if(any(sapply(tree$edge[tree$edge[,2]>N,2],
      "%in%",tree$edge[,1])==FALSE)) internal<-TRUE
    else internal<-FALSE
  }
  tree
}

The previous version ran into trouble when the number of remaining internal edges that the function tried to remove in a step is >1. This version is less ambitious & seems to work fine.

Dropping tips while retaining the ancestors of remaining extant tips as singleton nodes

Luke Mahler asked the following:

"Do you know of a way to drop a terminal branch from a phylogeny, yet preserve the node it came from as a singleton node? I initially thought drop.tip(trim.internal=F) would do this, but it does something a little different, apparently (it preserves internal branches that become tips by pruning, but not nodes that would become singleton nodes)."

In the simple case in which we just want to drop one tip, this is relatively straightforward. We just have to drop the corresponding row & element from tree$edge, tree$edge.length, and tree$tip.label, and then update our node & tip numbers in tree$edge to follow the "phylo" object convention. However, generalizing to drop an arbitrary number of tips (while retaining all ancestral nodes to extant tips, regardless of whether they now have one or multiple descendants) now becomes a little bit trickier. Here is my function for this:

drop.tip.singleton<-function(tree,tip){
  N<-length(tree$tip.label)
  m<-length(tip)
  ii<-sapply(tip,function(x,y) which(y==x),y=tree$tip.label)
  tree$tip.label<-tree$tip.label[-ii]
  ii<-sapply(ii,function(x,y) which(y==x),y=tree$edge[,2])
  tree$edge<-tree$edge[-ii,]
  tree$edge.length<-tree$edge.length[-ii]
  tree$edge[tree$edge<=N]<-
    as.integer(rank(tree$edge[tree$edge<=N]))
  tree$edge[tree$edge>N]<-tree$edge[tree$edge>N]-m
  N<-N-m
  if(any(sapply(tree$edge[tree$edge[,2]>N,2],"%in%",
    tree$edge[,1])==FALSE)) internal<-TRUE
  while(internal){
    ii<-which(sapply(tree$edge[,2],"%in%",c(1:N,
      tree$edge[,1]))==FALSE)
    nn<-sort(tree$edge[ii,2])
    tree$edge<-tree$edge[-ii,]
    tree$edge.length<-tree$edge.length[-ii]
    for(i in 1:length(nn)) tree$edge[tree$edge>nn[i]]<-
      tree$edge[tree$edge>nn[i]]-1
    tree$Nnode<-tree$Nnode-length(ii)
    if(any(sapply(tree$edge[tree$edge[,2]>N,2],
      "%in%",tree$edge[,1])==FALSE)) internal<-TRUE
    else internal<-FALSE
  }
  tree
}

Now try it:

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

> tip<-sample(LETTERS,10)
> tip
[1] "N" "M" "F" "I" "Z" "R" "P" "S" "G" "W"
> tt<-drop.tip.singleton(tree,tip)
> plotTree.singletons(tt)

This seems to be the correct result.