Inspect quality scores

Now that primers are removed, we can begin the DADA2 pipeline in full. We start with plotting quality profiles of our sequences to determine suitable quality filtering parameters. For the sake of speed we only show 2 samples here but feel free to look at all your samples (or more than 2 anyway).

plotQualityProfile(fwd_cut[1:2]) + ggtitle("Forward")

Everything looks pretty good here. The majority of our reads should be 230 bp based on the fact that we trimmed 20 bp primers off the beginning of the read. It looks like a very, very few did not get trimmed properly but from the red line across the bottom, which is the number of reads of that length, we can see that this is such a small percentage that we don’t need to worry about it.

plotQualityProfile(rev_cut[1:2]) + ggtitle("Reverse")

Reverse reads are lower quality as is always the case for Illumina data. However, it looks like the quality average is above 20 for most of the length of the read so we’ll go ahead with the analysis.

Filter reads

For the purposes of our analysis, we chose to conduct quality filtering with maximum expected errors of 2 for the forward sequence, 5 for the reverse, to truncate after a quality score of 2 or lower.

# Same as for the clippling we create an output directory to store the filtered files
filt_dir <- file.path(path, "filtered")
if (!dir.exists(filt_dir)) dir.create(filt_dir)

fwd_filt <- file.path(filt_dir, basename(fwd_files))
rev_filt <- file.path(filt_dir, basename(rev_files))

names(fwd_filt) <- samples
names(rev_filt) <- samples

filtered_out <- filterAndTrim(
  fwd = fwd_cut, 
  filt = fwd_filt,
  rev = rev_cut,
  filt.rev = rev_filt,
  maxEE = c(2, 5), 
  truncQ = 2, 
  rm.phix = TRUE, 
  compress = TRUE, 
  multithread = TRUE
  )  

head(filtered_out)

##                            reads.in reads.out
## G1_S1_L001_R1_001.fastq.gz    25276     24110
## G2_S2_L001_R1_001.fastq.gz    25908     24533

Learn errors

Here dada2 learns the error profile for your data using an interative approach. These error profiles are used in the next step to correct errors.

err_fwd <- learnErrors(fwd_filt, multithread = TRUE)

## 11163055 total bases in 48643 reads from 2 samples will be used for learning the error rates.

err_rev <- learnErrors(rev_filt, multithread = TRUE)

## 11118642 total bases in 48643 reads from 2 samples will be used for learning the error rates.

plotErrors(err_fwd, nominalQ = TRUE)

The plot shows how well the estimated error rates fit the observed rates. The black line should line up reasonably well with the black points. This looks reasonable so we continue on.

Denoising (aka error-correcting)

We can now denoise our sequences using the DADA algorithm. Note that older versions of dada2 required a dereplication step that is now done automatically by the dada function.

dada_fwd <- dada(fwd_filt, err = err_fwd, multithread = TRUE)

## Sample 1 - 24110 reads in 4919 unique sequences.
## Sample 2 - 24533 reads in 5351 unique sequences.

dada_rev <- dada(rev_filt, err = err_rev, multithread = TRUE)

## Sample 1 - 24110 reads in 15359 unique sequences.
## Sample 2 - 24533 reads in 15454 unique sequences.

Merge pairs

Now that we have accurate sequences we can merge the paired-end reads to reconstruct the full amplicon sequence. Note that most sequences should succesfully overlap, unless your trimming parameters were too aggresive (i.e. too much sequence was trimmed off the end) or you have amplicons that longer than the combined read length (>500 bp, in this case). In the latter case these reads are still useable but will need special handling (more details to come on those cases).

mergers <- mergePairs(
  dadaF = dada_fwd,
  dadaR = dada_rev,
  derepF = fwd_filt,
  derepR = rev_filt,
  maxMismatch = 1, 
  verbose=TRUE
)

Sequence table

Before we can assign taxonomy, we can construct a sequence table, which lists the number of each sequence present in each sample. We can see the dimensions of this sequence table, where the number of rows is the number of samples and the number of columns is the number of total unique sequences

seqtab <- makeSequenceTable(mergers)
dim(seqtab) 

## [1]  2 55

Remove chimeras

Don’t be alarmed here if many of your sequences are discarded - the number of chimeras varies but can be quite high in some cases.

seqtab_nochim <- removeBimeraDenovo(seqtab, method = "consensus", multithread = TRUE, verbose = TRUE)
dim(seqtab_nochim)

## [1]  2 19

Sequence length distribution

table(nchar(getSequences(seqtab_nochim)))

## 
## 273 274 281 283 287 288 289 291 
##   2   2   2   2   2   1   1   7

Number of reads at each step

This is an important step, particularly if you’re trying to track down where something went wrong.

# small function to get the number of sequences
getN <- function(x) sum(getUniques(x))

