Microbiome Data Basics

Interpreting observed counts

Each entry in a marker-gene count table is the number of retained reads assigned to one feature in one sample. The table reflects finite sampling and processing; it is not a census of microbial cells in the specimen.

Marker-gene terminology

Marker gene

A short DNA region used as a signpost for describing which microbes may be present.

Amplification

Making many copies of the chosen marker so the sequencing instrument can measure it.

Amplicon

The copied DNA fragment produced during amplification; this tutorial analyzes such fragments.

Taxonomy

Nested names used to classify organisms, from broad groups to narrower groups.

Gene-copy number

The number of copies of the marker gene in a genome; this can differ among microbes.

Phylogenetic tree

A branching diagram representing inferred evolutionary relatedness among features.

Sample

One collected specimen represented by one row in our teaching table.

Read

One sequencing observation retained after processing.

Feature

A unit produced by the pipeline, such as an ASV, OTU, taxon, gene, or pathway.

Count

The number of retained reads assigned to one feature in one sample.

ASV

An amplicon sequence variant: a sequence-level feature, not automatically a species.

The linked data objects

Many marker-gene microbiome projects use three linked tables. Sample IDs connect counts to sample metadata; feature IDs connect counts to taxonomy. Some projects also include representative sequences and a phylogenetic tree.

1. Sample-by-feature count table

Each cell is an observed count. The constructed table makes S2 exactly twice the counts of S1 and S4 exactly twice the counts of S3. B1 is a low-count negative control.

2. Sample metadata

Metadata describe samples: treatment, time, subject, site, batch, and similar variables. Batch is not a nuisance you may ignore; it can be confused with the scientific group if the study is poorly balanced.

3. Taxonomy table

A taxonomy table attaches hierarchical labels to feature IDs. Classification can stop at different levels or be uncertain. A phylum is a broader group than a genus. Keep the original ASV IDs; a missing genus label does not mean the feature count is zero.

Taxonomic labels are illustrative; no biological or functional inference is intended.

Validate IDs and calculate sequencing depth

Validate table structure and identifier alignment before analysis.

Import the three teaching tables

Download the files from the sidebar, preserve their names, and place them together in your project’s data folder.

# Run once if tidyverse is not already installed:
install.packages("tidyverse")

# Run at the start of each new session:
library(tidyverse)

counts <- read_csv(
  "data/microbiome_counts.csv",
  show_col_types = FALSE
)
metadata <- read_csv(
  "data/microbiome_metadata.csv",
  show_col_types = FALSE
)
taxonomy <- read_csv(
  "data/microbiome_taxonomy.csv",
  show_col_types = FALSE
)

count_matrix <- as.matrix(counts[, -1])
rownames(count_matrix) <- counts$sample_id

stopifnot(anyDuplicated(counts$sample_id) == 0)
stopifnot(anyDuplicated(metadata$sample_id) == 0)
stopifnot(anyDuplicated(taxonomy$feature_id) == 0)
stopifnot(setequal(counts$sample_id, metadata$sample_id))
stopifnot(setequal(colnames(count_matrix), taxonomy$feature_id))
stopifnot(!anyNA(count_matrix))
stopifnot(all(count_matrix >= 0))
stopifnot(all(count_matrix == floor(count_matrix)))

stopifnot() stops execution when an assumption is false. setequal() checks that the same IDs appear regardless of order. Before joining values, reorder explicitly and confirm identical() order.

In counts[, -1], the empty space means “all rows” and -1 removes the first column, which contains IDs rather than counts. as.matrix() then creates a numeric matrix; rownames() attaches the sample IDs to its rows.

Library size is a row sum

library_size <- rowSums(count_matrix)
library_size

B1 is a negative control with six retained reads; exclude it from biological summaries and evaluate it during contamination assessment.

Why an observed zero is difficult

A zero in the count table says only that no read was observed for that sample-feature combination. Different underlying processes can produce the same visible value.

Structural zero

The feature is genuinely absent from the target community under the scientific definition.

Sampling zero

The feature is present at a low level but was not captured by finite sampling or measurement.

NA

The value is unknown or unavailable. It is missing data, not an observed count of zero.

Zeros in the synthetic dataset

• ASV3 is zero in S1 and S2 but observed in S3 and S4. The table alone does not reveal why.
• ASV4 occurs only in negative control B1. That is a contamination warning to investigate, not automatic proof.
• ASV6 is observed in treatment samples only. This small descriptive dataset cannot establish a treatment effect.

