Microbiome data analysis with microViz

David Barnett

Part 3: PCA & Differential Abundance

Datasets

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 ]

mice

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

ibd

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 36349 taxa and 91 samples ]
sample_data() Sample Data:       [ 91 samples by 15 sample variables ]
tax_table()   Taxonomy Table:    [ 36349 taxa by 6 taxonomic ranks ]

We’ll primarily use (parts of) the shao19 data for demonstration. But during the exercises try exploring the other datasets too, if you want more practice.

More Dissimilarity

• Ecological & Phylogenetic dissimilarities / distances:
  • Binary Jaccard
  • Bray Curtis
  • UniFrac (unweighted/weighted/generalized)
  • Others…

• Euclidean distances
  • Widely used in data science
  • Useful for microbiome abundance data??

Euclidean distances

Pythagoras: \(c = \sqrt{a^2 + b^2}\)

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

Euclidean distance between samples p and q

dist(pq, method = "euclidean")

  p
q 5

• What are Euclidean distances? Here is a simple 2-dimensional example - two taxa.
• 2 samples: Pythagoras you probably remember from school \(c = \sqrt{a^2 + b^2}\) Distance in one dimension, square it, distance in another dimension, square it, add them together, and square root of the result.
• Euclidean distance is essentially just a generalisation to more dimensions, where each dimension is the counts of one taxon.
• Test it, indeed you get 5

Naïve Euclidean PCoA

• Sensitive to sparsity
• Excessive emphasis on high-abundance taxa

The PCoA looks weird!

We shouldn’t just use Euclidean distances on the raw data, why?

Sensitive to sparsity (double-zero problem). Excessive emphasis on high-abundance taxa. The PCoA looks weird! Most samples bunched in the middle with spindly projections..

Taxon Abundance transformations!

• Already used:
  • Compositional (proportions)
  • Binary (Presence or Absence)

• Log transformations
  • Log₂, Log₁₀, natural log?

• Centered Log Ratio - CLR

We already did two transformations with tax_transform(): binary (for Binary Jaccard distances) and compositional (for barplots).

Now we need log transformations, and the centered-log-ratio, CLR, transformation.

Transformations - none

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.

Transformations - none - all data!

Too much to show all the taxa (and samples, for this example)

Transformations - Compositions

Doesn’t actually help us with the problem of excess emphasis on high-abundance features

Transformations - Binary

Binary transformation is not really useful here as it discards all abundance info, but nice to show that even though we have picked the top 20 most abundant genera, there are still a lot of zeros. We need to deal with the zeros, because log(0) is undefined. The solution is to 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.

Transformations - log10(x+1)

What value should we use for the pseudo-count? One option is to just add 1.

Transformations - log10(x+“halfmin”)

Another popular option is to add half of the smallest observed real value (from across the whole dataset).

People fight about the best zero replacement strategy and some much more complicated methods exist: including generative statistical models. In practice: keep it simple and record your approach!

Log10-transformed Euclidean PCoA

Transformations - Centered Log Ratio

So I mentioned another transformation. The Centered-log-ratio transformation

The centered log-ratio (clr) transformation uses the geometric mean of the sample vector as the reference

Centered Log Ratio tranformation

• What?
  • Calculate geometric mean for each sample
  • Divide each taxon by the geometric mean of each sample
  • Take the natural log
• Why?
  • The compositional data problem

Compositional data problem

Gloor 2017 - Microbiome Datasets are Compositional: and this is not optional.

The total number of reads number of sequence counts is not a good measure of microbial biomass. It is an arbitrary sum, more reads only gets you higher precision and more reliable estimates of the proportions. I mentioned this briefly when we started making bar charts.

In this nice diagram from a paper by Prof Gregory Gloor, we see what this means for data interpretation. First part of B shows true abundance of 2 bacterial taxa in a simple ecosystem. Using sample 1 as the reference, if we compare this to sample 2, the true abundance of taxon B has decreased, and A is the same. In sample 3, A has increased, whilst B stays the same.

But from sequencing data we only know the relative abundance. So samples 2 and 3 look identical! That is the interpretation problem. But this also causes a statistical problem.

Gloor 2017 is a good read, including the extensive educational supplementary material!.

