dada2 Single-End Protocol - Boston_run1

Researchers/collaborators:

• Gabriella Edfeldt gabriella.edfeldt@ki.se
  
• Jiawu Xu jxu25@mgjh.harvard.edu
  
• Douglas S. Kwon dkwon@mgh.harvard.edu

References:

• Callahan, B. J., McMurdie, P. J., Rosen, M. J., Han, A. W., Johnson, A. J. A., & Holmes, S. P. (2016). DADA2: High-resolution sample inference from Illumina amplicon data. Nature Methods, 13(7), 581–3.
  • https://doi.org/10.1038/nmeth.3869
    
• McMurdie, P. J., & Holmes, S. (2013). phyloseq: an R package for reproducible interactive analysis and graphics of microbiome census data. PloS One, 8(4), e61217.
  • https://doi.org/10.1371/journal.pone.0061217
    
• DADA2 tutorial: http://benjjneb.github.io/dada2/tutorial.html
  
• DADA2 big data tutorial: http://benjjneb.github.io/dada2/bigdata.html
  
• Phyloseq: https://joey711.github.io/phyloseq/index.html
  
• Qiime 1 Mapping files: http://qiime.org/documentation/file_formats.html#mapping-file-overview

Expected output

1. A knitted R markdown document (default format is HTML, but can also select options including PDF or MSWord) that contains all comments, code, and figures produced. Key highlights include:
   • The document returns the names of any samples that failed to pass filter and were eliminated from the final output tables and phyloseq output.
     
   • The document calculates and returns cumulative processing time at various steps along the analysis.
     
2. Filtered reads with quality information for each sample stored in newly created subdirectory /Data/filtered.
   • If sequence files were originally named samplename.fastq, the corresponding filtered files will be named samplename.filt.fastq
     
3. A newly created subdirectory named /Output containing:
   • Plot of aggregated read quality scores from the {plot-unfiltered-run-quality} chunk
     
   • Plot(s) of read quality scores for individual sample(s) (if code for these plots is included in the {plot-unfiltered-run-quality} chunk)
     
   • Plot of error rates from the {learn-errors} chunk
     
   • Table of read counts and percentages retained after each sequence analysis step for each sample, saved in a .txt file in TSV format (can be opened in Excel)
     
   • Plot of read percentages retained after each sequence analysis step for each sample
     
4. A newly created subdirectory named /ps_objects containing a different phyloseq object for each taxonomic database used. Each phyloseq object is created in the {create-phyloseq} chunk and contains:
   • An otu_table (ASV table) based on amplicon sequence variants (ASVs) after dada sequence inference and chimera removal
     
   • A tax_table based on taxonomic assignments performed in the {assign-taxonomy} chunk
   • sample_data derived from the mapping file.
     
   • Phylogenetic trees are not constructed because depending on the type of sample sequenced, they may not be relevant (e.g. in vitro synthetic mixtures of laboratory bacterial strains for QC)
     
5. A newly created subdirectory named /RDS containing:
   • An RDS file containing the sequence table after chimera removal.
     
   • An RData workspaced imaging saving the final workspace for the analysis.

knitr::opts_knit$set(root.dir = getwd())

#Define variable with start time of running script.
start_time <- proc.time()

# Load libraries
library(ShortRead)
packageVersion("ShortRead")
## [1] '1.44.0'
library(dada2)
packageVersion("dada2")
## [1] '1.14.0'
library("phangorn")
packageVersion("phangorn")
## [1] '2.7.0'
library("phyloseq")
packageVersion("phyloseq")
## [1] '1.30.0'
library("ips")
packageVersion("ips")
## [1] '0.0.11'
library("tidyverse")
packageVersion("tidyverse")
## [1] '1.2.1'

# Set random seed for reproducibility purposes
set.seed(100)

if(params$run=="Boston_run1"){
  run.name <- "Boston_run1"
  run.date <-  "2019_02_25_v1" 
  mapping_file <- "Mapping_info_Boston_run1.csv"
  # Trim read lenght left and right
  FwdTrimLeft <- 10   
  FwdTrimRight <- 225  
}

if(params$run=="Boston_run2"){
  run.name <- "Boston_run2"
  run.date <-  "2019_05_25" 
  mapping_file <- "Mapping_info_Boston_run2.csv"
  # Trim read lenght left and right
  FwdTrimLeft <- 0   
  FwdTrimRight <- 252  
}

# Sequencing run date in format YYYY_MM_DD
print("Sequencing run date and version:")
run.date

## [1] "Sequencing run date and version:"
## [1] "2019_02_25_v1"

# if(params$run=="Boston_run1"){
#   url <- c(
#     "ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR919/001/SRR9198521/SRR9198521.fastq.gz",
#     "ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR919/007/SRR9198517/SRR9198517.fastq.gz",
#     "ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR919/008/SRR9198518/SRR9198518.fastq.gz"
#     ) }
# dir <- "../data/Fastq"
# dir.create(dir)
# purrr::walk(url, ~download.file(.x, file.path("../data/Fastq/", basename(.x)), method="auto"))
# system("gunzip ../data/Fastq/*.gz")
# 

if(params$run == "Boston_run1"){dir  <- "/Users/vilkal/Raw_data/Microbiome/Run1_Boston_CVLv3_March2019"}
if(params$run == "Boston_run2"){dir <- "/Users/vilkal/Raw_data/Microbiome/Run2_Boston_Tissuev3_June2019"}

files <- list.files(path = dir, pattern =".fastq$", full.names = T, recursive = T)

#Define and print path to working directory
pathwd <- getwd()
print("Working directory:")
pathwd