Relative abundance and compositionality

Relative abundance divides each feature count by its sample’s library size. It answers “What share of the observed reads belongs to this feature?” Every sample is rescaled to a total of 1, or 100%.

relative_matrix <- sweep(
  count_matrix,
  MARGIN = 1,
  STATS = rowSums(count_matrix),
  FUN = "/"
)

round(relative_matrix, 3)
rowSums(relative_matrix)

sweep() divides every row by its own row sum. Confirm that each nonempty row now sums to 1. Handle zero-library samples explicitly before division; otherwise they produce undefined values.

The compositional constraint

Because proportions sum to one, increasing one component’s share necessarily reduces at least one other share. Relative abundance therefore describes ratios among observed components, not absolute microbial load.

• Relative abundance is useful for descriptive plots.
• It does not recover absolute microbial load.
• “Twice the proportion” does not necessarily mean “twice as many cells.”
• Counts, proportions, and log-ratios are different representations and can answer different questions.
• Always retain the original count table.

Prevalence, abundance, and filtering

Prevalence is the proportion of samples in which a feature was detected. Abundance describes how much was observed. A feature can be widespread but low in abundance, or abundant in only a few samples.

counts_long <- counts |>
  pivot_longer(
    -sample_id,
    names_to = "feature_id",
    values_to = "count"
  )

biological_ids <- metadata |>
  filter(sample_type == "Biological") |>
  pull(sample_id)

feature_summary <- counts_long |>
  filter(sample_id %in% biological_ids) |>
  group_by(feature_id) |>
  summarise(
    detected = sum(count > 0),
    prevalence = mean(count > 0),
    total_reads = sum(count),
    .groups = "drop"
  )

feature_summary

Alpha diversity: within one sample

Alpha diversity summarizes within-sample feature variety. No single index represents all aspects of diversity; different indices emphasize different properties of the observed community.

Observed richness

Richness records which features were detected, but ignores how reads are distributed. It is sensitive to sequencing effort, detection, and filtering.

Shannon diversity

Shannon diversity reflects both richness and evenness. With the same richness, a sample whose reads are spread more evenly across features usually has a larger value than one dominated by a single feature. R’s log() uses the natural logarithm, so that is the convention used below.

bio_matrix <- count_matrix[
  rownames(count_matrix) %in% biological_ids,
  ,
  drop = FALSE
]

bio_relative <- sweep(
  bio_matrix,
  1,
  rowSums(bio_matrix),
  "/"
)

alpha_summary <- tibble(
  sample_id = rownames(bio_matrix),
  library_size = rowSums(bio_matrix),
  observed_richness = rowSums(bio_matrix > 0),
  shannon = apply(
    bio_relative,
    1,
    function(p) -sum(p[p > 0] * log(p[p > 0]))
  )
)

alpha_summary

Beta diversity, dissimilarity, and ordination

Beta-diversity analyses quantify differences among sample profiles. The primary object is a distance or dissimilarity matrix containing one value for each sample pair.

Terminology varies across fields, so do not report only “beta diversity.” Name the representation, exact distance or dissimilarity, feature set, and downstream method.

Bray–Curtis dissimilarity

• 0 means the two input profiles are identical.
• Values nearer 1 indicate greater dissimilarity under this definition.
• The answer depends on the values you put into the formula and which features you retained.

bray_curtis <- function(x, y) {
  sum(abs(x - y)) / sum(x + y)
}

bray_curtis(
  count_matrix["S1", ],
  count_matrix["S2", ]
)

bray_curtis(
  relative_matrix["S1", ],
  relative_matrix["S2", ]
)

Raw counts make S1 and S2 look different because S2 has twice the total count. Their relative profiles are identical, so relative-abundance Bray–Curtis is 0. Preprocessing is part of the comparison’s definition, not an invisible preliminary step.

In this matrix, S1 versus S3 is 0.718, whereas S1 versus S2 is 0. PCoA begins with a complete dissimilarity matrix such as this.

Principal coordinates analysis

Principal coordinates analysis (PCoA) places samples on a few axes so the displayed distances approximate the original pairwise distances. A two-dimensional display usually loses information, so report how much variation the displayed axes represent.

Microbiome analysis workflow

Begin with the study question, measurement process, and structure of each data object before selecting analytical or visualization methods.

