6.3.1 Principal Coordinates Analysis (PCoA)

We take as our first example, a reproduction of Figure 5 from the “Global Patterns” article(Caporaso 2011). The authors show a 3-dimensional representation of the first three axes of a Principal Coordinates Analysis (PCoA; This is also sometimes referred to as “Multi-Dimensional Scaling”, or “MDS”) performed on the unweighted-UniFrac distance using all of the available sequences (their approach included both 5’ and 3’ sequences). According to the authors, “the first axis [appears to be associated with a] host associated/free living [classification],” and similarly the third axis with “saline/nonsaline environment[s].”

The following reproduces the unweighted UniFrac distance calculation on the full dataset. Note that this calculation can take a long time because of the large number of OTUs. Parallelization is recommended for large datasets, typically if they are as large as Global Patterns, or larger. For details on parallelization, see the details section and examples in the UniFrac() documentation, and also the page dedicated to the topic on the phyloseq-demos site:

http://joey711.github.io/phyloseq-demo/unifrac.html

data(GlobalPatterns)

GPUF <- UniFrac(GlobalPatterns)

Load the pre-computed distance matrix, GPUF

load(system.file("doc", "Unweighted_UniFrac.RData", package="phyloseq"))

Calculate the PCoA on this distance matrix, GPUF.

GloPa.pcoa = ordinate(GlobalPatterns, method="PCoA", distance=GPUF)

Before we look at the results, let’s first investigate how much of the total distance structure we will capture in the first few axes. We can do this graphically with a “scree plot”, an ordered barplot of the relative fraction of the total eigenvalues associated with each axis.

plot_scree(GloPa.pcoa, "Scree plot for Global Patterns, UniFrac/PCoA")

Scree plot of the PCoA used to create Figure 5 from the “Global Patterns” article(Caporaso 2011). The first three axes represent 43% of the total variation in the distances. Interestingly, the fourth axis represents another 9%, and so may warrant exploration as well. A scree plot is an important tool for any ordination method, as the relative importance of axes can vary widely from one dataset to another.

Next, we will reproduce Figure 5 from the “Global Patterns” article(Caporaso 2011), but separating the three axes into 2 plots using plot_ordination().

(p12 <- plot_ordination(GlobalPatterns, GloPa.pcoa, "samples", color="SampleType") + 
  geom_point(size=5) + geom_path() + scale_colour_hue(guide = FALSE) )

(p13 <- plot_ordination(GlobalPatterns, GloPa.pcoa, "samples", axes=c(1, 3),
  color="SampleType") + geom_line() + geom_point(size=5) )

A reproduction in phyloseq / R of the main panel of Figure 5 from the “Global Patterns” article(Caporaso 2011), on two plots. The horizontal axis represents the first axis in the PCoA ordination, while the top and bottom vertical axes represent the second and third axes, respectively. Different points represent different samples within the dataset, and are shaded according to the environment category to which they belong. The color scheme is the default used by ggplot.

6.3.2 non-metric Multi-Dimensional Scaling (NMDS)

We repeat the previous example, but instead using non-metric multidimensional scaling (NMDS) limited to just two dimensions. This approach limits the amount of residual distance “not shown” in the first two (or three) axes, but forefeits some mathematical properties and does not always converge within the specified number of axes.

# (Re)load UniFrac distance matrix and GlobalPatterns data
data(GlobalPatterns)
load(system.file("doc", "Unweighted_UniFrac.RData", package="phyloseq"))
# perform NMDS, set to 2 axes
GP.NMDS <- ordinate(GlobalPatterns, "NMDS", GPUF)

