Evomics Workshop part 2

Microbiome analysis with microViz

Welcome to Part 2 👋

• Bar charts and diversity were relatively straightforward topics, conceptually.
• Let’s look at some more things we can do with microbiome composition data!

Topics for part 2:

1. Dissimilarity measures
2. Ordination: PCoA & PCA (with interactive data exploration)
3. Differential abundance testing

Further resources

Refer to the microViz website to see help pages for every function (as well as further tutorials).

---

Setup ⚙️

Load the R packages we will be using:

library(tidyverse)
library(seriation)
library(phyloseq)
library(microViz)
library(shiny)

Dissimilarity 💩 ↔︎ 💩 ?

• Ecosystem Diversity is sometimes referred to as “alpha” diversity or “within-sample” diversity
• Now we’re going to look at Dissimilarity between ecosystems
• Sometimes this is confusingly referred to as “beta diversity” analysis, or “between-sample” diversity

data("shao19") # this example data is built in to microViz

• For this part we’re going to use an infant gut microbiome shotgun metagenomics dataset.

shao19

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 819 taxa and 1644 samples ]
sample_data() Sample Data:       [ 1644 samples by 11 sample variables ]
tax_table()   Taxonomy Table:    [ 819 taxa by 6 taxonomic ranks ]
phy_tree()    Phylogenetic Tree: [ 819 tips and 818 internal nodes ]

• Do you remember how to examine a phyloseq object?
• Look at the rank names, sample data variables etc.

#
#
#
#

?shao19 # in the console, for info on the data source

---

A foreword on filtering

• You should check if any of your samples have a surprisingly low total number of (classified) reads.
• This can suggest that something went wrong in the lab (or during sample collection)
• The data from these samples might be unreliable.
• You might already do a check for total reads and remove poor quality samples during the fastq file processing.

shao19 %>%
  ps_mutate(reads = sample_sums(shao19)) %>%
  samdat_tbl() %>%
  ggplot(aes(x = reads)) +
  geom_freqpoly(bins = 500) +
  geom_rug(alpha = 0.5) +
  scale_x_log10(labels = scales::label_number()) +
  labs(x = "Number of classified reads", y = NULL) +
  theme_bw()

How many is enough? There is no easy answer!

These samples have great depth. There are a few with much less reads than the rest, and a few with under a million. You might consider dropping the samples with under a million reads, to see if it affects your results, but in this case we won’t.

But 100,000 is still a lot, compared to what older sequencing machines produced: 1000 reads might have been considered very good. So look at the distribution for your data, in case there are obvious outliers, and look at recent papers using a similar sequencing technique for what kind of threshold they used.

There might also be relevant information for the type of sequencer you used on e.g. Illumina website. e.g. for this type of sequencing Illumina suggests you should expect at least a million reads (and this is good for RNA seq analyses). https://support.illumina.com/bulletins/2017/04/considerations-for-rna-seq-read-length-and-coverage-.html

---

Dissimilarity measures

• What are we doing?
• Calculating the dissimilarities between pairs of microbiome samples
• We talked about these commonly-used dissimilarity measures in the lecture.
  • Binary Jaccard - presence-absence
  • Bray-Curtis - abundance-weighted
  • UniFrac distances (unweighted / weighted / generalised)
• To simplify and speed up the analyses, we’re going to take a smaller part of the dataset.
• We’ll only look at the 300 infant fecal samples from 4 days of age.

infants <- shao19 %>% ps_filter(family_role == "child", infant_age == 4)

• We’re going to filter out rare taxa quite strictly, for similar reasons
• But we won’t overwrite our smaller dataset: we’ll do the filtering per analysis
• You will find out more about the how and why of taxa filtering later

infants %>%
  tax_filter(min_prevalence = 2.5 / 100) %>%
  tax_agg(rank = "genus") %>%
  tax_transform("binary") %>% # converts counts to absence/presence: 0/1
  dist_calc(dist = "jaccard")

Proportional min_prevalence given: 0.025 --> min 8/306 samples.

psExtra object - a phyloseq object with extra slots:

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 35 taxa and 306 samples ]
sample_data() Sample Data:       [ 306 samples by 11 sample variables ]
tax_table()   Taxonomy Table:    [ 35 taxa by 5 taxonomic ranks ]

otu_get(counts = TRUE)       [ 35 taxa and 306 samples ]

psExtra info:
tax_agg = "genus" tax_trans = "binary" 

jaccard distance matrix of size 306 
0.6666667 0.7333333 0.9375 0.8125 0.6428571 ...

“Binary” Jaccard 010101

