Workflow part II:
ASV Inference and Taxon Assignment

---

Goals

Process metabarcoding reads, post-adapter removal, including:

• Evaluation of read quality
• Quality trimming and filtering
• Error correction, denoising
• ASV inference
• Read merging
• Taxonomy assignment

---

Before We Get Started

Notes

• This documented was adapted from Callahan et al. 2006 by Matthew Willman, with further edits by Soledad Benitez Ponce and Jelmer Poelstra.

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

Start an RStudio Server job at OSC

For more detailed instructions of the first steps, see this section from our intro to R session.

• Login to OSC at https://ondemand.osc.edu.

• Click on Interactive Apps (top bar) > RStudio Server (Owens and Pitzer).

• Fill out the form as shown here.

• Once your job has started, click Connect to RStudio Server.

• You should automatically be in your personal dir inside the dir /fs/project/PAS0471/workshops/2020-12_micro, with your RStudio Project open. You can see whether a Project is open and which one in the top-right of your screen:

  Here, the project jelmer is open. Your project name is also your username.

  If your Project isn’t open, click on the icon to open it:

  

• Now, click on the markdown directory in the Files pane, and open the file 07-ASV-inference.Rmd. That’s this file!

---

Step 1: Getting Started

Print some information to screen for when running this script as a job:

print('Starting ASV inference script...')
Sys.time()

cat('Working directory:', getwd(), '\n')

Set the number of cores

Most dada2 functions can use multiple cores. Because we are on a cluster and we have reserved only part of a node, auto-detection of the number of cores will not be appropriate (the program will detect more cores than we have available). Therefore, we should specify the appropriate number of cores in our function call below.

We will set the number of cores here, assumes you requested 4 cores in your job submission (change if needed):

n_cores <- 4

Install and load packages

To save time, we have already installed all the necessary R packages at OSC into a custom library. To add this library for the current R session:

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

Then, load the packages:

print('Loading packages...')

packages <- c('tidyverse', 'gridExtra', 'dada2',
              'phyloseq', 'DECIPHER', 'phangorn')
pacman::p_load(char = packages)

# If you wanted to install and load these packages yourself,
# just make sure you have the pacman package installed:
## install.packages('pacman')
# Then, the code above would work, as it will install pacakes as needed.

Set the file paths

We’ll set most of the file paths upfront, which will make it easier to change things or troubleshoot.

# Dir with input fastq files:
indir <- 'data/processed/fastq_trimmed'

# Dirs for output:
filter_dir <- 'data/processed/fastq_filtered'
outdir <- 'analysis/ASV_inference'

dir.create(filter_dir, showWarnings = FALSE, recursive = TRUE)
dir.create(outdir, showWarnings = FALSE, recursive = TRUE)

# Fasta file with training data:
# (Check for an up-to-date version at <https://benjjneb.github.io/dada2/training.html>)
tax_key <- 'data/ref/silva_nr99_v138_train_set.fa' 

# File with sample metadata:
metadata_file <- 'metadata/sample_meta.txt'

Assign fastq files to forward and reverse reads

We will assign the fastq files that we processed with cutadapt to two vectors: one with files with forward reads, and one with files with reverse reads. These files can be distinguished by having “R1” (forward) and “R2” (reverse) in their names.

fastqs_raw_F <- sort(list.files(indir, pattern = '_R1_001.fastq.gz', full.names = TRUE))
fastqs_raw_R <- sort(list.files(indir, pattern = '_R2_001.fastq.gz', full.names = TRUE))

print('First fastq files:')
head(fastqs_raw_F)

Check sample IDs

We’ll get the sample IDs from the fastq file names and from a file with metadata, and will check if they are the same. First we’ll prepare the metadata:

print('Load and prepare sample metadata...')

metadata_df <- read.table(file = metadata_file, sep = "\t", header = TRUE)

colnames(metadata_df)[1] <- 'SampleID'
rownames(metadata_df) <- metadata_df$SampleID

print('First fastq files:')
head(metadata_df)

Let’s compare the sample IDs from the metadata with the fastq filenames:

print('IDs from metadata:')
metadata_df$SampleID

print('Fastq file names:')
head(basename(fastqs_raw_F))    # basename() strips the dir name from the filename

To extract the sample IDs from the fastq file names, we remove everything after “-V4-V5” from the file names using the sub() function:

# sub() arguments: sub(pattern, replacement, vector)
sampleIDs <- sub("-V4-V5_.*", "", basename(fastqs_raw_F))

We can check whether the IDs from the fastq files and the metadata dataframe are the same:

