Continuous character methods

In this tutorial I'm going to quickly overview a range of plotting methods for phylogenies & comparative data that are implemented in the phytools package.

### Continuous character methods

The first plotting method I'll illustrate is called a 'traitgram' which is a projection of the tree into a space defined by trait values & time since the root. So, for example:

require(phytools)
## simulate a tree
tree<-pbtree(n=26,tip.label=LETTERS)
plotTree(tree)

## simulate data under Brownian motion
x<-fastBM(tree)
x
##          A          B          C          D          E          F
## -2.5058532 -2.5353028 -1.5198028 -3.3444568 -2.6098775 -2.6830252
##          G          H          I          J          K          L
## -2.7673502 -0.5024255 -0.4028755 -2.1903261 -1.3840878 -0.7009695
##          M          N          O          P          Q          R
## -1.6265088 -1.8593068 -2.7533331 -3.1032355 -3.1114163 -2.2118431
##          S          T          U          V          W          X
## -3.7306722 -2.7624590 -4.5961496 -2.5471489 -2.4448663 -2.7218879
##          Y          Z
##  0.7788399  1.6670663
## plot a simple traitgram

We can also plot a traitgram with the uncertainty about ancestral character states visualized using transparent probability density. This is implemented in the phytools function fancyTree. So, for instance:

## plot traitgram with 95% CI
## Computing density traitgram...

Next, we can explore the continuous character mapping function in phytools called contMap. Here is a quick demo using the same data:

## plot contMap
obj<-contMap(tree,x)

The function contMap returns an object of class "contMap" which we can then more easily replot using, for instance, different parameters:

## plot leftward
plot(obj,direction="leftwards")

## plot fan style tree
plot(obj,type="fan",lwd=5) ## for example

Here are some real data using anole.tre and svl.csv:

anole.tree
##
## Phylogenetic tree with 100 tips and 99 internal nodes.
##
## Tip labels:
##  ahli, allogus, rubribarbus, imias, sagrei, bremeri, ...
##
## Rooted; includes branch lengths.
svl<-as.matrix(svl)[,1]
svl
##            ahli         alayoni         alfaroi        aliniger
##        4.039125        3.815705        3.526655        4.036557
##        allisoni         allogus   altitudinalis         alumina
##        4.375390        4.040138        3.842994        3.588941
##       alutaceus     angusticeps     argenteolus     argillaceus
##        3.554891        3.788595        3.971307        3.757869
##         armouri   bahorucoensis        baleatus        baracoae
##        4.121684        3.827445        5.053056        5.042780
##       barahonae        barbatus        barbouri        bartschi
##        5.076958        5.003946        3.663932        4.280547
##         bremeri        breslini    brevirostris        caudalis
##        4.113371        4.051111        3.874155        3.911743
##       centralis  chamaeleonides    chlorocyanus     christophei
##        3.697941        5.042349        4.275448        3.884652
##       clivicola     coelestinus        confusus           cooki
##        3.758726        4.297965        3.938442        4.091535
##    cristatellus    cupeyalensis         cuvieri    cyanopleurus
##        4.189820        3.462014        4.875012        3.630161
##         cybotes     darlingtoni       distichus dolichocephalus
##        4.210982        4.302036        3.928796        3.908550
##       equestris      etheridgei   eugenegrahami       evermanni
##        5.113994        3.657991        4.128504        4.165605
##         fowleri         garmani         grahami           guafe
##        4.288780        4.769473        4.154274        3.877457
##       guamuhaya         guazuma       gundlachi       haetianus
##        5.036953        3.763884        4.188105        4.316542
##      hendersoni      homolechis           imias    inexpectatus
##        3.859835        4.032806        4.099687        3.537439
##       insolitus        isolepis           jubar           krugi
##        3.800471        3.657088        3.952605        3.886500
##      lineatopus   longitibialis        loysiana          lucius
##        4.128612        4.242103        3.701240        4.198915
##    luteogularis      macilentus        marcanoi          marron
##        5.101085        3.715765        4.079485        3.831810
##         mestrei       monticola          noblei        occultus
##        3.987147        3.770613        5.083473        3.663049
##         olssoni        opalinus      ophiolepis        oporinus
##        3.793899        3.838376        3.637962        3.845670
##        paternus        placidus       poncensis        porcatus
##        3.802961        3.773967        3.820378        4.258991
##        5.038034        3.799022        3.466860        3.901619
##      reconditus        ricordii     rubribarbus          sagrei
##        4.482607        5.013963        4.078469        4.067162
##    semilineatus        sheplani         shrevei      singularis
##        3.696631        3.682924        3.983003        4.057997
##      smallwoodi         strahmi       stratulus     valencienni
##        5.035096        4.274271        3.869881        4.321524
##       vanidicus    vermiculatus        websteri       whitemani
##        3.626206        4.802849        3.916546        4.097479
obj<-contMap(anole.tree,svl,fsize=c(0.6,1),outline=FALSE)