• Remember to run a “binary” transform on your data before computing “jaccard” distance.
• There is a quantitative form of the Jaccard distance, which is the default behaviour!
• But the qualitative (presence/absence) version is mostly used in microbial ecology.
• If you want an abundance-weighted ecological dissimilarity, use Bray-Curtis!

• We now have our pairwise dissimilarities! 🎉
• A distance matrix is attached as an extra part on the original phyloseq object

Dissimilarity or distance?

• These terms are often used interchangeably
• You will find dissimilarities in a distance matrix
• But if you want to be pedantic a true “distance metric” d, must satisfy 3 properties:
  1. Identity of indiscernibles: For any samples x and y, d(x, y) = 0 if and only if x = y
  2. Symmetry: For any samples x and y, d(x, y) = d(y, x)
  3. Triangle inequality: For any samples x, y, and z, d(x, z) ≤ d(x, y) + d(y, z).
• 3. can be interpreted as: “the direct path between two points must be at least as short as any detour”
• This is not true for e.g. Bray-Curtis, but in practice it is very rarely problematic.

• The object is now class “psExtra” (created by microViz)
• A psExtra also stores info about the aggregation and transformations you performed

distances <- infants %>%
  tax_filter(min_prevalence = 2.5 / 100, verbose = FALSE) %>%
  tax_agg(rank = "genus") %>%
  tax_transform("binary") %>%
  dist_calc(dist = "jaccard") %>%
  dist_get()

You can extract the distance matrix with dist_get.

as.matrix(distances)[1:5, 1:5]

            B01089_ba_4 B01190_ba_4 B01194_ba_4 B01196_ba_4 B01235_ba_4
B01089_ba_4   0.0000000   0.6666667   0.7333333   0.9375000   0.8125000
B01190_ba_4   0.6666667   0.0000000   0.7500000   0.9166667   0.8461538
B01194_ba_4   0.7333333   0.7500000   0.0000000   0.4615385   0.3076923
B01196_ba_4   0.9375000   0.9166667   0.4615385   0.0000000   0.4615385
B01235_ba_4   0.8125000   0.8461538   0.3076923   0.4615385   0.0000000

The Binary Jaccard dissimilarities range between 0 (identical) and 1 (no shared genera).

range(distances)

[1] 0 1

---

Ordination

• What can we do with these dissimilarities? 🤔
• We can make an ordination! 💡
• Conceptually, ordination refers to a process of ordering things (in our case: samples).
• Similar samples are placed closer to each other, and dissimilar samples are placed further away.

PCoA

Principal Co-ordinates Analysis is one kind of ordination.

• PCoA takes a distance matrix and finds new dimensions (a co-ordinate system, if you like)
• The new dimensions are created with the aim to preserve the original distances between samples
• And to capture the majority of this distance information in the first dimensions
• This makes it easier to visualize the patterns in your data, in 2D or 3D plots 👀

For more info, see “GUSTAME”

• There is helpful info about ordination methods including PCoA on the GUide to STatistical Analysis in Microbial Ecology website (GUSTA ME). https://sites.google.com/site/mb3gustame/dissimilarity-based-methods/principal-coordinates-analysis
• This website covers a lot of other topics too, which may be interesting for you to read at a later date if you’ll work on microbiome analysis.

infants %>%
  tax_filter(min_prevalence = 2.5 / 100, verbose = FALSE) %>%
  tax_transform(trans = "identity", rank = "genus") %>%
  dist_calc(dist = "bray") %>%
  ord_calc(method = "PCoA") %>%
  ord_plot(alpha = 0.6, size = 2) +
  theme_classic(12) +
  coord_fixed(0.7)

To get a little insight into what has happened here, we can colour each sample according to its dominant (most abundant) genus.

infants %>%
  tax_filter(min_prevalence = 2.5 / 100, verbose = FALSE) %>%
  ps_calc_dominant(rank = "genus", none = "Mixed", other = "Other") %>%
  tax_transform(trans = "identity", rank = "genus") %>%
  dist_calc(dist = "bray") %>%
  ord_calc(method = "PCoA") %>%
  ord_plot(color = "dominant_genus", alpha = 0.6, size = 2) +
  scale_color_brewer(name = "Dominant Genus", palette = "Dark2") +
  theme_classic(12) +
  coord_fixed(0.7)

Interactive ordination!

microViz provides a Shiny app ord_explore to interactively create and explore PCoA plots and other ordinations. See the code below to get started. But read the instructions first.

Instructions: a few things to try out

• Colour the samples using the variables in the sample data
• Select a few samples to view their composition on barplots!
• Change some ordination options:
  • Different rank of taxonomic aggregation
  • Different distances we’ve discussed
