Bimodality.Rmd

---
title: "Bimodality analysis"
author: "Leo Lahti, Sudarshan Shetty et al. `r Sys.Date()`"
bibliography: 
- bibliography.bib
output:
  BiocStyle::html_document:
    number_sections: no
    toc: yes
    toc_depth: 4
    toc_float: true
    self_contained: true
    thumbnails: true
    lightbox: true
    gallery: true
    use_bookdown: false
    highlight: haddock
---
<!--
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteIndexEntry{microbiome tutorial - stability}
  %\usepackage[utf8]{inputenc}
  %\VignetteEncoding{UTF-8}  
-->


Get example data - [HITChip Atlas of 130 genus-like taxa across 1006 healthy western adults](http://www.nature.com/ncomms/2014/140708/ncomms5344/full/ncomms5344.html). 

```{r bistability, message=FALSE}
# Load the example data
library(microbiome)
library(dplyr)
data(atlas1006)

# Rename the example data
pseq <- atlas1006

# Focus on specific DNA extraction method
pseq <- pseq %>% subset_samples(DNA_extraction_method == "r")

# Use relative abundances
pseq <- microbiome::transform(pseq, "compositional")

# Merge rare taxa to speed up examples
pseq <- aggregate_rare(pseq, level = "Genus", detection = .1/100, prevalence = 10/100)

# For cross-sectional analysis, include
# only the zero time point:
pseq0 <- subset_samples(pseq, time == 0)
```


# Bimodality indicators

Bimodality of the abundance distribution provides an indirect
indicator of bistability, although other explanations such as sampling
biases etc. should be controlled. Multiple bimodality scores are
available.


Multimodality score using [potential analysis with
bootstrap](http://www.nature.com/ncomms/2014/140708/ncomms5344/full/ncomms5344.html). Sarle's bimodality coefficient is available as well; and for classical test of unimodality, see the DIP test.

```{r bimodality2, message=FALSE, warning=FALSE}
# Bimodality is better estimated from log10 abundances
pseq0.clr <- microbiome::transform(pseq0, "clr")
bimodality <- bimodality(pseq0.clr, method = "potential_analysis", bs.iter = 20)
```

# Visualization   

**Visualize population densities for unimodal and bimodal groups**

```{r stability2, message=FALSE, warning=FALSE, fig.width=12, fig.height=5, out.width="500px"}
# Pick the most and least bimodal taxa as examples
unimodal  <- names(sort(bimodality))[[1]]
bimodal  <- rev(names(sort(bimodality)))[[1]]

# Visualize population frequencies at the baseline time point
library(ggplot2)
theme_set(theme_bw(20))
p1 <- plot_density(pseq0.clr, variable = unimodal) 
p2 <- plot_density(pseq0.clr, variable = bimodal) 
library(gridExtra)
library(ggplot2)
grid.arrange(p1, p2, nrow = 1)
```


## Variation lineplot and bimodality hotplot

Pick subset of the [HITChip Atlas data set](http://doi.org/10.5061/dryad.pk75d) and plot the subject abundance variation lineplot (**Variation tip plot**) and **Bimodality hotplot** for a given taxon as in [Lahti et al. 2014](http://www.nature.com/ncomms/2014/140708/ncomms5344/full/ncomms5344.html). The Dialister has bimodal population distribution and reduced temporal stability within subjects at intermediate abundances.

For examples on tipping point detection, see
[Stability](Stability). We set the tipping point manually in the
following example.

```{r stability-variationplot, message=FALSE, warning=FALSE, fig.show='hold', out.width="430px", eval=FALSE}
# Bimodality hotplot:
# Consider a unique sample from each subject: the baseline time point 
p <- hotplot(pseq0, tax, tipping.point = 0.005)
print(p)

# Visualize bimodality
pv <- plot_tipping(pseq, tax, tipping.point = 0.005)
print(pv)
```