We can also plot this using a circular or fan tree:

plot(obj,type="fan",outline=FALSE)

Note that the intense aliasing is just a byproduct of plotting to .png format. To get a high quality image we can plot to .pdf or other formats. For instance:

pdf(file="anole.contMap.pdf",width=10,height=10)
plot(obj,type="fan",outline=FALSE)
dev.off()
## png
##   2

The result can be viewed here.

Finally, a relatively simple new plotting method in phytools is the function plotTree.wBars. That function pretty much does what it sounds like it does:

plotTree.wBars(anole.tree,exp(svl),type="fan",scale=0.002)

It is not too difficult to combine this with a contMap plot. For example:

obj<-contMap(anole.tree,exp(svl),plot=FALSE)
plotTree.wBars(obj\$tree,exp(svl),method="plotSimmap",
tip.labels=TRUE,fsize=0.7,colors=obj\$cols,type="fan",scale=0.002)
x=0.9*par()\$usr[1],y=0.9*par()\$usr[3])

pdf(file="anole.plotTree.wBars.pdf",width=10,height=10)
obj<-contMap(anole.tree,exp(svl),plot=FALSE)
plotTree.wBars(obj\$tree,exp(svl),method="plotSimmap",
tip.labels=TRUE,fsize=0.7,colors=obj\$cols,type="fan",scale=0.002)
x=0.9*par()\$usr[1],y=0.9*par()\$usr[3])
dev.off()
## png
##   2

(A non-aliased version here.)

We can also do two dimensional visualizations of phylogenies in morphospace. The is called a 'phylomorphospace'. E.g.:

X<-fastBM(tree,nsim=3) ## simulate 3 characters
phylomorphospace(tree,X[,c(1,2)],xlab="trait 1",ylab="trait 2")

With phytools we can do static or animated version of the same. I'll just do a static version here, but animated is the default:

phylomorphospace3d(tree,X,method="static")

Try the animated version on your own.

Finally, it's possible to combine the continuous character mapping and 2D phylomorphospaces using a type of phylogenetic scatterplot. Here's a demo:

fancyTree(tree,type="scattergram",X=X)
## Computing multidimensional phylogenetic scatterplot matrix...

### Discrete character methods

Stochastic character maps are easy to plot using phytools. For instance:

data(anoletree)
plotSimmap(anoletree,type="fan")
## no colors provided. using the following legend:
##        CG        GB        TC        TG        Tr        Tw
##   "black"     "red"  "green3"    "blue"    "cyan" "magenta"
## setEnv=TRUE for this type is experimental. please be patient with bugs
## we can use any color scheme, but this is the default
cols<-setNames(palette()[1:length(unique(getStates(anoletree,
"tips")))],
sort(unique(getStates(anoletree,"tips"))))
y=0.9*par()\$usr[4],prompt=FALSE)

This is just one stochastically mapped tree. Let's perform the method of stochastic mapping & plot the posterior probabilities at internal nodes.

x<-getStates(anoletree,"tips")
trees<-make.simmap(anoletree,x,model="ER",nsim=100)
## make.simmap is sampling character histories conditioned on the transition matrix
## Q =
##             CG          GB          TC          TG          Tr          Tw
## CG -0.11570723  0.02314145  0.02314145  0.02314145  0.02314145  0.02314145
## GB  0.02314145 -0.11570723  0.02314145  0.02314145  0.02314145  0.02314145
## TC  0.02314145  0.02314145 -0.11570723  0.02314145  0.02314145  0.02314145
## TG  0.02314145  0.02314145  0.02314145 -0.11570723  0.02314145  0.02314145
## Tr  0.02314145  0.02314145  0.02314145  0.02314145 -0.11570723  0.02314145
## Tw  0.02314145  0.02314145  0.02314145  0.02314145  0.02314145 -0.11570723
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##        CG        GB        TC        TG        Tr        Tw
## 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667
## Done.
obj<-describe.simmap(trees,plot=FALSE)
obj
## 100 trees with a mapped discrete character with states:
##  CG, GB, TC, TG, Tr, Tw
##
## trees have 24.11 changes between states on average
##
## changes are of the following types:
##      CG,GB CG,TC CG,TG CG,Tr CG,Tw GB,CG GB,TC GB,TG GB,Tr GB,Tw TC,CG
## x->y  0.24  0.26  0.22  0.13  0.32  0.58  0.97  1.33  0.35  1.58  1.17
##      TC,GB TC,TG TC,Tr TC,Tw TG,CG TG,GB TG,TC TG,Tr TG,Tw Tr,CG Tr,GB
## x->y  0.55  0.32  0.69     1  1.17  3.58  2.02     1  2.07  0.16  0.22
##      Tr,TC Tr,TG Tr,Tw Tw,CG Tw,GB Tw,TC Tw,TG Tw,Tr
## x->y  0.18  0.15   0.3  0.55  0.83  1.24   0.5  0.43
##
## mean total time spent in each state is:
##              CG         GB         TC         TG          Tr         Tw
## raw  13.2553592 47.1436754 32.3207568 71.1973321 13.12972055 28.6201165
## prop  0.0644506  0.2292234  0.1571509  0.3461778  0.06383972  0.1391576
##        total
## raw  205.667
## prop   1.000
plot(obj,type="fan")
## setEnv=TRUE for this type is experimental. please be patient with bugs
y=0.9*par()\$usr[4],prompt=FALSE)

