Thursday, October 18, 2018

Plotting a "simmap" or "contMap" object in unrooted style

I haven't been posting to actively to this blog lately I have been busy with teaching and with my 'new' package learnPopGen as well as with the various shiny web interfaces that I have built for a number of the functions of the package.

Nonetheless, this morning for some reason I was thinking about the problem of plotting unrooted trees (which is actually much more complicated that plotting rooted trees in their various configurations, as discussed in Chapter 34 of Felsenstein 2004), and I suddenly realized that to plot trees with a mapped distrete character (objects of class "simmap") or trees with a mapped continuous character (objects of class "contMap") I could quite easily take advantage of ape::plot.phylo and simply project my tree on top of the coordinates of the unrooted tree plotted by plot.phylo in which this problem has already been solved.

This is what I mean.

library(phytools)
set.seed(777)
N<-16
tree<-rtree(n=N,tip.label=LETTERS[1:N])
Q<-matrix(c(-1,1,1,-1),2,2,dimnames=list(letters[1:2],letters[1:2]))
x<-sim.Mk(tree,Q)
x
## L J I H K B E M C N G D F A O P 
## a b b a a a a b b a b b b b b b 
## Levels: a b
map.tree<-make.simmap(tree,x)
## make.simmap is sampling character histories conditioned on the transition matrix
## 
## Q =
##           a         b
## a -1.038569  1.038569
## b  1.038569 -1.038569
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##   a   b 
## 0.5 0.5
## Done.
## here's where it gets interesting:
plot.phylo(map.tree,type="unrooted",lab4ut="axial",no.margin=TRUE,
    label.offset=0.05)
obj<-get("last_plot.phylo",envir=.PlotPhyloEnv)
tip.pos<-matrix(c(obj$xx[1:Ntip(map.tree)],obj$yy[1:Ntip(map.tree)]),
    Ntip(map.tree),2,dimnames=list(map.tree$tip.label,NULL))
node.pos<-matrix(c(obj$xx[1:map.tree$Nnode+Ntip(map.tree)],
    obj$yy[1:map.tree$Nnode+Ntip(map.tree)]),map.tree$Nnode,2,
    dimnames=list(1:map.tree$Nnode+Ntip(map.tree),NULL))
phylomorphospace(map.tree,tip.pos,A=node.pos,lwd=3,
    colors=setNames(c("blue","red"),sort(unique(getStates(map.tree)))),
    add=TRUE,ftype="off",node.size=c(0,1),node.by.map=TRUE)

plot of chunk unnamed-chunk-1

We can do something similar with a "contMap" object as follows:

set.seed(777)
y<-fastBM(tree)
cmap<-contMap(tree,y,plot=FALSE,res=1000)
par(lend=2)
plot.phylo(cmap$tree,type="unrooted",lab4ut="axial",no.margin=TRUE,
    label.offset=0.05,edge.width=7)
obj<-get("last_plot.phylo",envir=.PlotPhyloEnv)
tip.pos<-matrix(c(obj$xx[1:Ntip(cmap$tree)],obj$yy[1:Ntip(cmap$tree)]),
    Ntip(cmap$tree),2,dimnames=list(cmap$tree$tip.label,NULL))
node.pos<-matrix(c(obj$xx[1:cmap$tree$Nnode+Ntip(cmap$tree)],
    obj$yy[1:cmap$tree$Nnode+Ntip(cmap$tree)]),cmap$tree$Nnode,2,
    dimnames=list(1:cmap$tree$Nnode+Ntip(cmap$tree),NULL))
phylomorphospace(cmap$tree,tip.pos,A=node.pos,lwd=5,
    colors=cmap$cols,add=TRUE,ftype="off",
    node.size=c(0,0))
add.color.bar(h<-2,cmap$cols,title="trait value",
    lims=NULL,digits=3,direction="upwards",
    subtitle="",lwd=15,x=0.1,y=1.75,prompt=FALSE)
LWD<-diff(par()$usr[1:2])/dev.size("px")[1]
nticks<-6
Y<-cbind(seq(1.75,1.75+h,length.out=nticks),
    seq(1.75,1.75+h,length.out=nticks))
X<-cbind(rep(0.1+LWD*17/2,nticks),
    rep(0.1+LWD*17/2+0.02*h,nticks))
for(i in 1:nrow(Y)) lines(X[i,],Y[i,])
ticks<-seq(cmap$lims[1],cmap$lims[2],length.out=nticks)
text(x=X[,2],y=Y[,2],round(ticks,2),pos=4,cex=0.8)

plot of chunk unnamed-chunk-2

That's pretty cool (IMO).

Friday, September 28, 2018

New phytools version (0.6-60) on CRAN

In advance of a workshop that I will be teaching next week in Spain with Luke Harmon, I just submitted a new phytools version (0.6-60) to CRAN and overnight it was posted. (As of this morning, binaries have yet to be built but it can already be installed from source.)

