Dada2 Filter and Trim on Forward Reads

Contents

• 1 Introduction
• 2 Overview of the dada2 pipeline
• 3 Filter and Trim
• four Dereplicate
• 5 Learn the mistake rates
• 6 Infer sample composition
• 7 Merge forward/reverse reads
• eight Remove chimeras
• 9 A second sample
• ten Create sequence table

More detailed documentation is available at the DADA2 Domicile Folio. In particular, the online tutorial workflow is the virtually detailed and up-to-date demonstration of applying DADA2 to multi-sample amplicon datasets.

Introduction

The investigation of environmental microbial communities and microbiomes has been revolutionized past the development of high-throughput amplicon sequencing. In amplicon sequencing a particular genetic locus, for example the 16S rRNA gene in bacteria, is amplified from DNA extracted from the community of interest, and then sequenced on a side by side-generation sequencing platform. This technique removes the need to culture microbes in order to discover their presence, and cost-effectively provides a deep census of a microbial community.

Yet, the procedure of amplicon sequencing introduces errors into the sequencing data, and these errors severely complicate the interpretation of the results. DADA2 implements a novel algorithm that models the errors introduced during amplicon sequencing, and uses that error model to infer the truthful sample composition. DADA2 replaces the traditional "OTU-picking" pace in amplicon sequencing workflows, producing instead higher-resolution tables of amplicon sequence variants (ASVs).

As seen in the newspaper introducing DADA2 and in further benchmarking, the DADA2 method is more sensitive and specific than traditional OTU methods: DADA2 detects existent biological variation missed by OTU methods while outputting fewer spurious sequences. Some other contempo paper describes how replacing OTUs with ASVs improves the precision, comprehensiveness and reproducibility of marker-cistron data analysys.

Overview of the dada2 pipeline

The starting signal for the dada2 pipeline is a set of demultiplexed fastq files corresponding to the samples in your amplicon sequencing study. That is, dada2 expects there to be an private fastq file for each sample (or two fastq files, one forward and ane reverse, for each sample). Demultiplexing is often performed at the sequencing center, but if that has not been done at that place are a variety of tools practise accomplish this, including the popular QIIME python scripts split_libraries_fastq.py followed by split_sequence_file_on_sample_ids.py, and the utility idemp, amid others.

In addition, dada2 expects that there are no not-biological bases in the sequencing information. This may crave pre-processing with outside software if, as is relatively common, your PCR primers were included in the sequenced amplicon region.

One time demultiplexed fastq files without not-biological nucleotides are in hand, the dada2 pipeline proceeds every bit follows:

1. Filter and trim: filterAndTrim()
2. Dereplicate: derepFastq()
3. Acquire mistake rates: learnErrors()
4. Infer sample composition: dada()
5. Merge paired reads: mergePairs()
6. Make sequence table: makeSequenceTable()
7. Remove chimeras: removeBimeraDenovo()

The output of the dada2 pipeline is a feature table of amplicon sequence variants (an ASV table): A matrix with rows corresponding to samples and columns to ASVs, in which the value of each entry is the number of times that ASV was observed in that sample. This table is analogous to the traditional OTU table, except at higher resolution, i.e exact amplicon sequence variants rather than (usually 97%) clusters of sequencing reads.

