#análisis en base al tutorial de dada2 que está disponible en:
# https://benjjneb.github.io/dada2/tutorial.html
#setear la carpeta de trabajo (directory containing the fastq files)
setwd("C:/......................")
setwd("~/Library/CloudStorage/GoogleDrive-florbertoglio@gmail.com/Mi unidad/Curso genomica 2024") #ejemplo en mi caso
#instalar los siguientes paquetes:
install.packages("Rcpp") #en mi caso me pide que instale este paquete para luego poder instalar dada2.
#Esto depende de el sistema operativo de cada pc y de la version de R que se tenga
install.packages("dada2")
install.packages("phyloseq")
install.packages("Biostring")
install.packages("vegan")
library(Rcpp)
library(dada2)
library(phyloseq)
library(Biostrings)
library(vegan)

#### PROCESAMIENTO BIOINFORMÁTICO####
# crear la variabe "path" con el camino donde estan los fastq originales
path <-"~/Library/CloudStorage/GoogleDrive-florbertoglio@gmail.com/Mi unidad/Curso genomica 2024" 
list.files(path) #arroja la lista de documentos presentes en la carpeta especificada en el path
# Los archivos fastq forward y reverse tienen el formato: SAMPLENAME_R1.fastq y SAMPLENAME_R2.fastq
#R1 corresponde a los reads forward y R2 a los reverse
#para obtener listas coincidentes de los archivos fastq forward y reverse
fnFs <- sort(list.files(path, pattern="_R1.fastq", full.names = TRUE)) #genero una variable que contiene los nombres de los reads forward
fnRs <- sort(list.files(path, pattern="_R2.fastq", full.names = TRUE)) #genero una variable que contiene los nombres de los reads reverse

# Extraer los nombres de las muestras en la variable "sample.names" asumiendo que tienen el formato SAMPLENAME_XXX.fastq 
# Siempre chequear el nombre de los fastq antes de comenzar para asegurar que tengan dicho formato
sample.names <- sapply(strsplit(basename(fnFs), "_"), `[`, 1) #genero una variable con todos los nombres de los reads
View(sample.names)

# 1...........CHEQUEAR LA CALIDAD DE LOS READS...............................
plotQualityProfile(fnFs[1:2]) #para ver la calidad de los reads forward
plotQualityProfile(fnRs[1:2]) #para ver la calidad de los reads reverse

#2...............FILTRAR Y TRUNCAR (Filter and Trim).......................
#Crear subdirectorio, en este caso se va a llamar: filtered
#Assign the filenames for the filtered fastq.gz files
filtFs <- file.path(path, "filtered", paste0(sample.names, "_F_filt.fastq.gz")) #genero una variable que contiene los nombres de los reads forward filtrados que voy a generar mas abajo
filtRs <- file.path(path, "filtered", paste0(sample.names, "_R_filt.fastq.gz")) #mismo pero para los reads reverse
names(filtFs) <- sample.names
names(filtRs) <- sample.names
out <- filterAndTrim(fnFs, filtFs, fnRs, filtRs, truncLen=c(260,200),trimLeft=10,
                     maxN=0, maxEE=c(2,2), truncQ=2, rm.phix=TRUE,
                     compress=TRUE, multithread=TRUE)
# standard filtering parameters: maxN=0 (dada no acepta Ns), truncQ=2, rm.phix=TRUE (discard reads that match against the phiX genome), maxEE=2 
#el ouput se guardo en las variables creadas previamente (filtFs y filtRs)
#tambien generamos un objeto out el cual es una matriz que contiene los reads iniciales
#y los reads que quedaron luego de pasar el filtro
head(out) #muestra los primeros datos de dicha matriz

#3............... APRENDE LAS TASAS DE ERROR (learn the error rate): es el algoritmo que emplea dada2
#learnErrors: genera un modelo de error de nuestros datos
errF <- learnErrors(filtFs, multithread=FALSE)
errR <- learnErrors(filtRs, multithread=FALSE)
plotErrors(errF, nominalQ=TRUE) # In the plots, the black line is the error model, the dots are the actual errors
plotErrors(errR, nominalQ=TRUE) # In the plots, the black line is the error model, the dots are the actual errors