We can also overlay the posterior probabilities on a random stochastic map in a similar way:

plotSimmap(trees[[1]],type="fan")
## no colors provided. using the following legend:
##        CG        GB        TC        TG        Tr        Tw
##   "black"     "red"  "green3"    "blue"    "cyan" "magenta"
## setEnv=TRUE for this type is experimental. please be patient with bugs
nodelabels(pie=obj\$ace,piecol=cols,cex=0.5)
y=0.9*par()\$usr[4],prompt=FALSE)

For binary characters, we can also map the posterior density of stochastic maps onto the nodes & branches of a tree using the phytools function densityMap. Here's a demo of this using simulated data.

## first let's simulate data
Q<-matrix(c(-1,1,1,-1),2,2)
rownames(Q)<-colnames(Q)<-c(0,1)
tree<-sim.history(tree,Q)
## Done simulation(s).
x<-tree\$states
x
##   A   B   C   D   E   F   G   H   I   J   K   L   M   N   O   P   Q   R
## "1" "0" "0" "1" "1" "0" "0" "1" "0" "1" "0" "0" "0" "1" "1" "1" "1" "0"
##   S   T   U   V   W   X   Y   Z
## "0" "0" "0" "1" "1" "1" "0" "1"
## stochastic maps
trees<-make.simmap(tree,x,nsim=100)
## make.simmap is sampling character histories conditioned on the transition matrix
## Q =
##           0         1
## 0 -1.223814  1.223814
## 1  1.223814 -1.223814
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##   0   1
## 0.5 0.5
## Done.
## make densityMap
obj<-densityMap(trees,lwd=4,outline=TRUE)
## sorry - this might take a while; please be patient

### Combining discrete & continuous characters

We can, for instance, overlay a discrete character history (or stochastic map) onto a traitgram or phylomorphospace.

## simulate data with a high rate on some branches
x<-sim.rates(tree,setNames(c(1,10),c(0,1)))
## plot traitgram
phenogram(tree,x,colors=setNames(c("blue","red"),

## we can do the same with phylomorphospaces
X<-cbind(sim.rates(tree,setNames(c(1,10),c(0,1))),
sim.rates(tree,setNames(c(1,10),c(0,1))))
phylomorphospace(tree,X,colors=setNames(c("blue","red"),c(0,1)))

## it is even possible to overlay a posterior density from stochastic mapping
phylomorphospace(obj\$tree,X,colors=obj\$cols,
lwd=3,xlab="x",ylab="y")
prompt=FALSE,x=0.9*par()\$usr[1],y=0.9*par()\$usr[3])

### Other stuff

How about projecting a tree onto a geographic map? Here's an example using phytools function phylo.to.map:

# simulate tree & data
tree<-pbtree(n=26,scale=100)
tree\$tip.label<-LETTERS[26:1]
lat<-fastBM(tree,sig2=10,bounds=c(-90,90))
long<-fastBM(tree,sig2=80,bounds=c(-180,180))
# now plot projection
xx<-phylo.to.map(tree,cbind(lat,long),plot=FALSE)
## objective: 190
## objective: 190
## objective: 190
## objective: 188
## objective: 188
## objective: 188
## objective: 188
## objective: 186
## objective: 184
## objective: 184
## objective: 184
## objective: 130
## objective: 130
## objective: 130
## objective: 128
## objective: 128
## objective: 128
## objective: 104
## objective: 98
## objective: 98
## objective: 94
## objective: 92
## objective: 90
## objective: 90
## objective: 88
## tree points to map coordinates
plot(xx,type="phylogram",asp=1.3,mar=c(0.1,0.5,3.1,0.1))

## tree projected directly on map
plot(xx,type="direct",asp=1.3,mar=c(0.1,0.5,3.1,0.1))

That's it.

Written by Liam J. Revell. Last updated July 5, 2015