Workflow part III:
Analyzing ASV Data

---

Goals

Analyze microbiome experimental data as a phyloseq object - explore ecological metrics and identify differentially abundant taxa.

---

Before We Get Started

Notes

• This is a slightly modified version of a document created by Matthew Willman.

• To convert an Rmd (R Markdown) file to an R script, type the following in an R console: knitr::purl(input="<filename>.Rmd"). (You can download this Rmd file by clicking the Code button in the top-right of this page – but we will open it at OSC.)

• This document deliberately does not have any R output. In case you want to check your output against a reference, you can find a document with R output here.

Open this file in RStudio Server at OSC

• This assumes you still have an active RStudio Server job at OSC; if not see the instructions on the previous page to start one.

• In the Files pane, in the markdown directory, open the file 08-postASV-analysis.Rmd. That’s this file!

---

Getting Started

Install and Load Packages

Add our custom library again if you restarted your R session:

.libPaths(new = '/fs/project/PAS0471/.R/4.0/')

And load the necessary packages:

packages <- c('tidyverse', 'vegan', 'phyloseq', 'decontam', 'ape',
              'DESeq2', 'microbiome', 'metagenomeSeq', 'remotes',
              'ampvis2', 'breakaway')
pacman::p_load(char = packages)

# If you wanted to install and load these packages yourself,
# first make sure you have the pacman package installed:
## install.packages('pacman')
# Then, the code above would work except for the last two packages,
# which are available from Github only. To install these, you would run:
## remotes::install_github("MadsAlbertsen/ampvis2")
## remotes::install_github("adw96/breakaway")

Set Directories

# Dir with input files:
indir <- 'analysis/ASV_inference/'

# Dir for output:
outdir <- 'analysis/postASV_analysis/'

Load the Data

This session starts with a phyloseq object similar to the one generated using DADA2 in the previous session.

ps_raw <- readRDS(file.path(indir, 'ps_16S_V4_withtree.rds'))

# This would load the object from the previous session, which doesn't contain a tree:
# ps_raw <- readRDS(file.path(indir, 'ps_V4.rds'))

Our phyloseq object is made up of the following components (slots):

• otu_table
• sample_data
• tax_table
• phy_tree

Let’s have a look at the object:

ps_raw

# And the assigned taxon names:
taxa_names(ps_raw)[1:3]

As we can see, the taxon names are currently the associated sequences. We can create a new phyloseq component to store these sequences, then rename the ASVs so something shorter (“ASV_1”, “ASV_2”, etc.).

dna <- Biostrings::DNAStringSet(taxa_names(ps_raw))
names(dna) <- taxa_names(ps_raw)

# Merge sequence object into the phyloseq object:
ps_raw <- merge_phyloseq(ps_raw, dna)
ps_raw

# Rename ASVs:
taxa_names(ps_raw) <- paste("ASV", 1:ntaxa(ps_raw), sep = "_")
taxa_names(ps_raw)[1:3]

---

Filter taxa

Identify and Remove Contaminants

It is possible to introduce contaminating microbes during sample preparation. Before analyzing the data, we will identify and remove probable contaminants using the decontam package.

In this case, we will define a conataminant as an ASV whose abundance correlates with DNA concentration (post-PCR). Our assumption here is that each soil taxon’s abundance should be independant of DNA concentration. However, if we were to spike each sample with a contaminant, that contaminant would show a greater abundance in samples with lower DNA concentrations. To do this, we need DNA concentration measured after PCR and before pooling of the samples.

These data are in column 7 of the sample data component:

head(sample_data(ps_raw))

MCIC measured our DNA concentration by comparing band intensities to a single reference band. MCIC’s DNA concentration measurements range from 0 to 1 (the reference). Zero in this case is coded as NA. The decontam package can’t handle NA’s or 0’s, so we need to change them to a small value (e.g. 0.01).

# First, we identify the NAs in the DNA measurement column:
NAs <- is.na(sample_data(ps_raw)$PCR_prod_V4.V5)