# Define path to subdirectory containing demultiplexed forward-read fastq files
pathF <- file.path("../data/Fastq")
print("Directory containing raw demultiplexed .fastq sequence files:")
pathF

#Define path to directory containing 16S databases
print("Path to 16S databases:")
pathdb <- file.path("../Resources/RDP_database")
pathdb

#Define path to directory containing mapping file
print("Path to mapping file:")
pathMap <- file.path("../data/Mapping")
pathMap

#Create and define path to subdirectory in which filtered files will be stored
filtpathF <- file.path("../data/Fastq/filtered") 
print("Directory that will contain filtered .fastq files:")
filtpathF
dir.create(filtpathF)

#Create and define path to subdirectory in which to store output figures and tables
print("Directory in which to store output figures and tables:")
pathOut <- file.path("../results/dada2_output")
pathOut
dir.create(pathOut)

#Create and define path to subdirectory in which RDS files will be stored
print("Path to directory for saving RDS files:")
pathF.RDS <- file.path(pathOut, "RDS")
pathF.RDS
dir.create(pathF.RDS)

#Create and define path to subdirectory in which RDS files containing phyloseq objects will be stored
print("Path to sub-directory in which to save phyloseq objects:")
pathps0 <- file.path(pathOut, "ps_objects")
pathps0
dir.create(pathps0)

## [1] "Working directory:"
## [1] "/Users/vilkal/work/Brolidens_work/Projects/broliden_5325/reports"
## [1] "Directory containing raw demultiplexed .fastq sequence files:"
## [1] "../data/Fastq"
## [1] "Path to 16S databases:"
## [1] "../Resources/RDP_database"
## [1] "Path to mapping file:"
## [1] "../data/Mapping"
## [1] "Directory that will contain filtered .fastq files:"
## [1] "../data/Fastq/filtered"
## [1] "Directory in which to store output figures and tables:"
## [1] "../results/dada2_output"
## [1] "Path to directory for saving RDS files:"
## [1] "../results/dada2_output/RDS"
## [1] "Path to sub-directory in which to save phyloseq objects:"
## [1] "../results/dada2_output/ps_objects"

# File parsing
#Create vector of fastq files in pathF.directory
#If analyzing paired end reads, the search pattern below would need to be modified (see DADA2 tutorial)

print("Sequence files to be analyzed:")
fnFs <- sort(files)
head(fnFs)

#Extract sample names, assuming filenames have format: SAMPLENAME.fastq or SAMPLENAME.fastq.gz
sample.names <- sapply(strsplit(basename(fnFs), ".fast"), `[`, 1)
head(sample.names)

## [1] "Sequence files to be analyzed:"
## [1] "/Users/vilkal/Raw_data/Microbiome/Run1_Boston_CVLv3_March2019/Luminal_fastq_files/1872.022.rcbc214.2019.02.25.fastq"
## [2] "/Users/vilkal/Raw_data/Microbiome/Run1_Boston_CVLv3_March2019/Luminal_fastq_files/1907.032.rcbc226.2019.02.25.fastq"
## [3] "/Users/vilkal/Raw_data/Microbiome/Run1_Boston_CVLv3_March2019/Luminal_fastq_files/1954.043.rcbc238.2019.02.25.fastq"
## [4] "/Users/vilkal/Raw_data/Microbiome/Run1_Boston_CVLv3_March2019/Luminal_fastq_files/1996.054.rcbc250.2019.02.25.fastq"
## [5] "/Users/vilkal/Raw_data/Microbiome/Run1_Boston_CVLv3_March2019/Luminal_fastq_files/2051.064.rcbc262.2019.02.25.fastq"
## [6] "/Users/vilkal/Raw_data/Microbiome/Run1_Boston_CVLv3_March2019/Luminal_fastq_files/2140.110.rcbc320.2019.02.25.fastq"
## [1] "1872.022.rcbc214.2019.02.25" "1907.032.rcbc226.2019.02.25"
## [3] "1954.043.rcbc238.2019.02.25" "1996.054.rcbc250.2019.02.25"
## [5] "2051.064.rcbc262.2019.02.25" "2140.110.rcbc320.2019.02.25"

#Plot forward read quality aggregated across all forward read files
p.qual.f.1.agg <- plotQualityProfile(fnFs, aggregate = TRUE) + 
  ggtitle("Fwd read quality aggregated across samples")
p.qual.f.1.agg

#Save figure to file
ggsave(filename = file.path(pathOut, paste0("Read_quality_aggregate_fwd_", run.date, ".pdf")), 
       plot = p.qual.f.1.agg, device = "pdf", width = 8, height = 6, units = "in") 

#Calculate cumulative processing time
print("Cumulative processing time (seconds):")
proc.time() - start_time

## [1] "Cumulative processing time (seconds):"
##    user  system elapsed 
## 168.697   9.304 179.063

Define positions at which to trim filtered data

## ## ## ## ## ## ## ## ## ##
## USER INPUT REQUIRED HERE ##
## ## ## ## ## ## ## ## ## ##

# Define base positions at which to trim the filtered forward read sequences
# Decisions about where to trim should be based on the quality profile above.
# See DADA2 tutorial for additional information.
# The two variable defined here will then be supplied to the filterAndTrim function
  # FwdTrimLeft gives the posion from the left margin at which to trim
  # FwdTrimRight gives the position from the right margin at which to trim
# For example if FwdTrimLeft is 10 and FwdTrimRight is 240, sequence lengths will be 230 bp.

print(paste0("Trimming forward reads on left at base position ", FwdTrimLeft, 
       " and on right at base position ", FwdTrimRight))

