Sunday, May 7, 2023

Arc-style phylogeny plotting now in phytools

Recently, I described a clever trick to plot a phylogenetic tree in a semi-circular “arc” format. This type of plot is popular in the literature, but I wasn’t aware of an R phylogenetics plotting method that could make plots of this style.

I’m still working out the bugs, but I have now added this feature to various phytools plotting methods. Where available, we can access it by setting the argument type to type="arc".

Let’s see.

library(phytools)
packageVersion("phytools")
## [1] '1.8.6'

First, let’s plot a phylogeny using plotTree. This phylogeny of Mycalesina butterflies comes from Halali et al. (2020).

data(butterfly.tree)
butterfly.tree
## 
## Phylogenetic tree with 287 tips and 286 internal nodes.
## 
## Tip labels:
##   Bra_peitho_gigas, Bra_decira, Bra_simonsii, Bra_perspicua, Bra_phaea, Tel_adolphei, ...
## 
## Rooted; includes branch lengths.
plotTree(ladderize(butterfly.tree),type="arc",fsize=0.3,
  lwd=1,ftype="i",arc_height=1.5)

plot of chunk unnamed-chunk-8

That’s not bad. The argument arc_height is used to set the height of the arc compared to the total depth of the tree. arc_height=1.5 (in this case), makes the bottom of the arc 1.5 × the total tree depth in height above the origin.

Now, let’s try a stochastic character mapped tree. For this, I’ll use the anoletree object which is adapted from Mahler et al. (2010). It already includes a mapped character: ecomorph.

data(anoletree)
anoletree
## 
## Phylogenetic tree with 82 tips and 81 internal nodes.
## 
## Tip labels:
## 	ahli, allogus, rubribarbus, imias, sagrei, bremeri, ...
## 
## The tree includes a mapped, 6-state discrete character
## with states:
## 	CG, GB, TC, TG, Tr, Tw
## 
## Rooted; includes branch lengths.
cols<-setNames(c("green","#E4D96F","darkgreen",
    "brown","black","darkgrey"),
    c("CG","GB","TC","TG","Tr","Tw"))
cols
##          CG          GB          TC          TG          Tr          Tw 
##     "green"   "#E4D96F" "darkgreen"     "brown"     "black"  "darkgrey"
plot(anoletree,cols,type="arc",ftype="i",fsize=0.8,
  arc_height=1)
legend(x=0,y=1,names(cols),col=cols,pch=15,xjust=0.5,
  horiz=TRUE,cex=0.8,bty="n",pt.cex=1.5)

plot of chunk unnamed-chunk-11

Very nice.

Now, how about a continuous character map? In this case, I’ll graph the first phylogenetic PC from an analysis of various morphological traits of Anolis – again using data from Mahler et al. (2010).

data(anole.data)
anole_ppca<-phyl.pca(anoletree,anole.data)
anole_ppca
## Phylogenetic pca
## Standard deviations:
##        PC1        PC2        PC3        PC4        PC5        PC6 
## 0.33222579 0.09207406 0.05012082 0.04318123 0.02008655 0.01507125 
## Loads:
##            PC1         PC2         PC3         PC4         PC5         PC6
## SVL -0.9712234  0.16067288  0.01972570  0.14782215 -0.06211906  0.06935433
## HL  -0.9645111  0.16955087 -0.01203113  0.17994634  0.08064324 -0.04406887
## HLL -0.9814164 -0.02674808  0.10315533 -0.13790763  0.06887922  0.04126248
## FLL -0.9712265  0.17585694  0.10697935 -0.09105747 -0.06075142 -0.04864769
## LAM -0.7810052  0.37429334 -0.47398703 -0.15871456  0.00217418  0.00875408
## TL  -0.9014509 -0.42528918 -0.07614571  0.01709649 -0.01750404 -0.01088743

PC 1 is negative size, so let’s flip it.

pc1<--scores(anole_ppca)[,1]
pc1_cMap<-contMap(anoletree,pc1,plot=FALSE)
pc1_cMap<-setMap(pc1_cMap,viridisLite::viridis(n=10,direction=-1))

Since the biggest anoles should be the crown-giant ecomorphs, let’s add arrows that point to each of them in our tree.

ecomorph<-getStates(anoletree,"tips")
cg<-names(which(ecomorph=="CG"))
cg
##  [1] "garmani"      "baleatus"     "barahonae"    "ricordii"     "cuvieri"      "baracoae"     "noblei"      
##  [8] "smallwoodi"   "luteogularis" "equestris"
plot(pc1_cMap,type="arc",ftype="i",fsize=0.8,arc_height=1.5,
  legend=10,outline=FALSE,lwd=3,leg.txt="PC 1 (overall size)")
tips<-sapply(cg,function(x,y) which(y==x),y=anoletree$tip.label)
add.arrow(anoletree,tips,arrl=1,offset=1,lwd=2,
  col=palette()[4])

plot of chunk unnamed-chunk-16

Lastly, how about this?

plotTree(as.phylo(anoletree),offset=2,type="arc",lwd=1,
  fsize=0.8,ftype="i",arc_height=2.2)
tiplabels(pie=to.matrix(ecomorph,sort(unique(ecomorph))),
  piecol=cols,cex=0.2)
text(0,0.1,"fitted \"ARD\" model",adj=0.5)
ard_ecomorph<-fitMk(anoletree,ecomorph,model="ARD",
  opt.method="optimParallel",lik.func="pruning",
  pi="fitzjohn")
xy<-plot(as.Qmatrix(ard_ecomorph),show.zeros=FALSE,
  add=TRUE,ylim=c(-1.2,3),xlim=c(-2,2),spacer=0.2,
  text=FALSE,width=TRUE,max.lwd=5)
library(plotrix)
nulo<-mapply(draw.circle,xy$x,xy$y,col=cols,
    MoreArgs=list(radius=0.15))
text(xy$x,xy$y,xy$states,cex=0.8,col="white")

plot of chunk unnamed-chunk-17

That’s pretty cool, am I right?

No comments:

Post a Comment

Note: due to the very large amount of spam, all comments are now automatically submitted for moderation.