## Run 0 stress 0.1432774 
## Run 1 stress 0.1432625 
## ... New best solution
## ... Procrustes: rmse 0.003499558  max resid 0.01265589 
## Run 2 stress 0.2467057 
## Run 3 stress 0.2234687 
## Run 4 stress 0.1432774 
## ... Procrustes: rmse 0.003483686  max resid 0.01260292 
## Run 5 stress 0.1432774 
## ... Procrustes: rmse 0.003487121  max resid 0.01261713 
## Run 6 stress 0.1856306 
## Run 7 stress 0.1432625 
## ... New best solution
## ... Procrustes: rmse 1.405265e-05  max resid 3.151096e-05 
## ... Similar to previous best
## Run 8 stress 0.2574296 
## Run 9 stress 0.1809282 
## Run 10 stress 0.1809285 
## Run 11 stress 0.1432625 
## ... New best solution
## ... Procrustes: rmse 4.788746e-06  max resid 1.24343e-05 
## ... Similar to previous best
## Run 12 stress 0.1432774 
## ... Procrustes: rmse 0.003470474  max resid 0.01255052 
## Run 13 stress 0.1432625 
## ... Procrustes: rmse 1.152937e-05  max resid 3.31644e-05 
## ... Similar to previous best
## Run 14 stress 0.1432625 
## ... New best solution
## ... Procrustes: rmse 2.462092e-06  max resid 6.050872e-06 
## ... Similar to previous best
## Run 15 stress 0.1856315 
## Run 16 stress 0.1432625 
## ... Procrustes: rmse 9.826556e-06  max resid 3.260556e-05 
## ... Similar to previous best
## Run 17 stress 0.1809284 
## Run 18 stress 0.1856306 
## Run 19 stress 0.1432774 
## ... Procrustes: rmse 0.003490316  max resid 0.01261738 
## Run 20 stress 0.1432774 
## ... Procrustes: rmse 0.003471665  max resid 0.01255283 
## *** Solution reached

(p <- plot_ordination(GlobalPatterns, GP.NMDS, "samples", color="SampleType") +
  geom_line() + geom_point(size=5) )

An example exploratory ordination using non-metric multidimensional scaling (NMDS) on the unweighted UniFrac distance between samples of the “Global Patterns” dataset. Sample points are shaded by environment type, and connected by a line if they belong to the same type. Compare with Figure 5 from the “Global Patterns” article(Caporaso 2011).

The figure nicely shows the relative dissimilarities between microbial communities from different habitats. However, it fails to indicate what was different between the communities. For an ordination method that provides information on the taxa that explain differences between samples (or groups of samples), we use Correspondence Analysis.

6.3.3 Correspondence Analysis (CA)

In the following section we will show continue our exploration of the “GlobalPatterns” dataset using various features of an ordination method called Correspondence Analysis. We give special emphasis to exploratory interpretations using the biplot, because it provides additional information that is not available from PCoA or NMDS.

Let’s start by performing a Correspondence Analysis and investigating the scree plot. Both interestingly and challengingly, the scree plot suggests that the Global Patterns abundance data is quite high-dimensional, with the first two CA axes accounting for not quite 17% of the total (chi-square) variability. Note the absence of a steep decline in eigenvalue fraction as axis number increases. Each additional axis represents only marginally less variability than the previous. It is often more convenient if the first two (or three) axes account for most of the variability.

