Creating a phylogenetic tree in R using phytools::bind.tip in interactive mode

I just made a short video showing how to draw a tree “free-hand” in R using the interactive=TRUE mode of the phytools function bind.tip.

The video is here:

The syntax for this exercise (in case it cannot be easily read from the video) is as follows:

## first load our packages:
library(phytools)

## next create a two taxon tree that subtends the root:
tree<-pbtree(n=2,tip.label=c("monkey","human"))

## now interactively bind the addt'l tips we want onto this tree:
tips<-c("bonobo","chimp","gorilla","orangutan")
for(i in 1:length(tips)) tree<-bind.tip(tree,tips[i],interactive=TRUE)

## strip edge lengths (here they were arbitrary):
tree$edge.length<-NULL

## plot our tree:
roundPhylogram(tree)

## done.

The objected created in memory is just a regular "phylo" class object, so it can be plotted:

plotTree(tree,fsize=1.2)

written to file:

write.tree(tree)

## [1] "(monkey,(((human,(bonobo,chimp)),gorilla),orangutan));"

etc.

In a subsequent post I will demonstrate how to add edge lengths based on divergence dates associated with individual nodes of the tree.

plotTree.barplot with more user options, including stacked bars for multiple input traits

In response to a user comment I have added a feature to plotTree.barplot to permit multiple traits or values to be plotted for each species - either using stacked bars, or by plotting the bars side-by-side.

The update can be seen here. The problem turned out to be a little bit mroe complicated then it seemed it first, mainly because of a quirk in barplot that means that the object returned by boxplot (which I then use to correctly space the tips of the plotted tree) differs depending on whether the input data come in the form of a vector of single values or a matrix, &, if the latter, whether the argument beside=TRUE is selected or not.

Here's a quick & dirty demo. These data are (simulated to represent) relative frequencies - said dietary fraction of three different food types - so stacking of our bars is appropriate:

library(phytools)
packageVersion("phytools")

## [1] '0.5.70'

tree

## 
## Phylogenetic tree with 26 tips and 25 internal nodes.
## 
## Tip labels:
##  A, B, C, D, E, F, ...
## 
## Rooted; includes branch lengths.

X

##        plant  vertebrate invertebrate
## A 0.31086049 0.212175943   0.47696356
## B 0.15404215 0.440976975   0.40498087
## C 0.26982034 0.450834511   0.27934515
## D 0.12260206 0.358475976   0.51892197
## E 0.29979068 0.129743168   0.57046615
## F 0.20163877 0.440001107   0.35836012
## G 0.05465514 0.410528162   0.53481670
## H 0.47392116 0.121520364   0.40455848
## I 0.44344594 0.101251541   0.45530252
## J 0.24839375 0.268691107   0.48291514
## K 0.12309160 0.101699889   0.77520851
## L 0.24010676 0.065417274   0.69447597
## M 0.23424489 0.226431850   0.53932326
## N 0.40030662 0.057542453   0.54215092
## O 0.15332681 0.023255378   0.82341781
## P 0.44389924 0.491496743   0.06460401
## Q 0.27010640 0.292077796   0.43781581
## R 0.21643656 0.208244073   0.57531937
## S 0.19650159 0.317376110   0.48612230
## T 0.01261411 0.427646672   0.55973922
## U 0.32707264 0.164682573   0.50824478
## V 0.28931830 0.126987763   0.58369394
## W 0.08407750 0.252999348   0.66292315
## X 0.28631262 0.004105909   0.70958147
## Y 0.38118921 0.206806916   0.41200387
## Z 0.07941718 0.074216139   0.84636668

plotTree.barplot(tree,X)

Or, we can plot the bars side-by-side:

plotTree.barplot(tree,X,args.barplot=list(beside=TRUE,xlim=c(0,1),
    legend.text=TRUE,space=c(0,1.2),args.legend=list(x=1,y=17)))

Note that beside=FALSE is not a good option if our data have both positive & negative values as this can produce a weird result. For instance:

Y

##         [,1]        [,2]        [,3]
## A  0.7618018  2.59661203 -0.23590050
## B  0.8548063  1.40920654  0.73019040
## C -1.1970456  1.09259205  0.56974523
## D -0.4355176  2.31904920 -1.18347277
## E  2.2923226  2.35986558 -3.22201293
## F  2.5004466  3.98754763 -2.68168275
## G  2.7103706  1.94742983 -1.81626643
## H  1.4677625  1.90001387 -2.78720751
## I  1.5596683  4.57540030 -0.17596327
## J  1.5259543  5.22197986 -0.50073032
## K  1.4834068  4.99761663 -0.35070911
## L  1.3819468  4.84211948 -0.49581275
## M  1.4874512  4.99309837 -0.06583491
## N  0.3560632  3.50340705 -0.95027035
## O -0.2002478  3.80948866 -0.60808238
## P -0.9702761  3.93130784 -3.55745209
## Q -0.2638056  2.81862255 -2.47669725
## R  0.6922559  3.11378711 -3.53855007
## S -0.5630365  3.80538588 -2.87224419
## T  0.7185199  3.63159681 -0.96902150
## U  0.2464106  3.45631197 -0.05955653
## V -0.9173568  5.08132488 -3.08922774
## W -0.5006388  4.79906054 -3.26551534
## X  0.6836688  4.47601771 -3.28019095
## Y -0.5106440 -1.65927057 -0.89645983
## Z  0.3405579 -0.02014156 -0.76899774

plotTree.barplot(tree,Y)

