Friday, October 23, 2015

Node, edge, & tip labels for plotted cophylo objects

A recent blog comment asked if there was some way to 'maintain support values' (i.e., include node, edge, or tip labels) for a plotted "cophylo" object produced by phytools.

This is not quite as straightforward as it might seem because, unlike for a single plotted tree, there are two sets of node, two sets of tips, two sets of edges etc.

My solution was as follows:

First, I have the internally used tree plotting function, phylogram, create an environmental variable "last_plot.phylo", as do other tree plotting methods such as plot.phylo and plotSimmap.

Next, I had the S3 plot method for objects of class "cophylo" pull each of these two "last_plot.phylo" objects from their environment (one for each tree plotted) and create a new environmental variable, "last_plot.cophylo", in the same environment, consisting of a simple list of these two objects.

Finally, I created a series of wrappers around nodelabels, tiplabels, and edgelabels, which take as input which tree is to have labels added to it, and then sets the corresponding element of "last_plot.cophylo" to "last_plot.phylo" so that nodelabels etc. can be called internally.

Here's a demo using the following two trees which have matching tree labels. (Remember, we can use non-matching labels too if we are prepared to supply an association matrix.)

library(phytools)
packageVersion("phytools")
## [1] '0.5.4'
par(mfrow=c(1,2))
plotTree(t1)
plotTree(t2)

plot of chunk unnamed-chunk-1

Now, let's create an object of class "cophylo" and then try each of our new labeling methods:

obj<-cophylo(t1,t2)
## Rotating nodes to optimize matching...
## Done.
plot(obj)
nodelabels.cophylo() ## left side
nodelabels.cophylo(which="right")

plot of chunk unnamed-chunk-2

We can also change the values represented, of course. For instance:

library(stringr)
n.left<-vector()
for(i in 1:obj$trees[[1]]$Nnode){
    ii<-getDescendants(obj$trees[[1]],i+Ntip(obj$trees[[1]]))
    n.left[i]<-str_to_lower(paste(sort(obj$trees[[1]]$tip.label[ii[ii<=Ntip(obj$trees[[1]])]]),
        collapse=""))
}
n.right<-vector()
for(i in 1:obj$trees[[2]]$Nnode){
    ii<-getDescendants(obj$trees[[2]],i+Ntip(obj$trees[[2]]))
    n.right[i]<-str_to_lower(paste(sort(obj$trees[[2]]$tip.label[ii[ii<=Ntip(obj$trees[[2]])]]),
        collapse=""))
}
plot(obj)
nodelabels.cophylo(n.left,cex=0.9,bg="white",srt=-90)
nodelabels.cophylo(n.right,which="right",cex=0.9,bg="white",srt=-90)

plot of chunk unnamed-chunk-3

Functions to label edges & tips also work for this object class, for instance:

plot(obj)
## show rounded edge lengths:
edgelabels.cophylo(round(obj$trees[[1]]$edge.length,2),cex=0.7,
    frame="none",adj=c(0.5,1.2))
edgelabels.cophylo(round(obj$trees[[2]]$edge.length,2),cex=0.7,
    frame="none",adj=c(0.5,1.2),which="right")

plot of chunk unnamed-chunk-4

plot(obj,ftype="off")
## add tip labels
tiplabels.cophylo(pch=21,frame="none",bg="grey",cex=1.5)
tiplabels.cophylo(pch=21,frame="none",bg="grey",which="right",cex=1.5)

plot of chunk unnamed-chunk-5

We can use any of the options available from nodelabels and its sister functions. So, for instance, if we want to overlay reconstructed ancestral states for a discrete character, we can do so easily:

## tip data
x
##   C   E   H   L   J   B   A   D   K   G   F   I 
## "b" "a" "b" "b" "b" "b" "b" "a" "a" "a" "b" "a"
y
##   H   C   A   L   J   B   F   E   G   D   I   K 
## "a" "a" "b" "b" "a" "b" "a" "a" "a" "b" "b" "b"
fit.x<-ace(x,obj$trees[[1]],type="discrete")
fit.y<-ace(y,obj$trees[[2]],type="discrete")
## plot pies
plot(obj)
nodelabels.cophylo(pie=fit.x$lik.anc,piecol=c("grey","white"),cex=0.7)
nodelabels.cophylo(pie=fit.y$lik.anc,piecol=c("grey","white"),cex=0.7,
    which="right")
tiplabels.cophylo(pie=to.matrix(x[obj$trees[[1]]$tip.label],c("a","b")),
    piecol=c("grey","white"),cex=0.4)
tiplabels.cophylo(pie=to.matrix(y[obj$trees[[2]]$tip.label],c("a","b")),
    piecol=c("grey","white"),cex=0.4,which="right")

plot of chunk unnamed-chunk-6

Something very important to keep in mind is that the order of the tips, nodes, & edges in the object of class "cophylo" may be different than it was in the input trees, particulary if we set rotate=TRUE when we made the object.

That's it!

2 comments:

  1. This comment has been removed by the author.

    ReplyDelete
  2. I needed to use cophyl to do an awesome plot one nuclear data and one chloroplast data 135 tips. I started with RAXML bootstrap trees (bipartition files).
    #read in trees
    read.tree("nucTreeraxml")->nucTree
    read.tree("cTreeraxml")->cTree
    #I want the tips to line up and depth to be to be good enough to show nodes separately
    chronopl(nucTree, lambda=0.1)->nucTree
    chronopl(cTree, lambda=0.1)->cTree
    #plot cophylo with supports
    obj<-cophylo(cTree,nucTree) #default optimization worked best for me
    plot(obj, fsize=0.3)
    nodelabels.cophylo(obj$trees[[1]]$node.label, frame="none",cex=0.3)
    nodelabels.cophylo(obj$trees[[2]]$node.label, frame="none",cex=0.3, which="right")
    #this was pretty cool but I wanted to SET UP FOR COLOR coded dots at nodes as numbers are tedious to read
    obj$trees[[1]]$node.label<-as.numeric(obj$trees[[1]]$node.label) #oddly support labels were stored as characters so change to numeric
    obj$trees[[2]]$node.label<-as.numeric(obj$trees[[2]]$node.label)
    #set up empty object to store color names for nodes
    p1<-character(length(obj$trees[[1]]$node.label))
    p2<-character(length(obj$trees[[2]]$node.label))
    #set up ranges that correspond to node support levels
    p1[obj$trees[[1]]$node.label >= 95] <- "red"
    p1[obj$trees[[1]]$node.label < 95 & obj$trees[[1]]$node.label >= 75] <- "orange"
    p1[obj$trees[[1]]$node.label < 75] <- "grey"
    p1[is.na(obj$trees[[1]]$node.label)] <- "white"
    p2[obj$trees[[2]]$node.label >= 95] <- "red"
    p2[obj$trees[[2]]$node.label < 95 & obj$trees[[2]]$node.label >= 75] <- "orange"
    p2[obj$trees[[2]]$node.label < 75] <- "grey"
    p2[is.na(obj$trees[[2]]$node.label)] <- "white" #root node did not have a support value
    #check it worked
    obj$trees[[1]]$node.label
    p1
    obj$trees[[2]]$node.label
    p2
    #write to single page pdf
    pdf(file="support_cophylo_CP_NUC.pdf", 8.5, 11)
    plot(obj, fsize=0.3)
    nodelabels.cophylo(pch=21,bg=p1,cex=0.7)
    nodelabels.cophylo(pch=21,bg=p2,cex=0.7, which="right")
    dev.off()

    ReplyDelete