# Then, we replace NAs with 0.01:
sample_data(ps_raw)$PCR_prod_V4.V5[NAs] <- 0.01

We will use decontam::isContaminant(method = "frequency") to test each taxon against the null hypothesis: abundance is not associated with DNA concentration. Because the DNA concentration measurements are not very accurate, and we want to make sure to remove possible contaminants, we use a relaxed p-value (0.2) to reject the null hypothesis:

contam_df_freq <- isContaminant(ps_raw,
                                method = "frequency",
                                conc = "PCR_prod_V4.V5",
                                threshold = 0.2)
head(contam_df_freq)

As we can see, the top 3 ASVs are identified as probable contaminants. How many contaminants were identified (TRUE in the contaminant column) in total?

table(contam_df_freq$contaminant)

What are the most abundant contaminant taxa?

# Get abundance ranks for contaminant ASVs:
head(which(contam_df_freq$contaminant))

# Check which taxa are contaminants:
ps_contam <- prune_taxa(contam_df_freq$contaminant, ps_raw)
head(tax_table(ps_contam))

Let’s take a look at our tested associations for several taxa:

plot_frequency(ps_raw,
               taxa_names(ps_raw)[c(1, 2, 3, 4, 5, 11)],
               conc = "PCR_prod_V4.V5") +
  xlab("DNA Concentration (relative to single marker)")

As we can see, abundance of ASVs 1, 2, 3, and 11 are particularly high in samples with low DNA concentration.

Which samples have highest abundances of the contaminants?

asv1 <- otu_table(ps_raw)[, 1]
head(asv1[order(asv1, decreasing = TRUE)])

asv2 <- otu_table(ps_raw)[, 2]
head(asv2[order(asv2, decreasing = TRUE)])

How do these contaminants overlap with ASVs present in the negative control?

otu_table(ps_raw)[70:75, 1:11]

ASVs have high read counts in the extraction negative control. Thus, we can be fairly certain these are contaminants and not representative of our experimental samples.

Remove the contaminants that we identified:

ps_noncontam <- prune_taxa(!contam_df_freq$contaminant, ps_raw)

ps_noncontam

What proportion of our count data were removed as contaminants?

pre <- sum(sample_sums(ps_raw))
post <- sum(sample_sums(ps_noncontam))

(pre-post) / pre

Remove non-bacterial, non-archaeal ASVs

The V4 region of 16S rRNA is conserved within certain bacteria, archaea, chloroplasts, mitochondria, and eukaryotes. Since we are only interested in bacteria and archaea here, we’ll need to remove other taxa from our object.

# Create ps subsets for chloroplast, mitochondria, and eukaryotes:
chlr <- subset_taxa(ps_noncontam, Order == "Chloroplast")
mit <- subset_taxa(ps_noncontam, Family == "Mitochondria")
euk <- subset_taxa(ps_noncontam, Kingdom == "Eukaryota")

# Get all taxon IDs that we want to remove:
bad_taxa <- c(taxa_names(chlr), taxa_names(mit), taxa_names(euk))

# Get all taxon IDs that we want to keep:
all_taxa <- taxa_names(ps_noncontam)
good_taxa <- all_taxa[!(all_taxa %in% bad_taxa)]

# Subset phyloseq object:
ps_bac_arc <- prune_taxa(good_taxa, ps_noncontam)

What proportion of ASVs were kept?

pre <- sum(sample_sums(ps_noncontam))
post <- sum(sample_sums(ps_bac_arc))
post / pre

Finally, we’ll also remove non-experimental samples (e.g. negative controls). Note, though, that non-experimental samples are important for evaluating the quality of your reagents, consistency across sequencing runs, and contamination.

ps <- subset_samples(ps_bac_arc, Experiment == "HCC")

ps

---

Filter samples

After having filtered out unwanted taxa, we now need to filter out uninformative samples – those with low total taxon counts.

First, how many counts do we have for each sample?

sums <- sample_sums(ps)
sums[order(sums)]

# Lowest total count:
min(sums)

