Showing different orders or other taxonomic groups in different colors in a phylomorphospace plot

A few days ago, I received the following email about the phytools function phylomorphospace:

“I'm trying to make a phylomorphospace that shows me the Orders form the birds in different colors. I tried to follow the example of the sunfishes, where species are in red or in blue depending on the feeding mode. However, I do not know what is happening with my data. I'm sure that the Order column is a factor (with 21 different levels) and one color has been assigned to each level, but I don't know why the phylomorphospace is only showing one color (black). Is there a problem with the color argument in the phylomorphospace function?”

The quick answer to this is that if we want to use colors to show different groups or clades on our tree, we should first use paintSubTree to map these groups onto the phylogeny, and then we can use phylomorphospace to graph them.

To demonstrate this, I'll use the classic Garland et al. (1992) that comes loaded with the phytools package.

This consists of a phylogeny of mammals, along with a dataset for body mass and home range size for the different species in the tree. Conveniently, this phylogeny only includes representatives of three different mammalian orders (Carnivora, Perissodactyla, and Artiodactyla) so it should be pretty easy to map these three different orders on the tree.

First let's load our packages, tree, and dataset.

library(phytools)
data(mammal.tree)
data(mammal.data)

Now let's plot the tree as follows:

plotTree(mammal.tree,ftype="i",fsize=0.8)
nodelabels(cex=0.6,frame="circle",bg="white")

Maybe we can see from this plot that node 51 corresponds to the most recent common ancestor of all carnivorans, while the Perissodactyla and Artiodactyla are delineated by nodes 70 and 75 respectively. For a larger phylogenies with more orders we could use the function findMRCA to identify the MRCA based on a simple list of the members of each clade.

nodes<-c(51,70,75)
clades<-c("Carnivora","Perissodactyla","Artidactyla")

Next, let's map each of these three clades onto the tree using paintSubTree. paintSubTree will create an object of a new class: "simmap". In this case, however, we're not using this object type to show the stochastic history of a mapped discrete trait, but the known groups of the terminal taxa in the tree.

mammal.orders<-mammal.tree
for(i in 1:length(nodes)){
    mammal.orders<-paintSubTree(mammal.orders,node=nodes[i],
        state=clades[i],anc.state=if(i==1) "NA" else NULL,
        stem=FALSE)
}

Now we're ready to plot our tree to verify that we have correctly mapped our clades to each node.

plot(mammal.orders,colors<-setNames(c("grey",
    palette()[2:4]),c("NA",clades)),ftype="i",fsize=0.8)
nodelabels(cex=0.6,frame="circle",bg="white")

Terrific. This is exactly what we were going for.

Lastly, we plot our phylomorphospace as follows:

phylomorphospace(mammal.orders,log(mammal.data),colors=colors,
    ftype="off",bty="n",node.size=c(0,1),node.by.map=TRUE,
    xlab="log(body mass)",ylab="log(home range size)")
legend("topleft",names(colors)[2:4],pch=21,pt.bg=colors[2:4],
    bty="n",pt.cex=1.5)

Simulating different phylogenetic signal before & after a specified timepoint

A few days ago, an R-SIG-Phylo user asked:

“Is there any way to simulate the evolution of a trait with a given significant signal (say K=1) from the tree root to time x and no signal (K=0) from time X to the present?”

Well, there's actually no straightforward way to simulate a predetermined value of Blomberg et al.'s (2003); however, it is possible to simulate data in which evolution in one part of the tree is by a phylogenetic signal preserving process (e.g., Brownian motion), while evolution in a different part of the tree is by a phylogenetic signal destroying process.

The way I thought of to do this was using the phytools function treeSlice.

treeSlice cuts the tree at a user-specified point (it even has an interactive mode!) and then returns either the tree up to the point of the cut (for orientation="rootwards") or all the subtrees descended from the cut (for orientation="tipwards").

My idea was to simply simulate under Brownian motion on the rootward tree, pass the simulated “tip” values from this step as root states to all the descendant subtrees, and then simulate under a \(\lambda\) model (with \(\lambda\) = 0, in this case) for each subtree.

Let's try it.

## load packages & simulate tree
library(phytools)
tree<-pbtree(n=40,scale=100)
## get hidden lambdaTree function (this is just to simulate random
## values for the subtrees from time t
lt<-phytools:::lambdaTree
## set t (here I'm imagining the tree has a total length of 100 and
## I want to slice it at time = 70 units above the root
t<-70
## slice the tree rootward of time = t -- this generates an object
## of class "multiPhylo"
rootward<-treeSlice(tree,orientation="rootwards",slice=t)
## now slice the tree tipward of time t -- this generates an object
## of class "multiPhylo" with a length equal to the number of tips
## in rootward
tipward<-treeSlice(tree,orientation="tipwards",trivial=TRUE,
    slice=t)
## convert each subtree in tipward to a "phylo" object with an
## unbranching (singleton) root node. This is for simulation.
tipward<-lapply(tipward,rootedge.to.singleton)
## all of this is just plotting and can be ommitted
nn<-sum(sapply(tipward,Ntip)>1)
cc<-ceiling(sqrt(nn))
layout(matrix(c(rep(nn+1,cc),1:nn,rep(nn+2,cc^2-nn)),
    cc+1,cc,byrow=TRUE),heights=c(0.5,rep(1,cc)))
nulo<-lapply(tipward[sapply(tipward,Ntip)>1],plotTree,mar=rep(1.1,4),
    fsize=0.9)
plot(NA,xlim=c(-1,1),ylim=c(-1,1),axes=FALSE)
text(0,0,paste("Non-trivial subtrees after t =",t),cex=2,font=3)
## end plotting
## simulate the ancestral states for each subtree in tipward under
## pure Brownian motion
aa<-fastBM(rootward)
## create a vector of values 
xx<-vector()
## simulate on each subtree to populate it. a is the ancestral state
## for each simulation
for(i in 1:length(tipward))
    xx<-c(xx,fastBM(lt(tipward[[i]],0),a=aa[i]))

We can measure phylogenetic signal:

phylosig(tree,xx,method="lambda")

## 
## Phylogenetic signal lambda : 0.66126 
## logL(lambda) : -142.569

Let's plot our data simulated in this way using the phytools phenogram function:

phenogram(make.era.map(tree,c(0,t)),xx,ftype="off")
legend("topleft","switch in process",lty="dotted",bty="n")
clip(0,100,min(xx),max(xx))
abline(v=t,lty="dotted")

For fun, let's compare this to the opposite scenario – in which no phylogenetic signal is simulated through time = 70, and Brownian motion therafter….

t<-70
rootward<-treeSlice(tree,orientation="rootwards",slice=t)
tipward<-treeSlice(tree,orientation="tipwards",trivial=TRUE,
    slice=t)
tipward<-lapply(tipward,rootedge.to.singleton)
aa<-fastBM(lt(rootward,0))
xx<-vector()
for(i in 1:length(tipward))
    xx<-c(xx,fastBM(tipward[[i]],a=aa[i]))
phylosig(tree,xx,method="lambda")

## 
## Phylogenetic signal lambda : 1.00018 
## logL(lambda) : -123.278

phenogram(make.era.map(tree,c(0,t)),xx,ftype="off")
legend("topleft","switch in process",lty="dotted",bty="n")
clip(0,100,min(xx),max(xx))
abline(v=t,lty="dotted")