Saturday, December 24, 2016

Using phytools to plot different genera with different colors

A phytools user asks the following:

“I am following your blog. I have a phylogenetic tree. I would like to color the tree according to their genus name. Could you please help me, providing the code in R?”

There are in fact many ways to plot a tree with edges or clades in different colors. I'm going to give a demo using paintSubTree and plotSimmap in phytools as follows.

First, let's imagine the following species-level tree in which I have used the syntax Genus_species for all taxon labels.

library(phytools)
species.tree
## 
## Phylogenetic tree with 56 tips and 55 internal nodes.
## 
## Tip labels:
##  Rfom_nhijvq, Nvztxu_evipjo, Bswlbtch_azndgw, Ghsygvn_azkqlr, Qptbc_lwcqya, Qptbc_tpnozk, ...
## 
## Rooted; includes branch lengths.
plotTree(species.tree,ftype="i",fsize=0.7,color="black")

plot of chunk unnamed-chunk-1

First, let's identify all the genera:

genera<-sapply(species.tree$tip.label,function(x) strsplit(x,"_")[[1]][1])
genera<-sort(unique(genera))
genera
##  [1] "Bswlbtch" "Dxfnqt"   "Fyoed"    "Gdzxskm"  "Ghsygvn"  "Nvztxu"  
##  [7] "Pxvolmfp" "Qpebj"    "Qptbc"    "Rfom"     "Sehvikuc" "Sfownt"  
## [13] "Uspvk"    "Uvybs"    "Wpuqdc"   "Wqbhd"

Next, find the MRCA of each genus. In the event that the genus contains only one member, we can just paint the terminal edge:

for(i in 1:length(genera)){
    ii<-grep(genera[i],species.tree$tip.label)
    ca<-if(length(ii)>1) 
        findMRCA(species.tree,species.tree$tip.label[ii]) else ii
    species.tree<-paintSubTree(species.tree,ca,state=as.character(i),
        anc.state="0",stem=TRUE)
}

The following is a trick to remove map segments of zero length:

tol<-max(nodeHeights(species.tree))*1e-12
species.tree$maps<-lapply(species.tree$maps, function(x,tol) 
    if(length(x)>1) x[-which(x<tol)] else x,tol=tol)

Now we can set our colors & plot them:

cols<-setNames(c("grey",rainbow(length(genera))),
    0:length(genera))
plot(species.tree,fsize=0.7,colors=cols,ftype="i",split.vertical=TRUE,
    lwd=3,xlim=c(-24,120))
par(font=3)
add.simmap.legend(colors=setNames(cols[2:length(cols)],genera),fsize=0.8,
    prompt=FALSE,x=-24,y=24)

plot of chunk unnamed-chunk-5

That's basically the idea.

The tree for this example is obviously simulated. The following is the code that was used for simulation (taken from a previous post, here)::

library(phytools)
## first let's simulate our genus tree:
foo<-function() paste(sample(LETTERS,1),paste(sample(letters,
    round(runif(1,min=3,max=7))),collapse=""),sep="")
genera<-replicate(16,foo())
genus.tree<-pbtree(n=length(genera),tip.label=genera,scale=80)
genus.tree$edge.length[which(genus.tree$edge[,2]<=Ntip(genus.tree))]<-
    genus.tree$edge.length[which(genus.tree$edge[,2]<=Ntip(genus.tree))]+20
tips<-c()
for(i in 1:Ntip(genus.tree)){
    n.genus<-sample(1:6,1)
    if(n.genus>0) for(j in 1:n.genus) 
        tips<-c(tips,paste(genus.tree$tip.label[i],
        paste(sample(letters,6),collapse=""),sep="_"))
}
## add them:
species.tree<-genus.to.species.tree(genus.tree,tips)

Sunday, December 18, 2016

Wrap-up & t-shirt design for Bariloche (Argentina) Macroevolution course

I just got back from teaching a five-day workshop the Universidad Nacional del Comahue in San Carlos de Bariloche, Argentina.

It was tiring, but a lot of fun - and I hope it was as rewarding for the students as it was for me.

As has become a tradition in the Latin American workshops that I teach (e.g., 1, 2, 3, 4, etc.), we have a course t-shirt - and, not only that, the t-shirt must be designed in R.

This year's workshop in Bariloche was not spared. Firstly, the design:

library(phytools)
## Loading required package: ape
## Loading required package: maps
library(plotrix)
source("http://www.phytools.org/Bariloche2016/data/arctext.R")
set.seed(545)
tree<-pbtree(n=78,scale=0.8)
txt<-write.tree(tree)
txt<-strsplit(txt,";")[[1]]
txt<-paste("((Y:0.9,",txt,":0.1):0.1,Z:1);",sep="")
tree<-read.tree(text=txt)
tree<-phytools:::lambdaTree(tree,0.95)
tree<-make.era.map(tree,c(0,0.4,0.8))
par(bg="darkgrey")
par(ljoin=1,lmitre=30)
par(mar=rep(0,4))
col<-rgb(red=117/255,green=170/255,blue=219/255)
par(fg="transparent")
plot(tree,colors=setNames(rep("black",3),1:3),
    lwd=10,type="fan",part=0.5,ftype="off",ylim=c(-0.2,1.2),
    lend=1)
