Earlier today I tweeted about graphing discrete & continuous traits at the tips of a plotted tree in “arc” or “fan” style. The way I was doing this was just by added thickened lines that “grow” out of the tips of the tree.
A response from Ted Stankowich, however, has persuaded me to re-do this using polygons – even though the geometry is substantially more complicated. (It still requires nothing more than high school trig – but nonetheless!)
We've wrestled with whether to put the species names at the branch tips and the color rings outside of them. This plot is cool looking, but it seems like it's very hard to trace back the species name to its tip through the traits. Having white borders between markers might help? pic.twitter.com/uXdJNTh1wR
— Ted Stankowich 🦨🐆🦌🦔🦏🦡 (@CSULBMammalLab) March 27, 2024
The main advantages of polygons is that it allows us to incorporate “spacers” between adjacent species values, and easily widen the polygons as we go from the inside to the outside of our plotted traits.
OK. Here’s the code.
plotFanTree.wTraits<-function(tree,X,type=c("arc","fan"),
...){
X<-if(is.vector(X)) as.matrix(X[tree$tip.label]) else
X[tree$tip.label,]
h<-max(nodeHeights(tree))
d<-min(ncol(X)*0.07*h,h)
type<-type[1]
if(!(type%in%c("arc","fan"))) type<-"fan"
ftype<-if(hasArg(ftype)) list(...)$ftype else "i"
fsize<-if(hasArg(fsize)) list(...)$fsize else 0.5
part<-if(hasArg(part)) list(...)$part else
min(0.99,(Ntip(tree)-2)/Ntip(tree))
arc_height<-if(hasArg(arc_height)) list(...)$arc_height else
0.7
spacer<-if(hasArg(spacer)) list(...)$spacer else 0.025
spacer<-spacer*(2*pi*part/(Ntip(tree)-1))/2
xlim<-if(hasArg(xlim)) list(...)$xlim else NULL
ylim<-if(hasArg(ylim)) list(...)$ylim else NULL
if(hasArg(colors)) colors<-list(...)$colors
else {
colors<-list()
for(i in 1:ncol(X)){
if(is.numeric(X[,i])){
colors[[i]]<-setNames(hcl.colors(n=100),1:100)
} else {
if(!is.factor(X[,i])) X[,i]<-as.factor(X[,i])
colors[[i]]<-setNames(
palette.colors(n=length(levels(X[,i]))),
levels(X[,i]))
}
}
}
tt<-tree
tt$edge.length[which(tt$edge[,2]<=Ntip(tt))]<-
tt$edge.length[which(tt$edge[,2]<=Ntip(tt))]+d
plotTree(tt,type=type,ftype=ftype,fsize=fsize,
part=part,color="transparent",
arc_height=arc_height*h/max(nodeHeights(tt)),
xlim=xlim,ylim=ylim)
pp<-get("last_plot.phylo",envir=.PlotPhyloEnv)
outer_rad<-max(pp$xx)
plotTree(tree,type=type,ftype="off",part=part,
lwd=1,add=TRUE,xlim=pp$x.lim,ylim=pp$y.lim,
arc_height=arc_height,ftype="off")
pp<-get("last_plot.phylo",envir=.PlotPhyloEnv)
inner_rad<-max(pp$xx)
par(lend=3)
for(i in 1:ncol(X)){
if(is.numeric(X[,i])){
x_seq<-seq(min(X[,i]),max(X[,i]),length.out=100)
x_ind<-sapply(X[,i],function(x,y) which.min((x-y)^2),
y=x_seq)
colors[[i]]<-colorRampPalette(colors[[i]])(n=100)
cols<-colors[[i]][x_ind]
} else {
cols<-colors[[i]][X[tree$tip.label,i]]
}
for(j in 1:Ntip(tree)){
start<-if(pp$xx[j]>0)
(i-1)*(d/ncol(X))+(2/7)*(d/ncol(X)) else
-((i-1)*(d/ncol(X))+(2/7)*(d/ncol(X)))
end<-if(pp$xx[j]>0) i*d/ncol(X) else -i*d/ncol(X)
th<-atan(pp$yy[j]/pp$xx[j])
theta<-(2*pi*part/(Ntip(tree)-1))/2-spacer
sign<-if(pp$xx[j]>0) 1 else -1
H1<-(sign*inner_rad+start)/cos(theta)
H2<-(sign*inner_rad+end)/cos(theta)
th_up<-th+theta
th_down<-th-theta
x<-c(H1*cos(th_down),H2*cos(th_down),
H2*cos(th_up),H1*cos(th_up))
y<-c(H1*sin(th_down),H2*sin(th_down),
H2*sin(th_up),H1*sin(th_up))
polygon(x,y,col=cols[j],border=FALSE)
}
}
invisible(colors)
}
First, let’s apply it to the example from earlier today to see that it still works! Here, I’ll set the optional argument spacer = 0
so that we have no white space between species levels of each character.
library(phytools)
data("liolaemid.data")
data("liolaemid.tree")
liolaemid.data[liolaemid.tree$tip.label,]->liolaemid.data
colors<-list(
c("blue","white","red"),
terrain.colors(n=10),
setNames(c("#F0EAD6","#DF536B"),c("O","V")))
cols<-plotFanTree.wTraits(liolaemid.tree,
liolaemid.data[,3:1],lwd=12,colors=colors,ftype="off",
spacer=0)
legend(x=0,y=0.7*max(nodeHeights(liolaemid.tree)),
names(colors[[3]]),lwd=8,col=colors[[3]],
title="parity mode",bty="n",xjust=0.5,yjust=0.5)
add.color.bar(1.5*max(nodeHeights(liolaemid.tree)),cols[[2]],
title="maximum altitude (m)",
lims=range(liolaemid.data[,2]),digits=2,prompt=FALSE,
x=-0.75*max(nodeHeights(liolaemid.tree)),
y=0.2*max(nodeHeights(liolaemid.tree)),subtitle="",
lwd=8,outline=FALSE)
add.color.bar(1.5*max(nodeHeights(liolaemid.tree)),cols[[1]],
title="environmental temp.",
lims=range(liolaemid.data[,3]),digits=2,prompt=FALSE,
x=-0.75*max(nodeHeights(liolaemid.tree)),
y=-0.15*max(nodeHeights(liolaemid.tree)),subtitle="",
lwd=8,outline=FALSE)