We at present go through the pipeline on a highly simplified dataset of just one paired-end sample (we'll add a second later).

Filter and Trim

We'll commencement past getting the filenames of our example paired-terminate fastq files. Normally you lot volition ascertain these filenames directly, or read them out of a directory, but for this tutorial we're using files included with the packet, which we can identify via a particular function call:

            library(dada2); packageVersion("dada2")          

            ## [ane] 'i.22.0'          

            fnF1 <- system.file("extdata", "sam1F.fastq.gz", packet="dada2") fnR1 <- system.file("extdata", "sam1R.fastq.gz", package="dada2") filtF1 <- tempfile(fileext=".fastq.gz") filtR1 <- tempfile(fileext=".fastq.gz")          

Note that the dada2 package "speaks" the gzip format natively, therefore all fastq files tin remain in the infinite-saving gzip format throughout.

At present that we take the filenames, we're going to inspect the quality profile of our data:

            plotQualityProfile(fnF1) # Forward          

            ## Warning: `guides(<scale> = FALSE)` is deprecated. Please utilise `guides(<scale> = ## "none")` instead.          

            plotQualityProfile(fnR1) # Opposite          

            ## Warning: `guides(<scale> = Imitation)` is deprecated. Delight use `guides(<calibration> = ## "none")` instead.          

Later on inspecting the quality profiles, it is clear that the opposite read quality drops off more severely than in the forrad read. Thus we are going to trim the forward reads at position 240, and the reverse reads at position 200.

Filtering is an important step when dealing with sequence data, as low-quality sequences can incorporate unexpected and misleading errors. Trimming is also usually advised, every bit Illumina sequencing quality tends to drop off at the stop of reads, and the initial nucleotides can also exist problematic due to calibration issues:

            filterAndTrim(fwd=fnF1, filt=filtF1, rev=fnR1, filt.rev=filtR1,                   trimLeft=ten, truncLen=c(240, 200),                    maxN=0, maxEE=2,                   shrink=TRUE, verbose=TRUE)          

            ## Read in 1500 paired-sequences, output 890 (59.iii%) filtered paired-sequences.          

The filterAndTrim(...) function filters the frontward and reverse reads jointly, outputting only those pairs of reads that both pass the filter. In this part call we did four things: We removed the starting time trimLeft=10 nucleotides of each read. We truncated the frontward and reverse reads at truncLen=c(240, 200) nucleotides respectively. We filtered out all reads with more than than maxN=0 cryptic nucleotides. And nosotros filtered out all reads with more than two expected errors. The filtered output files were stored as gzipped fastq files (shrink=Truthful).

This represents a adequately standard set of filtering/trimming parameters. Still, it is ever worth evaluating whether the filtering and trimming parameters yous are using are appropriate for your data. Ane size does not fit all! (And are you sure y'all have removed your primers?)

An important consideration: If using paired-finish sequencing information, you lot must maintain a suitable overlap (>20nts) betwixt the forward and reverse reads after trimming! This is peculiarly important to go along in mind for mult-V-region amplicions (such equally V3-V4) in which there may exist relatively little overlap to begin with, and thus piddling read-truncation is possible if reads are to be merged later on.

Dereplicate

The next matter nosotros desire to practice is "dereplicate" the filtered fastq files. During dereplication, nosotros condense the data past collapsing together all reads that encode the same sequence, which significantly reduces after ciphering times:

            derepF1 <- derepFastq(filtF1, verbose=TRUE)          

            ## Dereplicating sequence entries in Fastq file: /tmp/Rtmp70itev/file391d6a52e201db.fastq.gz          

            ## Encountered 332 unique sequences from 890 full sequences read.          

            derepR1 <- derepFastq(filtR1, verbose=TRUE)          

            ## Dereplicating sequence entries in Fastq file: /tmp/Rtmp70itev/file391d6a226e435.fastq.gz          

            ## Encountered 485 unique sequences from 890 total sequences read.          

Dereplication is a mutual step in about all modern ASV inference (or OTU picking) pipelines, but a unique feature of derepFastq is that information technology maintains a summary of the quality information for each dereplicated sequence in $quals.

Larn the error rates

The dada algorithm uses a parametric model of the errors introduced by PCR amplification and sequencing. Those error parameters typically vary betwixt sequencing runs and PCR protocols, so our method provides a way to estimate those parameters from the data itself.

            errF <- learnErrors(derepF1, multithread=Imitation) # multithreading is bachelor on many functions          

            ## 204700 total bases in 890 reads from ane samples volition be used for learning the error rates.          

            errR <- learnErrors(derepR1, multithread=False)          

            ## 169100 total bases in 890 reads from 1 samples will be used for learning the error rates.          

We recommend using at least a subest of your data to learn your mistake rates for the nearly accurate sample inference.

Infer sample composition

The core method of the dada2 package is at the sample inference stage. The dada(...) function implements the algorithm described in our paper, and is simultaneously more than sensitive and more specific than whatever OTU algorithm we accept e'er tested.

Here we telephone call the dada(...) role, using the models of the error rates wer learned in the previous stride:

            dadaF1 <- dada(derepF1, err=errF, multithread=Fake)          

            ## Sample 1 - 890 reads in 332 unique sequences.          

            dadaR1 <- dada(derepR1, err=errR, multithread=FALSE)          

            ## Sample 1 - 890 reads in 485 unique sequences.          

            impress(dadaF1)          

            ## dada-class: object describing DADA2 denoising results ## 8 sequence variants were inferred from 332 input unique sequences. ## Cardinal parameters: OMEGA_A = 1e-twoscore, OMEGA_C = 1e-xl, BAND_SIZE = 16          

The dada(...) algorithm inferred viii sequence variants from the forward reads.

Merge forward/contrary reads

Nosotros've inferred the sample sequences in the frontward and reverse reads independently. Now information technology's time to merge those inferred sequences together, throwing out those pairs of reads that don't match:

            merger1 <- mergePairs(dadaF1, derepF1, dadaR1, derepR1, verbose=True)          

            ## 834 paired-reads (in 7 unique pairings) successfully merged out of 882 (in 10 pairings) input.          

The mergePairs(...) function returns a data.frame corresponding to each successfully merged unique sequence. The $forward and $reverse columns record which forrad and reverse sequence contributed to that merged sequence.

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 ane sequence is incompletely amplified, the incomplete amplicon primes the next amplification stride, yielding a spurious amplicon. The result is a sequence read which is half of ane sample sequence and one-half another.

We place and remove those sequence using the removeBimeraDenovo(...) part in the dada2 pipeline:

            merger1.nochim <- removeBimeraDenovo(merger1, multithread=FALSE, verbose=TRUE)          

            ## Identified 0 bimeras out of 7 input sequences.          

We now have a data.frame of merged, error-complimentary, non-chimeric, amplicon sequence variants!

A 2nd sample

In order to show an instance of making a sequence table, and to reiterate the workflow outlined above, we now process a second sample:

            # Assign filenames fnF2 <- organization.file("extdata", "sam2F.fastq.gz", package="dada2") fnR2 <- system.file("extdata", "sam2R.fastq.gz", package="dada2") filtF2 <- tempfile(fileext=".fastq.gz") filtR2 <- tempfile(fileext=".fastq.gz") # Filter and Trim filterAndTrim(fwd=fnF2, filt=filtF2, rev=fnR2, filt.rev=filtR2, maxN=0, trimLeft=ten, truncLen=c(240, 200), maxEE=2, compress=TRUE, verbose=TRUE)          

            ## Read in 1500 paired-sequences, output 858 (57.2%) filtered paired-sequences.          

            # Dereplicate derepF2 <- derepFastq(filtF2, verbose=TRUE)          

            ## Dereplicating sequence entries in Fastq file: /tmp/Rtmp70itev/file391d6a395a2793.fastq.gz          

            ## Encountered 305 unique sequences from 858 full sequences read.          

            derepR2 <- derepFastq(filtR2, verbose=True)          

            ## Dereplicating sequence entries in Fastq file: /tmp/Rtmp70itev/file391d6a1bee25fd.fastq.gz          

            ## Encountered 448 unique sequences from 858 total sequences read.          

            # Infer sample composition (using already learned mistake rates) dadaF2 <- dada(derepF2, err=errF, multithread=Simulated)          

            ## Sample 1 - 858 reads in 305 unique sequences.          

            dadaR2 <- dada(derepR2, err=errR, multithread=FALSE)          

            ## Sample 1 - 858 reads in 448 unique sequences.          

            # Merge merger2 <- mergePairs(dadaF2, derepF2, dadaR2, derepR2, verbose=True)          

            ## 731 paired-reads (in vi unique pairings) successfully merged out of 852 (in 10 pairings) input.          

With that second sample processed, we tin go ahead and create a sequence table.

Create sequence table

Finally we desire to combine the inferred samples into one unified table. For that purpose we use makeSequenceTable:

            seqtab <- makeSequenceTable(list(merger1, merger2)) seqtab.nochim <- removeBimeraDenovo(seqtab, verbose=TRUE)          

            ## Identified 0 bimeras out of 7 input sequences.          

            dim(seqtab.nochim)          

            ## [one] 2 7          

This is the final product of the dada2 pipeline, a matrix in which each row corresponds to a candy sample, and each column corresponds to an non-chimeric inferred sample sequence (a more precise analogue to the common "OTU tabular array"). From here we recommend proceeding forward with our friend the phyloseq package for further assay…

carawayuness2001.blogspot.com

Source: https://www.bioconductor.org/packages/release/bioc/vignettes/dada2/inst/doc/dada2-intro.html

Share this post