## [1] "Trimming forward reads on left at base position 10 and on right at base position 225"

Quality filter and trim

# Create vector of modified names of fastq files after filtering in
# which the modified file name changes from "samplename.fastq" 
# to "samplename.filt.fastq"

filtFs <- file.path(filtpathF, str_replace(basename(fnFs), "(.fastq)", ".filt\\1"))
head(filtFs)

# Filter and trim sequences listed in fnFs and store them in filtFs, 
# while storing summary in dataframe out

out <- filterAndTrim(fwd = fnFs, filt = filtFs, rev = NULL, filt.rev = NULL,
              trimLeft = FwdTrimLeft, truncLen = c(FwdTrimRight), maxEE = 2, truncQ = 11, maxN = 0, 
              rm.phix = TRUE, compress = TRUE, verbose = TRUE, multithread = TRUE)
head(out)
#truncate file extension (".fastq" or ".fastq.gz") from rownames
rownames(out) <- rownames(out) %>% str_replace(".fast.+", "")
head(out)

## [1] "../data/Fastq/filtered/1872.022.rcbc214.2019.02.25.filt.fastq"
## [2] "../data/Fastq/filtered/1907.032.rcbc226.2019.02.25.filt.fastq"
## [3] "../data/Fastq/filtered/1954.043.rcbc238.2019.02.25.filt.fastq"
## [4] "../data/Fastq/filtered/1996.054.rcbc250.2019.02.25.filt.fastq"
## [5] "../data/Fastq/filtered/2051.064.rcbc262.2019.02.25.filt.fastq"
## [6] "../data/Fastq/filtered/2140.110.rcbc320.2019.02.25.filt.fastq"
##                                   reads.in reads.out
## 1872.022.rcbc214.2019.02.25.fastq   110022     99554
## 1907.032.rcbc226.2019.02.25.fastq    67306     61170
## 1954.043.rcbc238.2019.02.25.fastq   114558    100639
## 1996.054.rcbc250.2019.02.25.fastq    73017     68506
## 2051.064.rcbc262.2019.02.25.fastq   157822    136551
## 2140.110.rcbc320.2019.02.25.fastq    32169     29611
##                             reads.in reads.out
## 1872.022.rcbc214.2019.02.25   110022     99554
## 1907.032.rcbc226.2019.02.25    67306     61170
## 1954.043.rcbc238.2019.02.25   114558    100639
## 1996.054.rcbc250.2019.02.25    73017     68506
## 2051.064.rcbc262.2019.02.25   157822    136551
## 2140.110.rcbc320.2019.02.25    32169     29611

# Calculate cumulative processing time
print("Cumulative processing time (seconds):")
proc.time() - start_time

## [1] "Cumulative processing time (seconds):"
##    user  system elapsed 
## 912.166  43.663 268.712

Sequence inference

## Set paramaters
# Create character of all files in directory filtpathF.1 (i.e. all filtered files)
# with names including string "fastq" (i.e. all *.fastq and *.fastq.gz files)
filtFs.pass <- list.files(filtpathF, pattern="fastq", full.names = TRUE)
head(filtFs.pass)

# Create character vector of sample names from the file names in filtFs.1
# This code assumes all filename = "samplename.filt.fastq"" or "samplename.filt.fastq.gz""
sample.names.filt.pass <- sapply(strsplit(basename(filtFs.pass), ".filt.fast"), `[`, 1) 
head(sample.names.filt.pass)

# Name the elements in filtFs according to the elements in sample.names.filt
names(filtFs.pass) <- sample.names.filt.pass

# Samples for which no sequences passed the filtering step
print("Samples for which no sequences passed filter criteria (i.e. samples excluded from filtered dataset):")
sample.names.filt.fail <- setdiff(sample.names, sample.names.filt.pass)
sample.names.filt.fail

## [1] "../data/Fastq/filtered/101u6866T.165.rcbc179.2019.06.10.filt.fastq"
## [2] "../data/Fastq/filtered/105u7143T.166.rcbc183.2019.06.10.filt.fastq"
## [3] "../data/Fastq/filtered/107u7152T.167.rcbc185.2019.06.10.filt.fastq"
## [4] "../data/Fastq/filtered/108u5343T.168.rcbc186.2019.06.10.filt.fastq"
## [5] "../data/Fastq/filtered/10u8310T.009.rcbc009.2019.06.10.filt.fastq" 
## [6] "../data/Fastq/filtered/111u7150T.169.rcbc189.2019.06.10.filt.fastq"
## [1] "101u6866T.165.rcbc179.2019.06.10" "105u7143T.166.rcbc183.2019.06.10"
## [3] "107u7152T.167.rcbc185.2019.06.10" "108u5343T.168.rcbc186.2019.06.10"
## [5] "10u8310T.009.rcbc009.2019.06.10"  "111u7150T.169.rcbc189.2019.06.10"
## [1] "Samples for which no sequences passed filter criteria (i.e. samples excluded from filtered dataset):"
## character(0)

# Set random seed for reproducibility purposes
set.seed(100)

# Learn forward read error rates
errF <- learnErrors(filtFs.pass, nread=1e6, multithread=TRUE, randomize = "TRUE")

## Warning in learnErrors(filtFs.pass, nread = 1e+06, multithread = TRUE, randomize
## = "TRUE"): The nreads parameter is DEPRECATED. Please update your code with the
## nbases parameter.

# Plot convergence of error rate computation
plot(dada2:::checkConvergence(errF), type = "o", col = "firebrick3", main = "Convergence")

# Plot calculated errors
p.err.F <- plotErrors(errF, nominalQ = TRUE)
p.err.F

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Transformation introduced infinite values in continuous y-axis

