Microbiome data science workflow

Standard components of a taxonomic profiling study

Monday, September 28, 2026

The data science workflow

Data science workflow: import, tidy, then an understand loop of transform, visualize and model, followed by communicate; programming surrounds all steps.

A typical workflow

Import Preprocess Alpha Beta Taxa Report

Key elements of a taxonomic profiling study
Mixed and matched in real case studies
Examples for each step in OMA

Your task

Build a clear, compact and reproducible report on your data.

📓

Reproducible

Quarto / R Markdown report that renders from raw data

🎯

Selective

Choose summaries and analyses that answer your question

📚

Build on OMA

Reuse material and examples from the OMA book

Setup

library(mia)
library(miaViz)
library(scater)
library(patchwork)
set.seed(1)  # reproducible permutation tests

# Example data set: skin microbiome with sample metadata
data("peerj13075", package = "mia")
tse <- peerj13075

# Relative abundances and CLR for downstream steps
tse <- transformAssay(tse, method = "relabundance")
tse <- transformAssay(tse, method = "clr", pseudocount = TRUE)

# Genus-level summary
tse_genus <- agglomerateByRank(tse, rank = "genus")

Alpha diversity

Alpha diversity: within-sample

Tasks

Estimate alpha diversity per sample
Draw a histogram
Compare two or more indices, visually and/or statistically

tse <- addAlpha(
    tse, index = c("shannon", "observed"))

p1 <- plotHistogram(
    tse, col.var = "shannon", bins = 15)
p2 <- plotColData(
    tse, x = "observed", y = "shannon")
p1 + p2

Tip

Do richness and evenness agree? Here Spearman \(\rho\) = 0.43.

Alpha diversity: within-sample

Beta diversity

Beta diversity: PCoA

Visualize community variation with PCoA, and check how the data transformation affects it:

Bray–Curtis

on compositional (relative abundance) data

Euclidean

on CLR-transformed data (Aitchison distance)

tse <- runMDS(tse, FUN = getDissimilarity, method = "bray",
              assay.type = "relabundance", name = "PCoA_BC")
tse <- runMDS(tse, FUN = getDissimilarity, method = "euclidean",
              assay.type = "clr", name = "PCoA_CLR")

plotReducedDim(tse, "PCoA_BC",  colour_by = "Geographical_location") +
plotReducedDim(tse, "PCoA_CLR", colour_by = "Geographical_location") +
    plot_layout(guides = "collect")

Beta diversity: PCoA

Community-level comparisons

Does community composition differ between groups?

PERMANOVA tests group differences in composition
Add covariates (e.g. age, sex) and see how the results change
Check homogeneity of dispersion

res <- getPERMANOVA(
    tse,
    assay.type = "relabundance",
    method = "bray",
    formula = x ~ Geographical_location + Gender + Age,
    test.homogeneity = TRUE)

res$permanova |>
    knitr::kable(digits = 3)

Community-level comparisons

	Df	SumOfSqs	R2	F	Pr(>F)
Geographical_location	2	1.956	0.093	2.887	0.001
Gender	1	0.238	0.011	0.702	0.698
Age	2	0.872	0.041	1.287	0.199
Residual	52	17.616	0.836	NA	NA
Total	57	21.079	1.000	NA	NA

Taxa-level analysis

Prevalence

What is the most prevalent genus in the data?

prev <- getPrevalence(
    tse_genus, assay.type = "relabundance",
    detection = 0, sort = TRUE)
round(head(prev, 4), 2)

 Staphylococcus   Paenibacillus        Bacillus Corynebacterium 
           1.00            1.00            1.00            0.93

prevalence = fraction of samples where a taxon is detected above a given abundance threshold

Core microbiota

Taxa above 0.1 % relative abundance in over 50 % of samples (Salonen et al. 2012). How sensitive is this to the thresholds?

core <- getPrevalent(
    tse_genus,
    assay.type = "relabundance",
    detection = 0.1 / 100,
    prevalence = 50 / 100)

plotAbundanceDensity(
    tse_genus[core, ],
    layout = "jitter",
    assay.type = "relabundance")

Differential abundance

Which genera are associated with a group difference (e.g. sex)?

library(DESeq2)

dds <- DESeqDataSet(tse_genus, design = ~ Gender)
dds <- DESeq(dds)
res <- results(dds)

top <- head(res[order(res$padj), ])
data.frame(log2FC = round(top$log2FoldChange, 2),
           padj = signif(top$padj, 2),
           row.names = rownames(top)) |>
    knitr::kable()

DESeq2 (Love, Huber, and Anders 2014); for alternatives and benchmarks see the OMA DAA chapter
Compare several methods: results can differ substantially

Differential abundance

	log2FC	padj
Clostridium	4.16	4.2e-06
Pseudomonas	5.06	4.2e-06
Finegoldia	5.93	2.1e-04
Orientia	4.23	6.1e-03
Kocuria	-4.51	9.7e-03
Amphritea	5.42	1.4e-02

Role of covariates

Key covariates (diet, medication, age, …) can change the interpretation.

Stool consistency, medication and diet explain a considerable part of population-level gut microbiome variation (Falony et al. 2016)
Include relevant covariates in PERMANOVA and differential abundance models
Compare results with and without covariates

Wrap-up

Checklist for your report

Alpha

Indices, histogram, comparison

Beta

PCoA (BC vs CLR), PERMANOVA with covariates

Taxa

Prevalence, core microbiota, differential abundance

Report

Clear, compact, reproducible

Resources

References

Falony, Gwen, Marie Joossens, Sara Vieira-Silva, Jun Wang, Youssef Darzi, Karoline Faust, Alexander Kurilshikov, et al. 2016. “Population-Level Analysis of Gut Microbiome Variation.” Science 352 (6285): 560–64. https://doi.org/10.1126/science.aad3503.

Love, Michael I., Wolfgang Huber, and Simon Anders. 2014. “Moderated Estimation of Fold Change and Dispersion for RNA-seq Data with DESeq2.” Genome Biology 15 (12): 550. https://doi.org/10.1186/s13059-014-0550-8.

Salonen, A., J. Salojärvi, L. Lahti, and W. M. de Vos. 2012. “The Adult Intestinal Core Microbiota Is Determined by Analysis Depth and Health Status.” Clinical Microbiology and Infection 18 (s4): 16–20. https://doi.org/10.1111/j.1469-0691.2012.03855.x.

Wickham, Hadley, Mine Çetinkaya-Rundel, and Garrett Grolemund. 2023. R for Data Science: Import, Tidy, Transform, Visualize, and Model Data. 2nd ed. O’Reilly Media. https://r4ds.hadley.nz.