---
title: "R DADA2 16S Pipeline Tutorial"
author: "John Quensen"
date: "2026-03-30"
output: pdf_document
editor_options: 
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

Load the dada2 library.  

```{r}
library(dada2)
library(tidyverse)
```

Set the path to the trimmed data. I will use it because there is no primer trimming included in this script.  

```{r}
setwd("//wsl.localhost/Ubuntu/home/johnq/novogene_new")
path <- "trimmed"
list.files(path)
```

Now we read in the names of the fastq files, and perform some string manipulation to get matched lists of the forward and reverse fastq files.  

```{r}
# Forward and reverse fastq file names have format: SAMPLENAME_R1.fastq and SAMPLENAME_R2.fastq
fnFs <- sort(list.files(path, pattern="_R1.fastq", full.names = TRUE))
fnRs <- sort(list.files(path, pattern="_R2.fastq", full.names = TRUE))
# Extract sample names, assuming file names have format: SAMPLENAME_XXX.fastq
sample.names <- sapply(strsplit(basename(fnFs), "_"), `[`, 1)
sample.names
```

# Inspect read quality profiles

First the forward reads.  

```{r}
png("forward_quality_scores.png", width = 900, height = 600)
plotQualityProfile(fnFs[1:2])
dev.off()
```

And then the same for the reverse reads.  

```{r}
png("reverse_quality_scores.png", width = 900, height = 600)
plotQualityProfile(fnRs[1:2])
dev.off()
```

# Filter and trim

Assign the file names for the filtered fastq.gz files.  

```{r}
# Place filtered files in filtered/ sub-directory
filtFs <- file.path(path, "filtered", paste0(sample.names, "_F_filt.fastq.gz"))
filtRs <- file.path(path, "filtered", paste0(sample.names, "_R_filt.fastq.gz"))
names(filtFs) <- sample.names
names(filtRs) <- sample.names
```

The standard filtering parameters are: maxN=0 (DADA2 requires no Ns), truncQ=2, rm.phix=TRUE and maxEE=2. The maxEE parameter sets the maximum number of “expected errors” allowed in a read, which is a better filter than simply averaging quality scores.  

I chose these parameters from the quality plots and the expected errors from USEARCH -fastx_info output. I am using the trimmed sequences because trimming is not included in this script.  

```{r}
out <- filterAndTrim(fnFs, filtFs, fnRs, filtRs, 
                     maxN=0, maxEE=c(2, 2), truncQ=2, rm.phix=TRUE,
              compress=TRUE, multithread=FALSE) # On Windows set multithread=FALSE
head(out)
```

# Learn error rates

If we used the "standard" method, not allowing for binned quality scores, we woud use the following script to learn the error rates. I do this now to generate plots of the transition error rates.  

```{r}
errF <- learnErrors(filtFs, multithread=TRUE)

errR <- learnErrors(filtRs, multithread=TRUE)

plt <- plotErrors(errF, nominalQ=TRUE)
plt2 <- plt +
  theme_bw() +
  theme(strip.background = element_rect(fill = "white"))
plt2

ggsave(filename = "R_dada2_novogene_default_forward_error_plot.png",
       width = 6, height = 4, units = "in", dpi = 300)

plt <- plotErrors(errR, nominalQ=TRUE)
plt2 <- plt +
  theme_bw() +
  theme(strip.background = element_rect(fill = "white"))
plt2

ggsave(filename = "R_dada2_novogene_default_reverse_error_plot.png",
       width = 6, height = 4, units = "in", dpi = 300)

```

Because the Q-scores are binned, we should instead use Callahan's makeBinnedQualErrFun to learn the transition error rates.  

```{r}
novaseqBinnedErrFun <- makeBinnedQualErrfun(c(2, 11, 25, 37))

errF <- learnErrors(filtFs, multithread=TRUE,
                    errorEstimationFunction = novaseqBinnedErrFun)

errR <- learnErrors(filtRs, multithread=TRUE,
                    errorEstimationFunction = novaseqBinnedErrFun)

png("novogene_R_forward_transition_error_rates.png", width = 1800, height = 1200)
plotErrors(errF, nominalQ=TRUE)
dev.off()

png("novogene_R_reverse_transition_error_rates.png", width = 1800, height = 1200)
plotErrors(errR, nominalQ=TRUE)
dev.off()

plt <- plotErrors(errF, nominalQ=TRUE)
plt2 <- plt +
  theme_bw() +
  theme(strip.background = element_rect(fill = "white"))
plt2
saveRDS(plt2, "novogene_R_forward_error_plot.rds")
ggsave(filename = "R_dada2_novogene_forward_error_plot.png",
       width = 6, height = 4, units = "in", dpi = 300)

plt <- plotErrors(errR, nominalQ=TRUE)
plt2 <- plt +
  theme_bw() +
  theme(strip.background = element_rect(fill = "white"))
plt2
saveRDS(plt2, "novogene_R_reverse_error_plot.rds")
ggsave(filename = "R_dada2_novogene_reverse_error_plot.png",
       width = 6, height = 4, units = "in", dpi = 300)
```

# Sample inference

```{r}
dadaFs <- dada(filtFs, err=errF, pool = "pseudo", multithread=TRUE)
dadaRs <- dada(filtRs, err=errR, pool = "pseudo",  multithread=TRUE)
dadaFs[[1]]
```

# Merge paired reads.  

```{r}
mergers <- mergePairs(dadaFs, filtFs, dadaRs, filtRs,
                      minOverlap = 10, 
                      verbose=TRUE, maxMismatch = 5,
                      trimOverhang = TRUE)
# Inspect the merger data.frame from the first sample
head(mergers[[1]])
```

# Construct sequence table

```{r}
seqtab <- makeSequenceTable(mergers)
dim(seqtab)
# Inspect distribution of sequence lengths
table(nchar(getSequences(seqtab)))
```

# Remove chimeras

```{r}
seqtab.nochim <- removeBimeraDenovo(seqtab,
                                    method="consensus",
                                    multithread=TRUE,
                                    verbose=TRUE)
dim(seqtab.nochim)
sum(seqtab.nochim)/sum(seqtab)
saveRDS(seqtab.nochim, file = "asv_table.rds")
```

# Track reads through pipeline  

```{r}
getN <- function(x) sum(getUniques(x))
track <- cbind(out, sapply(dadaFs, getN), sapply(dadaRs, getN), sapply(mergers, getN), rowSums(seqtab.nochim))
# If processing a single sample, remove the sapply calls: e.g. replace sapply(dadaFs, getN) with getN(dadaFs)
colnames(track) <- c("input", "filtered", "denoisedF", "denoisedR", "merged", "nonchim")
rownames(track) <- sample.names
track
class(track)
write.table(track, file = "R_mod_dada2_retention.tsv", sep = "\t")
```

Modify track to contain same columns as the QIIME2/DADA2 version of the retention table.

```{r}
retention <- track |> 
  as.data.frame() |> 
  rownames_to_column(var = "sample") |> 
  as_tibble() |> 
  rename(denoised = denoisedF) |> 
  dplyr::mutate(pc_non_chimeric = 100*(nonchim/input),
                pc_input_merged = 100 * merged / input,
                pc_input_filtered = 100 * filtered / input
                ) |> 
  dplyr::select(sample, input, filtered, pc_input_filtered, denoised, merged,
         pc_input_merged, nonchim, pc_non_chimeric)
retention
write_tsv(retention, file = "retention_R_dada2_binQ_merge_5mis_nostagger.tsv")
```

