By Kelly Shannon<\/p>\n\n\n\n
###This pipeline is adapted from the benjjneb dada2 pipeline found on GitHub: https:\/\/benjjneb.github.io\/dada2\/tutorial.html<\/a> and based on Callahan et al., 2016. DADA2: High Resolution Sample Inference from Illumina Amplicon Data. Nature Methods. As well as from one created by MK English in Dr. Ryan Mueller\u2019s lab at OSU. ##The pipeline will be performed with 14 total paired-end sequence fastq files (7 individual samples) as an example flowthrough. These are V4 hypervariable region amplicon sequences performed by the Oregon State University CGRB on an Illumina MiSeq NGS sequencer.<\/p>\n\n\n\n\n\n\n\n #Packages to load: I like to go ahead and load all of the packages I think that I will need before I start. If you do not yet have a package in your library, the command to install it is <install.packages(\u201cname_of_package\u201d). You may not need all of these but it doesn\u2019t hurt to load them regardless in case you do, since they are all commonly used packages to complement dada2.<\/p>\n\n\n\n #First things first, we want to set a path for R to find our paired-end reads, which are the fastq sequences that were recieved from the sequencing facility. This will be the same path as where they are located on your desktop, and can easily be copy and pasted from moving into the directory where they are located on the command line and copy pasting the output of . What I like to do is make a directory that contains only my fastq files of interest, that way the path to them will only display those files. We can store this path in a vector.<\/p>\n\n\n\n The output from this should show your fastq files, and maybe other directories that you have within this directory, that\u2019s fine too. Your files should end in \u201c.fastq\u201d<\/p>\n\n\n\n #We now want to create individual vectors containing just the forward and reverse reads with the <sort()> command. The pattern we put into this can really be anything, but it MUST match the common wording of your forward and reverse reads. For example, I had renamed mine before this pipeline and each forward read follows the pattern: \u201csamplename_F.fastq\u201d, where the sample name is unique to each sample, but they all end in \u201cF.fastq\u201d. Therefore, the pattern will be \u201cF.fastq\u201d; same deal for the reverse reads.<\/p>\n\n\n\n #Now that we have vectors made for our forward and reverse reads we can extract the sample names from them, again this is assuming that they have the file names following the structure of \u201csamplename_X.fastq\u201d.<\/p>\n\n\n\n The output of the vector should show only the names of your samples. #https:\/\/www.bioinformatics.babraham.ac.uk\/projects\/fastqc\/<\/a> #https:\/\/multiqc.info\/<\/a><\/p>\n\n\n\n #We\u2019ll continue with the DADA2 way of visualizing read qualities because the pipeline really won\u2019t change whether you are using the DADA2 <plotQualityProfile()> command, or multiqc output. This will output the total # of reads in each fastq file, as well as showing us where the reads drop off in quality. In general, the forward reads will be much better than the reverse, meaning the quality of the reads will not drop off until towards the end of the reads.<\/p>\n\n\n\n The mean quality score of all of the reads in each file is shown by the green line, so this is the one to pay attention to. You will generally want to trim off the last little section of nucleotides where the avg. read quality drops off, hence why we analyze these plots.<\/p>\n\n\n\n #First lets assign the file for what will be the filtered fastq.gz files, not that we will zip these so that they take up less space on our computer. The reason this is done is so that the output of the filtering command (the step after the following) can be assigned to file names, separate from our forward and reverse sequence vectors that we have already made.<\/p>\n\n\n\n That last step is so that we can still apply our sample names to the filtered reads.<\/p>\n\n\n\n #Okay now it\u2019s time to actually perform the filtering of our reads. We do this with the <filterAndTrim()> command. Now, there are tons of parameters that can be specified with this command and they need to be precise for your samples or else it may lead you to either drop too many reads (most likely outcome) and get too few eventual ASVs or keep too many reads and get spurious ASVs (happens less of the time). I\u2019d highly recommend reading the documentation of this command(I know scary right?) so that you can tailor it to your experiment, but this does take some practice and know-how.<\/p>\n\n\n\n A couple notes: the truncLen command should be set to where in your forward and reverse reads the avg. base quality steeply drops in quality. Depending on how strict you want to be, this means where the reads drop below ~25 on the y-axis (~99.5% confidence per nucleobase). For my values, my forward reads maintained high Phred Scores (above 30) all the way through so I truncated at 250 (where V4 reads should end). My reverse reads are of much lower quality, which is typical, so I truncated these at a lower value, where the Phred Scores drop below about 30. The values of these parameters are going to be highly specific to your own analysis, but you can also play it safe and go with the default values listed on the official DADA2 pipeline.<\/p>\n\n\n\n My reverse reads are pretty garbage so I will truncate those at 160. REMEMBER: there must be at least a 20 nt overlap between your forward and reverse reads for them to merge into contigs; obviously this won\u2019t be a problem for mine since I didn\u2019t truncate my forward reads, however.<\/p>\n\n\n\n We can view the output at the end so we can keep track of how many reads were filtered. This is a good way to learn how changing parameters influences the amount of outputted reads.<\/p>\n\n\n\n #It will speed things up down the line if we dereplicate our reads, meaning collapsing all identical reads.<\/p>\n\n\n\n #Now that we have filtered our reads, we want to learn the error rates of our dataset, which will be specific to the actual sequencing run. This command should therefore be run only on data that has all been sequenced on the same run (i.e. the same MiSeq run). This is looking for incorrect nucleotide substitutions as a result of sequencing errors, and will be performed on the forward and reverse reads separately. This command may take some time, depending on how many samples you have.<\/p>\n\n\n\n This will only use a subset of your reads to learn the error rates, so don\u2019t be alarmed if you don\u2019t see all of your reads being used for learning the error rates.<\/p>\n\n\n\n You can also visualize plots of your error rates with <plotErrors()>, which is a sanity check and should display a decrease in error frequency (y-axis) concurrent with the increase in per-base quality score (x-axis). This step is somewhat optional, but is a good check that the learnErrors command worked.<\/p>\n\n\n\n #We can now apply the \u201ccore sample inference algorithm\u201d to the filtered and dereplicated with the learned error rates. This is a pretty complicated algorithm that uses the error rates gathered from the command above to quantify the rate at which an amplicon read is produced from a sample sequence as a function of sequence composition and quality (don\u2019t worry if you don\u2019t understand what this means, I don\u2019t really either). The number of repeated observations of the sample sequence is used to calculate a p-value that the number of observed sequences is consistent with the error model. Low p-values mean that the number of sample sequences is too many than could be explained by sequencing errors, or in more simple terms, these are \u201ctrue biological sequences\u201d rather than spurious ones. The default p-value is 1e-40, which is REALLY low, but is the value recommended by the pipeline. You can change this, but it is not necessarily recommended to do so because you can be very sure that any reads that pass this p-value are true biological ones. Lowering this value could potentially allow spurious reads to pass through (but again, this depends on how strict you want to be in your analysis).<\/p>\n\n\n\n I\u2019ve changed up a couple major parameters above from the default pipeline. 1. I\u2019m keeping singletons because they make up such a large portion of my sequences, and I do not want to drop rare sequences. FAIR WARNING, this is controversial and not recommended by many because it is argued that singletons have higher chances of being erroneous than other sequences. I would like to maintain as many rare sequences as possible for my analysis so I am keeping them, but you may prefer to not take the risk and drop them. 2. I lowered the p-value, which dictates the threshold at which DADA2 calls sequences as being overabundant and uses them to base ASV assignments of other similar sequences. This is also not necessarily recommended because it can increase the chances of false positives, but again I\u2019d like to try and get my rarer sequences through. 3. I set pool equal to \u201cpseudo\u201d. The default is pool = FALSE, which will be much less computationally intensive, and will treat each sample independently. Pooling compares samples and takes this information in when performing denoising, which favors rare sequences that may have higher representation in more deeply sequenced samples than samples that have less sequences. Pseudo is sort of a hybrid between these two methods that is less computationally intensive than pool = TRUE. Phew, we got through that explanation.<\/p>\n\n\n\n #Now we are ready to merge our forward and reverse reads and make contigs! We\u2019ll also make a sequence table to display our formed contigs.<\/p>\n\n\n\n #Now we\u2019ll remove chimeric sequences from our contigs and track the amount of sequences we get at each contig size, and the amount that we lose from chimeras.<\/p>\n\n\n\n Find the percentage of chimeric sequences by dividing the total sequences by the chimeric ones. This percentage is generally below 1%, if it is much higher then this then spurious reads may have escaped through any one of the above filtering steps in the pipeline.<\/p>\n\n\n\n #Since these are 16S V4 amplicons, they should be around 250-255-ish base pairs in length. If there are sequences outside of this general region, as shown by the table() command, then they can be trimmed away. You can see how many unique sequences you have after all of this too.<\/p>\n\n\n\n Note that 250-255 is not a set-in-stone range of values. It\u2019s possible that your sequences may be outside of this range and this does not necessarily mean that they aren\u2019t true biological sequences (all kinds of things could cause this from the actual sequencing). However, for my data, this is the range that all of my contigs fall into.<\/p>\n\n\n\n If you were to view the chimera sequence table you would see the amount of reads, in each sample, that are attributed to each unique sequence, where the columns are the full ~250 bp sequences. The sequences in each of these columns are the \u201cASVs\u201d, which we will soon assign taxonomy to by aligning them to a reference database.<\/p>\n\n\n\n #Now this part is important. We want to know exactly how many reads, in each sample, were dropped from each step we have gone through, which we can track with the following commands\u2026<\/p>\n\n\n\nlibrary(BiocManager)\nlibrary(dplyr)\nlibrary(scales)\nlibrary(grid)\nlibrary(reshape2)\nlibrary(dada2)\nlibrary(tidyr)\nlibrary(stringr)\nlibrary(phyloseq); packageVersion(\"phyloseq\")\nlibrary(Biostrings); packageVersion(\"Biostrings\")\nlibrary(ggplot2); packageVersion(\"ggplot2\")\nlibrary(vegan)\nlibrary(lessR)\nlibrary(ape)\nlibrary(reshape2)\nlibrary(microbiome)<\/code><\/pre>\n\n\n\npath_to_files <- \"\/Users\/kellyshannon\/BIOINFORMATICS\/MiSeq_Beavers\/paired_reads\"\nlist.files(path_to_files)<\/code><\/pre>\n\n\n\nforward_seqs <- sort(list.files(path_to_files, pattern = \"_F.fastq\", full.names = TRUE))\nreverse_seqs <- sort(list.files(path_to_files, pattern = \"_R.fastq\", full.names = TRUE))<\/code><\/pre>\n\n\n\nsample_names <- sapply(strsplit(basename(forward_seqs), \"_\"), `[`, 1)\nsample_names<\/code><\/pre>\n\n\n\n
#We are now going to want to inspect the quality of our forward and reverse reads. Remember that this will be displaying the Phred Score on the y-axis (quality of the identity of nucleobases: https:\/\/en.wikipedia.org\/wiki\/Phred_quality_score<\/a>) and the length of the reads on the x-axis. Now, we can use DADA2 for this, but I think multiqc does a much better job of displaying read quality, and gives lots of additional info. Assuming your reads have been demultiplexed by the sequencing center (they will be generally from the OSU CGRB), they also may have come with their own FastQC reports, which can be put in directly to multiqc to generate an overall report. If not, you can download FastQC and MultiQC individually and use them on your sequences (links below). Esseentially, FastQC will generate a similar report to what we\u2019ll generate with DADA2, below, but then MultiQC can be used to compile all of your FastQC output into one beautiful file that can be visualized on your browser. It will also show you whether or not there are any adaptors in your reads, which will need to get trimmed away.<\/p>\n\n\n\nplotQualityProfile(forward_seqs)\nplotQualityProfile(reverse_seqs)<\/code><\/pre>\n\n\n\nfilt_forward_seqs <- file.path(path_to_files, \"filtered\", paste0(sample_names, \"_F_filt.fastq.gz\"))\nfilt_reverse_seqs <- file.path(path_to_files, \"filtered\", paste0(sample_names, \"_R_filt.fastq.gz\"))\nnames(filt_forward_seqs) <- sample_names\nnames(filt_reverse_seqs) <- sample_names<\/code><\/pre>\n\n\n\nfiltered_reads <- filterAndTrim(forward_seqs, filt_forward_seqs, reverse_seqs, filt_reverse_seqs, \ntruncLen = c(250,140),\nmaxN = 0, #This is getting rid of any ambiguous base calls. \nmaxEE = c(2,2), #Setting the amount of expected errors in reads, increasing this value lets more reads through. \ntruncQ = 2, #Truncates reads at the first sign of a quality score equal to the set value. \nrm.phix = TRUE, #Removes any phiX phage genomes, which are often used by sequencing facilities as controls. \ncompress = TRUE, multithread = TRUE)<\/code><\/pre>\n\n\n\nderep_forward <- derepFastq(filt_forward_seqs, verbose = TRUE)\nderep_reverse <- derepFastq(filt_reverse_seqs, verbose = TRUE)<\/code><\/pre>\n\n\n\nerrors_F <- learnErrors(filt_forward_seqs, multithread = TRUE)\nerrors_R <- learnErrors(filt_reverse_seqs, multithread = TRUE)<\/code><\/pre>\n\n\n\nsetDadaOpt(DETECT_SINGLETONS=TRUE, OMEGA_A=1e-30)\n\ndada_forward_seqs <- dada(derep_forward, err=errors_F, multithread=TRUE, pool = \"pseudo\")\ndada_reverse_seqs <- dada(derep_reverse, err=errors_R, multithread=TRUE, pool = \"pseudo\")<\/code><\/pre>\n\n\n\ncontigs <- mergePairs(dada_forward_seqs, derep_forward, dada_reverse_seqs, derep_reverse, verbose=TRUE)\nseqtable <- makeSequenceTable(contigs)<\/code><\/pre>\n\n\n\nseqtable_nochim <- removeBimeraDenovo(seqtable, method = \"consensus\", multithread=TRUE, verbose=TRUE)\ndim(seqtable_nochim)\ntable(nchar(getSequences(seqtable)))<\/code><\/pre>\n\n\n\nseqtable.2 <- seqtable[,nchar(colnames(seqtable)) %in% 250:257]\nseqtable <- seqtable.2\ntable(nchar(getSequences(seqtable)))\ndim(seqtable_nochim)<\/code><\/pre>\n\n\n\ngetN <- function(x) sum(getUniques(x))\ntrack_reads <- cbind(filtered_reads, sapply(dada_forward_seqs, getN), sapply(dada_reverse_seqs, getN), sapply(contigs, getN), rowSums(seqtable_nochim))\ncolnames(track_reads) <- c(\"input\", \"filtered\", \"denoisedF\", \"denoisedR\", \"merged\", \"nonchim\")\nrownames(track_reads) <- sample_names<\/code><\/pre>\n\n\n\n