How to plot tip & node labels without tiplabels or nodelabels

A recent tweet from Andrea Berardi asked how to adjust offset of pie graphs or circles plotted as tip labels in R.

Has anyone figured out how to plot a tree in R with tip circles & offset tip circles over? Been adjusting manually in Illustrator :( “offset” seems to adj text but not tips plot(tree, fsize=0.8, ftype=“i”, pie=tree$ace, colors=cols, piecol=cols, cex=0.5, offset=1) #phytools #ape pic.twitter.com/bT5mcTVWIx

— Andrea Berardi (@aeberardi) June 8, 2022

Lots of others responded with excellent suggestions, but I nonetheless thought it might be handy to provide another, alternative solution that doesn't use ape::nodelabels or ape::tiplabels at all!

You can find similar tips & tricks in my upcoming book with Luke Harmon.

First, I'm going to pull an example. These are data for diel activity pattern in primates from Kirk & Kay (2004) – that we also use in our book, by the way!

library(phytools)
## load tree
primate.tree<-read.tree(file="http://www.phytools.org/Rbook/3/primateEyes.phy")
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.

## load data
primate.data<-read.csv(file="http://www.phytools.org/Rbook/3/primateEyes.csv",
    row.names=1,stringsAsFactors=TRUE)
head(primate.data)

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

activity.pattern<-setNames(primate.data$Activity_pattern,rownames(primate.data))
head(activity.pattern)

## Allenopithecus_nigroviridis           Alouatta_palliata 
##                     Diurnal                     Diurnal 
##          Alouatta_seniculus           Aotus_trivirgatus 
##                     Diurnal                   Nocturnal 
##           Arctocebus_aureus     Arctocebus_calabarensis 
##                   Nocturnal                   Nocturnal 
## Levels: Cathemeral Diurnal Nocturnal

Neat.

(As an extra step here, which sometimes helps with plotting methods later, I'm going to run the function phytools::untangle. This has the effect of reordering the tip labels of the tree in what ape calls “cladewise” fashion, and can be pretty useful. This is not necessary here, because our tree was written in simple Newick format – but can be critical when we read trees from NEXUS format with a translation table!)

primate.tree<-untangle(primate.tree,"read.tree")

Next, I'm going to run stochastic character mapping, using the ARD model.

simmap.trees<-make.simmap(primate.tree,activity.pattern,model="ARD",
    pi="fitzjohn",nsim=100)

## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##              Cathemeral      Diurnal    Nocturnal
## Cathemeral -0.050529521  0.050529521  0.000000000
## Diurnal     0.002306398 -0.004367451  0.002061052
## Nocturnal   0.002586387  0.003565036 -0.006151422
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##  Cathemeral     Diurnal   Nocturnal 
## 0.003675779 0.021703763 0.974620458

## Done.

Here are my stochastic character mapped trees:

par(mfrow=c(10,10))
cols<-setNames(c("brown","goldenrod","black"),levels(activity.pattern))
plot(simmap.trees,cols,ftype="off",lwd=1)

par(mfrow=c(1,1))

Already, now let's compute the marginal frequencies at all the internal nodes of the tree using summary.

obj<-summary(simmap.trees)
plot(obj,colors=cols,ftype="i",fsize=0.7)
legend("bottomleft",levels(activity.pattern),pch=21,pt.bg=cols,
    pt.cex=2,bty="n")

I mean, this is OK, but we actually have limitless ability to add graphical elements to a plotted phylogeny in R, as long as we know where in our plot to add them!

To get this information, we want to pull the object last_plot.phylo from the R environment .PlotPhyloEnv. It doesn't matter if you don't really know what that means, or even what an environment is in R, you can still take advantage of this trick.

To start we'll re-plot our phylogeny, and then leave this graphical device open.

plotTree(primate.tree,ftype="i",fsize=0.5,lwd=1,offset=0.5)

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

This object that we just pulled, contains (among other things) the x and y coordinates of all the nodes and tips of the phylogeny!

str(pp)

## List of 20
##  $ type           : chr "phylogram"
##  $ use.edge.length: logi TRUE
##  $ node.pos       : num 1
##  $ show.tip.label : logi TRUE
##  $ show.node.label: logi FALSE
##  $ font           : num 3
##  $ cex            : num 0.5
##  $ adj            : num 0
##  $ srt            : num 0
##  $ no.margin      : logi FALSE
##  $ label.offset   : num 0.5
##  $ x.lim          : num [1:2] 0 90.5
##  $ y.lim          : num [1:2] 1 90
##  $ direction      : chr "rightwards"
##  $ 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             : num [1:179] 73 73 73 73 73 ...
##  $ yy             : num [1:179] 1 2 3 4 5 6 7 8 9 10 ...

To see that this is true, let's just plot simple points at all the nodes and tips as follows:

points(pp$xx,pp$yy)

Now, let's attempt the simple task of setting the colors of the points shown at the tips by the states of our discrete character for each tip.

h<-max(nodeHeights(primate.tree))
plotTree(primate.tree,ftype="i",fsize=0.5,lwd=1,offset=0.5)
points(pp$xx[1:Ntip(primate.tree)]+0.01*h,pp$yy[1:Ntip(primate.tree)],pch=16,
    col=cols[activity.pattern[primate.tree$tip.label]])

Let's add pie charts to the nodes. We're not going to do this, however, using ape::nodelabels. Instead, we'll use plotrix::floating.pies. plotrix is a dependency of phytools, so you already have it installed.

for(i in 1:primate.tree$Nnode){
    plotrix::floating.pie(pp$xx[i+Ntip(primate.tree)],
        pp$yy[i+Ntip(primate.tree)],radius=0.015*h,
        x=obj$ace[i,],col=cols,border="transparent")
}

Finally, add a legend.

legend("bottomleft",levels(activity.pattern),pch=16,col=cols,
    pt.cex=2,bty="n")

Just to see how this gives us more flexibility, here's an example in which I offset the tip labels, and then also adjust the internal node label size depending on whether or not there's uncertainty (max. PP < 0.99) for each state.

h<-max(nodeHeights(primate.tree))
plotTree(primate.tree,ftype="i",fsize=0.5,lwd=1,offset=0.5)
points(pp$xx[1:Ntip(primate.tree)]+0.01*h,pp$yy[1:Ntip(primate.tree)],pch=16,
    col=cols[activity.pattern[primate.tree$tip.label]])
for(i in 1:primate.tree$Nnode){
    plotrix::floating.pie(pp$xx[i+Ntip(primate.tree)],
        pp$yy[i+Ntip(primate.tree)],
        radius=if(max(obj$ace[i,]^2)<0.98) 0.02*h else 0.01*h,
        x=obj$ace[i,],col=cols,border="transparent")
}
legend("bottomleft",levels(activity.pattern),pch=16,col=cols,
    pt.cex=2,bty="n")

That's alright.