1 Introduction

The microbiome R package facilitates exploration and analysis of microbiome profiling data, in particular 16S taxonomic profiling.

This vignette provides a brief overview with example data sets from published microbiome profiling studies (Lahti et al. 2014, Lahti et al. (2013), O’Keefe et al. (2015)). A more comprehensive tutorial is available on-line.

Tools are provided for the manipulation, statistical analysis, and visualization of taxonomic profiling data. In addition to targeted case-control studies, the package facilitates scalable exploration of large population cohorts (Lahti et al. 2014). Whereas sample collections are rapidly accumulating for the human body and other environments, few general-purpose tools for targeted microbiome analysis are available in R. This package supports the independent phyloseq data format and expands the available toolkit in order to facilitate the standardization of the analyses and the development of best practices. See also the related PathoStat pipeline mare pipeline, phylofactor, and structSSI for additional 16S rRNA amplicon analysis tools in R. The aim is to complement the other available packages, but in some cases alternative solutions have been necessary in order to streamline the tools and to improve complementarity.

We welcome feedback, bug reports, and suggestions for new features from the user community via the issue tracker and pull requests. See the Github site for source code and other details. These R tools have been utilized in recent publications and in introductory courses (Salonen et al. 2014, Faust et al. (2015), Shetty et al. (2017)), and they are released under the Two-clause FreeBSD license.

Kindly cite the work as follows: “Leo Lahti et al. (Bioconductor, 2017). Tools for microbiome analysis in R. Microbiome package version . URL: (http://microbiome.github.io/microbiome)

2 Installation

To install microbiome package in R (Bioconductor development version), use

library(BiocManager)
if (!requireNamespace("BiocManager", quietly=TRUE))
    install.packages("BiocManager")
BiocManager::install(version = "devel")
BiocManager::install("microbiome")

Then load the package in R

library(microbiome)  

3 Data

The microbiome package relies on the independent phyloseq data format. This contains an OTU table (taxa abundances), sample metadata (age, BMI, sex, …), taxonomy table (mapping between OTUs and higher-level taxonomic classifications), and a phylogenetic tree (relations between the taxa).

# Data from 
# http://www.nature.com/ncomms/2014/140708/ncomms5344/full/ncomms5344.html
data(atlas1006) 
atlas1006

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 130 taxa and 1172 samples ]
## sample_data() Sample Data:       [ 1172 samples by 10 sample variables ]
## tax_table()   Taxonomy Table:    [ 130 taxa by 2 taxonomic ranks ]

data(dietswap) # Data from http://dx.doi.org/10.1038/ncomms7342
dietswap

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 130 taxa and 222 samples ]
## sample_data() Sample Data:       [ 222 samples by 8 sample variables ]
## tax_table()   Taxonomy Table:    [ 130 taxa by 2 taxonomic ranks ]

data(peerj32) # Data from https://peerj.com/articles/32/

# dietswap is a phyloseq object; see above
dietswap.compositional <- transform(dietswap, "compositional")

4 Ecosystem indices

g <- global(atlas1006, index = "gini")

# Estimate Shannon diversity and add it to the phyloseq object
sample_data(atlas1006)$diversity <- global(atlas1006, index = "shannon")[,1]

# Compare age and microbiome diversity
plot_regression(diversity ~ age, meta(atlas1006))

5 Core microbiota analysis

Population frequencies, or prevalence, of the taxonomic groups exceeding a given detection threshold can be calculated with

p <- prevalence(dietswap, detection = 0, sort = TRUE)

The core microbiota refers to the set of taxa that are detected in a remarkable fraction of the population above a given abundance threshold (see e.g. Jalanka-Tuovinen et al. 2011, Salonen et al. (2012)). The core subset can be sliced from a phyloseq object as follows.

# Taxa with over 50% prevance at .2% relative abundance
dietswap.core <- core(dietswap.compositional, 
                    detection = .2/100, prevalence = 50/100)
dietswap.core

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 39 taxa and 222 samples ]
## sample_data() Sample Data:       [ 222 samples by 8 sample variables ]
## tax_table()   Taxonomy Table:    [ 39 taxa by 2 taxonomic ranks ]

