Slanted phylogram (aka. cladogram) for "simmap" object class

I'm working on a slanted phylogram plotting method for objects of class "simmap" - that is, trees with a mapped discrete character.

I have basically been adapting the code of plotSimmap's internal function, plotPhylogram.

Note that I've called this function plotCladogram - to follow the convention in ape::plot.phylo - but, in reality, if the branch lengths contain information then it is properly a slanted phylogram & not a cladogram at all.

The design is such that the horizontal time spent in each color represents the relative branch length in that state along each edge.

It is a work in progress, but here is what I have so far:

plotCladogram<-function(tree,colors=NULL,fsize=1.0,ftype="reg",lwd=2,mar=NULL,
    add=FALSE,offset=NULL,direction="rightwards",xlim=NULL,ylim=NULL,
    nodes="intermediate",tips=NULL,lend=2,asp=NA,plot=TRUE){
    placement<-nodes
    ftype<-which(c("off","reg","b","i","bi")==ftype)-1
    if(!ftype) fsize=0 
    if(is.null(colors)){
        st<-sort(unique(unlist(sapply(tree$maps,names))))
        colors<-palette()[1:length(st)]
        names(colors)<-st
        if(length(st)>1){
            cat("no colors provided. using the following legend:\n")
            print(colors)
        }
    }
    # set offset fudge (empirically determined)
    offsetFudge<-1.37
    # reorder
    cw<-reorderSimmap(tree)
    pw<-reorderSimmap(tree,"postorder")
    # count nodes and tips
    n<-Ntip(cw)
    m<-cw$Nnode
    # Y coordinates for nodes
    Y<-matrix(NA,m+n,1)
    # first, assign y coordinates to all the tip nodes
    if(is.null(tips)) Y[cw$edge[cw$edge[,2]<=n,2]]<-1:n
    else Y[cw$edge[cw$edge[,2]<=n,2]]<-if(is.null(names(tips))) 
        tips[sapply(1:Ntip(cw),function(x,y) which(y==x),y=cw$edge[cw$edge[,2]<=n,2])]
        else tips[gsub(" ","_",cw$tip.label)]
    # get Y coordinates of the nodes
    nodes<-unique(pw$edge[,1])
    for(i in 1:m){
        if(placement=="intermediate"){ 
            desc<-pw$edge[which(pw$edge[,1]==nodes[i]),2]
            Y[nodes[i]]<-(min(Y[desc])+max(Y[desc]))/2
        } else if(placement=="centered"){
            desc<-getDescendants(pw,nodes[i])
            desc<-desc[desc<=Ntip(pw)]
            Y[nodes[i]]<-(min(Y[desc])+max(Y[desc]))/2
        } else if(placement=="weighted"){
            desc<-pw$edge[which(pw$edge[,1]==nodes[i]),2]
            n1<-desc[which(Y[desc]==min(Y[desc]))]
            n2<-desc[which(Y[desc]==max(Y[desc]))]
            v1<-pw$edge.length[which(pw$edge[,2]==n1)]
            v2<-pw$edge.length[which(pw$edge[,2]==n2)]
            Y[nodes[i]]<-((1/v1)*Y[n1]+(1/v2)*Y[n2])/(1/v1+1/v2)
        } else if(placement=="inner"){
            desc<-getDescendants(pw,nodes[i])
            desc<-desc[desc<=Ntip(pw)]
            mm<-which(abs(Y[desc]-median(Y[1:Ntip(pw)]))==min(abs(Y[desc]-
                median(Y[1:Ntip(pw)]))))
            if(length(mm>1)) mm<-mm[which(Y[desc][mm]==min(Y[desc][mm]))]
            Y[nodes[i]]<-Y[desc][mm]
        }
    }
    # compute node heights
    H<-nodeHeights(cw)
    # open plot
    par(mar=mar)
    if(is.null(offset)) offset<-0.2*lwd/3+0.2/3
    if(!add) plot.new()
    ###
    if(is.null(xlim)){
        pp<-par("pin")[1]
        sw<-fsize*(max(strwidth(cw$tip.label,units="inches")))+
            offsetFudge*fsize*strwidth("W",units="inches")
        alp<-optimize(function(a,H,sw,pp) (a*1.04*max(H)+sw-pp)^2,H=H,sw=sw,pp=pp,
            interval=c(0,1e6))$minimum
        xlim<-if(direction=="leftwards") c(min(H)-sw/alp,max(H)) else c(min(H),max(H)+sw/alp)
    }
    if(is.null(ylim)) ylim=range(Y)
    if(direction=="leftwards") H<-max(H)-H
    plot.window(xlim=xlim,ylim=ylim,asp=asp)
    if(plot){
        ####
        for(i in 1:nrow(cw$edge)){
            x<-H[i,1]
            y<-Y[cw$edge[i,1]]
            m<-(Y[cw$edge[i,2]]-Y[cw$edge[i,1]])/(H[i,2]-H[i,1])
            if(is.finite(m)){
                for(j in 1:length(cw$maps[[i]])){
                    if(direction=="leftwards")
                        lines(c(x,x-cw$maps[[i]][j]),c(y,y-cw$maps[[i]][j]*m),
                            col=colors[names(cw$maps[[i]])[j]],lwd=lwd,lend=lend)
                    else lines(c(x,x+cw$maps[[i]][j]),c(y,y+cw$maps[[i]][j]*m),
                        col=colors[names(cw$maps[[i]])[j]],lwd=lwd,lend=lend)
                    x<-x+if(direction=="leftwards") -cw$maps[[i]][j] else cw$maps[[i]][j]
                    y<-y+if(direction=="leftwards") -m*cw$maps[[i]][j] else m*cw$maps[[i]][j]
                    j<-j+1
                }
            } else {
                lines(rep(x,2),Y[cw$edge[i,]],col=colors[names(cw$maps[[i]])[1]],lwd=lwd,
                    lend=lend)
            }
        }
        if(direction=="leftwards") pos<-if(par()$usr[1]>par()$usr[2]) 4 else 2
        if(direction=="rightwards") pos<-if(par()$usr[1]>par()$usr[2]) 2 else 4
        for(i in 1:n) if(ftype) text(H[which(cw$edge[,2]==i),2],Y[i],cw$tip.label[i],pos=pos,
            offset=offset,cex=fsize,font=ftype)
    }
    PP<-list(type="phylogram",use.edge.length=TRUE,node.pos=1,
        show.tip.label=if(ftype) TRUE else FALSE,show.node.label=FALSE,
        font=ftype,cex=fsize,adj=0,srt=0,no.margin=FALSE,label.offset=offset,
        x.lim=xlim,y.lim=ylim,
        direction=direction,tip.color="black",Ntip=Ntip(cw),Nnode=cw$Nnode,
        edge=cw$edge,xx=sapply(1:(Ntip(cw)+cw$Nnode),
        function(x,y,z) y[match(x,z)],y=H,z=cw$edge),yy=Y[,1])
    assign("last_plot.phylo",PP,envir=.PlotPhyloEnv)
}

