`ro warning=FALSE, message=FALSE, comment=NA, cache=FALSE or`

---
name: phylogenetic-tree-balance
layout: post
title: Exploring phylogenetic tree balance metrics
date: 2012-10-10
author: Scott Chamberlain
tags: 
- R
- phylogenetic
- tree shape
---

I need to simulate balanced and unbalanced phylogenetic trees for some research I am doing.  In order to do this, I do rejection sampling: simulate a tree -> measure tree shape -> reject if not balanced or unbalanced __enough__.  But what is enough?  We need to define some cutoff value to determine what will be our set of balanced and unbalanced trees. 

### A function to calculate shape metrics
```{r theme } 
foo <- function(x, metric="colless") {
	if(metric == "colless"){
		xx <- as.treeshape(x) # convert to apTreeshape format
		colless(xx, "yule") # calculate colless' metric
	} else
		if(metric == "gamma"){
			gammaStat(x) 
		} else
			stop("metric should be one of colless or gamma")
}
````

### Simulate some trees
```{r simulatetrees } 
library(ape); library(phytools)

# lets simulate 1000 trees
numtrees <- 1000 

# simulate 500 pure-birth trees with 100 spp each
# ape = F makes it run faster
trees <- pbtree(n = 50, nsim = numtrees, ape = F) 
````

### Calculate Colless' shape metric on each tree
```{r collessprep }
library(plyr); library(apTreeshape)

# calculate metric for each tree
colless_df <- ldply(trees, foo, metric="colless") 
head(colless_df)

# Calculate the percent of trees that will 
# fall into the cutoff for balanced and unbalanced trees
col_percent_low <- round(length(colless_df[colless_df$V1 < -0.7, "V1"])/numtrees, 2)*100
col_percent_high <- round(length(colless_df[colless_df$V1 > 0.7, "V1"])/numtrees, 2)*100
````

### Create a distribution of the metric values
```{r fig.width=9, fig.height=5, collesshist }
library(ggplot2)

# plot histogram of distribution of values
ggplot(colless_df, aes(V1)) +  
	geom_histogram() + 
	theme_bw(base_size=18) + 
	scale_x_continuous(limits=c(-3,3), breaks=c(-3,-2,-1,0,1,2,3)) + 
	geom_vline(xintercept = -0.7, colour="red", linetype = "longdash") +
	geom_vline(xintercept = 0.7, colour="red", linetype = "longdash") +
	ggtitle(paste0("Distribution of Colless' metric for 1000 trees, cutoffs at -0.7 and 0.7 results in\n ", col_percent_low, "% (", numtrees*(col_percent_low/100), ") 'balanced' trees (left) and ", col_percent_low, "% (", numtrees*(col_percent_low/100), ") 'unbalanced' trees (right)")) +  
	labs(x = "Colless' Metric Value", y = "Number of phylogenetic trees") +
	theme(plot.title  = element_text(size = 16))

```

### And the same for Gamma stat, which measures the distribution of nodes in time. 
```{r fig.width=9, fig.height=5, gammaall }
gamma_df <- ldply(trees, foo, metric="gamma") # calculate metric for each tree
gam_percent_low <- round(length(gamma_df[gamma_df$V1 < -1, "V1"])/numtrees, 2)*100
gam_percent_high <- round(length(gamma_df[gamma_df$V1 > 1, "V1"])/numtrees, 2)*100
ggplot(gamma_df, aes(V1)) +  # plot histogram of distribution of values
	geom_histogram() + 
	theme_bw(base_size=18) + 
	scale_x_continuous(breaks=c(-3,-2,-1,0,1,2,3)) + 
	geom_vline(xintercept = -1, colour="red", linetype = "longdash") +
	geom_vline(xintercept = 1, colour="red", linetype = "longdash") +
	ggtitle(paste0("Distribution of Gamma metric for 1000 trees, cutoffs at -1 and 1 results in\n ", gam_percent_low, "% (", numtrees*(gam_percent_low/100), ") trees with deeper nodes (left) and\n ", gam_percent_high, "% (", numtrees*(gam_percent_high/100), ") trees with shallower nodes (right)")) +  
	labs(x = "Gamma Metric Value", y = "Number of phylogenetic trees") +
	theme(plot.title  = element_text(size = 16))
```

*********
#### Get the .Rmd file used to create this post [at my github account](https://github.com/sckott/sckott.github.io/blob/master/_drafts/2012-10-10-phylogenetic-tree-balance.Rmd) - or [.md file](https://github.com/sckott/sckott.github.io/tree/master/_posts/2012-10-10-phylogenetic-tree-balance.md).

#### Written in [Markdown](http://daringfireball.net/projects/markdown/), with help from [knitr](http://yihui.name/knitr/).