Constant sum, 1000 in this three taxon toy example. Causes taxa abundances to be appear correlated, on average. The counts exist in what is called the Aitchison simplex space. With 3 taxa, this space can be visualised on a ternary diagram, with 3 interconnected axes. We need to escape the simplex, back to unconstrained real space, this is what the CLR does.

Another slide from Prof Gloor. X is a sample vector where x1, x2 etc are the counts of each taxon in that sample.

After the CLR transformation, the abundance of each taxon is expressed relative to the other taxa in the sample. This transformation takes us back to real geometric space.

Confused? Don’t worry, if you’re interested check out that paper by Gloor. If not just know that you’ll often see microbiome data transformed in this way, and it makes the PCoA plots look better.

Aitchison PCoA (CLR + Euclidean)

PCA - where’s the “o” gone?

Principal Components Analysis (not Co-ordinates)

• Similar aim to PCoA (and name!)
• But doesn’t use distances!
• Rotation of original dimensions (taxa)
• Euclidean PCoA equivalent results to PCA

Principal Components Analysis - PCA

Looks exactly the same as the PCoA plot with CLR euclidean distances.

PCA provides (taxa) loadings!

How to interpret the taxa loading vectors? Cautiously. There are plenty of caveats and nuance to the interpretation of these plots, which are called PCA bi-plots, and you can read more about here.

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 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.

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.

PCA-sorted bar chart (iris plot)

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

Differential Abundance

• like RNA Seq. Differential Expression
• but for microbial taxon abundances

So far we have looked at patterns in microbiota data using multiple taxa at once. You should start microbiome analysis with multivariate exploration and testing, with PCA, PERMANOVA and alpha diversity, to establish if there are any interesting overall microbiome differences in your microbiota data. If you find evidence from multivariate analyses, you don’t need to stop there. It is of course possible and interesting to also look at one taxon at a time.

Differential Abundance - aim

• Is experimental condition X associated with relative abundance of bacterial genus Y?
• Is environmental exposure X associated with …?
• Is host characteristic X associated with… ?

Differential Abundance - methods

• Simple methods: Wilcoxon tests, Spearman corr., Linear Reg.
• CoDa methods: ALDEx2, ANCOM-II, ANCOM-BC
• RNA seq. methods: DESeq2, limma, edgeR
• Other: LEfSe, Beta-binomial reg., ZINB reg., ZIGDM reg., …!

Differential Abundance - which rank?

As well as choosing the statistical method, you also need to choose a rank at which to aggregate your taxa before differential abundance testing.

Differential Abundance with microViz

Keep it simple (-ish)

• Filter out rare taxa strictly (>5% prev or more)

• Linear regression models

  • Log2 transform of taxon proportions
  • (approach borrowed from Maaslin2)

• Model every taxon, at all ranks!

  • Visualize on a cladogram tree

Differential Abundance

Here is an example of visualizing the results of many statistical tests on a cladogram structure. The grey tree with labels on this side is just a key, so you know which taxa are which. Each node represents a taxon, the nodes close to the center are the phyla, and around the edge are the genera, in this case, which all the other ranks in between on the concentric circles (which are probably too faint to see on the screen)

The smaller plots with colours, show the results of many statistical tests. Each covariate in the statistical models gets its own tree. So in this example that’s only birth mode, and the age of the infant when the sample was taken. The colour of each node (each taxon) represents the direction and strength of that taxon’s association with that covariate. So the blue branch up there, representing Bacteroidetes phylum and genera inside that phylum, they’re negatively associated with C-section birth. You can also add markers for statistical significance, which in this case is unadjusted p-values < 0.05, although that should really be FDR corrected p-values!

Your turn

• Link to exercises guidance: https://david-barnett.github.io/evomics-material/exercises/microViz-3-transformations-PCA-DA-exercises

• Try out making heatmaps and PCA biplots, with transformed abundance data.

See Gloor 2017 Microbiome Datasets Are Compositional: And This Is Not Optional

• Try out differential abundance testing and creating taxonomic association tree plots

Nearing 2021 is a great intro to the Differential Abundance topic/problem: https://www.nature.com/articles/s41467-022-28034-z

• This is the last topic - thank you for your attention and effort :)