I haven't been posting to actively to this blog lately I have been busy with teaching and with my 'new' package learnPopGen as well as with the various shiny web interfaces that I have built for a number of the functions of the package.
Nonetheless, this morning for some reason I was thinking about the problem
of plotting unrooted trees (which is actually much more complicated that
plotting rooted trees in their various configurations, as discussed in
Chapter 34 of Felsenstein
2004),
and I suddenly realized that to plot trees with a mapped distrete character
(objects of class "simmap") or trees with a mapped continuous
character (objects of class "contMap") I could quite easily take
advantage of ape::plot.phylo and simply project my tree on top
of the coordinates of the unrooted tree plotted by plot.phylo in
which this problem has already been solved.
This is what I mean.
library(phytools)
set.seed(777)
N<-16
tree<-rtree(n=N,tip.label=LETTERS[1:N])
Q<-matrix(c(-1,1,1,-1),2,2,dimnames=list(letters[1:2],letters[1:2]))
x<-sim.Mk(tree,Q)
x
## L J I H K B E M C N G D F A O P
## a b b a a a a b b a b b b b b b
## Levels: a b
map.tree<-make.simmap(tree,x)
## make.simmap is sampling character histories conditioned on the transition matrix
##
## Q =
## a b
## a -1.038569 1.038569
## b 1.038569 -1.038569
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## a b
## 0.5 0.5
## Done.
## here's where it gets interesting:
plot.phylo(map.tree,type="unrooted",lab4ut="axial",no.margin=TRUE,
label.offset=0.05)
obj<-get("last_plot.phylo",envir=.PlotPhyloEnv)
tip.pos<-matrix(c(obj$xx[1:Ntip(map.tree)],obj$yy[1:Ntip(map.tree)]),
Ntip(map.tree),2,dimnames=list(map.tree$tip.label,NULL))
node.pos<-matrix(c(obj$xx[1:map.tree$Nnode+Ntip(map.tree)],
obj$yy[1:map.tree$Nnode+Ntip(map.tree)]),map.tree$Nnode,2,
dimnames=list(1:map.tree$Nnode+Ntip(map.tree),NULL))
phylomorphospace(map.tree,tip.pos,A=node.pos,lwd=3,
colors=setNames(c("blue","red"),sort(unique(getStates(map.tree)))),
add=TRUE,ftype="off",node.size=c(0,1),node.by.map=TRUE)
We can do something similar with a "contMap" object as follows:
set.seed(777)
y<-fastBM(tree)
cmap<-contMap(tree,y,plot=FALSE,res=1000)
par(lend=2)
plot.phylo(cmap$tree,type="unrooted",lab4ut="axial",no.margin=TRUE,
label.offset=0.05,edge.width=7)
obj<-get("last_plot.phylo",envir=.PlotPhyloEnv)
tip.pos<-matrix(c(obj$xx[1:Ntip(cmap$tree)],obj$yy[1:Ntip(cmap$tree)]),
Ntip(cmap$tree),2,dimnames=list(cmap$tree$tip.label,NULL))
node.pos<-matrix(c(obj$xx[1:cmap$tree$Nnode+Ntip(cmap$tree)],
obj$yy[1:cmap$tree$Nnode+Ntip(cmap$tree)]),cmap$tree$Nnode,2,
dimnames=list(1:cmap$tree$Nnode+Ntip(cmap$tree),NULL))
phylomorphospace(cmap$tree,tip.pos,A=node.pos,lwd=5,
colors=cmap$cols,add=TRUE,ftype="off",
node.size=c(0,0))
add.color.bar(h<-2,cmap$cols,title="trait value",
lims=NULL,digits=3,direction="upwards",
subtitle="",lwd=15,x=0.1,y=1.75,prompt=FALSE)
LWD<-diff(par()$usr[1:2])/dev.size("px")[1]
nticks<-6
Y<-cbind(seq(1.75,1.75+h,length.out=nticks),
seq(1.75,1.75+h,length.out=nticks))
X<-cbind(rep(0.1+LWD*17/2,nticks),
rep(0.1+LWD*17/2+0.02*h,nticks))
for(i in 1:nrow(Y)) lines(X[i,],Y[i,])
ticks<-seq(cmap$lims[1],cmap$lims[2],length.out=nticks)
text(x=X[,2],y=Y[,2],round(ticks,2),pos=4,cex=0.8)
That's pretty cool (IMO).