print('Are the sample IDs from the metadata and the fastq files the same?')
identical(sort(metadata_df$SampleID), sampleIDs)

print('Are any samples missing from the fastq files?')
setdiff(sort(metadata_df$SampleID), sampleIDs)

print('Are any samples missing from the metadata?')
setdiff(sampleIDs, sort(metadata_df$SampleID))

As it turns out, we don’t have sequences for three samples in the metadata.

---

Step 2: QC

Plot sequence quality data

DADA2 provide a function to plot the average base quality across sequence reads, plotQualityProfile(). You can generate and evaluate plots for each sample, e.g. the forward reads and reverse reads side-by-side like so:

plotQualityProfile(c(fastqs_raw_F[1], fastqs_raw_R[1]))

This code will generate a pdf file with plots for each sample:

pdf(file.path(outdir, 'error_profiles.pdf'))               # Open a pdf file
for (sample_idx in 1:length(fastqs_raw_F)) {               # Loop through samples
  print(plotQualityProfile(                                # Print plots into pdf
    c(fastqs_raw_F[sample_idx], fastqs_raw_R[sample_idx])) # F and R together
    )
}
dev.off()                                                  # Close the pdf file

    More on QC of fastq files

It is a good idea to run the fastqc program on your fastq files for more extensive QC. This is a stand-alone program that is easy to run from the command-line. When you have many samples, as is often the case, fastqc’s results can moreover be nicely summarized using [multiqc](https://multiqc.info/. In the interest of time, we skipped these steps during this workshop.

---

Step 3: Filtering and Quality Trimming

We will now perform quality filtering (removing poor-quality reads) and trimming (removing poor-quality bases) on the fastq files using DADA2’s filterAndTrim() function.

The filterAndTrim() function will write the filtered and trimmed reads to new fastq files. Therefore, we first define the file names for the new files:

fastqs_filt_F <- file.path(filter_dir, paste0(sampleIDs, '_F_filt.fastq'))
fastqs_filt_R <- file.path(filter_dir, paste0(sampleIDs, '_R_filt.fastq'))

The truncLen argument of filterAndTrim() defines the read lengths (for forward and reverse reads, respectively) beyond which additional bases should be removed, and these values should be based on the sequence quality visualized above. The trimming length can thus be different for forward and reverse reads, which is good because reverse reads are often of worse quality.

It is also suggested to trim the first 10 nucleotides of each read (trimLeft argument), since these positions are likely to contain errors.

maxEE is an important argument that will let DADA2 trim reads based on the maximum numbers of Expected Errors (EE) given the quality scores of the reads’ bases.

print('Filtering and Trimming...')
Sys.time()  # Print the time to keep track of running time for individual steps 
filter_results <-
  filterAndTrim(fastqs_raw_F, fastqs_filt_F,
                fastqs_raw_R, fastqs_filt_R,
                truncLen = c(250,210),
                trimLeft = 10,
                maxN = 0,
                maxEE = c(2,2),
                truncQ = 2,
                rm.phix = FALSE,
                multithread = n_cores, 
                compress = FALSE, verbose = TRUE) 
print('...Done!')
Sys.time()

head(filter_results)

---

Step 4: Dereplication & Error Training

Next, we want to “dereplicate” the filtered fastq files. During dereplication, we condense the data by collapsing together all reads that encode the same sequence, which significantly reduces later computation times.

fastqs_derep_F <- derepFastq(fastqs_filt_F, verbose = FALSE)
fastqs_derep_R <- derepFastq(fastqs_filt_R, verbose = FALSE)

names(fastqs_derep_F) <- sampleIDs
names(fastqs_derep_R) <- sampleIDs

The DADA2 algorithm makes use of a parametric error model (err) and every amplicon dataset has a different set of error rates. The learnErrors method learns this error model from the data, by alternating estimation of the error rates and inference of sample composition until they converge on a jointly consistent solution.

print('Learning errors...')
Sys.time()

err_F <- learnErrors(fastqs_derep_F, multithread = n_cores, verbose = TRUE)
err_R <- learnErrors(fastqs_derep_R, multithread = n_cores, verbose = TRUE)

print('...Done!')
Sys.time()

We’ll plot errors to verify that error rates have been reasonable well-estimated. Pay attention to the fit between observed error rates (points) and fitted error rates (lines):

plotErrors(err_F, nominalQ = TRUE)
plotErrors(err_R, nominalQ = TRUE)

---

Step 5: Infer ASVs

We will now run the core dada algorithm, which infers Amplicon Sequence Variants (ASVs) from the sequences.

This step is quite computationally intensive, and for this tutorial, we will therefore perform independent inference for each sample (pool = FALSE), which will keep the computation time down.

Pooling will increase computation time, especially if you have many samples, but will improve detection of rare variants seen once or twice in an individual sample, but many times across all samples. Therefore, for your own analysis, you will likely want to use pooling, though “pseudo-pooling” is also an option.

print('Inferring ASVs (running the dada algorithm)...')
Sys.time()

dada_Fs <- dada(fastqs_derep_F, err = err_F, pool = FALSE, multithread = n_cores)
dada_Rs <- dada(fastqs_derep_R, err = err_R, pool = FALSE, multithread = n_cores)

print('...Done.')
Sys.time()

Let’s inspect one of the resulting objects:

dada_Fs[[1]]

---

Step 6: Merge Read Pairs

In this step, we will first merge the forward and reverse read pairs: the fragment that we amplified with our primers was short enough to generate lots of overlap among the sequences from the two directions.

mergers <- mergePairs(dada_Fs, fastqs_derep_F,
                      dada_Rs, fastqs_derep_R,
                      verbose = TRUE)

Just like tables can be saved in R using write.table or write.csv, R objects can be saved using saveRDS. The resulting rds file can then be loaded into an R environment using readRDS. This is a convenient way to save R objects that require a lot of computation time.

We should not be needing the very large dereplicated sequence objects anymore, but to be able to quickly restart our analysis from a new R session if necessary, we now save these objects to rds files. And after that, we can safely remove these objects from our environment.

saveRDS(fastqs_derep_F, file = file.path(outdir, 'fastqs_derep_F.rds'))
saveRDS(fastqs_derep_R, file = file.path(outdir, 'fastqs_derep_R.rds'))

rm(fastqs_derep_F, fastqs_derep_R) # Remove objects from environment

---

Step 7: Construct a Sequence Table

Next, we construct an amplicon sequence variant table (ASV) table:

seqtab_all <- makeSequenceTable(mergers)

# The dimensions of the object are the nr of samples (rows) and the nr of ASVs (columns):
dim(seqtab_all)

Let’s inspect the distribution of sequence lengths:

table(nchar(getSequences(seqtab_all)))

# If you need to remove sequences of a particular length (e.g. too long):
# seqtab2 <- seqtab[, nchar(colnames(seqtab_all)) %in% seq(250,256)]

---

Step 8: Remove Chimeras

Now, we will remove chimeras. The dada algorithm models and removes substitution errors, but chimeras are another importance source of spurious sequences in amplicon sequencing. Chimeras are formed during PCR amplification. When one sequence is incompletely amplified, the incomplete amplicon primes the next amplification step, yielding a spurious amplicon. The result is a sequence read which is half of one sample sequence and half another.

Fortunately, the accuracy of the sequence variants after denoising makes identifying chimeras simpler than it is when dealing with fuzzy OTUs. Chimeric sequences are identified if they can be exactly reconstructed by combining a left-segment and a right-segment from two more abundant “parent” sequences.

seqtab <- removeBimeraDenovo(seqtab_all,
                             method = 'consensus',
                             multithread = n_cores,
                             verbose = TRUE)
ncol(seqtab)

# Proportion of retained sequences:
sum(seqtab) / sum(seqtab_all)

We will save the seqtab object as an rds file:

saveRDS(seqtab, file = file.path(outdir, 'seqtab_V4.rds'))

---

Step 9: Generate a Summary Table

In this step, we will generate a summary table of the number of sequences processed and outputs of different steps of the pipeline.

This information is generally used to further evaluate characteristics and quality of the run, sample-to-sample variation, and resulting sequencing depth for each sample.

To get started, we will define a function getN() that will get the number of unique reads for a sample. Then, we apply getN() to each element of the dada_Fs, dada_Rs, and mergers objects, which gives us vectors with the number of unique reads for each samples, during each of these steps:

getN <- function(x) {
  sum(getUniques(x))
}

denoised_F <- sapply(dada_Fs, getN)
denoised_R <- sapply(dada_Rs, getN)
merged <- sapply(mergers, getN)

We’ll join these vectors together with the “filter_results” dataframe, and the number of nonchimeric reads:

nreads_summary <- data.frame(filter_results,
                             denoised_F,
                             denoised_R,
                             merged,
                             nonchim = rowSums(seqtab),
                             row.names = sampleIDs)
colnames(nreads_summary)[1:2] <- c('input', 'filtered')

# Have a look at the first few rows:
head(nreads_summary)

Finally, we’ll write this table to file:

write.table(nreads_summary, file = file.path(outdir, 'nreads_summary.txt'),
            sep = "\t", quote = FALSE, row.names = TRUE)

---

Step 10: Assign Taxonomy to ASVs

Now, we will assign taxonomy to our ASVs.

Depending on the marker gene and the data, you will have to choose the appropriate reference file for this step. Several files have been formatted for taxonomy assignments in DADA2 pipeline and are available at the DADA2 website.

print('Assigning taxa to ASVs...')
Sys.time()

taxa <- assignTaxonomy(seqtab, tax_key, multithread = n_cores)
colnames(taxa) <- c('Kingdom', 'Phylum', 'Class', 'Order', 'Family', 'Genus')

print('...Done.')
Sys.time()

---

Step 11: Generate Output Files

In this last step, we will generate output files from the DADA2 outputs that are formatted for downstream analysis in phyloseq. First, we will write a fasta file with the final ASV sequences. (This fasta file can also be used for phylogenetic tree inference with different R packages.)

# Prepare sequences and headers:
asv_seqs <- colnames(seqtab)
asv_headers <- paste('>ASV', 1:ncol(seqtab), sep = '_')

# Interleave headers and sequences:
asv_fasta <- c(rbind(asv_headers, asv_seqs))

# Write fasta file:
write(asv_fasta, file = file.path(outdir, 'ASVs.fa'))

Now, we build the final phyloseq object. Notes:

• While we will not add a phylogenetic tree now, this can also be added to a phyloseq object.

• Our metadata dataframe contains three samples that we don’t have sequences for. However, this is not a problem: phyloseq will match the sample IDs in the metadata with those in the OTU table, and disregard IDs not present in the OTU table.

ps <- phyloseq(otu_table(seqtab, taxa_are_rows = FALSE),
               sample_data(metadata_df),
               tax_table(taxa))

# Saves the phyloseq object as an .rds file (which can be imported directly by phyloseq):
saveRDS(ps, file = file.path(outdir, 'ps_V4.rds'))

Report that we are done!

print('Done with ASV inference.')

---

Resources

• Callahan et al. 2006: “Bioconductor Workflow for Microbiome Data Analysis: from raw reads to community analyses”
• dada2 documentation and tutorials
• Taxonomic references for dada2
• cutadapt documentation and tutorials

---

Bonus: Phylogenetic Tree Estimation

A phylogenetic tree can be estimated for the sequence data you generated. Depending on the number of ASVs recovered and the phylogenetic tree algorithm of choice, this step could take several days. Simpler trees will be less computationally intensive. Depending on the marker gene you are working on, you may or may not choose to perform this step.

This step can be conducted after Step 10, and then the phylogeny can be included in the phyloseq object in Step 11.

#seqtab<- readRDS('seqtab_V4.rds')

seqs <- getSequences(seqtab)

# This propagates to the tip labels of the tree.
# At this stage ASV labels are full ASV sequence
names(seqs) <- seqs 
alignment <- AlignSeqs(DNAStringSet(seqs),
                       anchor = NA,
                       iterations = 5,
                       refinements = 5)

print('Computing pairwise distances from ASVs...')
Sys.time()
phang.align <- phyDat(as(alignment, 'matrix'), type = 'DNA')
dm <- dist.ml(phang.align)
treeNJ <- NJ(dm) # Note, tip order is not sequence order
fit = pml(treeNJ, data = phang.align)
print('...Done.')
Sys.time()

print('Fit GTR model...')
Sys.time()
fitGTR <- update(fit, k = 4, inv = 0.2)
print('...Done.')
Sys.time()

print('Computing likelihood of tree...')
Sys.time()
fitGTR <- optim.pml(fitGTR,
                    model = 'GTR',
                    optInv = TRUE,
                    optGamma = TRUE,
                    rearrangement = 'stochastic',
                    control = pml.control(trace = 0))
print('...Done'.)
Sys.time()

---

Bonus: Submit script as an OSC job

Extract R code from this document:

knitr::purl(input = 'markdown/07-ASV-inference.Rmd',
            output = 'scripts/02-ASV-inference.R')

Our script 02-ASV-inference.sh with SLURM directives that will submit the R script from the shell using Rscript:

#!/bin/bash
#SBATCH --nodes=1
#SBATCH --ntasks-per-node=4
#SBATCH --time=2:00:00
#SBATCH --account=PAS0471

module load gnu/9.1.0
module load mkl/2019.0.5
module load R/4.0.2

Rscript scripts/02-ASV-inference.R

Submit the script:

sbatch scripts/02-ASV-inference.sh