title: “Lydia Yanes, tutoriel dada2 pour ADM” output: github_document —
# je déclare l'entete du notebook Rmarkdown, définit le titre du document et le format de sortie, et genere un rendu qui sera compatible avec github, et qui pourra etre affiché directement dessus#je charge le package dada2 que je vais utiliser dans l'analyse des barcodes ARNr16s du tutoriel, et je verifie la version du package pour pouvoir reproduire les résultats lors de l'analyse des donnés de mon article :
library(dada2); packageVersion("dada2")## Loading required package: Rcpp
## [1] '1.28.0'#je définit le chemin vers le dossier où sont stockés mes fichier fastq, et je lis tout les fichiers qui sont dans ce dossier pour verifier qu'ils sont bien présents avant de continuer mon analyse :
path <- "~/dada2/data/MiSeq_SOP"
list.files(path)## [1] "F3D0_S188_L001_R1_001.fastq" "F3D0_S188_L001_R2_001.fastq"
## [3] "F3D1_S189_L001_R1_001.fastq" "F3D1_S189_L001_R2_001.fastq"
## [5] "F3D141_S207_L001_R1_001.fastq" "F3D141_S207_L001_R2_001.fastq"
## [7] "F3D142_S208_L001_R1_001.fastq" "F3D142_S208_L001_R2_001.fastq"
## [9] "F3D143_S209_L001_R1_001.fastq" "F3D143_S209_L001_R2_001.fastq"
## [11] "F3D144_S210_L001_R1_001.fastq" "F3D144_S210_L001_R2_001.fastq"
## [13] "F3D145_S211_L001_R1_001.fastq" "F3D145_S211_L001_R2_001.fastq"
## [15] "F3D146_S212_L001_R1_001.fastq" "F3D146_S212_L001_R2_001.fastq"
## [17] "F3D147_S213_L001_R1_001.fastq" "F3D147_S213_L001_R2_001.fastq"
## [19] "F3D148_S214_L001_R1_001.fastq" "F3D148_S214_L001_R2_001.fastq"
## [21] "F3D149_S215_L001_R1_001.fastq" "F3D149_S215_L001_R2_001.fastq"
## [23] "F3D150_S216_L001_R1_001.fastq" "F3D150_S216_L001_R2_001.fastq"
## [25] "F3D2_S190_L001_R1_001.fastq" "F3D2_S190_L001_R2_001.fastq"
## [27] "F3D3_S191_L001_R1_001.fastq" "F3D3_S191_L001_R2_001.fastq"
## [29] "F3D5_S193_L001_R1_001.fastq" "F3D5_S193_L001_R2_001.fastq"
## [31] "F3D6_S194_L001_R1_001.fastq" "F3D6_S194_L001_R2_001.fastq"
## [33] "F3D7_S195_L001_R1_001.fastq" "F3D7_S195_L001_R2_001.fastq"
## [35] "F3D8_S196_L001_R1_001.fastq" "F3D8_S196_L001_R2_001.fastq"
## [37] "F3D9_S197_L001_R1_001.fastq" "F3D9_S197_L001_R2_001.fastq"
## [39] "filtered" "HMP_MOCK.v35.fasta"
## [41] "Mock_S280_L001_R1_001.fastq" "Mock_S280_L001_R2_001.fastq"
## [43] "mouse.dpw.metadata" "mouse.time.design"
## [45] "stability.batch" "stability.files"#j'identifie les fichiers dans les deux lectures : forward (R1) et reverse (R2), comme les fragments sont sequencés sur les deux brins en illumina :
fnFs <- sort(list.files(path, pattern="_R1_001.fastq", full.names = TRUE))
fnRs <- sort(list.files(path, pattern="_R2_001.fastq", full.names = TRUE))
#j'extrait les noms des échantillons à partir des fichiers fastq, en supposant que le nom de l'échantillon est avant le premier "_" :
sample.names <- sapply(strsplit(basename(fnFs), "_"), `[`, 1)
# afin de permettre à dada2 d’associer correctement chaque paire de fichiers forward et reverse à l’échantillon correspondant et de suivre ces échantillons tout au long du pipeline# j'affiche un graphique de la qualité des deux premièrs reads pour visualiser la précision des bases et identifier les régions à tronquer afin d’améliorer l’analyse,
#ici les graphes sont representatifs de la qualité des deux premieres lectures forward (R1) :
plotQualityProfile(fnFs[1:2])
# de la meme maniere, j'affiche la qualité des deux premiers fichiers reverse R2 cette fois ci, afin avoir un aperçu rapide de la qualité avant de traiter tous les fichiers :
plotQualityProfile(fnRs[1:2])
# je prepare les chemins des fichiers filtrés avant d'executer le filtrage réel pour les deux lectures forward et reverse des sequences, puis je crée un chemin complet vers un sous dossier "filtered" où seront mit tout les fichiers filtrés et compréssés au format.gz :
filtFs <- file.path(path, "filtered", paste0(sample.names, "_F_filt.fastq.gz"))
filtRs <- file.path(path, "filtered", paste0(sample.names, "_R_filt.fastq.gz"))
#j'associe chaque fichier filtré et compréssé à son nom d’échantillon pour que dada2 puisse les traiter correctement plus tard :
names(filtFs) <- sample.names
names(filtRs) <- sample.names#je filtre et tronque les séquences afin de conserver uniquement les lectures de haute qualité de la maniere suivante :
#tronquer les lectures forward à 240 bases et les reverse à 160 bases pour éliminer les régions finales de faible qualité.
#supprimer toutes les séquences contenant des bases ambiguës (N).
#limiter le nombre d’erreurs attendues par lecture à 2 pour forward et 2 pour reverse.
#tronquer au premier point où la qualité est très faible (Q<2).
#supprimer les séquences issues du contrôle PhiX, un fragment d’ADN ajouté par Illumina pour calibrer le séquençage et vérifier le fonctionnement de la machine.
#les fichiers filtrés seront enregistrés compressés au format .gz.
# pas de parallélisation
out <- filterAndTrim(fnFs, filtFs, fnRs, filtRs, truncLen=c(240,160),
maxN=0, maxEE=c(2,2), truncQ=2, rm.phix=TRUE,
compress=TRUE, multithread=FALSE)
#j'affiche les premières lignes du tableau de filtrage, pour montrer le nombre de lectures conservées pour chaque échantillon après le tronquage et le filtrage :
head(out)## reads.in reads.out
## F3D0_S188_L001_R1_001.fastq 7793 7113
## F3D1_S189_L001_R1_001.fastq 5869 5299
## F3D141_S207_L001_R1_001.fastq 5958 5463
## F3D142_S208_L001_R1_001.fastq 3183 2914
## F3D143_S209_L001_R1_001.fastq 3178 2941
## F3D144_S210_L001_R1_001.fastq 4827 4312#j'estime les taux d'erreurs sur les sequences forward, dada2 ici va apprendre les models d'erreurs typiques introduite par le sequençage illumina, ce qui permettra par la suite de corriger ces sequences :
errF <- learnErrors(filtFs, multithread=TRUE)## 33514080 total bases in 139642 reads from 20 samples will be used for learning the error rates.# et le paramètre multithread=TRUE permet d’accélérer le calcul en utilisant plusieurs cœurs du processeur.#j'estime les taux d'erreurs mais pour les sequences reverse :
errR <- learnErrors(filtRs, multithread=TRUE)## 22342720 total bases in 139642 reads from 20 samples will be used for learning the error rates.#je visualise des erreurs apprises sur les lectures forward, en affichant un graphique montrant la probabilité d'erreur à chaque score de qualité :
plotErrors(errF, nominalQ=TRUE)## Warning: Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
#j'applique de l'algorithme dada pour corriger les erreurs sur les lectures forward, dada2 utilise le modèle d'erreurs appris (errF) :
dadaFs <- dada(filtFs, err=errF, multithread=TRUE)## Sample 1 - 7113 reads in 1979 unique sequences.
## Sample 2 - 5299 reads in 1639 unique sequences.
## Sample 3 - 5463 reads in 1477 unique sequences.
## Sample 4 - 2914 reads in 904 unique sequences.
## Sample 5 - 2941 reads in 939 unique sequences.
## Sample 6 - 4312 reads in 1267 unique sequences.
## Sample 7 - 6741 reads in 1756 unique sequences.
## Sample 8 - 4560 reads in 1438 unique sequences.
## Sample 9 - 15637 reads in 3590 unique sequences.
## Sample 10 - 11413 reads in 2762 unique sequences.
## Sample 11 - 12017 reads in 3021 unique sequences.
## Sample 12 - 5032 reads in 1566 unique sequences.
## Sample 13 - 18075 reads in 3707 unique sequences.
## Sample 14 - 6250 reads in 1479 unique sequences.
## Sample 15 - 4052 reads in 1195 unique sequences.
## Sample 16 - 7369 reads in 1832 unique sequences.
## Sample 17 - 4765 reads in 1183 unique sequences.
## Sample 18 - 4871 reads in 1382 unique sequences.
## Sample 19 - 6504 reads in 1709 unique sequences.
## Sample 20 - 4314 reads in 897 unique sequences.#j'applique l'algorithme dada pour corriger les erreurs sur les lectures reverse, dada2 utilise le modèle d'erreurs appris (errR) :
dadaRs <- dada(filtRs, err=errR, multithread=TRUE)## Sample 1 - 7113 reads in 1660 unique sequences.
## Sample 2 - 5299 reads in 1349 unique sequences.
## Sample 3 - 5463 reads in 1335 unique sequences.
## Sample 4 - 2914 reads in 853 unique sequences.
## Sample 5 - 2941 reads in 880 unique sequences.
## Sample 6 - 4312 reads in 1286 unique sequences.
## Sample 7 - 6741 reads in 1803 unique sequences.
## Sample 8 - 4560 reads in 1265 unique sequences.
## Sample 9 - 15637 reads in 3414 unique sequences.
## Sample 10 - 11413 reads in 2522 unique sequences.
## Sample 11 - 12017 reads in 2771 unique sequences.
## Sample 12 - 5032 reads in 1415 unique sequences.
## Sample 13 - 18075 reads in 3290 unique sequences.
## Sample 14 - 6250 reads in 1390 unique sequences.
## Sample 15 - 4052 reads in 1134 unique sequences.
## Sample 16 - 7369 reads in 1635 unique sequences.
## Sample 17 - 4765 reads in 1084 unique sequences.
## Sample 18 - 4871 reads in 1161 unique sequences.
## Sample 19 - 6504 reads in 1502 unique sequences.
## Sample 20 - 4314 reads in 732 unique sequences.# j'ffiche un résumé qui me permet de visualiser le nombre de séquences uniques (ASVs) identifiées dans ce premier échantillon (ça permets d'avoir le détail des séquences corrigées pour le premier échantillon forward) :
dadaFs[[1]]## dada-class: object describing DADA2 denoising results
## 128 sequence variants were inferred from 1979 input unique sequences.
## Key parameters: OMEGA_A = 1e-40, OMEGA_C = 1e-40, BAND_SIZE = 16#je fusionne les lectures forward et reverse, et combine les lectures correspondantes des deux directions pour reconstruire la séquence complète :
mergers <- mergePairs(dadaFs, filtFs, dadaRs, filtRs, verbose=TRUE)## 6540 paired-reads (in 107 unique pairings) successfully merged out of 6891 (in 197 pairings) input.
## 5028 paired-reads (in 101 unique pairings) successfully merged out of 5190 (in 157 pairings) input.
## 4986 paired-reads (in 81 unique pairings) successfully merged out of 5267 (in 166 pairings) input.
## 2595 paired-reads (in 52 unique pairings) successfully merged out of 2754 (in 108 pairings) input.
## 2553 paired-reads (in 60 unique pairings) successfully merged out of 2785 (in 119 pairings) input.
## 3646 paired-reads (in 55 unique pairings) successfully merged out of 4109 (in 157 pairings) input.
## 6079 paired-reads (in 81 unique pairings) successfully merged out of 6514 (in 198 pairings) input.
## 3968 paired-reads (in 91 unique pairings) successfully merged out of 4388 (in 187 pairings) input.
## 14233 paired-reads (in 143 unique pairings) successfully merged out of 15355 (in 352 pairings) input.
## 10528 paired-reads (in 120 unique pairings) successfully merged out of 11165 (in 278 pairings) input.
## 11154 paired-reads (in 137 unique pairings) successfully merged out of 11797 (in 298 pairings) input.
## 4349 paired-reads (in 85 unique pairings) successfully merged out of 4802 (in 179 pairings) input.
## 17431 paired-reads (in 153 unique pairings) successfully merged out of 17812 (in 272 pairings) input.
## 5850 paired-reads (in 81 unique pairings) successfully merged out of 6095 (in 159 pairings) input.
## 3716 paired-reads (in 86 unique pairings) successfully merged out of 3894 (in 147 pairings) input.
## 6865 paired-reads (in 99 unique pairings) successfully merged out of 7191 (in 187 pairings) input.
## 4426 paired-reads (in 67 unique pairings) successfully merged out of 4603 (in 127 pairings) input.
## 4576 paired-reads (in 101 unique pairings) successfully merged out of 4739 (in 174 pairings) input.
## 6092 paired-reads (in 109 unique pairings) successfully merged out of 6315 (in 173 pairings) input.
## 4269 paired-reads (in 20 unique pairings) successfully merged out of 4281 (in 28 pairings) input.## sequence
## 1 TACGGAGGATGCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGCAGGCGGAAGATCAAGTCAGCGGTAAAATTGAGAGGCTCAACCTCTTCGAGCCGTTGAAACTGGTTTTCTTGAGTGAGCGAGAAGTATGCGGAATGCGTGGTGTAGCGGTGAAATGCATAGATATCACGCAGAACTCCGATTGCGAAGGCAGCATACCGGCGCTCAACTGACGCTCATGCACGAAAGTGTGGGTATCGAACAGG
## 2 TACGGAGGATGCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGTAGGCGGCCTGCCAAGTCAGCGGTAAAATTGCGGGGCTCAACCCCGTACAGCCGTTGAAACTGCCGGGCTCGAGTGGGCGAGAAGTATGCGGAATGCGTGGTGTAGCGGTGAAATGCATAGATATCACGCAGAACCCCGATTGCGAAGGCAGCATACCGGCGCCCTACTGACGCTGAGGCACGAAAGTGCGGGGATCAAACAGG
## 3 TACGGAGGATGCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGTAGGCGGGCTGTTAAGTCAGCGGTCAAATGTCGGGGCTCAACCCCGGCCTGCCGTTGAAACTGGCGGCCTCGAGTGGGCGAGAAGTATGCGGAATGCGTGGTGTAGCGGTGAAATGCATAGATATCACGCAGAACTCCGATTGCGAAGGCAGCATACCGGCGCCCGACTGACGCTGAGGCACGAAAGCGTGGGTATCGAACAGG
## 4 TACGGAGGATGCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGTAGGCGGGCTTTTAAGTCAGCGGTAAAAATTCGGGGCTCAACCCCGTCCGGCCGTTGAAACTGGGGGCCTTGAGTGGGCGAGAAGAAGGCGGAATGCGTGGTGTAGCGGTGAAATGCATAGATATCACGCAGAACCCCGATTGCGAAGGCAGCCTTCCGGCGCCCTACTGACGCTGAGGCACGAAAGTGCGGGGATCGAACAGG
## 5 TACGGAGGATGCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGCAGGCGGACTCTCAAGTCAGCGGTCAAATCGCGGGGCTCAACCCCGTTCCGCCGTTGAAACTGGGAGCCTTGAGTGCGCGAGAAGTAGGCGGAATGCGTGGTGTAGCGGTGAAATGCATAGATATCACGCAGAACTCCGATTGCGAAGGCAGCCTACCGGCGCGCAACTGACGCTCATGCACGAAAGCGTGGGTATCGAACAGG
## 6 TACGGAGGATGCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGTAGGCGGGATGCCAAGTCAGCGGTAAAAAAGCGGTGCTCAACGCCGTCGAGCCGTTGAAACTGGCGTTCTTGAGTGGGCGAGAAGTATGCGGAATGCGTGGTGTAGCGGTGAAATGCATAGATATCACGCAGAACTCCGATTGCGAAGGCAGCATACCGGCGCCCTACTGACGCTGAGGCACGAAAGCGTGGGTATCGAACAGG
## abundance forward reverse nmatch nmismatch nindel prefer accept
## 1 579 1 1 148 0 0 1 TRUE
## 2 470 2 2 148 0 0 2 TRUE
## 3 449 3 4 148 0 0 1 TRUE
## 4 430 4 3 148 0 0 2 TRUE
## 5 345 5 6 148 0 0 1 TRUE
## 6 282 6 5 148 0 0 2 TRUE# je crée la table de séquences (ASV table), en prennant les séquences fusionnées et en construisant une matrice où chaque colonne correspond à une séquence unique (ASV : Amplicon Sequence Variant), chaque ligne correspond à un échantillon, et les valeurs indiquent combien de fois chaque séquence a été observée dans chaque échantillon :
seqtab <- makeSequenceTable(mergers)
# j'affichage les dimensions de la table, c'est à dire le nombre d'échantillons (lignes) et le nombre de séquences uniques (colonnes) :
dim(seqtab)## [1] 20 293# je vérifie la distribution des longueurs des séquences, pour voir si toutes les séquences ont la longueur attendue et s'il y a des anomalies avant l'étape de suppression des chimères :
table(nchar(getSequences(seqtab)))##
## 251 252 253 254 255
## 1 88 196 6 2# je supprime les séquences chimériques (generée lors de la PCR), et j'affiche les informations pendant l'exécution :
seqtab.nochim <- removeBimeraDenovo(seqtab, method="consensus", multithread=TRUE, verbose=TRUE)## Identified 61 bimeras out of 293 input sequences.# je vérifie le nombre de séquences conservées après suppression des chimères (les ASV) :
dim(seqtab.nochim)## [1] 20 232# je calcule la proportion des reads consérvés (donc aprés suppression des chimeres) par rapport au nombre total de reads (avant la suppression des chimeres):
sum(seqtab.nochim)/sum(seqtab)## [1] 0.9640374# je crée une fonction (getN) qui compte le nombre total de lectures uniques dans chaque étape du pipeline :
getN <- function(x) sum(getUniques(x))
# je construis un tableau (track) qui regroupe pour chaque échantillon le nombre de reads conservés après chaque phase : en entrée, après filtrage, après débruitage des reads forward et reverse, après fusion des paires, et après suppression des chimères :
track <- cbind(out, sapply(dadaFs, getN), sapply(dadaRs, getN), sapply(mergers, getN), rowSums(seqtab.nochim))
# je renomme les colonnes du tableau pour indiquer clairement les différentes étapes du traitement :
colnames(track) <- c("input", "filtered", "denoisedF", "denoisedR", "merged", "nonchim")
# j’associe chaque ligne du tableau au nom de l’échantillon correspondant :
rownames(track) <- sample.names
# j’affiche les premières lignes du tableau pour avoir un apperçu de la perte ou la conservation de lectures tout au long du processus :
head(track)## input filtered denoisedF denoisedR merged nonchim
## F3D0 7793 7113 6976 6979 6540 6528
## F3D1 5869 5299 5227 5239 5028 5017
## F3D141 5958 5463 5331 5357 4986 4863
## F3D142 3183 2914 2799 2830 2595 2521
## F3D143 3178 2941 2822 2868 2553 2519
## F3D144 4827 4312 4151 4228 3646 3507# j'assigne une taxonomie à chaque ASV de la table des ASVs, en utilisant SILVA comme base de reference, dada2 va comparer les sequences stockés dans le fichier data à cette base pour determiner à chaque taxon (jusqu'au niveau du genre) l'ASV appartient :
taxa <- assignTaxonomy(seqtab.nochim, "/home/rstudio/dada2/data/silva_nr99_v138.2_toGenus_trainset.fa.gz", multithread=TRUE)# je complète la taxonomie de chaque séquence pour obtenir une classification plus fine en ajoutant une assignation au niveau espèce pour chaque ASV déjà classée au niveau genre dans l’objet taxa, grace à la base de référence SILVA spécifique aux espèces, et le résultat final contient donc pour chaque ASV l’information du genre et, si possible, de l’espèce :
taxa <- addSpecies(taxa, "/home/rstudio/dada2/data/silva_v138.2_assignSpecies.fa.gz")# je crée une copie de l’objet taxa nommée taxa.print pour l’affichage uniquement
taxa.print <- taxa
# je supprime les noms de lignes (qui correspondent aux séquences ASV) afin de rendre le tableau plus lisible :
rownames(taxa.print) <- NULL
#j’affiche les premières lignes pour visualiser la taxonomie (de règne jusqu'à espèce) de quelques ASVs sans montrer les séquences elles-mêmes :
head(taxa.print)## Kingdom Phylum Class Order Family
## [1,] "Bacteria" "Bacteroidota" "Bacteroidia" "Bacteroidales" "Muribaculaceae"
## [2,] "Bacteria" "Bacteroidota" "Bacteroidia" "Bacteroidales" "Muribaculaceae"
## [3,] "Bacteria" "Bacteroidota" "Bacteroidia" "Bacteroidales" "Muribaculaceae"
## [4,] "Bacteria" "Bacteroidota" "Bacteroidia" "Bacteroidales" "Muribaculaceae"
## [5,] "Bacteria" "Bacteroidota" "Bacteroidia" "Bacteroidales" "Bacteroidaceae"
## [6,] "Bacteria" "Bacteroidota" "Bacteroidia" "Bacteroidales" "Muribaculaceae"
## Genus Species
## [1,] NA NA
## [2,] NA NA
## [3,] NA NA
## [4,] NA NA
## [5,] "Bacteroides" NA
## [6,] NA NA# je sélectionne les ASVs correspondant à l’échantillon controle "Mock" dans la matrice des ASVs :
unqs.mock <- seqtab.nochim["Mock",]
# je garde que les ASVs présents dans l'échantillon Mock (valeur >0) et je les trie par abondance décroissante :
unqs.mock <- sort(unqs.mock[unqs.mock>0], decreasing=TRUE)
# j'affiche dans la console le nombre de séquences uniques ASVs détectées par dada2 dans l'échantillon Mock :
cat("DADA2 inferred", length(unqs.mock), "sample sequences present in the Mock community.\n")## DADA2 inferred 20 sample sequences present in the Mock community.# je lis les séquences de référence de la Mock depuis le fichier FASTA (HMP_MOCK.v35.fasta) afin de connaître exactement quelles séquences devraient être présentes, et je les stocke dans mock.ref :
mock.ref <- getSequences(file.path(path, "HMP_MOCK.v35.fasta"))
# pour chaque ASV, je vérifie si elle correspond exactement à une séquence du fichier de référence, puis je compte le nombre total de correspondances exactes :
match.ref <- sum(sapply(names(unqs.mock), function(x) any(grepl(x, mock.ref))))
# j’affiche dans la console combien des ASVs détectées dans l’échantillon Mock sont des correspondances exactes avec les séquences attendues du fichier de référence, (cela permet d’évaluer la precision du pipeline dada2) :
cat("Of those,", sum(match.ref), "were exact matches to the expected reference sequences.\n")## Of those, 20 were exact matches to the expected reference sequences.