# Sample with lowest counts:
names(sums[order(sums)][1])

To remove uninformative samples, we will only keep those with over 1,000 counts. We will also remove a technical replicate, though note that technical replicates are important for quality control of runs and for understanding the reproducibility of the methodology.

# Remove samples with low counts:
ps <- subset_samples(ps, sample_sums(ps) > 1000)

# Remove the technical replicate:
ps <- subset_samples(ps, SampleID != "501S4")

Let’s take a look at our current phyloseq object:

sample_data(ps)[1:10]

head(t(otu_table(ps)[1:10]))

head(tax_table(ps))

---

Normalization

Normalization is currently a much-discussed issue of microbiome studies. Differences in read depth between samples often need to be corrected before analysis. Several normalization methods have been proposed, and no single method is perfect. It may be that the most appropriate method depends on the analysis.

In this tutorial, we will use:

• Proportion-normalized data to estimate ecological metrics.

• Stabilizing transformation normalized data to identify differentially abundant taxa.

However, for your reference, we will also present code for two other normalization methods: rarefaction and cumulative sum scaling.

For more details, read McMurdie and Holmes 2014 and McKnight et al. 2019.

---

Proportion Normalization

Proportion normalization involves dividing each OTU count by the total sum for each sample. The resulting count data will add up to 1 (100%) for each sample.

The microbiome::transform function can be used to easily normalize count data as proportions in a phyloseq object:

# Proportion normalization:
ps_prop <- transform(ps, "compositional")

# Have a look at the resulting OTU table:
otu_table(ps_prop)[1:5, 1:5]

# Check what the sums for each sample are now:
head(sample_sums(ps_prop))

---

Variance stabilizing transformation

This transformation method, available in DESeq2, normalizes data with respect to library size as well as dispersion.
Type ?vst or read Anders and Huber 2010 for more information.

# First we convert the phyloseq object to a DESeq object:
dds <- phyloseq_to_deseq2(ps, ~1)

# Then we estimate the size factors:
dds <- estimateSizeFactors(dds)

# Now can do the transformation:
dds_vst <- vst(dds, blind = TRUE)

# Have a look at the resulting counts:
vst_counts <- assay(dds_vst)
t(vst_counts)[1:5, 1:5]

---

Rarefaction

Rarefaction can be used to subset data such that the library depth is the same for each sample. Because sampling of the data is random, rarefaction can account for an effect of total read count on taxa richness. However, rarefaction is no longer considered to be a good way to normalize amplicon sequencing data (see McMurdie & Holmes 2014).

ps_rarefied <- rarefy_even_depth(ps,
                                 rngseed = 1,
                                 sample.size = min(sample_sums(ps)),
                                 replace = FALSE)

# Have a look at the result:
otu_table(ps_rarefied)[1:5, 1:5]

head(sample_sums(ps_rarefied))

Compare these sample sums to the non-rarified data sums:

head(sample_sums(ps))

---

Cumulative Sum Scaling

The metagenomeSeq Cumulative Sum Scaling (CSS) normalization is another option developed for microbiome data. For more information, read Paulson et al. 2013.

# Convert the phyloseq object to a metagenomeseq object:
mgs_css <- phyloseq_to_metagenomeSeq(ps)

# Perform the Cumulative Sum Scaling:
mgs_css <- cumNorm(mgs_css)

# Extract the counts and add them to a separate phyloseq object:
css_counts <- MRcounts(mgs_css, norm = TRUE)
ps_css <- ps
otu_table(ps_css) <- otu_table(t(css_counts), taxa_are_rows = FALSE)

Now lets compare the original data to the CSS normalized data:

otu_table(ps)[1:5, 1:5]
head(sample_sums(ps))

otu_table(ps_css)[1:5, 1:5]
head(sample_sums(ps_css))

---

Plot abundances

In this and following the examples we will use the proportions data.

# Subset to the T1 treatment
T1 <- subset_samples(ps_prop, Treatment == "T1")