# Save figure to file
ggsave(filename = file.path(pathOut, paste0("Error_rates_fwd_", run.date, ".pdf")),
       plot = p.err.F, device = "pdf", width = 8, height = 6, units = "in")

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Transformation introduced infinite values in continuous y-axis

## 246662803 total bases in 1131020 reads from 24 samples will be used for learning the error rates.

#Calculate cumulative processing time
print("Cumulative processing time (seconds):")
proc.time() - start_time

## [1] "Cumulative processing time (seconds):"
##     user   system  elapsed 
## 1677.122   56.091  416.242

# Dereplicate and apply error rate to resolve sequence variants

## ## ## ## ## ## ## ## ## ## ## ##
## This chunk gives two different options for analyzing the data, depending on the size of the dataset
## and availability of memory for analysis. Use only Option 1 or Option 2 for analysis, not both.
## Comment (or uncomment) the relevant option by highlighting and hitting <command> + <shift> + C
## ## ## ## ## ## ## ## ## ## ## ##

## ## ## ## ## ## ## ## ## ## ## ## ## ## ##
## OPTION 1: **************************** ##
## ## ## ## ## ## ## ## ## ## ## ## ## ## ##

### Code for smaller datasets (e.g. most MiSeq) for which memory is adequate for parallelization ###
### This code analyzes the samples in parallel, which is more efficient but more memory-intensive ###
### Derived from https://benjjneb.github.io/dada2/tutorial.html ###

# Dereplicate the filtered forward sequences by sample
derepFs <- derepFastq(filtFs.pass, verbose = FALSE)

# Set random seed for reproducibility purposes
set.seed(100)

# Run dada on the dereplicated forward sequences to resolve sequence variants
# based on the forward read error model calculated in errF
dadaFs <- dada(derepFs, err = errF, multithread = TRUE)


# ## ## ## ## ## ## ## ## ## ## ## ## ## ## ##
# ## OPTION 2: *************************** ##
# ## ## ## ## ## ## ## ## ## ## ## ## ## ## ##

# ### Code for very large (e.g. HiSeq) datasets if parallelization would exceed available memory ###
# ### This code analyzes the samples in series, which is less efficient but less memory-intensive ###
# ### Derived from https://benjjneb.github.io/dada2/bigdata.html ###
# 
# #Create lists derepFs and dadaFs of same length and element names as sample.names.filt.pass
# derepFs <- vector("list", length(sample.names.filt.pass))
# names(derepFs) <- sample.names.filt.pass
# dadaFs <- vector("list", length(sample.names.filt.pass))
# names(dadaFs) <- sample.names.filt.pass
# 
# #For each sequence file in filtFs.1, dereplicate the sequences from the file
# #then run dada to resolve sequence variants based on the error model errF.1
# #store results in the corresponding named element of singles.1
# for(sam in sample.names.filt.pass) {
#   cat("Processing:", sam, "\n")
#   derepFs[[sam]] <- derepFastq(filtFs.pass[[sam]])
#   # Set random seed for reproducibility purposes
#   set.seed(100)
#   dadaFs[[sam]] <- dada(derepFs[[sam]], err=errF, multithread=TRUE)
# }
# rm(sam)

