---
title: "TD script DADA2 du 05/12/2023"
output: html_document
date: "2023-11-22"
editor_options: 
  chunk_output_type: inline
---

```{r setup, include=FALSE}
setwd(Sys.getenv("HOME"))
knitr::opts_chunk$set(echo = TRUE)
```

#Packages
```{r}

```


#Library
```{r}
Start=Sys.time()
library(dada2)
library(phyloseq)
library(tidyverse)
library(magrittr)
library(ggpubr)
library(DT)
library(ade4)
library(doParallel)
library('readxl')
library(DECIPHER)
library(ggplot2)
```
#Start
##Import fastq files
```{r}
# setwd(Sys.getenv("HOME"))

path =  list.files("~/shared/LU3SV530 - Server/2019-11 Sci Rep Villette Scarcity Paper/", pattern = c("R1_001", "R2_001"))

```

##forward and reverse fastq
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}
#Obtenir les fichiers FASTQ "forward" (sens) et "reverse" (antisens)
# setwd(Sys.getenv("HOME"))

fnFs <- sort(list.files("~/shared/LU3SV530 - Server/2019-11 Sci Rep Villette Scarcity Paper/", pattern="_R1_001.fastq", full.names = TRUE))
fnRs <- sort(list.files("~/shared/LU3SV530 - Server/2019-11 Sci Rep Villette Scarcity Paper/", pattern="_R2_001.fastq", full.names = TRUE))

#Extraire un nom d'échantillon à partir du nom de fichier 

sample.names <- sapply(strsplit(basename(fnFs), "-"), `[`, 1)

```

```{r}
list_fastq= function(path, pattern, separator="_", level){

f=list.files(path)
# is= sum(grepl(pattern = pattern, f))
  if(length(pattern)==1){
   f= list.files(path, full.names = T)
   f= f[grepl(pattern = "fastq|fastq.gz", x = basename(f))]
   n= file.exists(f) %>% sum()
   message(paste(n, "fastq are single end or merged"))

   tmp <- str_split(basename(f), pattern = separator, simplify = T)[,level]

   return(list(fastq=f, names= tmp))
 } else {
   fwd=  list.files(path=path, pattern =  pattern[1], full.names = T)
   rv= list.files(path = path, pattern =  pattern[2],full.names = T)
   fwd= fwd[grepl(pattern = "fastq|fastq.gz", x = basename(fwd))]
   rv= rv[grepl(pattern = "fastq|fastq.gz", x = basename(rv))]

   n1= file.exists(fwd) %>% sum()
   n2= file.exists(rv) %>% sum()
   message(paste(n1, n2, "fastq are double end"))

   tmp1 <- str_split(basename(fwd), pattern=separator, simplify = T)[,level]
   tmp2 <- str_split(basename(rv), pattern=separator, simplify = T)[,level]

   return(list(fastq_fwd= fwd, fastq_rv= rv, names= tmp1))

 }

}
```


```{r}
# setwd(Sys.getenv("HOME"))

f_list = list_fastq("~/shared/LU3SV530 - Server/2019-11 Sci Rep Villette Scarcity Paper/", pattern = c("R1_001", "R2_001"), separator = "-")
```


#Inspection de la qualité des séquences
```{r}
#Pour les 10 premiers échantillons "forward" et "reverse" seulement
plotQualityProfile(fnFs[1:10])
plotQualityProfile(fnRs[1:10])

#Autre alternative : aggréger tous les fichiers fastq "forworde" ensemble, ainsi que les fichiers fastq "reverse".
p1 = plotQualityProfile(fnFs, aggregate = T)
p2 = plotQualityProfile(fnRs, aggregate = T)
ggarrange(p1, p2)
```

#Filtrage et tronquage

```{r}
filtFs = file.path(paste0(Sys.getenv("HOME"),"/filtered"), paste0(sample.names, "_F_filt.fastq.gz"))
filtRs = file.path(paste0(Sys.getenv("HOME"),"/filtered"), paste0(sample.names, "_R_filt.fastq.gz"))
names(filtFs) <- sample.names
names(filtRs) <- sample.names
```


```{r}
# setwd(Sys.getenv("HOME"))