• Copy the automatically generated code
  • Exit the app (press escape or click red button in R console!)
  • Paste and run the code to recreate the ordination plot
  • Customise the plot: change colour scheme, title, etc.
• Launch the app again with a different subset of the data
  • Practice using ps_filter etc.
  • e.g. plot the data of the mothers’ gut microbiomes!
  • compute and colour points by an alpha diversity measure?

Beware: some important notes on interactive analysis

• Unblock pop-ups: To allow the interactive analysis window to open in your browser, you may need to unblock pop-ups for your AMI instance address (check for messages about this after running the ord_explore command)
• Slow UniFrac: UniFrac distances can be quite slow (over a minute) to calculate!
  • Filter to fewer samples and fewer taxa to speed it up (Before launching the app)
• Many other distances: There are many distances available, feel free to try out ones we haven’t talked about
  • BUT:
    • You shouldn’t use a distance that you don’t understand in your actual work, even if the plot looks nice! 😉
    • A few of the distances might not work…
      • They are mostly implemented in the package vegan and I haven’t tested them all
      • Errors will appear in the RStudio R console
      • You can report to me any distances that don’t work (if you’re feeling helpful! 😇)
• Other ordination methods: There are other ordination methods available in ord_explore
  • Try out PCA, principal components analysis, which does NOT use distances
  • We will not cover constrained and conditioned ordinations
  • If you are interested in e.g. RDA, you can look this up later
  • See the Guide to Statistical Analysis in Microbial Ecology

# fire up the shiny app
# run these lines in your console (don't keep in script/notebook)
infants %>%
  tax_filter(min_prevalence = 2.5 / 100, verbose = FALSE) %>%
  # calculate new sample variables with dominant taxon (optional)
  ps_calc_dominant(rank = "genus", none = "Mixed", other = "Other") %>%
  # launch a Shiny app in your web browser!
  ord_explore()

# Again, with different options

# Run these lines in your console
shao19 %>%
  ps_filter(family_role == "mother") %>%
  tax_filter(min_prevalence = 2.5 / 100, verbose = FALSE) %>%
  # calculate a few sample variables for interest (optional)
  ps_calc_dominant(rank = "genus", none = "Mixed", other = "Other") %>%
  ps_calc_diversity(rank = "genus", index = "shannon") %>%
  ps_calc_richness(rank = "genus", index = "observed") %>%
  # launch a Shiny app in your web browser!
  ord_explore()

---

PERMANOVA

Permutational multivariate analysis of variance.

• ANOVA - analysis of variance (statistical modelling approach)
• Multivariate - more than one dependent variable (multiple taxa!)
• Permutational - statistical significance estimates obtained by shuffling the data many times

For more details on PERMANOVA

• See this excellent book chapter by Marti Anderson on PERMANOVA: https://onlinelibrary.wiley.com/doi/full/10.1002/9781118445112.stat07841
• Sometimes PERMANOVA is called NP-MANOVA (non-parametric MANOVA)
• e.g. on the GUide to STatistical Analysis in Microbial Ecology website.

• TLDR: Are those groups on the PCoA actually different??

infants %>%
  tax_filter(min_prevalence = 2.5 / 100, verbose = FALSE) %>%
  tax_agg(rank = "genus") %>%
  dist_calc(dist = "bray") %>%
  ord_calc(method = "PCoA") %>%
  ord_plot(alpha = 0.6, size = 2, color = "birth_mode") +
  theme_classic(12) +
  coord_fixed(0.7) +
  stat_ellipse(aes(color = birth_mode)) +
  scale_color_brewer(palette = "Set1")

infants %>%
  tax_filter(min_prevalence = 2.5 / 100, verbose = FALSE) %>%
  tax_agg(rank = "genus") %>%
  dist_calc(dist = "bray") %>%
  dist_permanova(variables = "birth_mode", n_perms = 99, seed = 123) %>%
  perm_get()

2024-12-17 12:06:59.922681 - Starting PERMANOVA with 99 perms with 1 processes

2024-12-17 12:07:00.255468 - Finished PERMANOVA

Permutation test for adonis under reduced model
Marginal effects of terms
Permutation: free
Number of permutations: 99

vegan::adonis2(formula = formula, data = metadata, permutations = n_perms, by = by, parallel = parall)
            Df SumOfSqs      R2      F Pr(>F)   
birth_mode   1   13.790 0.12366 42.898   0.01 **
Residual   304   97.727 0.87634                 
Total      305  111.518 1.00000                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Use more permutations for a more reliable p.value in your real work (slower)
# Set a random seed number for reproducibility of this stochastic method