# Agglomerate taxa to phylum level:
T1_phylum <- tax_glom(T1, taxrank = "Phylum", NArm = FALSE)

# Remove uncommon phyla:
T1_phylum <- subset_taxa(T1_phylum, taxa_sums(T1_phylum) > 0.1)

# Plot:
plot_bar(T1_phylum, fill = "Phylum") +
   facet_wrap(~Timepoint, scales = "free_x", nrow = 1)

---

Alpha diversity

Alpha diversity measures taxonomic diversity within a single population. Measures of alpha diversity include taxonomic richness (i.e. number of taxa), and indices which combine taxonomic richness with some measure of evenness.

Note: Many of these methods are greatly influenced by singleton data (i.e. taxa represented by a single count), meaning they may be unreliable if your data excludes singletons. DADA2 does not output singletons if it is run individually on each sample. Thus, if using dada2 to infer ASVs, it is advisable to use pool = TRUE when running the dada algorithm. This will allow calling per-sample singletons, but not per-study singletons.

Taxonomic richness: Rarefaction plot

Our goal in defining richness is to determine the number of unique taxa within each population (in our case, each plot). This is difficult, as our samples contain only a portion of the total richness within our plots. As sequence depth increases for each sample, so does the number of taxa. This can be illustrated with a rarefaction plot.

To make a rarefaction plot, we draw random samples from our data and count the number of unique ASVs as samples are drawn. The resulting rarefaction curve is expected to rise quickly then plateu as the most abundant taxa are represented. We can make a quick rarefaction curve plot directly from our phyloseq object of all samples using the vegan package:

rarecurve(otu_table(ps), step = 500, xlab = "Sample Size", ylab = "Taxa")

We can also split rarefaction curves by group using the ampvis2 package.

metadata <- data.frame(sample_data(ps), check.names = FALSE)

asv_table <- data.frame(t(otu_table(ps)),
                       tax_table(ps),
                       check.names = FALSE)

ampvis2 requires a Species column, which must be added to our data. If your table already has a Species column, skip this step:

asv_table$Species <- NA

Load an ampvis2 object and calculate rarefaction:

ps3 <- amp_load(asv_table, metadata)

rar_plot <- amp_rarecurve(ps3,
                          stepsize = 100,
                          facet_by = "Timepoint",
                          color_by = "Treatment") +
   ylab("Number of observed ASVs")

rar_plot

As we can see, as read depth increases for each sample, so does taxonomic richness. In all samples, increased sampling would result in increased richness. Thus, our observed sample richness is lower than the true population richness for every plot.

Further, sequence depth and total richness appear to be correlated across samples:

dp_by_asv <- data.frame(ct_sums = sample_sums(ps),
                        asvs = rowSums(otu_table(ps) != 0))

ggplot(data = dp_by_asv) +
  geom_point(aes(x =  ct_sums, y = asvs)) +
  labs(x = "Sequence depth", y= "ASV count") +
  theme_minimal()

We can test whether this correlation is significant:

cor <- cor.test(dp_by_asv$ct_sums, dp_by_asv$asvs, method = "spearman")
cor

Here we observe a significant correlation between read count and observed taxonomic richness.

Thus, a sample from a population with high actual diversity, but low read depth, could have low observed richness. Fortunately, there are methods to estimate the number of unobserved taxa. For an overview of these methods, read Bunge et al. 2014.

Frequency counts

To illustrate the difficulty of estimating population richness and importance of singletons to alpha diversity metrics, we’ll look at the frequency count distribution for our first sample: 101-S1.

We’ll count the number of singletons, doubletons, etc., then plot frequency as a function of count.

# Get the frequency count distribution for the 1st sample:
frequencytablelist <- build_frequency_count_tables(otu_table(ps))
freq_ct <- frequencytablelist[[1]]
colnames(freq_ct) <- c("Frequency", "Count")

# Plot:
ggplot(freq_ct, aes(Frequency, Count)) + 
  geom_point() + 
  scale_x_continuous(breaks = seq(0, 600, by = 50)) + 
  ggtitle("101-S1") +
  theme_minimal()