out <- filterAndTrim(fwd = fnFs, 
                     filt = filtFs, 
                     rev = fnRs, 
                     filt.rev = filtRs, 
                     trimLeft=20, 
                     truncLen = c(240,240), 
                     maxN=0, 
                     maxEE=c(2,2), 
                     truncQ = 2, 
                     rm.phix=TRUE, 
                     compress=TRUE, 
                     multithread=FALSE) 
                    
head(out)

```


# Proportion of sequences after filtering
```{r}
(mean(out[,2])/mean(out[,1]))*100
```

#Estimation du taux d'erreur
```{r}
errF <- learnErrors(filtFs, multithread=TRUE, MAX_CONSIST = 10, OMEGA_C = 0)
errR <- learnErrors(filtRs, multithread=TRUE, MAX_CONSIST = 25, OMEGA_C = 0)

p3=plotErrors(errF, nominalQ=TRUE)
p4=plotErrors(errR, nominalQ=TRUE)

ggarrange(p3,p4)
```

#Inférence des échantillons

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

#Déréplication
```{r}
derepFs <- derepFastq(filtFs)
names(derepFs) <- sample.names

derepRs <- derepFastq(filtRs)
names(derepRs) <- sample.names
```


#Fusion des séquences pairées

```{r}
mergers <- mergePairs(dadaFs, derepFs, dadaRs, derepRs, verbose=TRUE)
```


#Construct sequence table

```{r}
seqtab <- makeSequenceTable(mergers)

table(seqtab)

```


```{r}
dim(seqtab)
```



```{r}
table(nchar(getSequences(seqtab)))
```

#Inspecter la longueur des ASVs 
```{r}
hist(nchar(getSequences(seqtab)),xlab="Size", ylab="Frequency", main = "ASVs length") 
```


#Suppression des chimères
```{r}
seqtab.nochim <- removeBimeraDenovo(seqtab, method="consensus", multithread=TRUE, verbose=TRUE)
```

#Proportion kept after chimera removal
```{r}
sum(seqtab.nochim)/sum(seqtab)
```

#Tableau de suivi
```{r}
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)
```



#Assign taxonomy
```{r}
# setwd(Sys.getenv("HOME"))

taxa = assignTaxonomy(seqtab.nochim, "~/shared/LU3SV530 - Server/2019-11 Sci Rep Villette Scarcity Paper/Metadata/silva_nr_v138_train_set.fa.gz", multithread = 60, verbose = T, tryRC = T)
taxa= addSpecies(taxtab = taxa, "~/shared/LU3SV530 - Server/2019-11 Sci Rep Villette Scarcity Paper/Metadata/silva_species_assignment_v138.fa.gz", verbose = T, allowMultiple = T, tryRC = T )

taxa= as.data.frame(taxa)

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

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


#Créer l'objet Phyloseq
```{r}
# setwd(Sys.getenv("HOME"))

metadata_16S= read_excel(paste0(Sys.getenv("HOME"),"/shared/LU3SV530 - Server/2019-11 Sci Rep Villette Scarcity Paper/Metadata/2021-05 VILLETTE et al 16s Refinement Metadata.xlsx"), sheet=1, col_names = T) %>% as.data.frame()

metadata_16S=metadata_16S%>%
  dplyr::filter(SampleOrigin=="MCS Whole Cell" | SampleOrigin=="MCS Genomic")

metadata_16S=metadata_16S %>%
  separate(col= Filename, into=letters[1], sep="-") 

metadata_16S= metadata_16S[order(metadata_16S$a),]
rownames(metadata_16S)= metadata_16S$a

ps_16S= phyloseq(otu_table(seqtab.nochim, taxa_are_rows = F), 
                 sample_data(metadata_16S), 
                 tax_table(as.matrix(taxa)))

saveRDS(ps_16S, paste0(Sys.getenv("HOME"),"/Phyloseq_TD.rds"))
End=Sys.time()
End-Start
```