## Sample 1 - 3830 reads in 636 unique sequences.
## Sample 2 - 12506 reads in 1141 unique sequences.
## Sample 3 - 3266 reads in 494 unique sequences.
## Sample 4 - 2542 reads in 423 unique sequences.
## Sample 5 - 3337 reads in 430 unique sequences.
## Sample 6 - 5724 reads in 966 unique sequences.
## Sample 7 - 6767 reads in 819 unique sequences.
## Sample 8 - 22683 reads in 1061 unique sequences.
## Sample 9 - 2347 reads in 394 unique sequences.
## Sample 10 - 36268 reads in 3824 unique sequences.
## Sample 11 - 8755 reads in 1058 unique sequences.
## Sample 12 - 5053 reads in 536 unique sequences.
## Sample 13 - 1781 reads in 279 unique sequences.
## Sample 14 - 11438 reads in 761 unique sequences.
## Sample 15 - 99554 reads in 14313 unique sequences.
## Sample 16 - 84896 reads in 14497 unique sequences.
## Sample 17 - 425 reads in 92 unique sequences.
## Sample 18 - 61170 reads in 5667 unique sequences.
## Sample 19 - 100639 reads in 13896 unique sequences.
## Sample 20 - 68506 reads in 7558 unique sequences.
## Sample 21 - 5397 reads in 798 unique sequences.
## Sample 22 - 9035 reads in 779 unique sequences.
## Sample 23 - 136551 reads in 22538 unique sequences.
## Sample 24 - 2369 reads in 357 unique sequences.
## Sample 25 - 29611 reads in 3793 unique sequences.
## Sample 26 - 54119 reads in 6780 unique sequences.
## Sample 27 - 83040 reads in 9898 unique sequences.
## Sample 28 - 54658 reads in 9373 unique sequences.
## Sample 29 - 1982 reads in 404 unique sequences.
## Sample 30 - 14914 reads in 2086 unique sequences.
## Sample 31 - 1171 reads in 258 unique sequences.
## Sample 32 - 151505 reads in 22507 unique sequences.
## Sample 33 - 19758 reads in 3206 unique sequences.
## Sample 34 - 125170 reads in 20343 unique sequences.
## Sample 35 - 61002 reads in 12251 unique sequences.
## Sample 36 - 10503 reads in 772 unique sequences.
## Sample 37 - 16132 reads in 1372 unique sequences.
## Sample 38 - 52807 reads in 7331 unique sequences.
## Sample 39 - 72200 reads in 6623 unique sequences.
## Sample 40 - 53252 reads in 5892 unique sequences.
## Sample 41 - 6193 reads in 790 unique sequences.
## Sample 42 - 82649 reads in 12281 unique sequences.
## Sample 43 - 49580 reads in 5572 unique sequences.
## Sample 44 - 123635 reads in 22344 unique sequences.
## Sample 45 - 4898 reads in 544 unique sequences.
## Sample 46 - 46289 reads in 5440 unique sequences.
## Sample 47 - 28738 reads in 3620 unique sequences.
## Sample 48 - 66583 reads in 9588 unique sequences.
## Sample 49 - 68332 reads in 8997 unique sequences.
## Sample 50 - 102738 reads in 17395 unique sequences.
## Sample 51 - 26274 reads in 2806 unique sequences.
## Sample 52 - 204426 reads in 22584 unique sequences.
## Sample 53 - 58308 reads in 5431 unique sequences.
## Sample 54 - 3415 reads in 371 unique sequences.
## Sample 55 - 542 reads in 153 unique sequences.
## Sample 56 - 46684 reads in 6797 unique sequences.
## Sample 57 - 99987 reads in 13223 unique sequences.
## Sample 58 - 122217 reads in 20463 unique sequences.
## Sample 59 - 3234 reads in 380 unique sequences.
## Sample 60 - 77751 reads in 10614 unique sequences.
## Sample 61 - 19296 reads in 1230 unique sequences.
## Sample 62 - 66193 reads in 5949 unique sequences.
## Sample 63 - 5801 reads in 908 unique sequences.
## Sample 64 - 61538 reads in 6897 unique sequences.
## Sample 65 - 59419 reads in 5353 unique sequences.
## Sample 66 - 98522 reads in 12412 unique sequences.
## Sample 67 - 8386 reads in 1039 unique sequences.
## Sample 68 - 138784 reads in 25595 unique sequences.
## Sample 69 - 14956 reads in 1594 unique sequences.
## Sample 70 - 117485 reads in 19865 unique sequences.
## Sample 71 - 83965 reads in 8668 unique sequences.
## Sample 72 - 90306 reads in 6834 unique sequences.
## Sample 73 - 148625 reads in 21917 unique sequences.
## Sample 74 - 53315 reads in 5503 unique sequences.
## Sample 75 - 2638 reads in 373 unique sequences.
## Sample 76 - 182447 reads in 31774 unique sequences.
## Sample 77 - 60642 reads in 9587 unique sequences.
## Sample 78 - 54522 reads in 9191 unique sequences.
## Sample 79 - 54476 reads in 10831 unique sequences.
## Sample 80 - 4082 reads in 575 unique sequences.
## Sample 81 - 23956 reads in 2264 unique sequences.
## Sample 82 - 88668 reads in 13177 unique sequences.
## Sample 83 - 5259 reads in 867 unique sequences.
## Sample 84 - 23170 reads in 3003 unique sequences.
## Sample 85 - 55105 reads in 7408 unique sequences.
## Sample 86 - 13682 reads in 1271 unique sequences.
## Sample 87 - 11007 reads in 1757 unique sequences.
## Sample 88 - 7858 reads in 1006 unique sequences.
## Sample 89 - 42323 reads in 3821 unique sequences.
## Sample 90 - 41971 reads in 4779 unique sequences.
## Sample 91 - 129718 reads in 18484 unique sequences.
## Sample 92 - 154074 reads in 24765 unique sequences.
## Sample 93 - 53347 reads in 4218 unique sequences.
## Sample 94 - 1997 reads in 302 unique sequences.
## Sample 95 - 24712 reads in 2281 unique sequences.
## Sample 96 - 4628 reads in 580 unique sequences.
## Sample 97 - 3009 reads in 536 unique sequences.
## Sample 98 - 3072 reads in 544 unique sequences.
## Sample 99 - 4351 reads in 667 unique sequences.
## Sample 100 - 131987 reads in 20662 unique sequences.
## Sample 101 - 6628 reads in 1051 unique sequences.
## Sample 102 - 62350 reads in 11409 unique sequences.
## Sample 103 - 6949 reads in 934 unique sequences.
## Sample 104 - 119225 reads in 17305 unique sequences.
## Sample 105 - 48390 reads in 7883 unique sequences.
## Sample 106 - 2867 reads in 446 unique sequences.
## Sample 107 - 88113 reads in 8553 unique sequences.
## Sample 108 - 107725 reads in 15535 unique sequences.
## Sample 109 - 2022 reads in 370 unique sequences.
## Sample 110 - 130090 reads in 24860 unique sequences.
## Sample 111 - 113626 reads in 17016 unique sequences.
## Sample 112 - 47818 reads in 5341 unique sequences.
## Sample 113 - 91165 reads in 11363 unique sequences.
## Sample 114 - 102866 reads in 17551 unique sequences.
## Sample 115 - 53124 reads in 6810 unique sequences.
## Sample 116 - 6953 reads in 941 unique sequences.
## Sample 117 - 50081 reads in 3513 unique sequences.
## Sample 118 - 98203 reads in 17491 unique sequences.
## Sample 119 - 93145 reads in 11591 unique sequences.
## Sample 120 - 51403 reads in 10940 unique sequences.
## Sample 121 - 63737 reads in 5410 unique sequences.
## Sample 122 - 101586 reads in 22825 unique sequences.
## Sample 123 - 16904 reads in 2901 unique sequences.
## Sample 124 - 48191 reads in 4721 unique sequences.
## Sample 125 - 50825 reads in 5492 unique sequences.
## Sample 126 - 117705 reads in 16179 unique sequences.
## Sample 127 - 128179 reads in 21246 unique sequences.
## Sample 128 - 55943 reads in 3874 unique sequences.
## Sample 129 - 116646 reads in 18596 unique sequences.
## Sample 130 - 124978 reads in 21312 unique sequences.
## Sample 131 - 57007 reads in 5466 unique sequences.
## Sample 132 - 3774 reads in 584 unique sequences.
## Sample 133 - 52266 reads in 5671 unique sequences.
## Sample 134 - 468 reads in 104 unique sequences.
## Sample 135 - 32131 reads in 3863 unique sequences.
## Sample 136 - 56310 reads in 5237 unique sequences.
## Sample 137 - 875 reads in 186 unique sequences.
## Sample 138 - 108038 reads in 14585 unique sequences.
## Sample 139 - 82116 reads in 9330 unique sequences.
## Sample 140 - 12229 reads in 2108 unique sequences.
## Sample 141 - 12660 reads in 1227 unique sequences.
## Sample 142 - 6006 reads in 831 unique sequences.
## Sample 143 - 116130 reads in 18143 unique sequences.
## Sample 144 - 95444 reads in 11048 unique sequences.
## Sample 145 - 91957 reads in 11335 unique sequences.
## Sample 146 - 58194 reads in 8083 unique sequences.
## Sample 147 - 53233 reads in 4666 unique sequences.
## Sample 148 - 136166 reads in 23682 unique sequences.
## Sample 149 - 49833 reads in 3755 unique sequences.
## Sample 150 - 54888 reads in 8538 unique sequences.
## Sample 151 - 48866 reads in 4366 unique sequences.
## Sample 152 - 47998 reads in 3142 unique sequences.
## Sample 153 - 2317 reads in 217 unique sequences.
## Sample 154 - 127416 reads in 20523 unique sequences.
## Sample 155 - 159029 reads in 25178 unique sequences.
## Sample 156 - 3541 reads in 507 unique sequences.
## Sample 157 - 67434 reads in 6014 unique sequences.
## Sample 158 - 3079 reads in 496 unique sequences.
## Sample 159 - 55174 reads in 6339 unique sequences.
## Sample 160 - 63004 reads in 5951 unique sequences.
## Sample 161 - 1036 reads in 196 unique sequences.
## Sample 162 - 54539 reads in 5967 unique sequences.
## Sample 163 - 4846 reads in 643 unique sequences.
## Sample 164 - 7917 reads in 1017 unique sequences.
## Sample 165 - 79080 reads in 7725 unique sequences.
## Sample 166 - 30622 reads in 4218 unique sequences.
## Sample 167 - 11546 reads in 1622 unique sequences.
## Sample 168 - 57507 reads in 6934 unique sequences.
## Sample 169 - 69661 reads in 7285 unique sequences.
## Sample 170 - 10678 reads in 1328 unique sequences.
## Sample 171 - 2293 reads in 365 unique sequences.
## Sample 172 - 17662 reads in 2874 unique sequences.
## Sample 173 - 41767 reads in 3900 unique sequences.
## Sample 174 - 892 reads in 181 unique sequences.
## Sample 175 - 74259 reads in 7057 unique sequences.
## Sample 176 - 4501 reads in 638 unique sequences.
## Sample 177 - 43251 reads in 3353 unique sequences.
## Sample 178 - 134314 reads in 17056 unique sequences.
## Sample 179 - 3641 reads in 619 unique sequences.
## Sample 180 - 72443 reads in 8150 unique sequences.
## Sample 181 - 10813 reads in 1562 unique sequences.
## Sample 182 - 53295 reads in 5636 unique sequences.
## Sample 183 - 107956 reads in 13427 unique sequences.
## Sample 184 - 4160 reads in 679 unique sequences.
## Sample 185 - 12714 reads in 2106 unique sequences.
## Sample 186 - 8240 reads in 1216 unique sequences.
## Sample 187 - 7165 reads in 943 unique sequences.
## Sample 188 - 6769 reads in 1002 unique sequences.
## Sample 189 - 15919 reads in 2463 unique sequences.
## Sample 190 - 19986 reads in 2910 unique sequences.
## Sample 191 - 14453 reads in 2180 unique sequences.
## Sample 192 - 9983 reads in 1363 unique sequences.
## Sample 193 - 15411 reads in 2173 unique sequences.
## Sample 194 - 9294 reads in 1355 unique sequences.
## Sample 195 - 7557 reads in 1335 unique sequences.
## Sample 196 - 13601 reads in 1866 unique sequences.
## Sample 197 - 13172 reads in 1581 unique sequences.
## Sample 198 - 10927 reads in 1748 unique sequences.
## Sample 199 - 822 reads in 176 unique sequences.
## Sample 200 - 9734 reads in 1807 unique sequences.
## Sample 201 - 20985 reads in 1159 unique sequences.