We can try it as follows:

data(anoletree)
ss<-mapped.states(anoletree)
colors<-setNames(palette()[1:length(ss)],ss)
plotCladogram(anoletree,colors,fsize=0.7,mar=rep(0.1,4),ftype="i",
    ylim=c(-1,Ntip(anoletree)))
add.simmap.legend(colors=colors,prompt=FALSE,x=0,y=-1,vertical=FALSE)

Note that the line crossing is actually impossible to avoid for slanted phylograms (or, rather, it is impossible to guarantee there is no line crossing). ape::plot.phylo suffers the same problem.

plot.phylo(anoletree,type="cladogram",no.margin=TRUE,cex=0.7)

Finally, the aliasing that you see on the slanted edges can be avoided by exporting in vector format (such as PDF), which I always recommend for presentations & publications.

That's it for now.

Complete re-write of phytools densityTree function (& many new options)

Earlier today I posted about plotting multiple trees with conflicting topologies in the same plotting device.

There is theoretically a function in phytools called densityTree that automates at least some of these options; however in the process of preparing this post I discovered that the function appears to be a little bit broken.

In particular, this is what we get when we apply densityTree to the 9 trees from our previous example:

library(phytools)
trees

## 9 phylogenetic trees

densityTree(trees)

Oops! That's no good.

Anyhow, I'm going to try & fix that as well as add some additional features to the function, and in so doing I intend to employ some of the tricks that I used in today's post.

Here is the re-written function:

