Plotting an upwards facing traitgram using phytools

Earlier today I tweeted an illustration giving a vertical oriented traitgram (projection of the tree into a phenotype space) with a mapped discrete trait.

Fitting a multi-rate Brownian motion model with hidden states using the discrete approximation in #Rstats #phytools. 👉 https://t.co/Cj88iPdRjt pic.twitter.com/5TCOwEoE7L

— Liam Revell (@phytools_liam) June 12, 2024

In the popular phytools function phenogram the traitgram projection temporally runs from left to right in the plot. So how did I create this?

To see, let’s load the dataset from Kirk & Kay (2004) consisting of a phylogeny of primates and various phenotypic characteristics.

library(phytools)

data(primate.tree)
primate.tree

## 
## Phylogenetic tree with 90 tips and 89 internal nodes.
## 
## Tip labels:
##   Allenopithecus_nigroviridis, Cercopithecus_mitis, Cercopithecus_cephus, Cercopithecus_petaurista, Erythrocebus_patas, Chlorocebus_aethiops, ...
## 
## Rooted; includes branch lengths.

data(primate.data)
head(primate.data)

##                                  Group Skull_length Optic_foramen_area Orbit_area
## Allenopithecus_nigroviridis Anthropoid         98.5                7.0      298.7
## Alouatta_palliata           Anthropoid        109.8                5.3      382.3
## Alouatta_seniculus          Anthropoid        108.0                8.0      359.4
## Aotus_trivirgatus           Anthropoid         60.5                3.1      297.4
## Arctocebus_aureus            Prosimian         49.5                1.2      134.8
## Arctocebus_calabarensis      Prosimian         53.8                1.6      156.6
##                             Activity_pattern Activity_pattern_code
## Allenopithecus_nigroviridis          Diurnal                     0
## Alouatta_palliata                    Diurnal                     0
## Alouatta_seniculus                   Diurnal                     0
## Aotus_trivirgatus                  Nocturnal                     2
## Arctocebus_aureus                  Nocturnal                     2
## Arctocebus_calabarensis            Nocturnal                     2

Next, let’s pull out the log of skull length and create a traitgram for this character in the standard way.

lnSkull<-setNames(log(primate.data$Skull_length),
  rownames(primate.data))
phenogram(primate.tree,lnSkull,ftype="off")

Now, much like other phylogeny plotting functions in ape and phytools, phytools::phenogram exports an object called "last_plot.phylo" to the environment .PlotPhyloEnv so let’s retrieve it.

pp<-get("last_plot.phylo",envir=.PlotPhyloEnv)
str(pp)

## List of 20
##  $ type           : chr "phenogram"
##  $ use.edge.length: logi TRUE
##  $ node.pos       : num 1
##  $ show.tip.label : logi TRUE
##  $ show.node.label: logi FALSE
##  $ font           : num 0
##  $ cex            : num 0
##  $ adj            : num 0
##  $ srt            : NULL
##  $ no.margin      : logi FALSE
##  $ label.offset   : num 0.2
##  $ x.lim          : num [1:2] -2.92 75.92
##  $ y.lim          : num [1:2] 3.35 5.68
##  $ direction      : NULL
##  $ tip.color      : chr "black"
##  $ Ntip           : int 90
##  $ Nnode          : int 89
##  $ edge           : int [1:178, 1:2] 91 92 93 94 95 96 97 97 98 99 ...
##  $ xx             : Named num [1:179] 73 73 73 73 73 ...
##   ..- attr(*, "names")= chr [1:179] "1" "2" "3" "4" ...
##  $ yy             : Named num [1:179] 4.59 4.69 4.56 4.53 4.83 ...
##   ..- attr(*, "names")= chr [1:179] "1" "2" "3" "4" ...

This object contains all the coordinates of our plotted projection, including our x and y coordinates. To plot our traitgram vertically instead of horizontally, we just have to flip them!

plot.new()
par(mar=c(5.1,5.1,1.1,1.1))
plot.window(xlim=range(pp$yy),ylim=range(pp$xx))
for(i in 1:nrow(pp$edge)) lines(pp$yy[pp$edge[i,]],
  pp$xx[pp$edge[i,]])
axis(1,cex.axis=0.8)
axis(2,at=max(pp$xx)-seq(0,60,by=20),
  labels=seq(0,60,by=20),las=1,cex.axis=0.8)
title(xlab="log(skull length)",
  ylab="time before present")

But wait – in the example I included a mapped discrete character (actually, a set of simulated rate regimes). Our dataset contains a discrete trait, diel activity pattern, so let’s generate a map on the tree using simmap.

activity<-setNames(primate.data$Activity_pattern,
  rownames(primate.data))
er_fit<-fitMk(primate.tree,activity,model="ER")
smap_primate.tree<-simmap(er_fit,nsim=1)
smap_primate.tree

## 
## Phylogenetic tree with 90 tips and 89 internal nodes.
## 
## Tip labels:
## 	Allenopithecus_nigroviridis, Cercopithecus_mitis, Cercopithecus_cephus, Cercopithecus_petaurista, Erythrocebus_patas, Chlorocebus_aethiops, ...
## 
## The tree includes a mapped, 3-state discrete character
## with states:
## 	Cathemeral, Diurnal, Nocturnal
## 
## Rooted; includes branch lengths.

phenogram(primate.tree,lnSkull,
  spread.cost=c(1,0),ftype="i",fsize=0.4)

## Optimizing the positions of the tip labels...

OK, now let’s add some more features.

pp<-get("last_plot.phylo",envir=.PlotPhyloEnv)

tt<-map.to.singleton(smap_primate.tree)

cols<-setNames(hcl.colors(n=3)[c(2,3,1)],
  levels(activity))
plot.new()
par(mar=c(5.1,5.1,1.1,1.1),lend=1)
plot.window(xlim=c(25,exp(max(pp$yy))),ylim=pp$x.lim,
  log="x")
for(i in 1:nrow(pp$edge)) lines(exp(pp$yy[pp$edge[i,]]),
  pp$xx[pp$edge[i,]],
  col=make.transparent(cols[names(tt$edge.length)[i]],0.8),
  lwd=3)
axis(1,at=c(25,50,100,200),cex.axis=0.8)
axis(2,at=max(pp$xx)-seq(0,60,by=20),
  labels=seq(0,60,by=20),las=1,cex.axis=0.8)
title(xlab="skull length (mm)",
  ylab="time before present")
xx<-exp(seq(log(25),max(pp$yy),length.out=Ntip(tt)))
sorted<-exp(sort(lnSkull))
h<-max(nodeHeights(tt))
shift<-0.4*strwidth("W")*(diff(par()$usr[3:4])/
    diff(par()$usr[1:2]))
for(i in 1:Ntip(tt)){
  nn<-names(sorted)[i]
  text(x=xx[i],y=1.1*h,gsub("_"," ",nn),
    cex=0.4,font=3,pos=4,offset=0,srt=90)
  ntip<-which(tt$tip.label==nn)
  lines(x=c(exp(pp$yy[ntip]),xx[i]),
    y=c(pp$xx[ntip],1.1*h),lty="dotted",
    col="#36454F")
}
legend("bottomright",levels(activity),lwd=4,
  col=cols,bty="n")

That’s it, folks!