library(dada2)
library(ShortRead)
library(dplyr)
library(tidyr)
library(ggplot2)

## Clean the environment / remove everything
rm(list = ls())
# Path to the raw data
path <- "<path to raw data (FASTQ)>"
## List of the files present in the path
list.files(path)
## Path to the data once the primers were detected and removed
path_cut <- "<path to raw data which have been processed by cutadapt>"
if (!dir.exists(path_cut)) dir.create(path_cut)
## Path to the data after filtration of the reads
path_filt <- "<path to data after being processed by filterAndTrim command (reads filtering)>"
if (!dir.exists(path_filt)) dir.create(path_filt)
## Create a directory to keep a sum up of the different steps done
path_sum <- "<path to the folder which gathers overview of each bioinformatic steps>"
if (!dir.exists(path_sum)) dir.create(path_sum)
## Create a directory for final outputs of the pipeline
path_results <- "<path to the folder which gathers final outputs>"
if (!dir.exists(path_results)) dir.create(path_results)

## Load primers couples CYA359F Forward primer
Fwd <- "GGGGAATYTTCCGCAATGGG"
Fwd_RC <- dada2::rc(Fwd)
# CYA781Rd Reverse primer
Rv <- "GACTACWGGGGTATCTAATCCCWTT"
Rv_RC <- dada2::rc(Rv)

## Load primers couples 23S_UPA Forward primer
Fwd <- "ACAGWAAGACCCTATGAAGCTT"
Fwd_RC <- dada2::rc(Fwd)
# 23S_UPA Reverse primer
Rv <- "CCTGTTATCCCTAGAGTAACTT"
Rv_RC <- dada2::rc(Rv)

## Assign different names to forward and reverse corresponding FASTQ files and
## extract samples names
raw_F_reads <- sort(list.files(path, pattern = "_R1.fastq.gz", full.names = TRUE))
raw_R_reads <- sort(list.files(path, pattern = "_R2.fastq.gz", full.names = TRUE))
## Prepare arguments of the command Outputs (data with primers removed)
cut_F_reads <- file.path(path_cut, basename(raw_F_reads))
cut_R_reads <- file.path(path_cut, basename(raw_R_reads))
# Arguments for primers detection (see Cutadapt website)
R1_flags <- paste(paste("-g", Fwd, collapse = " "), paste("-a", Rv_RC, collapse = " "))
R2_flags <- paste(paste("-G", Rv, collapse = " "), paste("-A", Fwd_RC, collapse = " "))
# Cutadapt command
cutadapt <- "cutadapt"  # Path to the executable
for (i in seq_along(raw_F_reads)) {
    cat("Processing", "-----------", i, "/", length(raw_F_reads), "-----------\n")
    system2(cutadapt, args = c(R1_flags, R2_flags, "--discard-untrimmed", "--max-n 0",
        "-o", cut_F_reads[i], "-p", cut_R_reads[i], raw_F_reads[i], raw_R_reads[i]))
}

## Count number of reads before and after cutadapt command
out_R1 <- cbind(ShortRead::qa(raw_F_reads)[["readCounts"]][, "read", drop = FALSE],
    ShortRead::qa(cut_F_reads)[["readCounts"]][, "read", drop = FALSE])
out_R2 <- cbind(ShortRead::qa(raw_R_reads)[["readCounts"]][, "read", drop = FALSE],
    ShortRead::qa(cut_R_reads)[["readCounts"]][, "read", drop = FALSE])
## Create and fill table with first outputs
SynthTab <- rbind.data.frame(out_R1, out_R2)
colnames(SynthTab)[1:2] <- c("raw", "cut")
SynthTab$samples <- sapply(rownames(SynthTab), FUN = function(x) strsplit(x, split = "_")[[1]][1])
SynthTab$reads <- sapply(rownames(SynthTab), FUN = function(x) {
    substr(strsplit(x, split = ".fastq")[[1]][1], nchar(strsplit(x, split = ".fastq")[[1]][1]) -
        1, nchar(strsplit(x, split = ".fastq")[[1]][1]))
})
## Reorder
SynthTab <- SynthTab[, c(3, 4, 1, 2)]

## Have a quick look only at the 2 first ones to have an idea
plotQualityProfile(cut_F_reads[1:2])
plotQualityProfile(cut_R_reads[1:2])
## To get the quality profiles from all the reads (from all the samples) in PDF
## doc
pdf(file.path(path_sum, "Read_quality_profile_aggregated.pdf"))
p <- plotQualityProfile(sample(cut_F_reads, replace = FALSE, size = ifelse(length(cut_F_reads) <
    100, length(cut_F_reads), 100)), aggregate = TRUE)
p + ggplot2::labs(title = "Forward")
p <- plotQualityProfile(sample(cut_R_reads, replace = FALSE, size = ifelse(length(cut_R_reads) <
    100, length(cut_R_reads), 100)), aggregate = TRUE)
