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).

## No comments:

## Post a Comment

Note: due to the very large amount of spam, all comments are now automatically submitted for moderation.