New co-phylogenetic method for phytools

For a day or two, I've been playing around with a new co-phylogenetic plotting method.

Here's what it looks like right now:

cotangleplot<-function(tr1,tr2,type=c("cladogram","phylogram"),
    use.edge.length=TRUE,tangle=c("both","tree1","tree2"),...){
    tr1<-untangle(tr1,"read.tree")
    tr2<-untangle(tr2,"read.tree")
    type<-type[1]
    tangle<-tangle[1]
    if(!use.edge.length){
        tr1<-compute.brlen(tr1)
        tr2<-compute.brlen(tr2)
    }
    if(hasArg(layout)) layout<-list(...)$layout
    else layout<-c(0.45,0.1,0.45)
    if(hasArg(nodes)) nodes<-list(...)$nodes
    else nodes<-"centered"
    if(hasArg(lty)) lty<-list(...)$lty
    else lty<-if(type=="phylogram") c("dotted","solid") else
        if(type=="cladogram") "solid"
    if(type=="phylogram") if(length(lty)==1) lty<-rep(lty,2)
    if(hasArg(lwd)) lwd<-list(...)$lwd
    else lwd<-2
    if(hasArg(cex)) cex<-list(...)$cex
    else cex<-1
    if(hasArg(color)) color<-list(...)$color
    else color<-palette()[4]
    if(tangle=="both"){
        capture.output(tmp<-cophylo(tr1,tr2))
        tips.tr1<-setNames(1:Ntip(tr1),tmp$trees[[1]]$tip.label)
        tips.tr2<-setNames(1:Ntip(tr2),tmp$trees[[2]]$tip.label)
        tips<-sort(rowMeans(cbind(tips.tr1,tips.tr2[names(tips.tr1)])))
        tips[]<-1:length(tips)
    } else if(tangle=="tree1"){
        tips<-setNames(1:Ntip(tr2),tr2$tip.label)
    } else if(tangle=="tree2"){
        tips<-setNames(1:Ntip(tr1),tr1$tip.label)
    }
    layout(matrix(c(1,2,3),1,3),widths=layout)
    plotTree(tr1,color="transparent",ftype="off",tips=tips,type=type,
        nodes=nodes)
    h<-par()$usr[2]
    pp<-get("last_plot.phylo",envir=.PlotPhyloEnv)
    for(i in 1:Ntip(tr1)) lines(c(pp$xx[i],h),rep(pp$yy[i],2),
        lty="dotted")
    if(type=="phylogram"){
        for(i in 1:nrow(tr1$edge))
            lines(pp$xx[tr1$edge[i,]],rep(pp$yy[tr1$edge[i,2]],2),
                lwd=lwd,col=color,lty=lty[2])
        for(i in Ntip(tr1)+1:tr1$Nnode){
            dd<-Children(tr1,i)
            lines(rep(pp$xx[i],length(dd)),pp$yy[dd],lwd=lwd,
                col=color,lty=lty[1])   
        }
    } else if(type=="cladogram"){
        par(ljoin=2)
        for(i in 1:Ntip(tr1)){
            AA<-c(i,Ancestors(tr1,i))
            ii<-sapply(AA[1:(length(AA)-1)],function(x,y) 
                which(y==x), y=tr1$edge[,2])
            lines(c(pp$xx[tr1$edge[ii,2]],pp$xx[Ntip(tr1)+1]),
                c(pp$yy[tr1$edge[ii,2]],pp$yy[Ntip(tr1)+1]),
                lwd=lwd+2,
                col=if(par()$bg=="transparent") "white" else par()$bg,
                lty="solid")
            lines(c(pp$xx[tr1$edge[ii,2]],pp$xx[Ntip(tr1)+1]),
                c(pp$yy[tr1$edge[ii,2]],pp$yy[Ntip(tr1)+1]),
                lwd=lwd,col=color,lty=lty)
        }
    }   
    plot(NA,xlim=c(-1,1),ylim=pp$y.lim,axes=FALSE,xlab="",ylab="")
    text(rep(0,Ntip(tr1)),tips,gsub("_"," ",names(tips)),cex=cex,font=3)
    plotTree(tr2,color="transparent",ftype="off",direction="leftwards",
        tips=tips,type=type)
    h<-par()$usr[1]
    pp<-get("last_plot.phylo",envir=.PlotPhyloEnv)
    for(i in 1:Ntip(tr2)) lines(c(pp$xx[i],h),rep(pp$yy[i],2),
        lty="dotted")
    if(type=="phylogram"){
        for(i in 1:nrow(tr2$edge))
            lines(pp$xx[tr2$edge[i,]],rep(pp$yy[tr2$edge[i,2]],2),
                lwd=lwd,col=color,lty=lty[2])
        for(i in Ntip(tr2)+1:tr2$Nnode){
            dd<-Children(tr2,i)
            lines(rep(pp$xx[i],length(dd)),sort(pp$yy[dd]),
                lwd=lwd,col=color,lty=lty[1])
        }
    } else if(type=="cladogram"){
        for(i in 1:Ntip(tr2)){
            AA<-c(i,Ancestors(tr2,i))
            ii<-sapply(AA[1:(length(AA)-1)],function(x,y) 
                which(y==x), y=tr2$edge[,2])
            lines(c(pp$xx[tr2$edge[ii,2]],pp$xx[Ntip(tr2)+1]),
                c(pp$yy[tr2$edge[ii,2]],pp$yy[Ntip(tr2)+1]),
                lwd=lwd+2,
                col=if(par()$bg=="transparent") "white" else par()$bg,
                lty="solid")
            lines(c(pp$xx[tr2$edge[ii,2]],pp$xx[Ntip(tr2)+1]),
                c(pp$yy[tr2$edge[ii,2]],pp$yy[Ntip(tr2)+1]),
                lwd=lwd,col=color,lty=lty)
        }
    }   
}

Let's give it a try using a phylogeny for salamanders (from Highton & Larson 1979) – and then a “perturbed” two different ways using random NNIs (using phangorn::rNNI. In practice, though, we might use this method to compare (for example) trees estimated using different methodologies or data for the same set of taxa.

set.seed(32)
library(phytools)
library(phangorn)
data(salamanders)
tree1<-force.ultrametric(rNNI(salamanders,3),
    message=FALSE)
tree2<-force.ultrametric(rNNI(salamanders,3),
    message=FALSE)

First let's see what the two trees look like plotted side-by-side.

par(mfrow=c(1,2))
plotTree(tree1)
plotTree(tree2)

Now let's compare them using our cotangleplot method.

cotangleplot(tree1,tree2,lwd=4,tangle="tree2")

This function gives us some different options. For instance, we can choose to tangle tree 1, tree 2, or both. Let's try a different tree of mammals from Garland et al. (1992), but set tangle="both".

data(mammal.tree)
par(bg="black",fg="white")
tree1<-force.ultrametric(rNNI(mammal.tree,3),
    message=FALSE)
tree2<-force.ultrametric(rNNI(mammal.tree,3),
    message=FALSE)
cotangleplot(tree1,tree2,lwd=4,tangle="tree2",use.edge.length=FALSE,
    layout=c(0.43,0.14,0.43),cex=0.9)

Finally, here's a square phylogram style, in which I make the plotted edges semi-transparent.

cotangleplot(tree1,tree2,lwd=3,tangle="tree1",type="phylogram",
    color=make.transparent(palette()[4],0.75),layout=c(0.43,0.14,0.43))

OK, it's not amazing, but it's not bad.