It's been a while since my last CRAN update of phytools, so there is are a number of new features & functions in this version. As usual, all changes to phytools can be tracked via its GitHub page, but the following are nonetheless some highlights:

  1. A small but important fix to fitMk and make.simmap that affects some internals running matrix exponentiation which is otherwise sometimes very slow.

  2. A significant update to the function phyl.pairedttest for a phylogenetic paired t-test.

  3. Significant increases in user control of the S3 plotting method for "cophylo" objects (1, 2).

  4. A new object class and S3 methods for the multiple (partial) Mantel test function, multi.mantel.

  5. A new function, ctt, for visualizing the number of reconstructed changes on stochastic map trees through time (1, 2).

  6. A totally new method, implemented in fitmultiMk, for fitting a multi-rate Mk model in which the rate (or process) of character transition between states changes between different edges or clades in the tree (1, 2).

  7. A function gtt, along with print and plot methods, for plotting Pybus & Harvey's γ-statistic through time. (I actually forgot to put this in the namespace, but it can still be called using phytools:::gtt or by updating with the fix I just pushed to GitHub.)

  8. A totally new function, sim.multiMk, to simulate evolution by a variable-rate Mk process on the tree (1, 2).

  9. A small fix to the "static" method of phylomorphospace3d.

  10. A totally new plotting method (plotTree.datamatrix) to visualize a discrete character data matrix next to a tree (1, 2, 3, 4).

  11. Another new plotting function, but to animate the projection of a tree into morphospace called project.phylomorphospace (1, 2).

  12. Several different updates to the phytools function plotTree.wBars (1, 2).

  13. Some small updates to plotTree.barplot.

  14. And numerous other smaller things that are all catalogued on my GitHub.

Please check out the new version by updating from CRAN or GitHub & let me know if you run into any problems!

Tuesday, June 5, 2018

Preserving node & edge labels after optimizing node rotation in cophylo

A phytools user contacted me the other day about an issue in which her bootstrap values, stored as node labels in a Newick string, did not align with the correct edges when plotted using edgelabels.cophylo on a visualized "cophylo" object.

As it turns out, this had nothing to do with cophylo & everything to do with the fact that she wanted to plot node labels using edgelabels.

First, here is a demo showing that cophylo preserves the correct order of node labels.

A simulated tree with node labels:

library(phytools)
tree<-pbtree(n=26,tip.label=LETTERS)
tree$node.label<-letters[1:25]

Now I'll plot that tree labeling the nodes with their labels & node indices:

plotTree(tree)
nodelabels(tree$node.label,1:tree$Nnode+Ntip(tree))
nodelabels(adj=c(1.8,-1.8),frame="none",cex=0.5)

plot of chunk unnamed-chunk-2

Next, I'm going to randomly rotate nodes on this tree to produce two new trees as follows:

t1<-tree
for(i in 1:100) t1<-untangle(rotate(t1,sample(1:t1$Nnode+Ntip(t1),1)),
    'read.tree')
t2<-tree
for(i in 1:100) t2<-untangle(rotate(t2,sample(1:t2$Nnode+Ntip(t2),1)),
    'read.tree')

Now I'm ready to make my "cophylo" object. First, without node rotation to maximize matching:

obj<-cophylo(t1,t2,rotate=FALSE)
plot(obj)
nodelabels.cophylo(obj$trees[[1]]$node.label,1:obj$trees[[1]]$Nnode+
    Ntip(obj$trees[[1]]))
nodelabels.cophylo(adj=c(1.8,-1.8),frame="none",cex=0.5)
nodelabels.cophylo(obj$trees[[2]]$node.label,1:obj$trees[[2]]$Nnode+
    Ntip(obj$trees[[2]]),which="right")
nodelabels.cophylo(adj=c(-0.8,-1.8),frame="none",cex=0.5,which="right")

plot of chunk unnamed-chunk-4

We can see that the node indices (the little numbers that come from tree$edge have changed, but the node labels are still all associated with the right clades.

Now we can do the same thing, but in which we permit node rotation to be optimized as follows:

obj<-cophylo(t1,t2,rotate=TRUE)
## Rotating nodes to optimize matching...
## Done.
plot(obj)
nodelabels.cophylo(obj$trees[[1]]$node.label,1:obj$trees[[1]]$Nnode+
    Ntip(obj$trees[[1]]))
nodelabels.cophylo(adj=c(1.8,-1.8),frame="none",cex=0.5)
nodelabels.cophylo(obj$trees[[2]]$node.label,1:obj$trees[[2]]$Nnode+
    Ntip(obj$trees[[2]]),which="right")
nodelabels.cophylo(adj=c(-0.8,-1.8),frame="none",cex=0.5,which="right")

plot of chunk unnamed-chunk-5

So what was the problem with the edge labels in the user's inquiry? Well, basically, edgelabels (and edgelabels.cophylo, which just uses edgelabels internally) takes the labels in the order of the rows of tree$edge in the plotted phylogeny

  • not in the number order of the node indices, as in nodelabels and tree$node.label. Thus, we have to match the two in order to visualize our node labels along each preceding edge. This might be done as follows:
obj<-cophylo(t1,t2,rotate=TRUE)
## Rotating nodes to optimize matching...
## Done.
plot(obj)
edgelabels.cophylo(obj$trees[[1]]$node.label[2:obj$trees[[1]]$Nnode],
    edge=sapply(2:obj$trees[[1]]$Nnode+Ntip(obj$trees[[1]]),
    function(n,e) which(e==n),e=obj$trees[[1]]$edge[,2]),
    frame="none",adj=c(0.5,1))
edgelabels.cophylo(obj$trees[[2]]$node.label[2:obj$trees[[2]]$Nnode],
    edge=sapply(2:obj$trees[[2]]$Nnode+Ntip(obj$trees[[2]]),
    function(n,e) which(e==n),e=obj$trees[[2]]$edge[,2]),
    frame="none",adj=c(0.5,1),which="right")

plot of chunk unnamed-chunk-6

Neat.