You can see from the model output that the p value, Pr(>F) is below 0.05. So there is good statistical evidence that the bacterial gut microbiota composition of C-section delivered infants has a different composition than vaginally delivered infants at 4 days of age.

Reporting PCoA and PERMANOVA methods

• Your methodological choices matter, you should report what you did:
  • any relevant rare taxon filtering thresholds
  • the taxonomic rank of aggregation
  • the dissimilarity measure used to compute the pairwise distances

It’s probably a good idea to decide on a couple of appropriate distance measures up front for these tests, and report both (at least in supplementary material), as the choice of distance measure can affect results and conclusions!

Covariate-adjusted PERMANOVA

You can also adjust for covariates in PERMANOVA, and often should, depending on your study design. Let’s fit a more complex model, adjusting for infant sex, birth weight, and the total number of assigned reads.

infants %>%
  tax_filter(min_prevalence = 2.5 / 100, verbose = FALSE) %>%
  tax_agg(rank = "genus") %>%
  dist_calc(dist = "bray") %>%
  dist_permanova(
    variables = c("birth_mode", "sex", "birth_weight", "number_reads"),
    n_perms = 99, seed = 111
  ) %>%
  perm_get()

Dropping samples with missings: 15

2024-12-17 12:07:00.391675 - Starting PERMANOVA with 99 perms with 1 processes

2024-12-17 12:07:01.301649 - Finished PERMANOVA

Permutation test for adonis under reduced model
Marginal effects of terms
Permutation: free
Number of permutations: 99

vegan::adonis2(formula = formula, data = metadata, permutations = n_perms, by = by, parallel = parall)
              Df SumOfSqs      R2       F Pr(>F)   
birth_mode     1   10.794 0.10163 34.8045   0.01 **
sex            1    0.280 0.00264  0.9031   0.43   
birth_weight   1    0.565 0.00532  1.8215   0.06 . 
number_reads   1    2.873 0.02705  9.2656   0.01 **
Residual     286   88.696 0.83509                  
Total        290  106.211 1.00000                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Use more permutations for a more reliable p.value in your real work (slower)
# Set a random seed number for reproducibility of this stochastic method

PCA

• Principal Components Analysis.
• For practical purposes, PCA is quite similar to Principal Co-ordinates Analysis.
• In fact, PCA produces equivalent results to PCoA with Euclidean distances.

“Help, what are Euclidean distances?”

• Euclidean distances are essentially a generalization of Pythagoras’ theorem to more dimensions.
• In our data every taxon is a feature, a dimension, on which we calculate Euclidean distances.

Pythagoras’ theorem:

\[c = \sqrt{a^2 + b^2}\]

Euclidean distance:

\[d\left(p, q\right) = \sqrt{\sum _{i=1}^{n_{taxa}} \left( p_{i}-q_{i}\right)^2 }\]

• Distance \(d\) between samples \(p\) and \(q\), with \(n\) taxa.

Euclidean PCoA

infants %>%
  tax_agg(rank = "genus") %>%
  dist_calc(dist = "euclidean") %>%
  ord_calc(method = "PCoA") %>%
  ord_plot(alpha = 0.6, size = 2) +
  geom_rug(alpha = 0.1)

PCA on counts

infants %>%
  tax_agg(rank = "genus") %>%
  ord_calc(method = "PCA") %>%
  ord_plot(alpha = 0.6, size = 2) +
  geom_rug(alpha = 0.1)

Problems with PCA (or PCoA with Euclidean Distances) on microbiome data

• These plots look weird! most samples bunch in the middle, with spindly projections..
• Sensitive to sparsity (double-zero problem) –> filter rare taxa
• Excessive emphasis on high-abundance taxa –> log transform features first

---

Log transformations, and CLR

• First let’s look at the abundance again, this time with heatmaps.
• Each column is a sample (from an infant), and each row is a taxon.

infants %>%
  tax_sort(by = sum, at = "genus", trans = "compositional", tree_warn = FALSE) %>%
  tax_transform(trans = "compositional", rank = "genus") %>%
  comp_heatmap(samples = 1:20, taxa = 1:20, name = "Proportions", tax_seriation = "Identity")

• Even though we have picked the top 20 most abundant genera, there are still a lot of zeros
• Problem: We need to deal with the zeros, because log(0) is undefined.
• Solution: add a small amount to every value (or just every zero), before applying the log transformation.
• This small value is often called a pseudo-count.

What value should we use for the pseudo-count?

• One easy option is to just add a count of 1
• Another popular option is to add half of the smallest observed real value (from across the whole dataset)
• In general, for zero replacement, keep it simple and record your approach