1. Write the scientific question in one sentence.
2. Identify the count table’s row and column orientation.
3. Match sample and feature IDs across objects.
4. Inspect metadata, controls, batch, missingness, and possible confounders.
5. Calculate library sizes and examine their distribution.
6. Summarize prevalence and abundance without silently discarding features.
7. Choose and document filtering based on the question and design.
8. Choose a representation and distance appropriate to the target quantity.
9. Calculate within-sample or between-sample summaries.
10. Interpret results together with study design and uncertainty.

Build a long table for plotting

count_feature_ids <- setdiff(names(counts), "sample_id")

unmatched_samples <- anti_join(
  counts |> distinct(sample_id),
  metadata |> distinct(sample_id),
  by = "sample_id"
)
unmatched_features <- anti_join(
  tibble(feature_id = count_feature_ids),
  taxonomy |> distinct(feature_id),
  by = "feature_id"
)

stopifnot(nrow(unmatched_samples) == 0)
stopifnot(nrow(unmatched_features) == 0)

plot_data <- counts |>
  pivot_longer(
    -sample_id,
    names_to = "feature_id",
    values_to = "count"
  ) |>
  left_join(metadata, by = "sample_id") |>
  left_join(taxonomy, by = "feature_id") |>
  group_by(sample_id) |>
  mutate(
    library_size = sum(count),
    relative_abundance = if_else(
      library_size > 0,
      count / library_size,
      NA_real_
    )
  ) |>
  ungroup()

stopifnot(
  nrow(plot_data) == nrow(counts) * (ncol(counts) - 1)
)

The two anti_join() checks detect count-table IDs with no matching metadata or taxonomy record. The final row-count check detects unexpected row multiplication; ncol(counts) - 1 is the number of feature columns after excluding sample_id. Never treat a successful join as proof that labels are correct.

Plot observed composition

plot_data |>
  filter(sample_type == "Biological") |>
  ggplot(
    aes(
      x = sample_id,
      y = relative_abundance,
      fill = feature_id
    )
  ) +
  geom_col(width = 0.78) +
  facet_grid(~ group, scales = "free_x", space = "free_x") +
  scale_y_continuous(
    labels = function(x) paste0(x * 100, "%")
  ) +
  labs(
    title = "Observed composition by sample",
    subtitle = "Synthetic tutorial data",
    x = NULL,
    y = "Relative abundance",
    fill = "Feature"
  ) +
  theme_minimal(base_size = 12)

Stacked bars reveal dominant features and sample composition. They are weak for comparing thin segments across many samples. For larger data, show only a declared set of features, combine the remainder as “Other,” and preserve the full table for analysis.

Exercises and reference

Glossary

Abundance

How much of a feature was observed under a stated representation.

Alpha diversity

Diversity summarized within one sample.

Beta diversity

Variation or dissimilarity among samples; always state the definition used.

Composition

Parts represented relative to a constrained whole.

Feature

The unit represented by a count-table row or column, such as an ASV or pathway.

Library size

Total retained reads observed for one sample.

Metadata

Variables describing samples, subjects, collection, or processing.

Ordination

A lower-dimensional display of multivariate relationships.

Prevalence

Proportion of samples in which a feature was detected.

Relative abundance

A feature count divided by the sample’s total count.

Sampling zero

A possibly present feature was not observed in finite sampling.

Structural zero

A feature is genuinely absent from the target community.

Taxonomy

Hierarchical labels assigned to microbial features.

References and further reading

• Callahan et al. (2016) — exact sequence inference and ASVs with DADA2.
• McMurdie and Holmes (2013) — organizing linked microbiome data objects in phyloseq.
• Vandeputte et al. (2017) — relative versus quantitative microbiome profiling.
• Kaul et al. (2017) — microbiome data with excess and structural zeros.
• Schmidt et al. (2022) — sampling, structural, and missing zeros in diversity analysis.
• Gloor et al. (2017) — why microbiome datasets are compositional.
• Weiss et al. (2017) — how normalization choices interact with data characteristics.
• McMurdie and Holmes (2014) — tradeoffs in normalizing microbiome count data.
• Whittaker (1960) — foundational alpha- and beta-diversity concepts.
• Bray and Curtis (1957) — the dissimilarity used in the worked example.
• Gower (1966) — principal coordinates analysis.
• Aitchison (1982) — statistical analysis of compositional data.
• Lozupone and Knight (2005) — UniFrac and phylogenetic comparison.