# Calculate cumulative processing time
print("Cumulative processing time (seconds):")

## [1] "Cumulative processing time (seconds):"

proc.time() - start_time

##     user   system  elapsed 
## 2620.622   69.918  667.978

Construct sequence table and remove chimeras

# Make sequence table and review dimensions
seqtab.prechimeraremoval <- makeSequenceTable(dadaFs)
print("Sequence table dimensions (samples, resolved sequences) before chimera removal:")
dim(seqtab.prechimeraremoval)

# Inspect distribution of sequence lengths
print("Distribution of sequence lengths before chimera removal:")
table(nchar(getSequences(seqtab.prechimeraremoval)))

# Remove chimeras
seqtab.nochim <- removeBimeraDenovo(seqtab.prechimeraremoval,
                                    method = "consensus",
                                    multithread = TRUE,
                                    verbose = TRUE)

## Identified 6287 bimeras out of 8386 input sequences.

print("Sequence table dimensions (samples, resolved sequences) after chimera removal:")
dim(seqtab.nochim)

# Inspect distribution of sequence lengths
print("Distribution of sequence lengths after chimera removal:")
table(nchar(getSequences(seqtab.nochim)))
print("Proportion of reads eliminated by chimera removal:")
sum(seqtab.nochim) / sum(seqtab.prechimeraremoval)