If one is determined to use a barplot for these kind of data, I would recommend beside=TRUE which produces a more sensible result:

plotTree.barplot(tree,Y,args.barplot=list(beside=TRUE,space=c(0,1.2)))

That's it.

Overlaying a contMap style continuous character map on a dotTree plot

A phytools user recently posted the following question:

“Is there a blog post / function capable of combining the functionality of contmap and dottree? I'm looking to plot a continuous character with a color ramp on the tree with a discrete character dot next to the tips of the tree. If I there is a post on this I have missed it in my search and would appreciate being pointed in the right direction if it is out there.”

The answer is that this is pretty straightforward to do - though it does take a little investigating of the internals of each function.

Here my approach will be to compute a "contMap" object & then overlay it on top of a plotted dotTree.

We can start by imagining the following data:

library(phytools)
tree

## 
## Phylogenetic tree with 26 tips and 25 internal nodes.
## 
## Tip labels:
##  A, B, C, D, E, F, ...
## 
## Rooted; includes branch lengths.

x

## A B C D E F G H I J K L M N O P Q R S T U V W X Y Z 
## c c c b a c c c c c c b a a b b b b b c c c b a a a 
## Levels: a b c

y

##       A       B       C       D       E       F       G       H       I 
## -0.7592 -0.9114 -1.1249  0.0006 -0.8630  2.5855  2.8233  2.4789  2.3854 
##       J       K       L       M       N       O       P       Q       R 
## -2.5600 -2.3400 -2.3566 -1.8551 -1.6360 -0.4659 -3.1787 -2.8704 -0.0257 
##       S       T       U       V       W       X       Y       Z 
##  0.3514  0.0436  0.9419 -0.2918 -0.6332 -0.5828 -0.3969 -2.1122

Next, we compute our "contMap" object, but we don't plot it:

obj<-contMap(tree,y,plot=FALSE)
obj

## Object of class "contMap" containing:
## 
## (1) A phylogenetic tree with 26 tips and 25 internal nodes.
## 
## (2) A mapped continuous trait on the range (-3.1787, 2.8233).

Finally, we plot our discrete character using dotTree after which we overlay our contMap style plot. This is the tricky part because I had to look up & duplicate the internal plotting parameters of dotTree to ensure that the two trees where directly overlain:

cols<-setNames(c("black","red","blue"),c("a","b","c"))
dotTree(tree,x,colors=cols,lwd=7)
par(fg="transparent")
plot(obj$tree,add=TRUE,lwd=5,colors=obj$cols,
    ylim=c(-1/25*Ntip(tree),Ntip(tree)),offset=1.7)
par(fg="black")
add.color.bar(0.3*max(nodeHeights(tree)),obj$cols,title="trait value",
    lims=obj$lims,digits=3,prompt=FALSE,x=0.3*max(nodeHeights(tree)),
    y=0.4*(1+par()$usr[3]),lwd=4,fsize=1,subtitle="")

That's it!

Plotting terminal edges of the tree different colors depending (or not) on a discrete character using phytools

I recently received the following request:

“I would love to know if there is any way of extracting and colouring only a certain set of the terminal edges of the phylogeny, based on the tip labels. For example, using edge.color within plot.phylo based on a vector of tips. Do you know of any code that could be used to do this.”

This is pretty straightforward to do using phytools. Note that it can also be done using plot.phylo in the ape package, but I'm going to focus on the phytools way for obvious reasons.

I'll image that I'm painting the terminal edges leading to each tip with different colors based on the state of a discrete trait - but we could modify our technique for any arbitrary grouping of tips.

Here's the tree & data:

library(phytools)
tree

## 
## Phylogenetic tree with 26 tips and 25 internal nodes.
## 
## Tip labels:
##  A, B, C, D, E, F, ...
## 
## Rooted; includes branch lengths.

x

## A B C D E F G H I J K L M N O P Q R S T U V W X Y Z 
## c a c a a c c a a a b c c a c b a a b c c b c a a b 
## Levels: a b c

Now, let's say we want to color the edges blue & red respectively if the corresponding edge is in state b or c, respectively, and leave it black if the tip is in state a.

## identify the tips in states b or c
b<-names(x)[x=="b"]
b

## [1] "K" "P" "S" "V" "Z"

c<-names(x)[x=="c"]
c

##  [1] "A" "C" "F" "G" "L" "M" "O" "T" "U" "W"

## paint the edges
tt<-paintBranches(tree,edge=sapply(b,match,tree$tip.label),
    state="b",anc.state="a")
tt<-paintBranches(tt,edge=sapply(c,match,tree$tip.label),
    state="c")

Now we can plot our tree:

tt

## 
## Phylogenetic tree with 26 tips and 25 internal nodes.
## 
## Tip labels:
##  A, B, C, D, E, F, ...
## 
## The tree includes a mapped, 3-state discrete character with states:
##  a, b, c
## 
## Rooted; includes branch lengths.

cols<-setNames(c("black","blue","red"),c("a","b","c"))
plot(tt,colors=cols,lwd=4,split.vertical=TRUE,ftype="i")

Let's check this against a plot generated using dotTree:

dotTree(tree,x,colors=cols)

That's about it.

Tree & data for this demo were simulated as follows:

tree<-pbtree(n=26,tip.label=LETTERS)
Q<-matrix(c(-2,1,1,1,-2,1,1,1,-2),3,3)
rownames(Q)<-colnames(Q)<-letters[1:3]
x<-as.factor(sim.history(tree,Q)$states)