dada-pipeline.R

#!/usr/bin/env Rscript

library(dada2)
library(ggplot2)

pdf('results.pdf')

path <- "reads/"

fns <- list.files(path)
fns

fastqs <- fns[grepl(".fastq$", fns)]
fastqs <- sort(fastqs) # Sort ensures forward/reverse reads are in same order

fnFs <- fastqs[grepl("_R1", fastqs)] # Just the forward read files
fnRs <- fastqs[grepl("_R2", fastqs)] # Just the reverse read files

# Get sample names, assuming files named as so: SAMPLENAME_XXX.fastq
sample.names <- sapply(strsplit(fnFs, "_"), `[`, 1)

# Specify the full path to the fnFs and fnRs
fnFs <- file.path(path, fnFs)
fnRs <- file.path(path, fnRs)

plotQualityProfile(fnFs[[1]])
plotQualityProfile(fnRs[[1]])

filt_path <- file.path(path, "filtered")
if(!file_test("-d", filt_path)) dir.create(filt_path)

filtFs <- file.path(filt_path, paste0(sample.names, "_F_filt.fastq.gz"))
filtRs <- file.path(filt_path, paste0(sample.names, "_R_filt.fastq.gz"))

# Filter
for(i in seq_along(fnFs)) {
      fastqPairedFilter(
                        c(fnFs[i],
                          fnRs[i]),
                        c(filtFs[i],
                          filtRs[i]),
                        truncLen=c(240,160),
                        trimLeft=20, # trim adapter sequences (first ~20 bases)
                        maxN=0,
                        maxEE=c(2,2),
                        truncQ=2,
                        rm.phix=TRUE,
                        compress=TRUE,
                        verbose=TRUE)
}

derepFs <- list(derepFastq(filtFs, verbose=TRUE))
derepRs <- list(derepFastq(filtRs, verbose=TRUE))

# Name the derep-class objects by the sample names
names(derepFs) <- sample.names
names(derepRs) <- sample.names

# learn the error rates
dadaFs.lrn <- list(dada(derepFs, err=NULL, selfConsist = TRUE, multithread=TRUE))
errF <- dadaFs.lrn[[1]]$err_out

plotErrors(dadaFs.lrn[[1]], nominalQ=TRUE)

dadaRs.lrn <- list(dada(derepRs, err=NULL, selfConsist = TRUE, multithread=TRUE))
errR <- dadaRs.lrn[[1]]$err_out

plotErrors(dadaFs.lrn[[1]], nominalQ=TRUE)

# sample inference
dadaFs <- dada(derepFs, err=errF, multithread=TRUE)
dadaRs <- dada(derepRs, err=errR, multithread=TRUE)

# merge paired reads
mergers <- mergePairs(dadaFs, derepFs, dadaRs, derepRs, verbose=TRUE)

# write mergers.. we need this for loading into dada2 later
saveRDS(mergers, file='mergers.RData')

# Inspect the merger data.frame from the first sample

head(mergers[[1]])

# construct the sequence table
seqtab <- makeSequenceTable(mergers)

# Inspect distribution of sequence lengths
table(nchar(getSequences(seqtab)))

# remove chimeras
seqtab.nochim <- removeBimeraDenovo(seqtab, verbose=TRUE)
dim(seqtab.nochim)

# fraction of chimeras
sum(seqtab.nochim)/sum(seqtab)

# assign taxonomy
taxa <- assignSpecies(seqtab.nochim, "rdp_species_assignment_14.fa.gz")

message("saving plots to results.pdf")
dev.off()

message("saving sequence table to seqtab.csv")
write.csv(seqtab.nochim, "seqtab.csv")

message("saving results to taxa.csv")
write.csv(taxa, "taxa.csv")