## [1] "Sequence table dimensions (samples, resolved sequences) before chimera removal:"
## [1]  201 8386
## [1] "Distribution of sequence lengths before chimera removal:"
## 
##  215  252 
## 6510 1876 
## [1] "Sequence table dimensions (samples, resolved sequences) after chimera removal:"
## [1]  201 2099
## [1] "Distribution of sequence lengths after chimera removal:"
## 
##  215  252 
## 1171  928 
## [1] "Proportion of reads eliminated by chimera removal:"
## [1] 0.9265224

# Save seqtab.nochim RDS file
print("Save seqtab.nochim as RDS file named:")
f.seqtab.nochim.RDS <- file.path(pathF.RDS, paste0("MiSeq_", run.date, "_preprocessing_single_nochim.RDS"))
f.seqtab.nochim.RDS
saveRDS(seqtab.nochim, f.seqtab.nochim.RDS)

## [1] "Save seqtab.nochim as RDS file named:"
## [1] "../results/dada2_output/RDS/MiSeq_2019_02_25_v1_preprocessing_single_nochim.RDS"

# Calculate cumulative processing time
print("Cumulative processing time (seconds):")
proc.time() - start_time

## [1] "Cumulative processing time (seconds):"
##     user   system  elapsed 
## 2667.893   70.258  672.825

Track reads through pipeline

The code in the chunk below is substantially altered from the code suggested in the DADA2 tutorial because
the tutorial code can lead to errors if any samples are eliminated at a step (e.g. if no reads pass filter).

# Function to calculate reads for each element in dada output
getN <- function(x) sum(getUniques(x))

# Read count remaining at each step of processing. Would need additional columns if analyzing
# paired-end data.

# Create interim dataframes listing the sequence count after each step of processing per sample
print("Read count remaining at each step of processing:")

## [1] "Read count remaining at each step of processing:"

outdf <- out %>% 
  as.data.frame() %>%
  tibble::rownames_to_column() %>% 
  dplyr::rename(Sample = "rowname", input = "reads.in", filtered = "reads.out")

dadacountdf <- sapply(dadaFs, getN) %>% 
  as.data.frame() %>% 
  tibble::rownames_to_column() %>%
  dplyr::rename(Sample = "rowname", denoisedF = ".")

nonchimdf <- rowSums(seqtab.nochim) %>% 
  as.data.frame() %>% 
  tibble::rownames_to_column() %>%
  dplyr::rename(Sample = "rowname", nonchim = ".")

# Join dataframes by Sample column to track counts for each sample at each step of processing
track <- dplyr::left_join(outdf, dadacountdf, by = c("Sample" = "Sample")) %>%
  dplyr::left_join(nonchimdf, by = c("Sample" = "Sample"))
head(track)

##                        Sample  input filtered denoisedF nonchim
## 1 1872.022.rcbc214.2019.02.25 110022    99554     99301   90227
## 2 1907.032.rcbc226.2019.02.25  67306    61170     61091   60477
## 3 1954.043.rcbc238.2019.02.25 114558   100639    100248   93753
## 4 1996.054.rcbc250.2019.02.25  73017    68506     68342   65224
## 5 2051.064.rcbc262.2019.02.25 157822   136551    135715  117825
## 6 2140.110.rcbc320.2019.02.25  32169    29611     29529   29474

# Percentage of original unfiltered reads remaining at each step of processing
print("Percentage of original unfiltered foward reads remaining at each step of processing:")

## [1] "Percentage of original unfiltered foward reads remaining at each step of processing:"

track.percent <- dplyr::mutate_at(track, vars(-matches("Sample")), funs(. * 100 / track[["input"]]))

## Warning: `funs()` is deprecated as of dplyr 0.8.0.
## Please use a list of either functions or lambdas: 
## 
##   # Simple named list: 
##   list(mean = mean, median = median)
## 
##   # Auto named with `tibble::lst()`: 
##   tibble::lst(mean, median)
## 
##   # Using lambdas
##   list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.

track.percent <- rename_at(track.percent, vars(-matches("Sample")), funs(paste0(., "_percent")))
head(track.percent)

##                        Sample input_percent filtered_percent denoisedF_percent
## 1 1872.022.rcbc214.2019.02.25           100         90.48554          90.25559
## 2 1907.032.rcbc226.2019.02.25           100         90.88343          90.76605
## 3 1954.043.rcbc238.2019.02.25           100         87.84982          87.50851
## 4 1996.054.rcbc250.2019.02.25           100         93.82199          93.59738
## 5 2051.064.rcbc262.2019.02.25           100         86.52216          85.99245
## 6 2140.110.rcbc320.2019.02.25           100         92.04825          91.79334
##   nonchim_percent
## 1        82.00814
## 2        89.85380
## 3        81.83889
## 4        89.32714
## 5        74.65689
## 6        91.62237

# Create a "long" form version of track.percent for use with plotting
track.percent.long <- track.percent %>% 
  tidyr::gather(key = "Step", value = "Percent", -Sample) 