#4......INFERENCIA DE LAS MUESTRAS (SAMPLE INFERENCE).......#
#A partir del modelo de error generado infiere las verdaderas secuencias biológicas de los datos  
dadaFs <- dada(filtFs, err=errF, multithread=TRUE) 
dadaRs <- dada(filtRs, err=errR, multithread=TRUE)

#5.....FUSIÓN DE AMBOS READS EN UNA MISMA SECUENCIA (MERGE PAIRED READS).......#
#en este paso se fusionan ambos reads (mediante una region overlap) en un único read final que abarca la longitud del amplicon original 
#esto permite corregir errores y determinar con una mayor precisión la cobertura de lectura
mergers <- mergePairs(dadaFs, filtFs, dadaRs, filtRs, verbose=TRUE)

# Chequear el merger data.frame para la primer muestra
head(mergers[[1]])

#6....GENERAR UNA TABLA DE SECUENCIAS (ASV)....(CONTRUCT SEQUENCE TABLES).........#
seqtab <- makeSequenceTable(mergers) # The sequences being tabled vary in length.
dim(seqtab)
colnames(seqtab)
table(nchar(getSequences(seqtab))) # Inspect distribution of sequence lengths
#Dado que los primers usados para secuenciar la región V4 fueron 806R y 515F: 806-515=291 pb (largo de region secuenciada),
#saco secuencias que sean mas cortas que 250 (no hay mas cortas) y mas largas que 300.
seqtab2 <- seqtab[,nchar(colnames(seqtab)) %in% 250:300] 
dim(seqtab2) #dim(seqtab2) < dim(seqtab)
table(nchar(getSequences(seqtab2))) # Inspect distribution of sequence lengths
saveRDS(seqtab2, "seqtab.rds") #guarda la tabla en el directorio de trabajo

#7....IDENTIFICACIÓN Y REMOCIÓN DE QUIMERAS (REMOVE CHIMERAS).......#
seqtab.nochim <- removeBimeraDenovo(seqtab2, method="consensus", multithread=FALSE, verbose=TRUE) 
dim(seqtab.nochim) 
sum(seqtab.nochim)/sum(seqtab2) #proporción de secuencias con las que se quedó luego de eliminar las quimeras: eliminó un 4% 
1-(sum(seqtab.nochim)/sum(seqtab2)) #proporción de secuencias que son quimeras
1- (dim(seqtab.nochim)[2]/dim(seqtab2)[2]) #proporción de ASVs que son quimeras
saveRDS(seqtab.nochim, "seqtab.nochim.rds")

#8......RESUMEN DE LOS RECUENTOS EN TODOS LOS PASOS.....#
getN <- function(x) sum(getUniques(x))
track <- cbind(out, sapply(dadaFs, getN), sapply(dadaRs, getN), sapply(mergers, getN), rowSums(seqtab.nochim))
colnames(track) <- c("input", "filtered", "denoisedF", "denoisedR", "merged", "nonchim")
rownames(track) <- sample.names
head(track)
write.table(track, "Track")


#9...ASIGNACIÓN DE TAXONOMÍA (ASSIGN TAXONOMY)..........
#Descargar la base de datos de SILVA para el 16S y dejarla fija en una carpeta, por ej en descargas en mi caso
#página para descargarla: https://zenodo.org/record/4587955#.YUNGny2xBQI
#El nombre del archivo a descargar es silva_nr99_v138.1_train_set.fa.gz (pueden haber versiones más recientes ya que se actualiza con frecuencia)
taxa <- assignTaxonomy(seqtab.nochim, "~/Downloads/silva_nr99_v138.1_train_set.fa.gz", multithread=FALSE, tryRC=TRUE) #tryRC: the reverse-complement of each sequences
#Chequear algunos asignaciones taxonómicas:
taxa.print <- taxa # Removing sequence rownames for display only
rownames(taxa.print) <- NULL
head(taxa.print)
taxa.print <- taxa # Removing sequence rownames for display only
rownames(taxa.print) <- NULL
head(taxa.print)
View(taxa.print)