16S Microbiome Downstream Analysis in R:
from QIIME2 output to publication figures.

You have run your samples through QIIME2. The denoising is done, the feature table exists, taxonomy has been assigned. You are now looking at hundreds of ASVs across dozens of samples — and the question that always follows is the same: what do I do with this now, and what are these numbers actually telling me about biology?

This is a guide through the downstream data processing using R — not a step-by-step tutorial, but the conceptual map of the analysis steps.

Taxonomic composition in phyloseq: phylum-level profiling with ggplot2

Every investigation starts with a survey of the scene. Before any statistical test, the first question is compositional: which microbial phyla dominate, and does the pattern vary across groups? If it does not, everything downstream will be harder to interpret. If it does — as it almost always does across distinct body sites or conditions — you already have your first biological result.

The standard approach aggregates ASVs to phylum level, normalizes each sample to relative abundance, and displays the result as stacked bars. This visualization is intended to provide an initial overview of the community structure, allowing you to spot broad patterns between groups before moving to more granular statistical testing.

The signal is immediate: the gut is defined by Firmicutes and Bacteroidetes, a ratio shaped by diet, age, and antibiotic history. The palms carry more Actinobacteria — a phylum that includes skin-adapted genera like Corynebacterium. The tongue occupies its own distinct space, with Fusobacteria visible only there.

This composition plot sets the expectation for everything that follows. Groups that look different here should show statistically significant differences in diversity analyses. Groups that look similar here should make you question what you are really comparing.

Alpha diversity in R: Shannon index, Simpson, and Observed features with Kruskal-Wallis

Alpha diversity answers a deceptively simple question: how diverse is the microbial community within a single sample? The answer changes depending on what aspect of diversity you are measuring.

Running multiple indices in parallel gives a more complete picture. Shannon and Simpson capture different aspects of the same community: if they tell the same story, your result is robust. If they diverge, it suggests the difference between groups is driven by rare taxa (Shannon is more sensitive to rare taxa than Simpson).

Rarefaction in R: how to standardize sequencing depth before diversity analysis

Before calculating any diversity metric, sequencing depth must be standardized. If one sample has several times the reads of another, comparing their Observed features is meaningless — you will always find more species when you look harder, regardless of the true community diversity. The first step is to look at how reads are distributed across your samples.

The histogram immediately shows the problem: depth varies nearly tenfold across samples. Rarefaction solves this by randomly subsampling every sample down to a common depth — here, 90% of the minimum, with set.seed(42) for reproducibility. Samples below the threshold are excluded. This is conservative, easy to justify in a Methods section, and retains nearly all samples in a well-designed study.

But rarefying to an arbitrary depth is only justified if that depth is sufficient to capture the real community diversity. The rarefaction curve tells you whether it is.

A curve that reaches a clear plateau means you have sequenced deeply enough — additional reads would not reveal meaningfully more taxa. In this dataset, all groups plateau before 2,000 reads, which validates the rarefaction depth choice and confirms that diversity comparisons across body sites are on solid ground.

The right statistical test for alpha diversity

Alpha diversity distributions are rarely normal. They are bounded, often skewed, and heteroskedastic across groups. The appropriate test is non-parametric: Kruskal-Wallis for the global comparison, Dunn's post-hoc test for pairwise comparisons, and Benjamini-Hochberg FDR correction throughout. The significance stars in the plot above reflect p-values after that correction — not raw p-values, which would dramatically overstate the evidence.

Beta diversity in R: Bray-Curtis, Jaccard, and UniFrac with phyloseq and vegan

Where alpha diversity asks "how diverse is this sample?", beta diversity asks "how different is this sample from that one?" The answer is a distance matrix — one number for every pair of samples — and the choice of distance metric changes what biological signal you capture.

Using Bray-Curtis and Weighted UniFrac together gives you a fuller picture. Bray-Curtis captures abundance-weighted compositional differences; Weighted UniFrac embeds those differences in an evolutionary context. When both metrics show the same pattern, the result is robust.

NMDS ordination in R: visualizing microbiome beta diversity with phyloseq

The full distance matrix — all pairwise distances in this dataset — cannot be shown directly. NMDS (Non-metric Multidimensional Scaling) compresses it into two dimensions while preserving the rank order of distances. Points that cluster together have similar communities; points that are far apart do not. The stress value tells you how faithfully the 2D representation captures the true distances: below 0.10 is good.

The NMDS reveals something biologically meaningful beyond the simple group separation: the gut cluster is tight — gut communities are consistent across samples and time points. The palm sites overlap with each other but not with gut or tongue — left palm and right palm share a similar skin environment. Tongue samples cluster separately from all three, reflecting the unique oral microbiome ecology.