We can see there are ~3,200 taxa represented as singletons, and ~1,600 taxa as doubletons. This is a high proportion of the total number taxa are observed in this sample:

# Total number of taxa:
sum(freq_ct$Count)

# Proportion of singletons
freq_ct$Count[1] / sum(freq_ct$Count)

How many more taxa might we observe if we increased our sampling?

Diversity Indices

We can use several indices to explore how evenly species are distributed within each sample. Note that ‘Observed’ is the sample taxonomic richness. The other three indices combine richness and abundance data for each taxon.

plot_richness(ps,
              x = "Timepoint",
              shape = "Treatment",
              color = "Treatment", 
              measures = c("Observed", "Shannon", "Simpson", "InvSimpson"))

Notice that several T2:S3 samples have low measures, suggesting they have low diversity (i.e. abundances are dominated by a few taxa) compared to the other samples.

---

Beta diversity

Beta diversity measures differences in microbial compositions between populations. If available (use ps_16S_V4_withtree.rds), we can use a phylogenetic tree and otu table to estimate the UniFrac for each sample pair.

The result is a distance matrix with height and width equal to the number of samples. UniFrac may be weighted or unweighted. Selection of weighted or unweighted will depend on the question you want to answer. Unweighted UniFrac uses only presence or absence of a taxon from the otu table (i.e. 1000 will be treated the same as 1). It is appropriate if you want to test qualitative taxa differences between samples. Weighted UniFrac uses presence as well as quantity data (i.e. 1 < 1000).

We are interested to know what changes to taxon abundance may be driven by our treatments, so we will use weighted UniFrac.

Ordination

First we’ll make a weighted UniFrac distance matrix. For UniFrac, you require a tree as part of your phyloseq object.

set.seed(100)
wUF <- phyloseq::distance(ps_prop, method = "wunifrac")

as.matrix(wUF)[1:6, 1:6]

# Note: If your *phyloseq* object lacks a phylogenetic tree,
# you could use a dissimilarity index that does not require one, such as Bray-Curtis:
# bray <- vegan::vegdist(otu_table(ps_prop), method = "bray")
# as.matrix(bray)[1:6,1:6]

Then perform Principal Coordinate Analysis (PCoA) on your favorite distance matrix. In this case, we are using weighted UniFrac.

do <- ordinate(ps_prop, method = "PCoA", distance = wUF)

baseplot <- plot_ordination(ps_prop, do,
                            type = "samples",
                            color = "Treatment",
                            shape = "Timepoint") +
  ggtitle("HCC ordination, all samples\nmethod=PCoA, distance=weighted UniFrac") +
  theme_bw()

# Single plot with all samples:
baseplot +
  geom_point(size = 2)
  
# Idem, but with sample names on the plot: 
baseplot +
  geom_point(size = 1) + 
  geom_label(aes(label = SampleID), size = 3.5, color = "black")

# Facet by time point:
baseplot +
  geom_point(size = 2) + 
  ggtitle("HCC ordination, all samples\nmethod=PCoA, distance=weighted UniFrac") +
  facet_wrap(~Timepoint, 2)

We can see some grouping by timepoint. Depending on your experiment and questions, you might want to consider other ordination or clustering approaches.

Permutational analysis of variance (PERMANOVA)

Finally, we’ll test whether there is a significant ecological level treatment effect. PERMANOVA sounds fancy, but it is just an ANOVA performed using permutations. Permutations are used to determine how data may appear if there is no treatment effect and group differences are due to random chance. Observed data are then compared to the randomized data to calculate a p-value.

Treatments are independent within Timepoints, so we will first subset the data by Timepoint:

ps_S1 <- subset_samples(ps_prop, Timepoint == "S1")
ps_S3 <- subset_samples(ps_prop, Timepoint == "S3")
ps_S4 <- subset_samples(ps_prop, Timepoint == "S4")

Subset the metadata and estimate UniFrac for S3:

