Update to demarcating changes on a plotted "simmap" style tree using markChanges

I posted recently about a new phytools function to mark reconstructed changes on the tree from stochastic mapping. I will be adding this function to phytools; however I have made a couple of changes. Foremost among these is that the function now returns (invisibly) a matrix containing the x & y coordinates of the plotted changes on the tree, as well as the type of each change (e.g., a->b, b->c, etc. We can also turn off plotting so that only this matrix is computed & returned. This theoretically enables us to overlay any symbol that we want on the tree to mark changes.

Here's a quick demo of the function as presently written:

library(phytools)
data(anoletree)
## function to mark changes
markChanges<-function(tree,colors=NULL,cex=1,lwd=2,plot=TRUE){
    states<-sort(unique(getStates(tree)))
    if(is.null(colors)) colors<-setNames(palette()[1:length(states)],
        states)
    obj<-get("last_plot.phylo",envir=.PlotPhyloEnv)
    nc<-sapply(tree$maps,length)-1
    ii<-which(nc>0)
    nc<-nc[ii]
    xx<-yy<-vector()
    for(i in 1:length(ii)){
        for(j in 1:nc[i]){
            ss<-names(tree$maps[[ii[i]]])[j+1]
            mm<-tree$edge[ii[i],1]
            dd<-tree$edge[ii[i],2]
            x<-rep(obj$xx[mm]+cumsum(tree$maps[[ii[i]]])[j],2)
            y<-c(obj$yy[dd]-0.5*mean(strheight(LETTERS)*cex),
                obj$yy[dd]+0.5*mean(strheight(LETTERS)*cex))
            if(plot) lines(x,y,lwd=lwd,col=colors[ss],lend=2)
            xx<-c(xx,setNames(x[1],
                paste(names(tree$maps[[ii[i]]])[j:(j+1)],
                collapse="->")))
            yy<-c(yy,mean(y))
        }
    }
    XY<-cbind(xx,yy)
    colnames(XY)<-c("x","y")
    invisible(XY)
}
plotSimmap(anoletree,ftype="i",fsize=0.6,ylim=c(-1,Ntip(anoletree)))

## no colors provided. using the following legend:
##        CG        GB        TC        TG        Tr        Tw 
##   "black"     "red"  "green3"    "blue"    "cyan" "magenta"

add.simmap.legend(x=0,y=-1,colors=
    setNames(palette()[1:length(unique(getStates(anoletree,"tips")))],
    sort(unique(getStates(anoletree,"tips")))),prompt=FALSE,vertical=FALSE)
obj<-markChanges(anoletree,lwd=2,cex=0.6)

Now obj is a matrix containing the x & y coordinates of each change as recorded in this particular stochastic map of (here) ecomorph state on the tree. So, in this case:

obj

##                x        y
## TG->GB 4.8320281  8.00000
## TG->TC 3.9596523 15.25000
## TC->CG 4.9301502 14.00000
## TG->Tw 3.5809367 17.00000
## TG->TC 3.8002031 19.50000
## TG->GB 3.8962124 21.50000
## TG->GB 4.2372427 24.00000
## TG->Tr 2.0060270 29.93750
## TG->GB 2.9155935 33.25000
## TG->Tw 3.9550415 35.00000
## TG->CG 1.6044174 47.12500
## TG->Tw 0.4992279 59.88281
## Tw->TC 2.3342945 52.68750
## TC->Tw 3.1128996 54.00000
## Tw->Tr 5.0515733 54.00000
## Tw->TG 2.2167389 64.03125
## TG->GB 2.8428529 64.03125
## TG->CG 4.2711895 71.93750
## TG->GB 1.9843757 75.87500
## GB->CG 3.6874052 77.00000
## CG->Tw 5.6935007 77.00000
## TG->TC 2.3486612 80.12500
## TG->TC 1.5791488 82.00000
## TC->Tw 4.2869186 82.00000

We can use different symbols to plot these changes if we want to, for instance, simply:

plotSimmap(anoletree,ftype="i",fsize=0.6,ylim=c(-1,Ntip(anoletree)))

## no colors provided. using the following legend:
##        CG        GB        TC        TG        Tr        Tw 
##   "black"     "red"  "green3"    "blue"    "cyan" "magenta"

obj<-markChanges(anoletree,plot=FALSE)
add.simmap.legend(x=0,y=-1,colors=
    setNames(palette()[1:length(unique(getStates(anoletree,"tips")))],
    sort(unique(getStates(anoletree,"tips")))),prompt=FALSE,vertical=FALSE)
points(obj,pch=19)

## or
plotSimmap(anoletree,ftype="i",fsize=0.6,ylim=c(-1,Ntip(anoletree)))

## no colors provided. using the following legend:
##        CG        GB        TC        TG        Tr        Tw 
##   "black"     "red"  "green3"    "blue"    "cyan" "magenta"

obj<-markChanges(anoletree,plot=FALSE)
add.simmap.legend(x=0,y=-1,colors=
    setNames(palette()[1:length(unique(getStates(anoletree,"tips")))],
    sort(unique(getStates(anoletree,"tips")))),prompt=FALSE,vertical=FALSE)
points(obj,pch=8)

etc.

That's it for now.