Beta Diversity, UniFrac, and Ordination

Why beta diversity

Alpha diversity asks how rich or even a single sample is. Beta diversity asks a different question: how different are two communities from each other. You quantify that with a pairwise dissimilarity, which most metrics scale from 0 (identical composition) to 1 (nothing in common), computed for every pair of samples. Those pairwise values stack into a square distance matrix, which is the object that all downstream ordination and statistics consume.

Four workhorse metrics

Bray-Curtis uses relative abundances and is the classic count-based dissimilarity. Jaccard is its presence/absence cousin, comparing only which features are shared. UniFrac goes further by using a phylogenetic tree: it measures the fraction of branch length unique to each sample, so evolutionary relatedness counts. Unweighted UniFrac looks at presence/absence of branches, while weighted UniFrac scales each branch by abundance. Because UniFrac needs a tree, you must build one from your representative sequences before you can compute it, using qiime phylogeny or FastTree.

# Align, mask, build, and root a tree in one QIIME 2 pipeline step
qiime phylogeny align-to-tree-mafft-fasttree \
  --i-sequences rep-seqs.qza \
  --o-alignment aligned-rep-seqs.qza \
  --o-masked-alignment masked-aligned-rep-seqs.qza \
  --o-tree unrooted-tree.qza \
  --o-rooted-tree rooted-tree.qza

# Standalone equivalent outside QIIME 2: FastTree on a masked alignment
# (midpoint-root the result before computing UniFrac)
# FastTree -nt -gtr masked-alignment.fasta > unrooted-tree.nwk

# Rarefies every sample to --p-sampling-depth, then computes the standard
# alpha-diversity vectors plus all four beta-diversity distance matrices,
# their PCoA ordinations, and Emperor 3D plots in one command.
# Pick the depth from your feature-table summary (table.qzv).
qiime diversity core-metrics-phylogenetic \
  --i-phylogeny rooted-tree.qza \
  --i-table table.qza \
  --p-sampling-depth 10000 \
  --m-metadata-file metadata.tsv \
  --output-dir core-metrics-results

Saved FeatureTable[Frequency] to: core-metrics-results/rarefied_table.qza
Saved SampleData[AlphaDiversity] to: core-metrics-results/faith_pd_vector.qza
Saved SampleData[AlphaDiversity] to: core-metrics-results/observed_features_vector.qza
Saved SampleData[AlphaDiversity] to: core-metrics-results/shannon_vector.qza
Saved SampleData[AlphaDiversity] to: core-metrics-results/evenness_vector.qza
Saved DistanceMatrix to: core-metrics-results/unweighted_unifrac_distance_matrix.qza
Saved DistanceMatrix to: core-metrics-results/weighted_unifrac_distance_matrix.qza
Saved DistanceMatrix to: core-metrics-results/jaccard_distance_matrix.qza
Saved DistanceMatrix to: core-metrics-results/bray_curtis_distance_matrix.qza
Saved PCoAResults to: core-metrics-results/unweighted_unifrac_pcoa_results.qza
Saved PCoAResults to: core-metrics-results/weighted_unifrac_pcoa_results.qza
Saved PCoAResults to: core-metrics-results/jaccard_pcoa_results.qza
Saved PCoAResults to: core-metrics-results/bray_curtis_pcoa_results.qza
Saved Visualization to: core-metrics-results/unweighted_unifrac_emperor.qzv
Saved Visualization to: core-metrics-results/weighted_unifrac_emperor.qzv
Saved Visualization to: core-metrics-results/jaccard_emperor.qzv
Saved Visualization to: core-metrics-results/bray_curtis_emperor.qzv

library(phyloseq)
library(vegan)
library(ggplot2)

# ps: a phyloseq object holding counts, sample_data, and the rooted tree.
# Reproduce Emperor with ggplot via a PCoA ordination:
ord <- ordinate(ps, method = "PCoA", distance = "wunifrac")
plot_ordination(ps, ord, color = "body_site") +
  stat_ellipse() + theme_bw()

# PERMANOVA: do group centroids differ?
wuf  <- phyloseq::distance(ps, method = "wunifrac")
meta <- as(sample_data(ps), "data.frame")
set.seed(42)
adonis2(wuf ~ body_site, data = meta, permutations = 999)

# Dispersion: is within-group spread unequal? (a PERMANOVA assumption)
bd <- betadisper(wuf, meta$body_site)
permutest(bd, permutations = 999)

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

adonis2(formula = wuf ~ body_site, data = meta, permutations = 999)
          Df SumOfSqs      R2      F Pr(>F)    
body_site  3   6.9346 0.3572 8.5209  0.001 ***
Residual  46  12.4788 0.6428                  
Total     49  19.4134 1.0000                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Metric choice is a biological choice, not a formality. Jaccard and unweighted UniFrac only ask which taxa are present, so they amplify rare and low-abundance lineages and are sensitive to sequencing depth and contamination. Bray-Curtis and weighted UniFrac weight by abundance, so they track shifts in the dominant members of a community. UniFrac additionally credits how related the differing taxa are on the tree, so two samples sharing no exact features but many close relatives look more similar than they would under Bray-Curtis. A separation that is significant under unweighted UniFrac but vanishes under weighted UniFrac usually points to turnover among rare lineages rather than a reshuffling of abundant taxa. Always confirm a significant adonis2 result with betadisper, because PERMANOVA can be driven by unequal within-group spread rather than a true difference in centroid location.