Centered Log Ratio transformation:

Remember, Microbiome Datasets Are Compositional: And This Is Not Optional.

More details on the “CoDa” problem:

• The sequencing data gives us relative abundances, not absolute abundances.
• The total number of reads sequenced per sample is an arbitrary total.

This leads to two main types of problem:

• Interpretation caveats: see differential abundance section later
• Statistical issues: taxon abundances are not independent, and may appear negatively correlated
• These issues are worse with simpler ecosystems

Example: If one taxon blooms, the relative abundance of all other taxa will appear to decrease, even if they did not.

There is the same problem in theory with RNAseq data, but I suspect it is less bothersome because there are many more competing “species” of RNA transcript than there are bacterial species in even a very complex microbiome. The centered-log-ratio transformation (along with some other similar ratio transformations) are claimed to help with the statistical issues by transforming the abundances from the simplex to the real space.

TL;DR - the CLR transformation is useful for compositional microbiome data.

• Practically, the CLR transformation involves finding the geometric mean of each sample
• Then dividing abundance of each taxon in that sample by this geometric mean
• Finally, you take the natural log of these ratios

infants %>%
  tax_sort(by = sum, at = "genus", trans = "compositional", tree_warn = FALSE) %>%
  tax_agg(rank = "genus") %>%
  tax_transform(trans = "clr", zero_replace = "halfmin", chain = TRUE) %>%
  comp_heatmap(
    samples = 1:20, taxa = 1:20, grid_lwd = 2, name = "CLR",
    colors = heat_palette(sym = TRUE),
    tax_seriation = "Identity"
  )

PCA on CLR-transformed taxa

infants %>%
  tax_filter(min_prevalence = 2.5 / 100, verbose = FALSE) %>%
  tax_transform(rank = "genus", trans = "clr", zero_replace = "halfmin") %>%
  ord_calc(method = "PCA") %>%
  ord_plot(alpha = 0.6, size = 2, color = "birth_mode") +
  theme_classic(12) +
  coord_fixed(0.7)

• After the CLR transformation, the plot looks better
• We can see a pattern where the gut microbiomes of infants cluster by birth mode

So why is PCA interesting for us?

• Principal components are built directly from a (linear) combination of the original features.

• That means we know how much each taxon contributes to each PC dimension

• We can plot this information (loadings) as arrows, alongside the sample points

pca <- infants %>%
  tax_filter(min_prevalence = 2.5 / 100, verbose = FALSE) %>%
  tax_transform(rank = "genus", trans = "clr", zero_replace = "halfmin") %>%
  ord_calc(method = "PCA") %>%
  ord_plot(
    alpha = 0.6, size = 2, color = "birth_mode",
    plot_taxa = 1:6, tax_vec_length = 0.275,
    tax_lab_style = tax_lab_style(
      type = "text", max_angle = 90, aspect_ratio = 0.7,
      size = 3, fontface = "bold"
    ),
  ) +
  theme_classic(12) +
  coord_fixed(0.7, clip = "off")
pca

Interestingly, samples on the right of the plot (which tend to be vaginally-delivered infants) seem to have relatively more Bifidobacterium, Bacteroides and Escherichia, whilst the C-section born infants have relatively more Veillonella.

Wait, how to interpret these taxa loadings?

Cautiously

• There are caveats and nuance to the interpretation of these plots, which are called PCA bi-plots
• You can read more here: https://sites.google.com/site/mb3gustame/indirect-gradient-analysis/principal-components-analysis

In general:

The relative length and direction of an arrow indicates how much that taxon contributes to the variation on each visible PC axis, e.g. Variation in Enterococcus abundance contributes quite a lot to variation along the PC2 axis.

This allows you to infer that samples positioned at the bottom of the plot will tend to have higher relative abundance of Enterococcus than samples at the top of the plot.

(Side note, Phocaeicola were considered part of Bacteroides until this year!)

Fancy circular bar charts?

We can make another kind of barplot, using the PCA information to order our samples in a circular layout.

iris <- infants %>%
  tax_filter(min_prevalence = 2.5 / 100, verbose = FALSE) %>%
  tax_transform(rank = "genus", trans = "clr", zero_replace = "halfmin") %>%
  ord_calc(method = "PCA") %>%
  ord_plot_iris(
    tax_level = "genus", n_taxa = 12, other = "Other",
    anno_colour = "birth_mode",
    anno_colour_style = list(alpha = 0.6, size = 0.6, show.legend = FALSE)
  )

patchwork::wrap_plots(pca, iris, nrow = 1, guides = "collect")