metadata <- as(sample_data(ps_S3), "data.frame")
unifrac_dist <- UniFrac(ps_S3, weighted = TRUE)

Then we’ll use vegan::adonis2 to perform the PERMANOVA.

set.seed(100)  # Set the seed to make the analysis reproducible

permanova <- adonis2(unifrac_dist ~ Treatment,
                     data = metadata,
                     permutations = 10000)

permanova

The low p-value suggests that within the S3 timepoint (post-termination), there is a significant treatment effect.

---

Differential taxon abundance (DESeq)

Here, we will use DESeq2 to identify differentially abundant taxa between time points and treatments. DESeq2 was designed to analyze RNAseq datasets, which are similar to OTU/ASV data sets in that both handle large, sparse contingency tables generated from Illumina sequencing data. For more details, see McMurdie and Holmes 2014.

First, we’ll convert our non-normalized count data to a DESeq object.

dds <- phyloseq_to_deseq2(ps, ~Timepoint + Treatment)

dds
colData(dds)

There are two ways to analyze interaction effects using DESeq2. The first is to fit a multivariate model (e.g. ~A+B+A:B) and explore the model coefficients. The second is to fit a univariate model and explore pairwise contrasts. Here, we will group our two factors (Treatment and Timepoint) and use the latter approach.

# Prepare the data:
dds$Group <- factor(paste0(dds$Treatment, dds$Timepoint))
design(dds) <- formula(~Group)
dds <- dds[, -50]                    # Remove tech_rep sample
dds$Group <- droplevels(dds$Group)   # Drop tech_rep level from design matrix

# Perform the differential abundance analysis:
dds <- DESeq(dds)

results(dds)

Let’s explore the contrast, “T1S4 : T4S4” (no cc, seedling stage : late termination, seedling stage):

res <- results(dds, contrast = c("Group", "T1S4", "T4S4"))
summary(res)

# Number of taxa with significantly higher abundance in T4S4:
sum(res$padj < 0.1 & res$log2FoldChange > 0, na.rm = TRUE)

# Number of taxa with significantly lower abundance in T4S4:
sum(res$padj < 0.1 & res$log2FoldChange < 0, na.rm = TRUE)

We can sort results by p-value and look at the top ASVs:

res.order <- res[order(res$pvalue), ]
head(res.order)[c(1:2, 6)]

And plot counts for the first in the ordered list ASV:

plotCounts(dds, gene = res.order@rownames[1], intgroup = "Group")

We can also plot counts by sample_ID:

plotCounts(dds, gene = "ASV_48", intgroup = "Group")

The plot can also be saved and customized with ggplot:

d <- plotCounts(dds, gene = "ASV_48", intgroup = "Group", returnData = TRUE)

ggplot(d, aes(x = Group, y = count)) +
  geom_point(position = position_jitter(w = 0.1, h = 0), size = 2) +
  scale_y_log10() +
  ggtitle(expression(paste("ASV 48, log"[10], "-scale"))) +
  theme_bw()

We can also subset our differentially abundant taxa to keep those with an adjusted p-value < 0.05:

sig <- subset(res.order, padj <= 0.05)
head(sig)

What taxa are associated with significantly differentially abundant ASVs?

sig_asv <- rownames(sig)
sig_taxa <- tax_table(ps)[sig_asv,]
head(sig_taxa)

sig_write <- cbind(sig, sig_taxa)
write.csv(sig_write, file.path(outdir, 'DATaxa_T1S1_T1S4.csv'))

Report that we are done!

print('Done with ASV inference.')

---

Resources

• decontam tutorial

• vegan tutorial

• DESeq2 tutorial

• Intro course on microbiome analysis

• UniFrac publication

• Additional UniFrac publication

• Articles on data normalization:

  • Waste not, want not: why rarifying microbiome data is inadmissible
  • Microbiome datasets are compositional: and this is not optional
  • Methods for normalizing microbiome data: an ecological perspective
  • Differential abundance analysis for microbial marker-gene surveys