## Exercise 15: Plotting phylogenies & comparative data

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.

### 1. 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:

``````## first load phytools
library(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
## -0.009685797  0.509191927  0.619909980 -0.382375084 -0.759424172
##            F            G            H            I            J
## -1.017865681 -1.423530384 -1.390218699  0.071477307  0.016275700
##            K            L            M            N            O
##  0.093906234 -0.333070445 -1.584253105  0.197771724  0.809517227
##            P            Q            R            S            T
##  0.948086563 -0.238809121  0.205425189 -0.165778994 -0.225704617
##            U            V            W            X            Y
## -0.084803920  1.624216958  1.905642115  1.253396087  1.421246681
##            Z
##  1.527713040
``````
``````## plot a simple traitgram
``````

``````## change the spread costs
``````

We can visualize the states just at the tips of the tree using circles:

``````dotTree(tree,x)
``````

``````## or
dotTree(tree,x,standardize=TRUE)
``````

If we have multiple characters we have other options still:

``````X<-fastBM(tree,nsim=10)
dotTree(tree,X)
``````

``````## or
dotTree(tree,X,standardize=TRUE)
``````

``````## or
phylo.heatmap(tree,X)
``````

``````## or
phylo.heatmap(tree,X,standardize=TRUE)
``````

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:

``````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 (-1.584253, 1.905642).
``````
``````## plot leftward
plot(obj,direction="leftwards")
``````

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

We can also change the color map, to invert it, or to set it to a gray scale:

``````obj<-setMap(obj,invert=TRUE)
plot(obj)
``````

``````obj<-setMap(obj,colors=c("white","black"))
plot(obj)
``````

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

``````anole.tree<-read.tree("Anolis.tre")
anole.tree
``````
``````##
## Phylogenetic tree with 100 tips and 99 internal nodes.
##
## Tip labels:
##  ahli, allogus, rubribarbus, imias, sagrei, bremeri, ...
##
## Rooted; includes branch lengths.
``````
``````svl<-read.csv("svl.csv",header=TRUE,row.names=1)
svl<-setNames(svl[,1],rownames(svl))
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.5,1),outline=FALSE)
``````

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

``````plot(obj,type="fan",outline=FALSE,fsize=c(0.7,1),legend=0.7)
``````

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, another relatively simple 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)
``````

``````## or with tip labels:
plotTree.wBars(anole.tree,exp(svl),type="fan",scale=0.002,tip.labels=TRUE,
fsize=0.7)
``````

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])
``````

Or, conversely, to color just the bars!

``````## first, let's invert the color map
obj<-setMap(obj,invert=TRUE)
ii<-sapply(1:Ntip(obj\$tree), function(x,tree) which(tree\$edge[,2]==x),
tree=obj\$tree)
cols<-setNames(obj\$cols[sapply(ii,function(x,maps)
names(maps[[x]])[length(maps[[x]])],
maps=obj\$tree\$maps)],obj\$tree\$tip.label)
## I'm going to change the scale too
plotTree.wBars(obj\$tree,exp(svl),method="plotTree",
tip.labels=TRUE,fsize=0.7,col=cols,type="fan",scale=0.004)
x=0.9*par()\$usr[1],y=0.9*par()\$usr[3])
``````

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 computer:

``````phylomorphospace3d(tree,X)
``````

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...
``````

### 2. Discrete character methods

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

``````data(anoletree)
plot(anoletree,type="fan",fsize=0.8,lwd=3)
``````
``````## no colors provided. using the following legend:
##        CG        GB        TC        TG        Tr        Tw
##   "black"     "red"  "green3"    "blue"    "cyan" "magenta"
``````
``````## 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,fsize=0.9)
``````

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<-summary(trees,plot=FALSE)
obj
``````
``````## 100 trees with a mapped discrete character with states:
##  CG, GB, TC, TG, Tr, Tw
##
## trees have 23.99 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.31  0.23  0.26  0.13  0.36  0.63  0.65  1.32   0.4  1.41  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.41  0.29  0.61  0.85  0.98  3.28  1.96  1.08  1.76  0.18  0.17
##      Tr,TC Tr,TG Tr,Tw Tw,CG Tw,GB Tw,TC Tw,TG Tw,Tr
## x->y  0.29  0.12  0.26  0.61  1.37  1.66  0.72  0.52
##
## mean total time spent in each state is:
##               CG         GB         TC         TG          Tr         Tw
## raw  13.46162668 46.0141794 31.1704542 69.9454383 13.03485493 32.0404070
## prop  0.06545352  0.2237315  0.1515579  0.3400908  0.06337846  0.1557878
##        total
## raw  205.667
## prop   1.000
``````
``````plot(obj,type="fan",fsize=0.8,cex=c(0.5,0.3))
y=0.9*par()\$usr[4],prompt=FALSE,fsize=0.9)
``````

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

