This tutorial is aimed at being a walkthrough of the DADA2 pipeline. It uses the data of the now famous MiSeq SOP by the Mothur authors but analyses the data using DADA2.
DADA2 is a relatively new method to analyse amplicon data which uses exact variants instead of OTUs.
The advantages of the DADA2 method is described in the paper
There are two ways to follow this tutorial: you can copy and paste all the codes blocks below in R directly, or you can download this document in the Rmarkdown format and execute the cells.
Install and Load Packages¶
First install DADA2 and other necessary packages
source('https://bioconductor.org/biocLite.R') biocLite('dada2') biocLite('phyloseq') biocLite('DECIPHER') install.packages('ggplot2') install.packages('phangorn')
Now load the packages and verify you have the correct DADA2 version
library(dada2) library(ggplot2) library(phyloseq) library(phangorn) library(DECIPHER) packageVersion('dada2')
Download the Data¶
You will also need to download the data, as well as the SILVA database
If you are following the tutorial on the website, the following block of commands has to be executed outside of R. If you run this tutorial with the R notebook, you can simply execute the cell block
wget http://www.mothur.org/w/images/d/d6/MiSeqSOPData.zip unzip MiSeqSOPData.zip rm -r __MACOSX/ cd MiSeq_SOP wget https://zenodo.org/record/824551/files/silva_nr_v128_train_set.fa.gz wget https://zenodo.org/record/824551/files/silva_species_assignment_v128.fa.gz cd ..
Back in R, check that you have downloaded the data
path <- 'MiSeq_SOP' list.files(path)
Filtering and Trimming¶
First we create two lists with the sorted name of the reads: one for the forward reads, one for the reverse reads
raw_forward <- sort(list.files(path, pattern="_R1_001.fastq", full.names=TRUE)) raw_reverse <- sort(list.files(path, pattern="_R2_001.fastq", full.names=TRUE)) # we also need the sample names sample_names <- sapply(strsplit(basename(raw_forward), "_"), `[`, # extracts the first element of a subset 1)
then we visualise the quality of our reads
What do you think of the read quality?
The forward reads are good quality (although dropping a bit at the end as usual) while the reverse are way worse.
Based on these profiles, we will truncate the forward reads at position 240 and the reverse reads at position 160 where the quality distribution crashes.
in this tutorial we perform the trimming using DADA2's own functions. If you wish to do it outside of DADA2, you can refer to the Quality Control tutorial
Dada2 requires us to define the name of our output files
# place filtered files in filtered/ subdirectory filtered_path <- file.path(path, "filtered") filtered_forward <- file.path(filtered_path, paste0(sample_names, "_R1_trimmed.fastq.gz")) filtered_reverse <- file.path(filtered_path, paste0(sample_names, "_R2_trimmed.fastq.gz"))
We’ll use standard filtering parameters:
maxN=0 (DADA22 requires no Ns),
The maxEE parameter sets the maximum number of “expected errors” allowed in a read, which according to the USEARCH authors is a better filter than simply averaging quality scores.
out <- filterAndTrim(raw_forward, filtered_forward, raw_reverse, filtered_reverse, truncLen=c(240,160), maxN=0, maxEE=c(2,2), truncQ=2, rm.phix=TRUE, compress=TRUE, multithread=TRUE) head(out)
Learn the Error Rates¶
The DADA2 algorithm depends on a parametric error model and every amplicon dataset has a slightly different error rate.
learnErrors of Dada2 learns the error model from the data and will help DADA2 to fits its method to your data
errors_forward <- learnErrors(filtered_forward, multithread=TRUE) errors_reverse <- learnErrors(filtered_reverse, multithread=TRUE)
then we visualise the estimated error rates
plotErrors(errors_forward, nominalQ=TRUE) + theme_minimal()
Do you think the error model fits your data correctly?
From the Dada2 documentation:
Dereplication combines all identical sequencing reads into into “unique sequences” with a corresponding “abundance”: the number of reads with that unique sequence. Dereplication substantially reduces computation time by eliminating redundant comparisons.
derep_forward <- derepFastq(filtered_forward, verbose=TRUE) derep_reverse <- derepFastq(filtered_reverse, verbose=TRUE) # name the derep-class objects by the sample names names(derep_forward) <- sample_names names(derep_reverse) <- sample_names
We are now ready to apply the core sequence-variant inference algorithm to the dereplicated data.
dada_forward <- dada(derep_forward, err=errors_forward, multithread=TRUE) dada_reverse <- dada(derep_reverse, err=errors_reverse, multithread=TRUE) # inspect the dada-class object dada_forward[]
The DADA2 algorithm inferred 128 real sequence variants from the 1979 unique sequences in the first sample.
Merge Paired-end Reads¶
Now that the reads are trimmed, dereplicated and error-corrected we can merge them together
merged_reads <- mergePairs(dada_forward, derep_forward, dada_reverse, derep_reverse, verbose=TRUE) # inspect the merger data.frame from the first sample head(merged_reads[])
Construct Sequence Table¶
We can now construct a sequence table of our mouse samples, a higher-resolution version of the OTU table produced by traditional methods.
seq_table <- makeSequenceTable(merged_reads) dim(seq_table) # inspect distribution of sequence lengths table(nchar(getSequences(seq_table)))
dada method used earlier removes substitutions and indel errors but chimeras remain.
We remove the chimeras with
seq_table_nochim <- removeBimeraDenovo(seq_table, method='consensus', multithread=TRUE, verbose=TRUE) dim(seq_table_nochim) # which percentage of our reads did we keep? sum(seq_table_nochim) / sum(seq_table)
As a final check of our progress, we’ll look at the number of reads that made it through each step in the pipeline
get_n <- function(x) sum(getUniques(x)) track <- cbind(out, sapply(dada_forward, get_n), sapply(merged_reads, get_n), rowSums(seq_table), rowSums(seq_table_nochim)) colnames(track) <- c('input', 'filtered', 'denoised', 'merged', 'tabled', 'nonchim') rownames(track) <- sample_names head(track)
We kept the majority of our reads!
Now we assign taxonomy to our sequences using the SILVA database
taxa <- assignTaxonomy(seq_table_nochim, 'MiSeq_SOP/silva_nr_v128_train_set.fa.gz', multithread=TRUE) taxa <- addSpecies(taxa, 'MiSeq_SOP/silva_species_assignment_v128.fa.gz')
for inspecting the classification
taxa_print <- taxa # removing sequence rownames for display only rownames(taxa_print) <- NULL head(taxa_print)
DADA2 is reference-free so we have to build the tree ourselves
We first align our sequences
sequences <- getSequences(seq_table) names(sequences) <- sequences # this propagates to the tip labels of the tree alignment <- AlignSeqs(DNAStringSet(sequences), anchor=NA)
Then we build a neighbour-joining tree then fit a maximum likelihood tree using the neighbour-joining tree as a starting point
phang_align <- phyDat(as(alignment, 'matrix'), type='DNA') dm <- dist.ml(phang_align) treeNJ <- NJ(dm) # note, tip order != sequence order fit = pml(treeNJ, data=phang_align) ## negative edges length changed to 0! fitGTR <- update(fit, k=4, inv=0.2) fitGTR <- optim.pml(fitGTR, model='GTR', optInv=TRUE, optGamma=TRUE, rearrangement = 'stochastic', control = pml.control(trace = 0)) detach('package:phangorn', unload=TRUE)
First load the metadata
sample_data <- read.table( 'https://hadrieng.github.io/tutorials/data/16S_metadata.txt', header=TRUE, row.names="sample_name")
We can now construct a phyloseq object from our output and newly created metadata
physeq <- phyloseq(otu_table(seq_table_nochim, taxa_are_rows=FALSE), sample_data(sample_data), tax_table(taxa), phy_tree(fitGTR$tree)) # remove mock sample physeq <- prune_samples(sample_names(physeq) != 'Mock', physeq) physeq
Let's look at the alpha diversity of our samples
plot_richness(physeq, x='day', measures=c('Shannon', 'Fisher'), color='when') + theme_minimal()
No obvious differences. Let's look at ordination methods (beta diversity)
We can perform an MDS with euclidean distance (mathematically equivalent to a PCA)
ord <- ordinate(physeq, 'MDS', 'euclidean') plot_ordination(physeq, ord, type='samples', color='when', title='PCA of the samples from the MiSeq SOP') + theme_minimal()
now with the Bray-Curtis distance
ord <- ordinate(physeq, 'NMDS', 'bray') plot_ordination(physeq, ord, type='samples', color='when', title='PCA of the samples from the MiSeq SOP') + theme_minimal()
There we can see a clear difference between our samples.
Let us take a look a the distribution of the most abundant families
top20 <- names(sort(taxa_sums(physeq), decreasing=TRUE))[1:20] physeq_top20 <- transform_sample_counts(physeq, function(OTU) OTU/sum(OTU)) physeq_top20 <- prune_taxa(top20, physeq_top20) plot_bar(physeq_top20, x='day', fill='Family') + facet_wrap(~when, scales='free_x') + theme_minimal()
We can place them in a tree
bacteroidetes <- subset_taxa(physeq, Phylum %in% c('Bacteroidetes')) plot_tree(bacteroidetes, ladderize='left', size='abundance', color='when', label.tips='Family')