---
title: "Plotting trees with ggfree"
output: github_document
vignette: >
  %\VignetteIndexEntry{ggfree-trees}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_knit$set(global.par = TRUE)
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/README-",
  out.width = "500px",
  dpi = 150,
  fig.align = "center",
  fig.width=4,
  fig.height=4,
  out.width=400
)
```



The [ggtree](https://github.com/YuLab-SMU/ggtree) package by Guangchuang Yu is an excellent resource for drawing information-rich visualizations of phylogenetic trees.
(See https://guangchuangyu.github.io/software/ggtree/gallery/ for a gallery of figures generated using `ggtree`.)
`ggtree` employs the same grammar as `ggplot`, so writing code looks quite similar:
```{r}
# load the package without all the preamble
suppressMessages(require(ggtree))

# generate a random tree
set.seed(1999); phy <- rtree(50)

# draw the tree gg-style!
ggtree(phy) + theme_tree2() + geom_tiplab(size=3)
```

Since I am accustomed to base R, however, I find this syntax unfamiliar and frustrating, and was motivated to start implementing similar plot functions in `ggfree`.
Of course there are several other packages in R that provide plotting functions for phylogenetic tree objects.
For example, the [ape](https://cran.r-project.org/web/packages/ape/index.html) package &mdash; which defines the `phylo` [S3 object](http://adv-r.had.co.nz/S3.html) class that many packages for phylogenetics will accept as inputs &mdash; adapts the generic S3 `plot` function for trees (called `plot.phylo`).
However, it is not straight-forward to annotate the resulting plot with additional information, such as highlighting a subset of internal nodes.
This is because there are two things happening when we call `plot.phylo`: 

1. the method maps the nodes and edges (branches) of the tree to the `x` and `y` coordinates according to one of several layout algorithms, and;

2. the method draws the tree on these coordinates in the active plot device.

Consequently, the `x` and `y` coordinates are not readily available for the user to add additional information to the plot.
(I'll demonstrate why it's nice to have access to these data in subsequent sections.)

```{r echo=FALSE}
par(family="Helvetica")
```

`ggfree` separates these tasks into different functions.
Generating the coordinates for a given layout algorithm is handled by the function `tree.layout`, which returns an [S3 object](http://adv-r.had.co.nz/S3.html) that I call `phyloLayout`.
```{r}
require(ggfree)
L <- tree.layout(phy)
plot(L)
axis(side=1)  # add a horizontal axis like ggtree example
```

There are clearly some differences here.
For instance, the order of tips is different because `ggtree` has automatically ladderized the phylogeny by default.
We can amend this using the `ape:ladderize` function:
```{r}
phy2 <- ladderize(phy, right=F)
plot(tree.layout(phy2), col='black', lwd=1, cex=0.75, mar=c(2,0,0,0))
axis(side=1, cex.axis=0.75, mgp=c(3,0.5,0))
```
Note that I've also added some other basic R graphical parameters to make the plot more similar to the `ggtree` example.
A lot of this code is tucked into the `ggtree` function 


### `phyloLayout` objects

Let's have a look at what's getting returned when we call `tree.layout`:
```{r}
L <- tree.layout(phy2)
L  # calls print.phyloLayout function
```

Here we can see that a `phyloLayout` object contains two data frames for the nodes and edges of the tree, respectively.
(Only the first six rows of each are being displayed because this `print` function is calling `head` for each data frame.)
The `nodes` data frame contains x and y coordinates, and the `edges` data frame contains two sets of x and y coordinates.
In addition, any node- or edge-level attributes that you happened to attach to the original `phylo` object are automatically carried over to these data frames.
These are the main outputs of the layout algorithm.

### Layout algorithms

`ggfree` provides three popular layout algorithms for phylogenetic trees.
We specify which layout algorithm we want to use by passing a `type` argument to the function `tree.layout`:

1. **Rectangular/slanted trees.**  The above trees have been plotted using a rectangular layout.
The horizontal lines (branches) represent the amount of time or expected amount of evolution along each lineage.
By convention, time "flows" forward from the left to right of the figure.
The vertical lines are only used to connect linegaes to their common ancestors and to separate lineages in the vertical dimension - their lengths do not represent any information.

   A slanted tree uses the same coordinate information as a rectangular layout, except that lineages are connected directly to their ancestral nodes without vertical line segments:
   ```{r}
   plot(tree.layout(phy, type='s', unscaled=TRUE), cex=0.5, mar=c(0,0,0,3))
   ```

2. **Circular (radial) trees.**  This is essentially the rectangular layout algorithm in a polar coordinate system.  Time flows outward from the origin (centre) of the plot.
   ```{r}
   plot(tree.layout(phy, type='o'))
   ```

3. **Unrooted trees.**  Like many other R packages that provide plot functions for trees, `ggtree` implements [Joe Felsenstein's](https://en.wikipedia.org/wiki/Joseph_Felsenstein) equal-angle layout algorithm for unrooted trees.
Note that I've turned off tip labels because they become cluttered in this layout:
   ```{r}
   plot(tree.layout(phy, type='u'), label='n', mar=c(0,0,0,0))
   ```


## Annotating nodes



## Example: Influenza tree from `ggtree`

To demonstrate how we can reproduce plots generated by `ggtree` with base R graphics, I've taken the example data sets from the `ggtree` distribution (https://github.com/YuLab-SMU/ggtree/blob/master/inst/examples/) and reprocessed them using base R, with the exception of `read.nexus` from the `ape` package:

```R
# read maximum clade credibility tree
flu <- read.nexus('MCC_FluA_H3.tree')