Similar functions are available also for rare and variable taxa, see the microbiome tutorial for a full description.

To visualize the core as in (Shetty et al. 2017), use

library(ggplot2, quiet = TRUE)
p <- plot_core(transform(dietswap.core, "compositional"), 
    plot.type = "heatmap", 
    colours = gray(seq(0,1,length=5)),
    prevalences = seq(.05, 1, .05), 
    detections = 10^seq(log10(1e-3), log10(.2), length = 10), 
    horizontal = TRUE) +
    xlab("Detection Threshold (Relative Abundance (%))") 
print(p)    

6 Microbiome composition

Composition heatmap: Z-transformed taxon abundances

tmp <- plot_composition(dietswap.core, plot.type = "heatmap", transform = "Z", 
            mar = c(6, 13, 1, 1), sample.sort = "nationality")

Composition barplot

plot_composition(transform(dietswap.core, "compositional"), 
    plot.type = "barplot", sample.sort = "neatmap")

Visualize sample similarities, or the microbiome landscape (Shetty et al. 2017) on an ordination map based on various ordination methods (PCoA, NMDS, RDA etc) that are available through various other R packages.

# Project the samples with the given method and dissimilarity measure. 
# Ordinate the data; note that some ordinations are sensitive to random seed
# "quiet" is used to suppress intermediate outputs
set.seed(423542)

# Get ordination
x <- dietswap.core
quiet(x.ord <- ordinate(x, method = "NMDS", distance = "bray"))
# Pick the projected data (first two columns + metadata)
quiet(proj <- phyloseq::plot_ordination(x, x.ord, justDF=TRUE))
# Rename the projection axes
names(proj)[1:2] <- paste("Comp", 1:2, sep=".")

p <- plot_landscape(proj[, 1:2], col = proj$nationality, legend = TRUE)
print(p)

7 Association heatmaps

Let us cross-correlate example data containing microbiome profiling and blood serum lipids study (Lahti et al. 2013).

# Define data sets to cross-correlate
data(peerj32)
x <- log10(peerj32$microbes)   # OTU Log10 (44 samples x 130 genera)
y <- as.matrix(peerj32$lipids) # Lipids (44 samples x 389 lipids)

# Cross correlate data sets and return a table
correlation.table <- associate(x, y, 
    method = "spearman", mode = "table", p.adj.threshold = 0.05, n.signif = 1)

kable(head(correlation.table))

Let us rearrange the data and plot the heatmap and highlight significant correlations with stars as in (Lahti et al. 2013).

heat(correlation.table, "X1", "X2", fill = "Correlation", 
    star = "p.adj", p.adj.threshold = 0.05) 

8 Bistability and tipping elements

The microbiome package provides tools to quantify stability and bimodality in abundance data, following (Lahti et al. 2014). Let use Dialister as an example. The tipping point is here set manually but it can also be estimated from the data.

# Use relative abundances and plot subject variation 
# during the follow-up period
pseq <- transform(atlas1006, "compositional")
plot_tipping(pseq, "Dialister", tipping.point = 0.004)

# Show the bimodal population distribution
hotplot(pseq, "Dialister", tipping.point = 0.004)

9 Further reading

The on-line tutorial provides many additional tools and examples, with more thorough descriptions. This package and tutorials are work in progress. We welcome feedback, for instance via issue Tracker, pull requests, or via Gitter.

10 Acknowledgements

Thanks to all contributors. Financial support has been provided by Academy of Finland (grants 256950 and 295741), University of Turku, Department of Mathematics and Statistics. In addition, the work has been supported by VIB lab for Bioinformatics and (eco-)systems biology, VIB/KULeuven, Belgium, Molecular Ecology group, Laboratory of Microbiology, Wageningen University, Netherlands, and Department of Veterinary Bioscience, University of Helsinki, Finland. This work relies on the independent phyloseq package and data structures for R-based microbiome analysis developed by Paul McMurdie and Susan Holmes. This work also utilizes a number of independent R extensions, including dplyr (Wickham and Francois 2016), ggplot2 (Wickham 2009), phyloseq (McMurdie and Holmes 2013), and vegan (Oksanen et al. 2015).