PCoA ordination: Weighted UniFrac and variance explained in phyloseq

PCoA (Principal Coordinates Analysis) preserves actual distances, not just their rank order. The axes are interpretable — they represent the directions of maximum variance in the distance matrix — and the percentage of variance explained by each axis is shown directly on the plot. When applied to Bray-Curtis, PCoA and NMDS should tell the same story. If they do not, investigate why.

PERMANOVA and betadisper in R: testing beta diversity significance with vegan adonis2

An ordination plot shows that groups look separated. PERMANOVA tests whether that separation is real. It asks: if you randomly shuffled group labels 999 times, how often would you see a separation at least as strong as what you actually observed?

The result is an R² (what fraction of total community variance is explained by your grouping variable) and a permutation-based p-value. Body site explains a large fraction of the variance in community composition — a large and robust effect.

Differential abundance analysis in R: Wilcoxon, Kruskal-Wallis, and FDR correction for 16S data

Beta diversity tells you the communities are different. Differential abundance testing names the suspects. It is the transition from the general — "the groups differ" — to the specific: "here is the genus that distinguishes them, and here is the evidence."

For each genus, test whether its relative abundance differs significantly across groups. Two groups: Wilcoxon rank-sum test. Three or more groups: Kruskal-Wallis, followed by Dunn's post-hoc for pairwise comparisons. Throughout: Benjamini-Hochberg FDR correction. When you run one test per genus, the probability of finding at least one false positive by chance approaches certainty. Correction is not optional — it is the difference between a finding and a noise artifact.

The compositionality problem

Microbiome data is compositional. You are measuring relative proportions, not absolute abundances. If Bacteroides increases in relative abundance, every other genus mathematically decreases — even if absolute cell counts are unchanged. Standard non-parametric tests were not designed for this constraint and can generate spurious associations when compositional effects are strong.

Wilcoxon/Kruskal-Wallis on relative abundances is standard practice in published microbiome research and appropriate for exploratory analysis. For studies where compositionality is the central concern, ALDEx2 or ANCOM-BC provide more rigorous alternatives. Use these results as hypothesis-generating — be explicit about this in your Methods section.

Microbiome metadata correlation in R: Spearman with CLR transformation and FDR correction

The final layer of analysis asks whether specific taxa correlate with continuous metadata variables — clinical measurements, days since treatment, age, BMI, or any numeric variable in your sample data. This is where microbiome composition meets the biology you actually care about.

The approach: CLR (Centered Log-Ratio) transformation of genus-level counts, then Spearman rank correlation between each taxon and each metadata variable. Spearman is the right choice — it makes no distributional assumptions, handles the zero-inflated nature of microbiome data, and is robust to outliers. Global FDR correction across all taxon × variable pairs controls the false discovery rate for the entire correlation matrix simultaneously.

Blue cells: as the metadata variable increases, genus abundance decreases. Red cells: positive correlation. Intensity reflects the strength of association. Blank cells (grey background) did not survive FDR correction — they are excluded rather than shown as weak or ambiguous signals. This is a deliberate design choice: a correlation heatmap that shows everything is a map of noise.

Identifying consensus key taxa

After differential abundance testing and correlation analysis, the workflow synthesizes results into a consensus view: which genera appear in both analyses, and which show the strongest combined signal? This is not a formal statistical step — it is a guided look at which taxa are most robustly associated with your study variables across multiple lines of evidence.

16S microbiome analysis R script: complete downstream workflow for QIIME2 data

Every analysis in this guide — composition, alpha diversity, beta diversity with PERMANOVA, differential abundance with FDR, correlation heatmap, and the consensus synthesis — is implemented in a single structured R script organized into 13 labeled sections.

You open the script in RStudio. You set three variables in the configuration block — your folder path, your sample ID column, your grouping variable. You run one section at a time, check the output, and move on. The script handles everything else: package installation, data import from QIIME2 or CSV, filtering, rarefaction, all analyses, all figures, and two text files with a complete Methods section and software citations ready for journal submission.

Methods reference

McMurdie & Holmes (2013) phyloseq, PLoS ONE 8(4):e61217 · Oksanen et al. (2024) vegan · Kassambara (2023) rstatix · Barnett et al. (2021) microViz, JOSS 6(63):3201 · Revelle (2024) psych · Wickham et al. (2016) ggplot2 · Gu et al. (2016) ComplexHeatmap, Bioinformatics 32:2847 · Caporaso et al. (2011) Moving Pictures, Genome Biology 12:R50.