# extract dN/dS ratios from text file
mlc <- readLines('mlc.txt')
nwk <- mlc[ which(grepl("w ratios", mlc)) + 1 ]
nwk <- gsub("#", ":", nwk)
nwk <- gsub(" ", "", nwk)
dnds <- read.tree(text=nwk)

# annotate flu tree with dN/dS outcomes
flu <- ladderize(flu, right=F)
dnds <- ladderize(dnds, right=F)
flu$dnds <- dnds$edge.length

# export the result for loading with `ggfree`
save(flu, file="flu.RData")
```

Here is the example code from `ggtree` to draw the phylogeny annotated with host species and genotype (with a heatmap):
```R
cols <- scale_color(merged_tree, "dN_vs_dS", low="#0072B2", high="#D55E00",
interval=seq(0, 1.5, length.out=100))
p <- ggtree(merged_tree, size=.8, mrsd="2013-01-01", ndigits = 2, color=cols)
p <- add_colorbar(p, cols, font.size=3)

## add annotation of clade posterior and colored in darkgreen
p <- p + geom_text(aes(label=posterior), color="darkgreen", vjust=-.1,
hjust=-.03, size=1.8)

## add annotation of amino acid substitution inferred by joint probabilities
p <- p + geom_text(aes(x=branch, label=joint_AA_subs), vjust=-.03, size=1.8)

## get the tip labels
tip <- get.tree(merged_tree)$tip.label

## create a new data.frame containing information of host species.
host <- rep("Human", length(tip))
host[grep("Swine", tip)] <- "Swine"
host.df <- data.frame(taxa=tip, host=factor(host))

## use %<+% operator to attach the information to the tree view
p <- p %<+% host.df

## add circles in tips and tip labels align to the right hand side.
## after the attachment via %<+% operator,
## we can use 'host' information to color circles and labels of tips.
p <- p + geom_tippoint(aes(color=host), size=2) +
geom_tiplab(aes(color=host), align=TRUE, size=3, linesize=.3)

## rescale the color to blue (for human) and red (for swine)
p <- p + scale_color_manual(values = c("#377EB8", "#E41A1C"), guide='none')

## setting the coordinate system
p <- p + theme_tree2() + scale_y_continuous(expand=c(0, 0.6)) + xlab("Time")

## setting size of axis text and title.
p <- p + theme(axis.text.x=element_text(size=10), axis.title.x =
element_text(size=12))
p <- gheatmap(p, genotype, width=.4, offset=7, colnames=F) %>% scale_x_ggtree
p <- p + scale_fill_brewer(palette="Set2")
p <- p+theme(legend.text=element_text(size=8), legend.key.height=unit(.5,
"cm"),
legend.key.width=unit(.4, "cm"),
legend.position=c(.13, y=.945))
```
Note that I've omitted some of the code for loading and reformatting the data for brevity and for a more equitable comparison to `ggfree` syntax.


```{r dpi=300, fig.width=10, fig.height=10, out.width="80%"}
# now let's try the same thing with ggfree
L <- tree.layout(flu, 'r')
plot(L, cex=0.6, type='n', mar=c(3,1,0,20), label='n')