``````plot(trees[[1]],type="fan",fsize=0.8)
``````
``````## no colors provided. using the following legend:
##        CG        GB        TC        TG        Tr        Tw
##   "black"     "red"  "green3"    "blue"    "cyan" "magenta"
``````
``````nodelabels(pie=obj\$ace,piecol=cols,cex=0.5)
y=0.9*par()\$usr[4],prompt=FALSE,fsize=0.9)
``````

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 using the elopomorph data from the ancestral state reconstruction part.

In case you no longer have them, these data are here:

``````tree<-read.tree("elopomorph.tre")
tree
``````
``````##
## Phylogenetic tree with 61 tips and 60 internal nodes.
##
## Tip labels:
##  Moringua_edwardsi, Kaupichthys_nuchalis, Gorgasia_taiwanensis, Heteroconger_hassi, Venefica_proboscidea, Anguilla_rostrata, ...
##
## Rooted; includes branch lengths.
``````
``````X<-read.csv("elopomorph.csv",row.names=1)
fmode<-setNames(X[,1],rownames(X))
## stochastic maps
trees<-make.simmap(tree,fmode,nsim=100)
``````
``````## make.simmap is sampling character histories conditioned on the transition matrix
##
## Q =
##                bite     suction
## bite    -0.01582783  0.01582783
## suction  0.01582783 -0.01582783
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##    bite suction
##     0.5     0.5
``````
``````## Done.
``````
``````## make densityMap
obj<-densityMap(trees,plot=FALSE)
``````
``````## sorry - this might take a while; please be patient
``````
``````plot(obj,lwd=4,outline=TRUE,states=c("suction","bite"),
fsize=c(0.7,0.9))
``````

### 3. Combining discrete & continuous characters

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

``````## plot traitgram   with single stochastic map
bsize<-setNames(X[,2],rownames(X))
phenogram(trees[[1]],bsize,colors=
setNames(c("blue","red"),c("suction","bite")),
ftype="i")
``````
``````## Optimizing the positions of the tip labels...
``````

``````## or with our posterior density
ylab="body size")
``````
``````## Optimizing the positions of the tip labels...
``````
``````add.color.bar(50,obj\$cols,title="PP(state=\"bite\")",subtitle="",
prompt=FALSE,x=0,y=240,fsize=0.8)
``````

``````## or with our posterior probabilities
obj<-summary(trees)
phenogram(tree,bsize,fsize=0.7,ftype="i",ylab="body size")
``````
``````## Optimizing the positions of the tip labels...
``````
``````nodelabels(pie=obj\$ace,piecol=setNames(c("red","blue"),colnames(obj\$ace)),
cex=0.6)
tiplabels(pie=to.matrix(fmode[tree\$tip.label],colnames(obj\$ace)),
piecol=setNames(c("red","blue"),colnames(obj\$ace)),cex=0.4)
``````

It's also possible to do the same kind of thing with `phylomorphospace`.

### 4. Other methods

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

``````# simulate tree & data
set.seed(1)
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: 126
## objective: 126
## objective: 126
``````
``````## objective: 106
## objective: 106
``````
``````## objective: 100
## objective: 100
## objective: 100
## objective: 100
## objective: 100
## objective: 100
## objective: 100
``````
``````## objective: 98
## objective: 98
``````
``````## objective: 96
## objective: 96
``````
``````## objective: 94
## objective: 94
## objective: 94
``````
``````## objective: 84
``````
``````## objective: 76
## objective: 76
``````
``````## objective: 72
## objective: 72
``````
``````## objective: 70
``````
``````## 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))
``````

Co-phylogenetic plotting? Let's try it:

``````t1<-rtree(n=26,tip.label=LETTERS)
t2<-rtree(n=26,tip.label=LETTERS)
obj<-cophylo(t1,t2)
``````
``````## Rotating nodes to optimize matching...
## Done.
``````
``````plot(obj)
``````

``````## or non-matching:
t2<-rtree(n=40)
assoc<-cbind(t1\$tip.label,sample(t2\$tip.label,26,replace=TRUE))
assoc
``````
``````##       [,1] [,2]
##  [1,] "C"  "t29"
##  [2,] "W"  "t8"
##  [3,] "J"  "t12"
##  [4,] "K"  "t36"
##  [5,] "P"  "t3"
##  [6,] "F"  "t29"
##  [7,] "U"  "t39"
##  [8,] "N"  "t20"
##  [9,] "S"  "t26"
## [10,] "A"  "t14"
## [11,] "I"  "t10"
## [12,] "X"  "t9"
## [13,] "R"  "t28"
## [14,] "V"  "t34"
## [15,] "D"  "t14"
## [16,] "E"  "t19"
## [17,] "B"  "t13"
## [18,] "O"  "t4"
## [19,] "T"  "t35"
## [20,] "Z"  "t29"
## [21,] "G"  "t38"
## [22,] "Y"  "t29"
## [23,] "M"  "t6"
## [24,] "Q"  "t7"
## [25,] "L"  "t35"
## [26,] "H"  "t2"
``````
``````obj<-cophylo(t1,t2,assoc=assoc)
``````
``````## Rotating nodes to optimize matching...
## Done.
``````
``````plot(obj)
``````