par(fg="black")
plot(tree,colors=setNames(c(col,"white",col),1:3),
    lwd=8,type="fan",part=0.5,ftype="off",ylim=c(-0.2,1.2),
    lend=1,add=TRUE)
text(mean(par()$usr[1:2]),-0.15,"Bariloche, Argentina 2016",
    col="black",cex=2.7,font=2)
chars<-strsplit("Latin American Macroevolution Workshop","")[[1]]
stretch<-rep(1.1,length(chars))
arctext("Latin American Macroevolution Workshop",center=c(0,0),radius=1.1,
    cex=3.5,font=2,stretch=stretch)

plot of chunk unnamed-chunk-1

(High resolution PDF with Times font for the text here.)

The only small wrinkle is that I tweaked the code of arctext in the plotrix package to get the text to look the way I wanted it to.

Finally, but most importantly, the students (& me) sporting our t-shirts:

(Click on image for larger version.)

Many thanks to course organizers Lina Moreno Azócar & Félix Cruz, who were wonderful hosts, helped put together a terrific course, and who also had the wonderful patience to drive me back & forth across the border to Chile (no small task) so I didn't have to fly through Buenos Aires!

Until next time Argentina!

Wednesday, December 7, 2016

Coloring tip nodes in phylomorphospace3d

I was just corresponding with a phytools user & essentially the lack of an argument node.col in the function phylomorphospace3d was broached.

Why this option doesn't exist, I couldn't tell you; however it is relatively easy to circumvent because the function invisibly returns a list of functions which can be used to plot points, lines, planes, etc., back in the original (simulated) space.

Here's a quick example, using simulated data:

library(phytools)
tree ## our tree
## 
## Phylogenetic tree with 26 tips and 25 internal nodes.
## 
## Tip labels:
##  A, B, C, D, E, F, ...
## 
## Rooted; includes branch lengths.
X ## our data
##          [,1]         [,2]        [,3]
## A -2.84172235  1.148037248  1.64782403
## B -2.05052328  1.431708869  1.15202298
## C -1.24686168  1.946696167  1.54897262
## D -0.52310503  0.573110584  0.46436125
## E -0.55578283  0.813448906  0.40254794
## F -0.51150759  1.733066759  1.49429840
## G  0.53936198  1.816801753  1.86050506
## H  0.01234429  2.302339855  2.63078614
## I -0.01684737  1.810997083  2.22857705
## J  1.32317051  2.492348110  0.55660832
## K  0.39032083  1.138053408 -0.42215854
## L  1.35224610  0.701773929 -0.73233992
## M  0.06642055  2.393143196  1.67769111
## N  0.53170462  2.387858206  1.07395647
## O  1.33893651  1.869361233  1.18002858
## P  0.16788769  1.722551818  0.03908209
## Q  0.72826260  0.591128118  1.23662094
## R  0.79698877  1.310841641  1.38182879
## S  0.08976056  1.727227579  2.10565084
## T  1.37506375 -0.551729936 -1.73123720
## U -1.09244120  0.960226117 -1.01694632
## V -0.75692309 -0.002878853 -1.02230472
## W -0.53535698 -0.168545803 -1.33846396
## X -0.26869690 -0.241183840 -1.14188545
## Y  0.19315024  0.282576551 -0.89620379
## Z -0.56460363 -0.613852876  0.98093216
y ## our character we want to map on the tips
##   A   B   C   D   E   F   G   H   I   J   K   L   M   N   O   P   Q   R 
## "c" "c" "b" "b" "b" "a" "a" "b" "b" "c" "c" "b" "c" "c" "c" "a" "c" "b" 
##   S   T   U   V   W   X   Y   Z 
## "b" "b" "c" "b" "c" "c" "c" "c"
tip.col<-y
tip.col[tip.col=="a"]<-"red"
tip.col[tip.col=="b"]<-"blue"
tip.col[tip.col=="c"]<-"yellow"
obj<-phylomorphospace3d(tree,X,method="static")
obj$points3d(X,cex=1.4,pch=21,bg=tip.col[rownames(X)])

plot of chunk unnamed-chunk-1

That's it.

The tree & data were simulated as follows:

X<-fastBM(tree,nsim=3)
Q<-matrix(c(-1,1,0,1,-2,1,0,1,-1),3,3)
rownames(Q)<-colnames(Q)<-c("a","b","c")
y<-sim.history(tree,Q)$states