dada2_workflow.Rmd

---
title: "DADA2 tutorial"
output: html_notebook
---

Example of DADA2 tutorial can be found from [here](https://benjjneb.github.io/dada2/tutorial.html).
This example include only the code. The pipeline ends to _TreeSummarizedExperiment_.

We recommend to use the most recent reference database. You can find it from
[here](https://benjjneb.github.io/dada2/training.html).

```{r dada2_1}
# GETTING READY

library(dada2); packageVersion("dada2")
```

```{r dada2_2}
path <- "data/MiSeq_SOP" # CHANGE ME to the directory containing the fastq files after unzipping.
list.files(path)
```

```{r dada2_3}
# Forward and reverse fastq filenames have format: SAMPLENAME_R1_001.fastq and SAMPLENAME_R2_001.fastq
fnFs <- sort(list.files(path, pattern="_R1_001.fastq", full.names = TRUE))
fnRs <- sort(list.files(path, pattern="_R2_001.fastq", full.names = TRUE))
# Extract sample names, assuming filenames have format: SAMPLENAME_XXX.fastq
sample.names <- sapply(strsplit(basename(fnFs), "_"), `[`, 1)
```

```{r dada2_4}
# INSPECT READ QUALITY PROFILES

plotQualityProfile(fnFs[1:2])
```

```{r dada2_5}
plotQualityProfile(fnRs[1:2])
```

```{r dada2_6}
# FILTER AND TRIM

# Place filtered files in filtered/ subdirectory
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
```

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

```{r dada2_8}
# LEARN THE ERROR RATES

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

```{r dada2_9}
errR <- learnErrors(filtRs, multithread=TRUE)
```

```{r dada2_10}
plotErrors(errF, nominalQ=TRUE)
```

```{r dada2_11}
# SAMPLE INFERENCE

dadaFs <- dada(filtFs, err=errF, multithread=TRUE)
```

```{r dada2_12}
dadaRs <- dada(filtRs, err=errR, multithread=TRUE)
```

```{r dada2_13}
dadaFs[[1]]
```

```{r dada2_14}
# MERGE PAIRED READS

mergers <- mergePairs(dadaFs, filtFs, dadaRs, filtRs, verbose=TRUE)
# Inspect the merger data.frame from the first sample
head(mergers[[1]])
```

```{r dada2_15}
# CONSTRUCT SEQUENCE TABLE

seqtab <- makeSequenceTable(mergers)
dim(seqtab)
```

```{r dada2_16}
# Inspect distribution of sequence lengths
table(nchar(getSequences(seqtab)))
```

```{r dada2_17}
# REMOVE CHIMERAS

seqtab.nochim <- removeBimeraDenovo(seqtab, method="consensus", multithread=TRUE, verbose=TRUE)
dim(seqtab.nochim)
```

```{r dada2_18}
sum(seqtab.nochim)/sum(seqtab)
```

```{r dada2_19}
# TRACK READS THROUGH THE PIPELINE

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
head(track)
```

```{r dada_20}
# ASSIGN TAXONOMY

# You can find updated reference databases from 
# https://benjjneb.github.io/dada2/training.html

# Change path
taxa <- assignTaxonomy(seqtab.nochim, "data/silva_nr_v132_train_set.fa.gz", multithread=TRUE)
```

```{r dada2_21}
# Change path
taxa <- addSpecies(taxa, "data/silva_species_assignment_v132.fa.gz")
```

```{r dada2_22}
taxa.print <- taxa # Removing sequence rownames for display only
rownames(taxa.print) <- NULL
head(taxa.print)
```

```{r dada2_23}
# EVALUATE ACCURACY

unqs.mock <- seqtab.nochim["Mock",]
unqs.mock <- sort(unqs.mock[unqs.mock>0], decreasing=TRUE) # Drop ASVs absent in the Mock
cat("DADA2 inferred", length(unqs.mock), "sample sequences present in the Mock community.\n")
```

```{r dada2_24}
mock.ref <- getSequences(file.path(path, "HMP_MOCK.v35.fasta"))
match.ref <- sum(sapply(names(unqs.mock), function(x) any(grepl(x, mock.ref))))
cat("Of those,", sum(match.ref), "were exact matches to the expected reference sequences.\n")
```

## Create a TreeSE data object

```{r dada2_25}
library(mia)
library(ggplot2)

if( !require("BiocManager") ){
    install.packages("BiocManager")
    library("BiocManager")
}

if( !require("Biostrings") ){
    BiocManager::install("Biostrings")
    library("Biostrings")
}
library(Biostrings)
```

```{r dada2_26}
samples.out <- rownames(seqtab.nochim)
subject <- sapply(strsplit(samples.out, "D"), `[`, 1)
gender <- substr(subject,1,1)
subject <- substr(subject,2,999)
day <- as.integer(sapply(strsplit(samples.out, "D"), `[`, 2))
samdf <- data.frame(Subject=subject, Gender=gender, Day=day)
samdf$When <- "Early"
samdf$When[samdf$Day>100] <- "Late"
rownames(samdf) <- samples.out
```

```{r dada2_27}
tse <- TreeSummarizedExperiment(assays = SimpleList(counts = t(seqtab.nochim)),
                                colData = DataFrame(samdf),
                                rowData = DataFrame(taxa)
                                )

# Remove mock sample like it is also done in DADA2 pipeline tutorial
tse <- tse[ , colnames(tse) != "mock"]
```

```{r dada2_38}
dna <- Biostrings::DNAStringSet( rownames(tse) )
referenceSeq(tse) <- dna
rownames(tse) <- paste0("ASV", seq( nrow(tse) ))
tse
```