p + ggplot2::labs(title = "Reverse")
dev.off()

## Extract sample names, assuming filenames have format: SAMPLENAME_XXX.fastq
sample_names <- sapply(strsplit(basename(raw_F_reads), "_"), function(x) x[1])
## Create names for output files
filt_F_reads <- file.path(path_filt, basename(raw_F_reads))
filt_R_reads <- file.path(path_filt, basename(raw_R_reads))
## Assign names to samples
names(filt_F_reads) <- sample_names
names(filt_R_reads) <- sample_names
## To merge our reads together with a sufficient overlapping length ()
out_2 <- filterAndTrim(cut_F_reads, filt_F_reads, cut_R_reads, filt_R_reads, truncLen = c(220,
    220), maxN = 0, maxEE = c(2, 2), truncQ = 2, rm.phix = TRUE, compress = TRUE,
    multithread = TRUE)
outFilt <- rbind(ShortRead::qa(filt_F_reads)[["readCounts"]][, "read", drop = FALSE],
    ShortRead::qa(filt_R_reads)[["readCounts"]][, "read", drop = FALSE])
## Add info to the synthesis table and reformat it
SynthTab <- merge(SynthTab, outFilt, by = 0)
rownames(SynthTab) <- SynthTab[, 1]
SynthTab <- SynthTab[, -1]
colnames(SynthTab)[5] <- "filtered"

## set seed to ensure that randomized steps are replicable
set.seed(100)
## Learn forward error rates
error_F <- learnErrors(filt_F_reads, multithread = TRUE, randomize = TRUE)
## Learn reverse error rates
error_R <- learnErrors(filt_R_reads, multithread = TRUE, randomize = TRUE)
## To save the error rates estimations in a PDF doc
pdf(file.path(path_sum, "Error_rates_learning.pdf"))
p <- plotErrors(error_F, nominalQ = TRUE)
p + ggplot2::labs(title = "Error Forward")
p <- plotErrors(error_R, nominalQ = TRUE)
p + ggplot2::labs(title = "Error Reverse")
dev.off()

## Here we work with list, with samples as elements, containing their relative
## outputs
dada_F_reads <- list()
dada_R_reads <- list()
merged_reads <- list()
SynthTab$uniq <- NA
for (r in 1:length(sample_names)) {
    cat("Processing (", r, "/", length(sample_names), ") :", sample_names[r], "\n")
    # Dereplication step
    derepF_tmp <- derepFastq(filt_F_reads[r], verbose = 1)
    derepR_tmp <- derepFastq(filt_R_reads[r], verbose = 1)
    # Sample inference (denoising)
    dada_F_reads[[sample_names[r]]] <- dada(derepF_tmp, err = error_F, multithread = TRUE,
        verbose = 1)
    dada_R_reads[[sample_names[r]]] <- dada(derepR_tmp, err = error_R, multithread = TRUE,
        verbose = 1)
    # Merging the reads together
    merged_reads[[sample_names[r]]] <- mergePairs(dadaF = dada_F_reads[[sample_names[r]]],
        derepF = derepF_tmp, dadaR = dada_R_reads[[sample_names[r]]], derepR = derepR_tmp,
        verbose = 1)
    # To fill synthesis table
    SynthTab[which(grepl(sample_names[r], SynthTab$samples) & SynthTab$reads == "R1"),
        "uniq"] <- length(derepF_tmp$uniques)
    SynthTab[which(grepl(sample_names[r], SynthTab$samples) & SynthTab$reads == "R2"),
        "uniq"] <- length(derepR_tmp$uniques)
}
## Keep ASV with a given overlapping length
for (i in 1:length(merged_reads)) {
    merged_reads[[i]] <- merged_reads[[i]] %>%
        filter(nmatch > 35 & nmatch < 45)
}
## Create from the list the sequence table
seqtab <- makeSequenceTable(merged_reads)
## Add infos to the synthesis table and reformat it
SynthTab$merged <- NA
for (i in sample_names) {
    SynthTab[which(grepl(i, SynthTab$samples)), "merged"] <- nrow(merged_reads[[i]])
}

seqtab.nochim <- removeBimeraDenovo(seqtab, method = "consensus", multithread = TRUE,
    verbose = TRUE)
## never stop being curious
dim(seqtab.nochim)  # number of samples * number of ASV (- chimera)
table(nchar(getSequences(seqtab)))  # length distribution of the ASV (- chimera)
sum(seqtab.nochim)/sum(seqtab) * 100  # % of kept ASV to have an idea of the chimera proportion
## Add infos to the synthesis table and reformat it
SynthTab$no_chimera <- NA
for (i in sample_names) {
    SynthTab[which(grepl(i, SynthTab$samples)), "no_chimera"] <- length(which(seqtab.nochim[i,
        ] > 0))
}