# Modify track.percent.long so that Step column is a factor with levels in the correct order
track.percent.long[["Step"]] <- track.percent.long[["Step"]] %>% forcats::fct_relevel(colnames(track.percent[-1]))

# Plot percentages of input forward reads remaining after each step on a per sample basis.
p.track <- ggplot(track.percent.long, aes(x = Step, y = Percent)) + 
  geom_line(aes(group = Sample, colour = Sample)) +
  geom_point() +
  labs(title = "Percent of input forward reads remaining\nafter each processing step per sample",
       subtitle = paste0("Fwd read left trim position: ", FwdTrimLeft, "\nFwd read right trim position: ", FwdTrimRight),
       x = "Processing Step", y = "% of Input Reads Remaining") +
  theme(plot.title = element_text(hjust = 0.5, size = 18), 
        plot.subtitle = element_text(hjust = 0.5, size = 16),
        legend.position = "none",
        axis.text = element_text(size = 12), axis.title = element_text(size = 16))
p.track

# Save figure to file
ggsave(filename = file.path(pathOut, paste0("read_track_fwd_L", FwdTrimLeft, "_R", FwdTrimRight, "_", run.date, ".pdf")), 
       plot = p.track, device = "pdf", width = 8, height = 6, units = "in")

# Merge track and track.percent dataframes
track.merged <- dplyr::left_join(track, track.percent, by = c("Sample" = "Sample"))
head(track.merged)

##                        Sample  input filtered denoisedF nonchim input_percent
## 1 1872.022.rcbc214.2019.02.25 110022    99554     99301   90227           100
## 2 1907.032.rcbc226.2019.02.25  67306    61170     61091   60477           100
## 3 1954.043.rcbc238.2019.02.25 114558   100639    100248   93753           100
## 4 1996.054.rcbc250.2019.02.25  73017    68506     68342   65224           100
## 5 2051.064.rcbc262.2019.02.25 157822   136551    135715  117825           100
## 6 2140.110.rcbc320.2019.02.25  32169    29611     29529   29474           100
##   filtered_percent denoisedF_percent nonchim_percent
## 1         90.48554          90.25559        82.00814
## 2         90.88343          90.76605        89.85380
## 3         87.84982          87.50851        81.83889
## 4         93.82199          93.59738        89.32714
## 5         86.52216          85.99245        74.65689
## 6         92.04825          91.79334        91.62237

# Save log file of tracking reads to file.
print("Save track.merged dataframe tracking read counts as TSV file named:")

## [1] "Save track.merged dataframe tracking read counts as TSV file named:"

f.track.merged <- file.path(pathOut, paste0("MiSeq_", run.date, "_fwd_read_tracking_log.txt"))
head(f.track.merged)

## [1] "../results/dada2_output/MiSeq_2019_02_25_v1_fwd_read_tracking_log.txt"

readr::write_tsv(track.merged, f.track.merged)

Assign taxonomy

# RDP
# Record start time for this taxonomy assignment step
tax_start_time = proc.time()
print(paste0("RDP training database: ", "rdp_train_set_16.fa.gz"))
taxa.rdp <- assignTaxonomy(seqtab.nochim, file.path(pathdb, "rdp_train_set_16.fa.gz"), multithread = TRUE) 

# Add RDP species assignment
print(paste0("RDP species assignment database: ", "rdp_species_assignment_16.fa.gz"))
taxa.rdp.plus <- addSpecies(taxa.rdp, file.path(pathdb, "rdp_species_assignment_16.fa.gz"))
writeLines("RDP taxonomy assignment step processing time (seconds):")
proc.time() - tax_start_time

## [1] "RDP training database: rdp_train_set_16.fa.gz"
## [1] "RDP species assignment database: rdp_species_assignment_16.fa.gz"
## RDP taxonomy assignment step processing time (seconds):
##    user  system elapsed 
## 943.440   9.687 115.800

# Calculate cumulative processing time
print("Cumulative processing time (seconds):")
proc.time() - start_time

## [1] "Cumulative processing time (seconds):"
##     user   system  elapsed 
## 3612.008   80.130  789.518

Generate PhyloSeq Objects and final save

## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   .default = col_character(),
##   PCR_temp = col_double(),
##   primer_number = col_double()
## )
## ℹ Use `spec()` for the full column specifications.

## [1] "Sequencing run mapping file:"
## [1] "../data/Mapping/Mapping_info_Boston_run1.csv"
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 2099 taxa and 109 samples ]
## sample_data() Sample Data:       [ 109 samples by 23 sample variables ]
## tax_table()   Taxonomy Table:    [ 2099 taxa by 7 taxonomic ranks ]
##  [1] "Firmicutes"                "Actinobacteria"           
##  [3] "Fusobacteria"              "Bacteroidetes"            
##  [5] "Proteobacteria"            "Tenericutes"              
##  [7] NA                          "Deinococcus-Thermus"      
##  [9] "Cyanobacteria/Chloroplast" "Spirochaetes"             
## [11] "Chloroflexi"               "SR1"                      
## [13] "Chlamydiae"                "Synergistetes"            
## [15] "Euryarchaeota"             "Verrucomicrobia"          
## [17] "Planctomycetes"            "candidate_division_WPS-1" 
## [1] "Save R image (workspace) as .RData file named:"
## [1] "../results/dada2_output/RDS/MiSeq_2019_02_25_v1_preprocessing_image.RData"

# Calculate total processing time
print("Total processing time (seconds):")
proc.time() - start_time

## [1] "Total processing time (seconds):"
##     user   system  elapsed 
## 3686.430   80.553  864.476