First, let’s severely subset the number of species for the sake of run-time.\footnote{This is for illustration purposes only, do not repeat unless you are very sure you have a good reason for doing this.

data(GlobalPatterns)
# Take a subset of the GP dataset, top 200 species
topsp <- names(sort(taxa_sums(GlobalPatterns), TRUE)[1:200])
GP    <- prune_taxa(topsp, GlobalPatterns)
# Subset further to top 5 phyla, among the top 200 OTUs.
top5ph <- sort(tapply(taxa_sums(GP), tax_table(GP)[, "Phylum"], sum), decreasing=TRUE)[1:5]
GP     <- subset_taxa(GP, Phylum %in% names(top5ph))
# Re-add human variable to sample data:
sample_data(GP)$human <- factor(human)

Now perform the correspondence analysis.

# Now perform a unconstrained correspondence analysis
gpca  <- ordinate(GP, "CCA")
# Scree plot
plot_scree(gpca, "Scree Plot for Global Patterns Correspondence Analysis")

The correspondence analysis (CA) scree plot of the “Global Patterns” dataset.

Now let’s investigate how the samples behave on the first few CA axes.

(p12 <- plot_ordination(GP, gpca, "samples", color="SampleType") + 
  geom_line() + geom_point(size=5) )

(p34 <- plot_ordination(GP, gpca, "samples", axes=c(3, 4), color="SampleType") + 
  geom_line() + geom_point(size=5) )

First 4 axes of Correspondence Analysis (CA) of the “Global Patterns” dataset (Caporaso 2011).

A clear feature of these plots is that the feces and mock communities cluster tightly together, far away from all other samples on the first axis (CA1). The skin and tongue samples separate similarly, but on the second axis. Taken together, it appears that the first two axes are best explained by the separation of human-associated “environments” from the other non-human environments in the dataset, with a secondary separation of tongue and skin samples from feces.

We will now investigate further this top-level structure of the data, using an additional feature of correspondence analysis that allows us to compare the relative contributions of individual taxa on the same graphical space: the “biplot”. However, because we just displayed the position of samples in the ordination and there are often many thousands of OTUs, we will focus on creating an interpretable plot of the OTUs. For creating graphics that combine the two plots, try the "biplot" or "split" option for type in plot_ordination().

p1  <- plot_ordination(GP, gpca, "species", color="Phylum")
(p1 <- ggplot(p1$data, p1$mapping) + geom_point(size=5, alpha=0.5) + 
  facet_wrap(~Phylum) +  scale_colour_hue(guide = FALSE) )

Species plot of the “Global Patterns” correspondence analysis first two axes, with each phylum on a different panel (“facet”). Only the most abundant 5 phyla among the most abundant 200 taxa (cumulative, all samples) are included. Arbitrary reduction, for computational efficiency of example.

Let’s try drawing the figure again, only this time summarizing the species points as a 2D density estimate, without any individual points.

(p3 <- ggplot(p1$data, p1$mapping) + geom_density2d() +
  facet_wrap(~Phylum) +  scale_colour_hue(guide = FALSE) )

Redrawn figure, which is severely overplotted, as a 2-dimensional species-density topographic map, faceted in the same way.

These figures reveal some useful patterns and interesting outliers, but what if we want a complete summary of how each phylum is represented along each axis? The following code is a way to show this using boxplots, while still avoiding the occlusion problem (points layered on top of each other), and also conveying some useful information about the pattern of taxa that contribute to the separation of human-associated samples from the other sample types. It re-uses the data that was stored in the ggplot2 plot object created in the previous figure, p, and creates a new boxlot graphical summary of the positions of each phylum.

library("reshape2")
# Melt the species-data.frame, DF, to facet each CA axis separately
mdf <- melt(p1$data[, c("CA1", "CA2", "Phylum", "Family", "Genus")], 
            id=c("Phylum", "Family", "Genus") )
# Select some special outliers for labelling
LF <- subset(mdf, variable=="CA2" & value < -1.0)
# build plot: boxplot summaries of each CA-axis, with labels
p <- ggplot(mdf, aes(Phylum, value, color=Phylum)) + 
  geom_boxplot() + 
  facet_wrap(~variable, 2) + 
  scale_colour_hue(guide = FALSE) +
  theme_bw() + 
  theme( axis.text.x = element_text(angle = -90, vjust = 0.5) )
# Add the text label layer, and render ggplot graphic
(p <- p + geom_text(data=subset(LF, !is.na(Family)),
  mapping = aes(Phylum, value+0.1, color=Phylum, label=Family), 
  vjust=0,
  size=2))

Boxplot of taxa (species in this case) of the “Global Patterns” CA first two axes, shaded/separated by phylum. Through this approach it is much easier to see particular species that cluster unusually relative to the rest of their phylum, for example the Bacteroidetes species (Prevotellaceae family) that is positioned most in the negative CA2 direction toward the Tongue/Skin samples.

One way to relate some of the high-level patterns we observed from correspondence analysis is to directly visualize the abundances in an organized, quantitative way, to see if this does in fact support / explain the human/environment microbiome differences. Here is an example using the plot_bar function described in an earlier section.

plot_bar(GP, x="human", fill="SampleType", facet_grid= ~ Phylum)

Phylum-level comparison of relative abundance of taxa in samples that are from human microbiomes (or not).

In this figure we’ve used the threshold parameter to omit all but phyla accounting for the top 90% of phyla in any one sample. Some patterns emerging from this display appear to be: (1) Cyanobacteria, Actinobacteria appear under-represented in human samples; (2) conversely, Firmicutes appear over-represented in human samples; (3) Acidobacteria, Verrucomicrobia appear over-represented in the fecal samples; (4) the only Crenarchaeota were observed in the Mock sample, which is not really a community but a simulated community used as a control. These are not hugely surprising based on previous biological observations from the field, but it is hopefully useful code that can be applied on other datasets that you might have.

6.3.4 Double Principle Coordinate Analysis (DPCoA)

Here is a quick example illustrating the use of Double Principal Coordinate Analysis (DPCoA~\cite{Pavoine2004523), using the using the ordinate() function in phyloseq, as well as the “biplot” option for `plot_ordination(). For a description that includes an applied example using the “enterotype” dataset and comparison with UniFrac/PCoA, see Fukuyama et al~\cite{fukuyama2012com.

# Perform ordination
GP.dpcoa <- ordinate(GP, "DPCoA") 
# Generate default ordination bi-plot
pdpcoa <- 
  plot_ordination(
    physeq = GP, 
    ordination = GP.dpcoa, 
    type="biplot",
    color="SampleType", 
    shape="Phylum")
# Adjust the shape scale manually 
# to make taxa hollow and samples filled (advanced)
shape.fac <- pdpcoa$data$Phylum
man.shapes <- c(19, 21:25)
names(man.shapes) <- c("Samples", levels(shape.fac)[levels(shape.fac)!="Samples"])
p2dpcoa <- pdpcoa + scale_shape_manual(values=man.shapes)
p2dpcoa

A biplot representation of a Double Principal Coordinate Analysis (DPCoA), on a simplified version of the “Global Patterns” dataset with only the most abundant 200 OTUs included.

# Show just Samples or just Taxa
plot_ordination(GP, GP.dpcoa, type="taxa", shape="Phylum")

plot_ordination(GP, GP.dpcoa, type="samples", color="SampleType")

# Split
plot_ordination(GP, GP.dpcoa, type="split",
                color="SampleType", shape="Phylum") +
  ggplot2::scale_colour_discrete()

6.4.1 distance(): Central Distance Function

Many comparisons of microbiome samples, including the graphical model and the PCoA analysis, require a calculation for the relative dissimilarity/distance between one microbial community and another. The phyloseq-package provides a general “wrapper” function for calculating ecological distance matrices between the samples in an experiment.

distance() currently supports 43 method options, as well as user-provided arbitrary methods via an interface to vegan’s designdist() function. Currrently only sample-wise distances are supported (the type argument), but eventually species-wise (OTU-wise) distances will be supported as well. In addition to supporting any of the method options to the three main distance functions of the vegan-package~\cite{veganpkg — including the 14 distances of the vegdist() function and all 24 ecological distances reviewed in Koleff et al. 2003~\cite{Koleff2003JAE coded by betadiver — distance() will eventually support many of the distance methods in the ade4-package~\cite{Dray:2007vs, and others. It also supports the Fast UniFrac distance function (Section~\ref{sec:unifrac) included in phyloseq as native R code, and a wrapper for retreiving the sample-distances from Double Principal Coordinate Analysis (DPCoA).

The function takes a phyloseq-class object and an argument indicating the distance type; and it returns a `dist-class distance matrix.

data(esophagus)
distance(esophagus, "bray") 

##           B         C
## C 0.4061135          
## D 0.4976303 0.5907173

distance(esophagus, "wunifrac") # weighted UniFrac

##           B         C
## C 0.2035424          
## D 0.2603371 0.2477016

distance(esophagus, "jaccard") # vegdist jaccard

##           B         C
## C 0.5776398          
## D 0.6645570 0.7427056

distance(esophagus, "g") # betadiver method option "g"

##           B         C
## C 0.6136364          
## D 0.6250000 0.6078431

6.4.2 UniFrac and weighted UniFrac

UniFrac is a recently-defined~\cite{Lozupone:2005gn and popular distance metric to summarize the difference between pairs of ecological communities. All UniFrac variants use a phylogenetic tree of the relationship among taxa as central information to calculating the distance between two samples/communities. An unweighted UniFrac distance matrix only considers the presence/absence of taxa, while weighted UniFrac accounts for the relative abundance of taxa as well as their phylogenetic distance. Prior to phyloseq, a non-parallelized, non-Fast implementation of the unweighted UniFrac was available in \R{ packages (picante::unifrac~\cite{Kembel:2010ft). In the phyloseq package we provide optionally-parallelized implementations of Fast UniFrac~\cite{Hamady:2009fk (both weighted and unweighted, with plans for additional UniFrac variants), all of which return a sample-wise distance matrix from any `phyloseq-class object that contains a phylogenetic tree component.

The following is an example calculating the UniFrac distance (both weighted and unweighted) matrix using the “esophagus” example dataset:

data(esophagus)
distance(esophagus, "wUniFrac")
distance(esophagus, "uUniFrac")

See the phyloseq demo page about fast parallel UniFrac.