PartIIphyloseq.rnw

\subsection*{phyloseq}

\BioCpkg{phyloseq}\cite{phyloseqplosone}
is an R package to import, store, analyze, and graphically display
microbiome census data that has already been
organized into counts of sequence features for each sample.
These sequence features can be RSVs, OTUs, 
or other analogous features.
It is most useful when there is also associated
sample data,
phylogeny,
and/or taxonomic assignment for each feature/taxa.
\BioCpkg{phyloseq} leverages and builds upon
many of the tools available in R
for ecology and phylogenetic analysis (\citeCRANpkg{vegan}, \citeCRANpkg{ade4}, \citeCRANpkg{ape}),
while also using the powerful and extensible 
graphics framework, ggplot2~\citeCRANpkg{ggplot2},
to easily produce publication-quality graphics of complex phylogenetic data.
The {\tt phyloseq} package uses a system of S4 data classes
to store related data
as a single, self-consistent, self-describing experiment-level object,
making it easier to share data and reproduce analyses.
In general, phyloseq seeks to facilitate the use of R
for efficient interactive and reproducible analysis
of microbiome count data, 
jointly with important sample or taxonomic covariates.

\subsection*{Further documentation}

This tutorial shows a useful example workflow,
but many more analyses are available to you in phyloseq, and R in general,
than can fit in a single workflow.
The \href{http://joey711.github.io/phyloseq/}{phyloseq home page}
is a good place to begin browsing additional phyloseq documentation,
as are the three vignettes included within the package,
and linked directly at
\href{http://bioconductor.org/packages/release/bioc/html/phyloseq.html}{the phyloseq release page on Bioconductor}.


%\subsection*{Load Packages}

<<joeysetupknitr, cache=FALSE, include=FALSE>>=
##For later version
library("BiocStyle")
@
<<xtrasetupknitr, cache=FALSE, echo=FALSE>>=
options(wdith=88)
options(digits=3)
opts_chunk$set(cache = FALSE, fig.path="phyloseqfigure/",
               fig.width = 6, fig.height = 3,
               warning = FALSE, fig.align = "center", message = FALSE)
@
<<joeysetup, echo=FALSE>>=
library("ggplot2")
theme_set(theme_bw())
min_theme <- theme_update(panel.border = element_blank(),
                          panel.grid = element_blank(),
                          axis.ticks = element_blank(),
                          legend.title = element_text(size = 8),
                          legend.text = element_text(size = 6),
                          axis.text = element_text(size = 6),
                          axis.title = element_text(size = 8),
                          strip.background = element_blank(),
                          strip.text = element_text(size = 8),
                          legend.key = element_blank())

@


\subsection*{Loading Data}
Many use cases result in the need to import and combine different
data into a phyloseq class object,
this can be done using th \Rfunction{import\_biom}
function to read recent QIIME format files, older files
can still be imported with \Rfunction{import\_qiime}.
More complete details can be found on the
\href{https://www.bioconductor.org/packages/release/bioc/vignettes/phyloseq/inst/doc/phyloseq-FAQ.html}{phyloseq FAQ page}.

In the previous section the results of \BioCpkg{dada2} sequence processing
were organized into a phyloseq object.
This object was also saved in R-native serialized RDS format.
We will re-load this here for completeness as the initial object \Robject{p0}.

<<>>=
library("phyloseq")
library("gridExtra")
ps = readRDS("data/ps.rds")
ps
@


\subsubsection*{Shiny-phyloseq}

It can be beneficial to start the data exploration
process interactively, this often saves time
in detecting outliers and specific features of the data.
\href{http://joey711.github.io/shiny-phyloseq/}{Shiny-phyloseq}
\cite{mcmurdie2015}
is an interactive web application that provides a graphical user interface to
the phyloseq package.
The object just loaded into the R session in this workflow
is  suitable for this graphical interaction with Shiny-phyloseq.


\subsection*{Filtering}

\BioCpkg{phyloseq}
provides useful tools for filtering, subsetting, and agglomerating taxa --
a task that is often appropriate or even necessary
for effective analysis of microbiome count data.
In this subsection, we graphically explore the prevalence of taxa
in the example dataset,
and demonstrate how this can be used as a filtering criteria.
One of the reasons to filter in this way
is to avoid spending much time analyzing taxa
that were seen only rarely among  samples.
This also turns out to be a useful filter of noise
(taxa that are actually just artifacts of the data collection process),
a step that should probably be considered essential for datasets constructed
via heuristic OTU-clustering methods,
which are notoriously prone to generating spurious taxa.

\subsubsection*{Taxonomic Filtering}

In many biological settings, the set of all organisms from all samples are well-represented in the available taxonomic reference database.
When (and only when) this is the case, it is reasonable or even advisable to filter taxonomic features
for which a high-rank taxonomy could not be assigned.
Such ambiguous features in this setting are almost always sequence artifacts that don't exist in nature.
It should be obvious that such a filter is not appropriate for samples from poorly characterized or novel specimens, at least until the possibility of taxonomic novelty can be satisfactorily rejected.
Phylum is a useful taxonomic rank to consider using for this purpose, but others may work effectively for your data.

To begin, create a table of read counts for each Phylum present in the dataset.

<<taxfilter0>>=
# Show available ranks in the dataset
rank_names(ps)
# Create table, number of features for each phyla
table(tax_table(ps)[, "Phylum"], exclude = NULL)
@

This shows a few phyla for which only one feature was observed.
Those may be worth filtering, and we'll check that next.
First, notice that in this case, six features were annotated with a Phylum of NA.
These features are probably artifacts in a dataset like this,
and should be removed.

The following ensures that features with ambiguous phylum annotation are also removed.
Note the flexibility in defining strings that should be considered ambiguous annotation.

<<removeNAphyla>>=
ps0 <- subset_taxa(ps, !is.na(Phylum) & !Phylum %in% c("", "uncharacterized"))
@

A useful next step is to explore feature \emph{prevalence} in the dataset,
which we will define here as
the number of samples in which a taxon appears at least once.

<<prevfilter0>>=
# Compute prevalence of each feature, store as data.frame
prevdf = apply(X = otu_table(ps0),
               MARGIN = ifelse(taxa_are_rows(ps0), yes = 1, no = 2),
               FUN = function(x){sum(x > 0)})
# Add taxonomy and total read counts to this data.frame
prevdf = data.frame(Prevalence = prevdf,
                    TotalAbundance = taxa_sums(ps0),
                    tax_table(ps0))
@


Are there phyla that are comprised of mostly low-prevalence features?
Compute the total and average prevalences of the features in each phylum.

<<>>=
plyr::ddply(prevdf, "Phylum", function(df1){cbind(mean(df1$Prevalence),sum(df1$Prevalence))})
@

Deinococcus-Thermus appeared in just over one percent of samples,
and Fusobacteria appeared in just 2 samples total.
In some cases it might be worthwhile to explore
these two phyla in more detail despite this
(though probably not Fusobacteria's two samples).
For the purposes of this example, though,
they will be filtered from the dataset.


<<taxfilter>>=
# Define phyla to filter
filterPhyla = c("Fusobacteria", "Deinococcus-Thermus")
# Filter entries with unidentified Phylum.
ps1 = subset_taxa(ps0, !Phylum %in% filterPhyla)
ps1
@


\subsubsection*{Prevalence Filtering}

The previous filtering steps are considered \emph{supervised},
because they relied on prior information
that is external to this experiment
(a taxonomic reference database).
This next filtering step is completely \emph{unsupervised},
relying only on the data in this experiment,
and a parameter that we will choose after exploring the data.
Thus, this filtering step can be applied even in settings
where taxonomic annotation is unavailable or unreliable.

First, explore the relationship of
prevalence and total read count
for each feature.
Sometimes this reveals outliers that should probably be removed,
and also provides insight into the ranges of either feature
that might be useful.
This aspect depends quite a lot on the experimental design
and goals of the downstream inference,
so keep these in mind.
It may even be the case that different types of downstream inference
require different choices here.
There is no reason to expect ahead of time that
a single filtering workflow is appropriate for all analysis.



<<plotprevalence, fig.show="hide",fig.width=7, fig.height=4>>=
# Subset to the remaining phyla
prevdf1 = subset(prevdf, Phylum %in% get_taxa_unique(ps1, "Phylum"))
ggplot(prevdf1, aes(TotalAbundance, Prevalence / nsamples(ps0),color=Phylum)) +
  # Include a guess for parameter
  geom_hline(yintercept = 0.05, alpha = 0.5, linetype = 2) +  geom_point(size = 2, alpha = 0.7) +
  scale_x_log10() +  xlab("Total Abundance") + ylab("Prevalence [Frac. Samples]") +
  facet_wrap(~Phylum) + theme(legend.position="none")
@
\begin{figure}
\includegraphics[width=\linewidth]{phyloseqfigure/plotprevalence-1}
\caption{Taxa prevalence versus total counts. Each point is a different taxa. Exploration of the data in this way is often useful for selecting filtering parameters, like the minimum prevalence criteria we will used to filter the data above.}
\end{figure}
Sometimes a natural separation in the dataset reveals itself,
or at least, a conservative choice that is in a stable region
for which small changes to the choice would have
minor or no effect on the biological interpreation (stability).
Here no natural separation is immediately evident,
but it looks like we might reasonably define a prevalence threshold
in a range of zero to ten percent or so.
Take care that this choice does not introduce bias
into a downstream analysis of association of differential abundance.

The following uses five percent of all samples
as the prevalence threshold.

<<prevalencefilter>>=
# Define prevalence threshold as 5% of total samples
prevalenceThreshold = 0.05 * nsamples(ps0)
prevalenceThreshold
# Execute prevalence filter, using `prune_taxa()` function
keepTaxa = rownames(prevdf1)[(prevdf1$Prevalence >= prevalenceThreshold)]
ps2 = prune_taxa(keepTaxa, ps0)
@


\subsection*{Agglomerate taxa}

When there is known to be a lot of
species or sub-species functional redundancy
in a microbial community,
it might be useful to agglomerate the data features
corresponding to closely related taxa.
Ideally we would know the functional redundancies perfectly ahead of time,
in which case we would agglomerate taxa using those defined relationships
and the \Rfunction{merge\_taxa()} function in phyloseq.
That kind of exquisite functional data is usually not available,
and different pairs of microbes will have different sets of overlapping functions,
complicating the matter of defining appropriate grouping criteria.

While not necessarily the most useful or functionally-accurate
criteria for grouping microbial features
(sometimes far from accurate),
taxonomic agglomeration has the advantage of being much easier to define ahead of time.
This is because taxonomies are usually defined
with a comparatively simple tree-like graph structure
that has a fixed number of internal nodes, called ``ranks''.
This structure is simple enough for the phyloseq package
to represent taxonomies as table of taxonomy labels.
Taxonomic agglomeration groups all the ``leaves'' in the hierarchy
that descend from the user-prescribed agglomerating rank, this is sometimes called
`glomming'.

The following example code shows how one
would combine all features that descend from the same genus.

<<taxglom>>=
# How many genera would be present after filtering?
length(get_taxa_unique(ps2, taxonomic.rank = "Genus"))
ps3 = tax_glom(ps2, "Genus", NArm = TRUE)
@

If taxonomy is not available or not reliable,
tree-based agglomeration is a "taxonomy-free"
alternative to combine data features corresponding to closely-related taxa.
In this case, rather than taxonomic rank, the user specifies
a tree height corresponding to the phylogenetic distance between features
that should define their grouping.
This is very similar to ``OTU Clustering'',
except that in many OTU Clustering algorithms the sequence distance being used
does not have the same (or any) evolutionary definition.

<<tipglom>>=
h1 = 0.4
ps4 = tip_glom(ps2, h = h1)
@

Here phyloseq's \Robject{plot\_tree()} function
compare the original unfiltered data,
the tree after taxonoic agglomeration,
and the tree after phylogenetic agglomeration.
These are stored as separate plot objects,
then rendered together in one combined graphic
using \Robject{gridExtra::grid.arrange}.

<<plotglomprep>>=
multiPlotTitleTextSize = 8
p2tree = plot_tree(ps2, method = "treeonly",
                   ladderize = "left",
                   title = "Before Agglomeration") +
  theme(plot.title = element_text(size = multiPlotTitleTextSize))
p3tree = plot_tree(ps3, method = "treeonly",
                   ladderize = "left", title = "By Genus") +
  theme(plot.title = element_text(size = multiPlotTitleTextSize))
p4tree = plot_tree(ps4, method = "treeonly",
                   ladderize = "left", title = "By Height") +
  theme(plot.title = element_text(size = multiPlotTitleTextSize))
@

\begin{figure}[H]
<<plotglomtree>>=
# group plots together
grid.arrange(nrow = 1, p2tree, p3tree, p4tree)
@
\caption{The original tree (left), taxonomic agglomeration at Genus rank (middle), phylogeentic agglomeration at a fixed distance of 0.4 (right).}
\label{fig:glomsummary}
\end{figure}


\subsection*{Abundance value transformation}

It is usually necessary to transform microbiome count data
to account for differences in library size, variance, scale, etc.
The phyloseq package provides a flexible interface for defining
new functions to accomplish these transformations of the abundance values
via \Robject{transform\_sample\_counts()}.
The first argument to this function is the phyloseq object you want to transform,
and the second argument is an R function that defines the transformation.
The R function is applied sample-wise, expecting that the first unnamed argument
is a vector of taxa counts in the same order as the phyloseq object.
Additional arguments are passed on to the function
specified in the second argument,
providing an explicit means to include
pre-computed values, previously defined parameters/thresholds,
or any other object that might be appropriate
for computing the transformed values of interest.

This example begins by defining a custom plot function,
\Robject{plot\_abundance()},
that uses phyloseq's \Robject{psmelt()} function
to define a relative abundance graphic.
We will use this to compare differences
in scale and distribution of the abundance values
in our phyloseq object before and after transformation.

<<abundancetransformation>>=
plot_abundance = function(physeq,title = "",
                          Facet = "Order", Color = "Phylum"){
  # Arbitrary subset, based on Phylum, for plotting
  p1f = subset_taxa(physeq, Phylum %in% c("Firmicutes"))
  mphyseq = psmelt(p1f)
  mphyseq <- subset(mphyseq, Abundance > 0)
  ggplot(data = mphyseq, mapping = aes_string(x = "sex",y = "Abundance",
                              color = Color, fill = Color)) +
    geom_violin(fill = NA) +
    geom_point(size = 1, alpha = 0.3,
               position = position_jitter(width = 0.3)) +
    facet_wrap(facets = Facet) + scale_y_log10()+
    theme(legend.position="none")
}
@

The transformation in this case converts the counts
from each sample into their frequencies,
often referred to as \emph{proportions}
or \emph{relative abundances}.
This function is so simple that it is easiest
to define it within the function call to
\Robject{transform\_sample\_counts()}.

<<abundancetransformation2>>=
# Transform to relative abundance. Save as new object.
ps3ra = transform_sample_counts(ps3, function(x){x / sum(x)})
@

Now plot the abundance values before and after transformation.

\begin{figure}[H]
<<abundancetransformation3, fig.height=4, fig.width=8.5>>=
plotBefore = plot_abundance(ps3,"")
plotAfter = plot_abundance(ps3ra,"")
# Combine each plot into one graphic.
grid.arrange(nrow = 2,  plotBefore, plotAfter)
@
\caption{
Comparison of original abundances (top panel)
and relative abundances (lower).
}
\end{figure}


\subsection*{Subset by taxonomy}

Notice on the previous plot that
\emph{Lactobacillales} appears to be a taxonomic Order
with bimodal abundance profile in the data.
We can check for a taxonomic explanation of this pattern
by plotting just that taxonomic subset of the data.
For this, we subset with the \Rfunction{subset\_taxa()} function,
and then specify a more precise taxonomic rank
to the \Robject{Facet} argument of the \Robject{plot\_abundance} function
that we defined above.

\begin{figure}[H]
<<subsettaxa>>=
psOrd = subset_taxa(ps3ra, Order == "Lactobacillales")
plot_abundance(psOrd, Facet = "Genus", Color = NULL)
@
\caption{Violin plot of the relative abundances of Lactobacillales taxonomic Order,
grouped by host sex and genera.
Here it is clear that the apparent biomodal distribution
of Lactobacillales on the previous plot
was the result of a mixture of two different genera,
with the typical \emph{Lactobacillus} relative abundance
much larger than \emph{Streptococcus}.}
\end{figure}