## function to be used internally
rescaleTree<-function(tree,scale){
    tree$edge.length<-tree$edge.length/max(nodeHeights(tree))*scale
    tree
}
## densityTree function
densityTree<-function(trees,colors="blue",alpha=NULL,method="plotTree",
    fix.depth=FALSE,use.edge.length=TRUE,compute.consensus=TRUE,
    use.gradient=FALSE,show.axis=TRUE,...){
    N<-length(trees)
    if(any(sapply(trees,function(x) is.null(x$edge.length)))) 
        use.edge.length<-FALSE
    if(!use.edge.length) trees<-lapply(trees,compute.brlen)
    h<-sapply(trees,function(x) max(nodeHeights(x)))
    if(fix.depth){
        if(method=="plotTree"){
            trees<-lapply(trees,rescaleTree,mean(h))
            class(trees)<-"multiPhylo"
        } else if(method=="plotSimmap"){ 
            trees<-rescaleSimmap(trees,depth=mean(h))
            print(class(trees))
            print(class(trees[[1]]))
        }
        h<-sapply(trees,function(x) max(nodeHeights(x)))
    }
    tips<-setNames(1:Ntip(trees[[1]]), 
        if(compute.consensus) untangle(consensus(trees),
            "read.tree")$tip.label 
        else trees[[1]]$tip.label)
    if(is.null(alpha)) alpha<-max(c(1/N,0.01))
    args<-list(...)
    args$direction<-"leftwards"
    args$tips<-tips
    args$add<-FALSE
    if(is.null(args$nodes)) args$nodes<-"inner"
    if(is.null(args$mar)) args$mar<-if(show.axis) c(4.1,1.1,1.1,1.1) else rep(1.1,4)
    if(is.null(args$ftype)) args$ftype<-"i"
    if(!use.gradient){
        plotTree(trees[[which(h==max(h))[1]]],direction="leftwards",mar=args$mar,
            plot=FALSE)
        args$xlim<-get("last_plot.phylo",envir=.PlotPhyloEnv)$x.lim[2:1]
        if(method=="plotTree"){
            args$color<-make.transparent(colors[1],alpha)
            for(i in 1:length(trees)){
                args$tree<-trees[[i]]
                do.call(plotTree,args)
                if(i==1){ 
                    if(show.axis) axis(1)
                    args$ftype<-"off"
                    args$add<-TRUE
                }
            }
        } else if(method=="plotSimmap"){
            states<-sort(unique(as.vector(mapped.states(trees))))
            print(states)
            print(colors)
            if(length(colors)!=length(states)){
                colors<-setNames(c("grey",palette()[2:length(states)]),
                    states)
            }
            colors<-sapply(colors,make.transparent,alpha=alpha)
            args$colors<-colors
            for(i in 1:length(trees)){
                args$tree<-trees[[i]]
                do.call(plotSimmap,args)
                if(i==1){ 
                    if(show.axis) axis(1)
                    args$ftype<-"off"
                    args$add<-TRUE
                }
            }
        }
    } else if(use.gradient){
        nulo<-capture.output(rf<-multiRF(trees),type="message")
        mds<-cmdscale(rf,k=1)[,1]
        trees<-trees[order(mds)]
        args$ylim<-c(0,Ntip(trees[[1]])+1)
        plotTree(trees[[which(h==max(h))[1]]],direction="leftwards",mar=args$mar,
            ylim=args$ylim,plot=FALSE)
        args$xlim<-get("last_plot.phylo",envir=.PlotPhyloEnv)$x.lim[2:1]
        colors<-sapply(rainbow(n=length(trees)),make.transparent,alpha=alpha)
        ftype<-args$ftype
        for(i in 1:length(trees)){
            y.shift<-(i-median(1:length(trees)))/length(trees)/2
            args$tree<-trees[[i]]
            args$tips<-tips+y.shift
            args$color<-colors[i]
            args$ftype<-if(i==floor(median(1:length(trees)))) ftype else "off"
            do.call(plotTree,args)
            if(i==1){ 
                if(show.axis) axis(1)
                args$ftype<-"off"
                args$add<-TRUE
            }
        }
    }
}

We can try it out using various options. Note, to follow along you will have to update phytools from GitHub.

First, the basic method:

densityTree(trees)

or fix the total depth to be a constant value:

densityTree(trees,color="red",fix.depth=TRUE)

We can also use a color gradient to be able to better distinguish the different topologies. I demonstrated this option today. I also added the feature that the method first sorts the trees (in a way) by computing the Robinson-Foulds distances between topologies & the uses MDS to sort the trees into one dimension. That means that plotted trees with similar gradient colors tend to have more similar topologies!

densityTree(trees,use.gradient=TRUE,fix.depth=TRUE,alpha=0.2)

## Some trees are rooted. Unrooting all trees.

This function also allows a mapped discrete character by changing the argument method to plotSimmap. Let's try it:

mtrees

## 90 phylogenetic trees with mapped discrete characters

densityTree(mtrees,colors=setNames(c("blue","red"),letters[1:2]),method="plotSimmap",
    nodes="intermediate")

Lastly, let's try it on a larger set of slightly bigger trees to see how it works:

tt

## 99 phylogenetic trees

densityTree(tt)

densityTree(tt,use.gradient=TRUE,alpha=0.05,fix.depth=TRUE)

## Some trees are rooted. Unrooting all trees.

Cool.