Get example data - HITChip Atlas of 130 genus-like taxa across 1006 healthy western adults.

# 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. Sarle’s bimodality coefficient is available as well; and for classical test of unimodality, see the DIP test.

# 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

# 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 and plot the subject abundance variation lineplot (Variation tip plot) and Bimodality hotplot for a given taxon as in Lahti et al. 2014. The Dialister has bimodal population distribution and reduced temporal stability within subjects at intermediate abundances.

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

# 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)