# map dN/dS values to colours
pal <- colorRampPalette(c('#0072B2', '#D55E00'))(20)
breaks <- c(seq(0, 1.5, length.out=19), 1000)
col <- pal[as.integer(cut(L$edges$dnds, breaks=breaks))]

# draw labels
host <- ifelse(grepl('Swine', L$nodes$label), '#E41A1C', '#377EB8')
text(L, align=TRUE, cex=0.75, col=host)

# draw tree
lines(L, col=col, lwd=3)
axis(side=1, at=seq(0, 20, 5), labels=seq(1990, 2010, 5), line=-2)

# draw points on tips
points(L, pch=20, col=host, cex=ifelse(L$nodes$n.tips==0, 1.5, 0))

# map colors from edges to nodes
index <- L$edges$child[L$edges$isTip] 
draw.guidelines(L, col=host[index])

# load genotype data (example from ggtree)
path <- system.file("extdata/Genotype.txt", package="ggfree")
geno <- read.table(path, header=T, sep='\t', na.strings='')
geno <- geno[match(flu$tip.label, row.names(geno)), ]

# draw boxes
require(RColorBrewer)
col <- brewer.pal(3, 'Set2')
image(L, geno, xlim=c(30, 37), col=col, cex.axis=0.75, line=-2)
```

This end result is pretty similar to the `ggtree` example (Figure 5 from [this paper](https://doi.org/10.1111/2041-210X.12628)), except that we're not embedding small images and not annotating the internal nodes with amino acid substitutions (which I find rather messy and uninformative!).


## Example: a bird tree

The `ape` package has several examples of phylogenetic trees that relates birds at different taxonomic levels.
I'm going to use one that relates birds at the level of families.
```{r fig.width=10, fig.height=10, out.width="50%"}
# load the example tree
data("bird.families")
plot(tree.layout(bird.families, type='o'), cex.lab=0.8)
```

I'd like to annotate this tree with the number of named species per family.
Since these data are not provided with this tree object, I had to search around for some supplementary data.
I settled on the dataset published by the non-governmental conservation organization [BirdLife International](http://datazone.birdlife.org/home), where I retrieved the ZIP archive including an Excel spreadsheet that I converted to a CSV and processed in R.
The resulting data frame covers substantially more bird families than our phylogeny:
```{r}
path <- system.file("extdata/birdlife.csv", package='ggfree')
birds <- read.csv(path, row.names=1)
head(birds)
nrow(birds)
```

Let's try to match up family names:
```{r}
index <- match(bird.families$tip.label, birds$Family)
# some family names failed to match
bird.families$tip.label[is.na(index)]
# There's only 13 of these missing, so I'm willing to manually restore them.
# I *think* Dacelonidae is the tree kingfishers subfamily...
missing <- data.frame(
  Family=c("Dendrocygnidae", "Bucorvidae", "Rhinopomastidae", "Dacelonidae", "Cerylidae", "Centropidae", "Coccyzidae", "Crotophagidae", "Neomorphidae", "Batrachostomidae", "Eurostopodidae", "Chionididae", "Eopsaltriidae"),
  Count=c(8, 2, 3, 70, 9, 10, 13, 4, 6, 5, 3, 2, 44)
)
birds <- rbind(birds, missing)
index <- match(bird.families$tip.label, birds$Family)
```

Now we're ready to use this information to annotate the tree: 
```{r fig.width=10, fig.height=10, out.width="100%"}
# the count distribution is more even if we do a log-transform
bins <- as.integer(cut(log(birds$Count[index]), breaks=8))
require(RColorBrewer)
pal <- brewer.pal(9, 'Blues')[2:9]

L <- tree.layout(bird.families, type='o')
plot(L, cex.lab=0.7, offset=2, mar=rep(5,4), col='chocolate')
image(L, z=as.matrix(bins), xlim=c(28.5,30), col=pal)
```