track <- cbind(
  filtered_out, 
  sapply(dada_fwd, getN), 
  sapply(dada_rev, getN), 
  sapply(mergers, getN), 
  rowSums(seqtab_nochim)
)

colnames(track) <- c("raw", "filtered", "denoised_fwd", "denoised_rev", "merged", "no_chim")
rownames(track) <- samples  
head(track)

##      raw filtered denoised_fwd denoised_rev merged no_chim
## G1 25276    24110        24074        21483  21254   20304
## G2 25908    24533        24479        21994  21681   19947

Training the classifier

The IDTAXA classification method relies on a training set that contains sequences that are representative of known taxa. Although there exists trained classifiers, we will be using the ITS2 Nematode database as our classifier. Since it is not trained, we will have to train the classifier using LearnTaxa, which is also available as a DECIPHER package. The train parameter of LearnTaxa contains the sequences included in the classifier as a DNAStringSet. Each sequence is formatted such that the names of each taxonomic level is separated by a semicolon. Once the sequence headers of the classifier are parsed, the taxonomy parameter will contain the taxonomic assignment, separated by a semicolon, for each sequence in train.

train <- readDNAStringSet("~/Nematode_ITS2_1.0.0_idtaxa.fasta")
tax <- read_tsv("~/Nematode_ITS2_1.0.0_idtaxa.tax")


# Train the classifier
trainingSet <- LearnTaxa(train, names(train), tax)

## ================================================================================
## 
## Time difference of 131.39 secs

Running IdTaxa

To use IdTaxa, we will need to convert the unique sequence variants that we want to classify into a DNAStringSet object. Now that we have our trained classifier (i.e. training set) and the sequences we want to classify, we can run IdTaxa. You can change various parameters of IdTaxa depending on your data. We recommend reviewing the [IdTaxa documentation] (https://rdrr.io/bioc/DECIPHER/man/IdTaxa.html) for further information. Here, we have specified several parameters:

• strand = "both" : by setting the strand parameter to “both”, each sequence variant is classified using both the forward and reverse complement orientation. The sequence varient will be classified as the result with the highest confidence.
• threshold = 60 : setting the threshold to 60 specifies when the taxonomic classifications are truncated. A lower threshold usually results in higher taxonomic level classifications, but lower accuracy (i.e. confidence). A higher threshold usually results in lower taxonomic level classifications, but with higher accuracy.
• bootstraps = 100 : this specifies the amount of times bootstrap replicates are performed for each sequence.
• processors = NULL : automatically uses all available processors.
• verbose = TRUE : displays progress.
• type = "extended" : by setting the type to “extended”, the output for each sequence variant will contain the taxonomic classification, the names of each taxonomic rank, and the confidence of the assignment.

dna <- DNAStringSet(getSequences(seqtab_nochim))

idtaxa <- IdTaxa(dna, 
                 trainingSet, 
                 strand = "both", 
                 threshold = 60, 
                 bootstraps = 100, 
                 processors = NULL, 
                 verbose = TRUE, 
                 type = "extended")

## ================================================================================
## 
## Time difference of 0.56 secs

Although it contains all the necessary information, the output of IdTaxa is not as intuitive as we would like. It outputs a list of the sequence variants. Within this list, each element will have its taxonomic classification and the confidence for each taxon. To facilitate easier manipulation and comprehension of the data, we will convert this list into a matrix, with the rows as the sequence variants and the columns as their taxonomic classifications. We also drop the “Root” column.

taxid <- t(sapply(idtaxa, function(x) setNames(x$taxon, x$rank)))[, -1]

Over to phyloseq

The individual pieces of the analysis can be save to text or excel files but if you’ve made it this far then why not continue on in R! Phyloseq is a powerful toolkit for working with these types of data in R and is the recommended next step for analysis. More tutorials on downstream analysis to come!

library(phyloseq)
# This is just a placeholder - normally this would be a dataframe containing all your sample inforamtion,
# with rownames matching the sample names
samp_data <-  data.frame(
  row.names = samples,
  sample = samples
)

# We need better tables for our sequences than the actual sequence which is the dada2 default
asvs <- paste0("ASV_", 1:length(dna))
rownames(taxid) <- asvs
colnames(seqtab_nochim) <- asvs
names(dna) <- asvs

# Now we put everything into one phyloseq object (even the sequences) so it is easy to use
physeq = phyloseq(
  otu_table(seqtab_nochim, taxa_are_rows = FALSE),
  tax_table(taxid),
  sample_data(samp_data),
  dna
)

physeq

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 19 taxa and 2 samples ]
## sample_data() Sample Data:       [ 2 samples by 1 sample variables ]
## tax_table()   Taxonomy Table:    [ 19 taxa by 8 taxonomic ranks ]
## refseq()      DNAStringSet:      [ 19 reference sequences ]