DADA2 pipeline - Français
Alexis Carteron & Simon Morvan
2019-05-10 13:52:03
Introduction
DADA2 est un pipeline bio-informatique (Callahan et al., 2016). Il consiste en une série d’étapes permettant de filtrer les séquences brutes obtenues grâce au séquençage Illumina. La dernière étape vise à obtenir la taxonomie des séquences ayant été filtrées en vue d’étudier la communauté microbienne.
DADA2 a deux particularités qui le distingue des autres pipeline courament utilisés. D’une part, il va procéder à modélisation de l’erreur dûe au séquençage ce qui est censé permettre des distinguer les séquences mutantes, des séquences érronées. D’autre part, contrairement à d’autres pipelines comme QIIME ou Mothur, DADA2 ne regroupe pas les séquences similaires à 97% en Unités Taxonomiques Opérationelles (OTUs). Ses Variants de Séquences d’Amplicons (ASVs) ne subissent pas de regroupement si les séqences ne sont pas identiques à 100%. Voir figure ci-dessous.
A l’origine construit pour les séquences de gène marqueur 16S (Bacteria), nous l’utiliserons avec des séquences de gène marqueur ITS (Fungi) provenant d’un séquençage en paire Illumina MiSEQ 2x300 paires de bases. Afin d’accélérer l’exécution de chaques étapes, nous avons sous-échantilloner aléatoirement un jeu de données afin d’avoir 1000 séquences par échantillons.
Enfin, Redde Caesari quae sunt Caesaris : ce tutoriel s’est largement inspiré du propre tutoriel de DADA2.
De manière générale avant de débuter ce pipeline, il faut prendre quelques précautions:
- Les échantillons doivent être démultiplexés chaque échantillon doit avoir son propre fichier fastq.
- En cas de séquençage en paire, les séquences sens et anti-sens doivent être dans deux fichier fastq distincts et être dans le même ordre dans les deux fichiers.
- Les nucléotides qui ne font pas partie de l’amplicon (amorces, adaptateurs, bar-code) doivent avoir été retirées. Dans le cas contraire, ils devront être à l’étape de filtrage.
- La plupart des fonctions présentées ont une option de multithreading qui permet d’accélerer les temps de calcul en accédant à plusieurs prorcesseurs. Il suffit d’indiquer multithread = TRUE pour l’activer. Attention, cette option ne marche pas sous Windows.
Cette figure extraite de Hugerth et Andersson, 2017 illustre la différence théorique entre OTUs et ASVs. Chaque couleur représente une clade. Les étoiles jaunes indiquent des mutations, les étoiles rouges indiquent des erreurs d’amplification ou de séquençage. La taille de l’espace entre les séquences indique leur regroupement.
(A) OTUs regroupés à 100 % d’identité.
La moindre variation de séquences provoque la création d’un nouveau groupe. Les séquences mutantes et les séquences erronées sont traitées de la même manière.
(B) OTUs regroupés à 97 % d’identité.
Un regroupement plus large permet de ne plus considérer les séquences erronées, cependant les séquences mutantes seront également regroupées dans le groupe consensus.
(C) ASVs
L’apprentissage du taux d’erreur permet théoriquement de regrouper les séquences erronées avec les séquences consensus. En revanche, les séquences mutantes sont considérées à part entière.
Commençons !
Nous allons tout d’abord chargé la librairie DADA2. Vous devriez avoir la denière version: packageVersion('dada2').
Puis nous allons créer une variable (path) indiquant le chemin qui permettra d’accéder aux objets dont nous allons avoir besoin.
library(dada2); packageVersion("dada2")## Loading required package: Rcpp## [1] '1.11.1'path <- "data/ITS_sub/"Vérifions ce qu’il y a au bout du chemin…
list.files(path)## [1] "cutadapt" "filtered_pairedend" "filtN"
## [4] "S10_R1.fastq.gz" "S10_R2.fastq.gz" "S11_R1.fastq.gz"
## [7] "S11_R2.fastq.gz" "S12_R1.fastq.gz" "S12_R2.fastq.gz"
## [10] "S13_R1.fastq.gz" "S13_R2.fastq.gz" "S14_R1.fastq.gz"
## [13] "S14_R2.fastq.gz" "S15_R1.fastq.gz" "S15_R2.fastq.gz"
## [16] "S16_R1.fastq.gz" "S16_R2.fastq.gz" "S17_R1.fastq.gz"
## [19] "S17_R2.fastq.gz" "S18_R1.fastq.gz" "S18_R2.fastq.gz"
## [22] "S1_R1.fastq.gz" "S1_R2.fastq.gz" "S20_R1.fastq.gz"
## [25] "S20_R2.fastq.gz" "S21_R1.fastq.gz" "S21_R2.fastq.gz"
## [28] "S22_R1.fastq.gz" "S22_R2.fastq.gz" "S23_R1.fastq.gz"
## [31] "S23_R2.fastq.gz" "S24_R1.fastq.gz" "S24_R2.fastq.gz"
## [34] "S25_R1.fastq.gz" "S25_R2.fastq.gz" "S26_R1.fastq.gz"
## [37] "S26_R2.fastq.gz" "S27_R1.fastq.gz" "S27_R2.fastq.gz"
## [40] "S28_R1.fastq.gz" "S28_R2.fastq.gz" "S29_R1.fastq.gz"
## [43] "S29_R2.fastq.gz" "S2_R1.fastq.gz" "S2_R2.fastq.gz"
## [46] "S30_R1.fastq.gz" "S30_R2.fastq.gz" "S31_R1.fastq.gz"
## [49] "S31_R2.fastq.gz" "S33_R1.fastq.gz" "S33_R2.fastq.gz"
## [52] "S34_R1.fastq.gz" "S34_R2.fastq.gz" "S35_R1.fastq.gz"
## [55] "S35_R2.fastq.gz" "S36_R1.fastq.gz" "S36_R2.fastq.gz"
## [58] "S37_R1.fastq.gz" "S37_R2.fastq.gz" "S39_R1.fastq.gz"
## [61] "S39_R2.fastq.gz" "S3_R1.fastq.gz" "S3_R2.fastq.gz"
## [64] "S41_R1.fastq.gz" "S41_R2.fastq.gz" "S42_R1.fastq.gz"
## [67] "S42_R2.fastq.gz" "S43_R1.fastq.gz" "S43_R2.fastq.gz"
## [70] "S44_R1.fastq.gz" "S44_R2.fastq.gz" "S45_R1.fastq.gz"
## [73] "S45_R2.fastq.gz" "S46_R1.fastq.gz" "S46_R2.fastq.gz"
## [76] "S47_R1.fastq.gz" "S47_R2.fastq.gz" "S48_R1.fastq.gz"
## [79] "S48_R2.fastq.gz" "S49_R1.fastq.gz" "S49_R2.fastq.gz"
## [82] "S4_R1.fastq.gz" "S4_R2.fastq.gz" "S50_R1.fastq.gz"
## [85] "S50_R2.fastq.gz" "S51_R1.fastq.gz" "S51_R2.fastq.gz"
## [88] "S52_R1.fastq.gz" "S52_R2.fastq.gz" "S53_R1.fastq.gz"
## [91] "S53_R2.fastq.gz" "S54_R1.fastq.gz" "S54_R2.fastq.gz"
## [94] "S55_R1.fastq.gz" "S55_R2.fastq.gz" "S56_R1.fastq.gz"
## [97] "S56_R2.fastq.gz" "S58_R1.fastq.gz" "S58_R2.fastq.gz"
## [100] "S59_R1.fastq.gz" "S59_R2.fastq.gz" "S5_R1.fastq.gz"
## [103] "S5_R2.fastq.gz" "S60_R1.fastq.gz" "S60_R2.fastq.gz"
## [106] "S61_R1.fastq.gz" "S61_R2.fastq.gz" "S62_R1.fastq.gz"
## [109] "S62_R2.fastq.gz" "S63_R1.fastq.gz" "S63_R2.fastq.gz"
## [112] "S64_R1.fastq.gz" "S64_R2.fastq.gz" "S65_R1.fastq.gz"
## [115] "S65_R2.fastq.gz" "S66_R1.fastq.gz" "S66_R2.fastq.gz"
## [118] "S67_R1.fastq.gz" "S67_R2.fastq.gz" "S68_R1.fastq.gz"
## [121] "S68_R2.fastq.gz" "S69_R1.fastq.gz" "S69_R2.fastq.gz"
## [124] "S6_R1.fastq.gz" "S6_R2.fastq.gz" "S72_R1.fastq.gz"
## [127] "S72_R2.fastq.gz" "S73_R1.fastq.gz" "S73_R2.fastq.gz"
## [130] "S74_R1.fastq.gz" "S74_R2.fastq.gz" "S75_R1.fastq.gz"
## [133] "S75_R2.fastq.gz" "S76_R1.fastq.gz" "S76_R2.fastq.gz"
## [136] "S77_R1.fastq.gz" "S77_R2.fastq.gz" "S78_R1.fastq.gz"
## [139] "S78_R2.fastq.gz" "S79_R1.fastq.gz" "S79_R2.fastq.gz"
## [142] "S7_R1.fastq.gz" "S7_R2.fastq.gz" "S80_R1.fastq.gz"
## [145] "S80_R2.fastq.gz" "S81_R1.fastq.gz" "S81_R2.fastq.gz"
## [148] "S8_R1.fastq.gz" "S8_R2.fastq.gz" "S9_R1.fastq.gz"
## [151] "S9_R2.fastq.gz"Vous devriez voir les noms des fichiers fastq.
Nous allons maintenant lire les noms des fichiers fastq, et manipuler leur chaine de charactères variables pour obtenir une liste des fastq sens et antisens. La fonction sort permet d’obtenir le même ordre entre les fastq sens et antisens.
fnFs <- sort(list.files(path, pattern="_R1.fastq"))
fnRs <- sort(list.files(path, pattern="_R2.fastq"))Etant donné que les paires de fichiers fastq sens et antisens appartiennent au même échantillon, nous allons créer une variable qui extrait le nom de cet échantillon. Dans ce cas, nous partons du principe que les noms des fichiers fastq ont un format: SAMPLENAME_XXX.fastq.
sample.names <- sapply(strsplit(fnFs, "_"), `[`, 1)
sample.names ## [1] "S10" "S11" "S12" "S13" "S14" "S15" "S16" "S17" "S18" "S1" "S20"
## [12] "S21" "S22" "S23" "S24" "S25" "S26" "S27" "S28" "S29" "S2" "S30"
## [23] "S31" "S33" "S34" "S35" "S36" "S37" "S39" "S3" "S41" "S42" "S43"
## [34] "S44" "S45" "S46" "S47" "S48" "S49" "S4" "S50" "S51" "S52" "S53"
## [45] "S54" "S55" "S56" "S58" "S59" "S5" "S60" "S61" "S62" "S63" "S64"
## [56] "S65" "S66" "S67" "S68" "S69" "S6" "S72" "S73" "S74" "S75" "S76"
## [67] "S77" "S78" "S79" "S7" "S80" "S81" "S8" "S9"Nous allons maintenant préciser le chemin d’accès aux objets fnFs et fnRS.
fnFs <- file.path(path, fnFs)
fnRs <- file.path(path, fnRs)Profil de qualité
Cette première étape permet de visualiser la qualité des séquences grâce au Q score associé à chaque nucléotide.
plotQualityProfile(fnFs[1]) # 1er echantillon sens
plotQualityProfile(fnRs[1]) # 1er echantillon anti-sens
Dans ces figures, la médianne est en vert, et les quartiles en orange pointillé. Ici, nous avons choisi de visualiser le premier échantillon sens (fnFs[1]) et antisens (fnRs[1]), mais il est possible de visualiser plusieurs graphiques en même temps (fnFs[x:y]) ou les aggréger comme ce qui suit.
plotQualityProfile(fnFs, aggregate = TRUE)
plotQualityProfile(fnRs, aggregate = TRUE)
L’analyse de ces graphiques nous permet de choisir les paramètres de filtrage et de rognage de la prochaine étape. En effet, l’indice de Q score nous renseigne sur la précision du séquençage (voir tableau ci-dessous).
| Q score | Precision |
|---|---|
| 10 | 90 % |
| 20 | 99 % |
| 30 | 99.9 % |
| 40 | 99.99 % |
Un autre outil plus complet pour évaluer la qualité des séquence: FastQC.
Filtrage et tronquage
Tout d’abord nous allons créer un dossier (filtered_pairedend) et des objets (filtFs et filtRs) pour stoquer les séquences filtrées.
filt_path <- file.path(path, "filtered_pairedend")
filtFs <- file.path(filt_path, paste0(sample.names, "_F_filt.fastq.gz"))
filtRs <- file.path(filt_path, paste0(sample.names, "_R_filt.fastq.gz"))Procédons avec la fonction filterAndTrim, sa sortie sera stocké dans l’objet out.
out <- filterAndTrim(fnFs, filtFs, fnRs, filtRs,
truncQ=6,
truncLen = c(280,280),
trimLeft=c(18,20),
maxEE=c(2,2))
#multithread=TRUE)En premier lieu, la fonction a besoin des séquences non filtrées (fnFs et FnRs) ainsi que les noms des objets des séquences filtrées (filtFs et filtRs). Plusieurs paramètres peuvent ensuite être modifiés à notre guise:
- truncQ : définit un l’indice Q score minimale. A la première instance d’un score de qualité inférieur ou égal à truncQ, la séquence est tronquée.
- truncLen : définit à quelle longueur les séquences vont être tronquées. Les séquences plus courtes que la longueur sont éliminées.
- trimLeft : définit la longueur que l’on coupe du côté 5’ des séquences. Celà permet d’enlever les amorces si ça n’a pas été fait préalablement.
(Amorces:
ITS3_KYO2: GATGAAGAACGYAGYRAA = 18bp
ITS4: TCCTCCGCTTATTGATATGC = 20b) - maxEE : définit le nombre maximum d’ “erreurs attendues” autorisées dans une lecture. Ce filtre se base sur l’indice Q score. Plus on augmente le chiffre, moins on est strict.
D’autres paramètres peuvent également être modifiés, ils sont accessibles à la page d’aide de la fonction: ?filterAndTrim.
Visualisation du filtrage / Filtering visualization
pourc <- cbind((out[,2]/out[,1])*100) # Pourcentage des séquences filtrées / séquences non-filtrées
pourc_disc <- cbind(out, pourc) # Combine l'objet out et pourc
pourc_disc ## reads.in reads.out
## S10_R1.fastq.gz 1000 366 36.6
## S11_R1.fastq.gz 1000 432 43.2
## S12_R1.fastq.gz 1000 546 54.6
## S13_R1.fastq.gz 1000 417 41.7
## S14_R1.fastq.gz 1000 229 22.9
## S15_R1.fastq.gz 1000 512 51.2
## S16_R1.fastq.gz 1000 481 48.1
## S17_R1.fastq.gz 1000 634 63.4
## S18_R1.fastq.gz 1000 607 60.7
## S1_R1.fastq.gz 1000 595 59.5
## S20_R1.fastq.gz 1000 614 61.4
## S21_R1.fastq.gz 1000 598 59.8
## S22_R1.fastq.gz 1000 539 53.9
## S23_R1.fastq.gz 1000 561 56.1
## S24_R1.fastq.gz 1000 599 59.9
## S25_R1.fastq.gz 1000 598 59.8
## S26_R1.fastq.gz 1000 567 56.7
## S27_R1.fastq.gz 1000 650 65.0
## S28_R1.fastq.gz 1000 607 60.7
## S29_R1.fastq.gz 1000 526 52.6
## S2_R1.fastq.gz 1000 482 48.2
## S30_R1.fastq.gz 1000 551 55.1
## S31_R1.fastq.gz 1000 571 57.1
## S33_R1.fastq.gz 1000 482 48.2
## S34_R1.fastq.gz 1000 644 64.4
## S35_R1.fastq.gz 1000 558 55.8
## S36_R1.fastq.gz 1000 453 45.3
## S37_R1.fastq.gz 1000 446 44.6
## S39_R1.fastq.gz 1000 521 52.1
## S3_R1.fastq.gz 1000 556 55.6
## S41_R1.fastq.gz 1000 517 51.7
## S42_R1.fastq.gz 1000 447 44.7
## S43_R1.fastq.gz 1000 500 50.0
## S44_R1.fastq.gz 1000 369 36.9
## S45_R1.fastq.gz 1000 359 35.9
## S46_R1.fastq.gz 1000 540 54.0
## S47_R1.fastq.gz 1000 501 50.1
## S48_R1.fastq.gz 1000 391 39.1
## S49_R1.fastq.gz 1000 516 51.6
## S4_R1.fastq.gz 1000 425 42.5
## S50_R1.fastq.gz 1000 422 42.2
## S51_R1.fastq.gz 1000 567 56.7
## S52_R1.fastq.gz 1000 521 52.1
## S53_R1.fastq.gz 1000 599 59.9
## S54_R1.fastq.gz 1000 472 47.2
## S55_R1.fastq.gz 1000 571 57.1
## S56_R1.fastq.gz 1000 416 41.6
## S58_R1.fastq.gz 1000 525 52.5
## S59_R1.fastq.gz 1000 559 55.9
## S5_R1.fastq.gz 1000 387 38.7
## S60_R1.fastq.gz 1000 593 59.3
## S61_R1.fastq.gz 1000 537 53.7
## S62_R1.fastq.gz 1000 470 47.0
## S63_R1.fastq.gz 1000 644 64.4
## S64_R1.fastq.gz 1000 595 59.5
## S65_R1.fastq.gz 1000 561 56.1
## S66_R1.fastq.gz 1000 551 55.1
## S67_R1.fastq.gz 1000 544 54.4
## S68_R1.fastq.gz 1000 532 53.2
## S69_R1.fastq.gz 1000 372 37.2
## S6_R1.fastq.gz 1000 619 61.9
## S72_R1.fastq.gz 1000 558 55.8
## S73_R1.fastq.gz 1000 581 58.1
## S74_R1.fastq.gz 1000 431 43.1
## S75_R1.fastq.gz 1000 356 35.6
## S76_R1.fastq.gz 1000 403 40.3
## S77_R1.fastq.gz 1000 387 38.7
## S78_R1.fastq.gz 1000 400 40.0
## S79_R1.fastq.gz 1000 560 56.0
## S7_R1.fastq.gz 1000 521 52.1
## S80_R1.fastq.gz 1000 457 45.7
## S81_R1.fastq.gz 1000 594 59.4
## S8_R1.fastq.gz 1000 423 42.3
## S9_R1.fastq.gz 1000 492 49.2(mean(out[,2])/mean(out[,1]))*100 # Pourcentage moyen## [1] 50.98243Près de la moitié des séquences n’ont pas passé nos paramètres de filtrage!
CHALLENGE
Tracer le profil de qualité du premier échantillon une fois filtré et le comparer à son profil de qualité non-filtré. Utiliser la fonction plotQualityProfile.
D’autres outils de filtrage existent tel que : Trimmomatic.
Apprentissage des taux d’erreur
Cet étape consiste à estimer le taux d’erreur de séquençage. Son but est de permettre de différencier les séquences mutantes et les séquences érronées.Le modèle d’erreur est calculé en alternant l’estimation des taux d’erreur et l’inférence de la composition de l’échantillon jusqu’à ce qu’ils convergent vers une solution cohérente.
errF <- learnErrors(filtFs)## 9884474 total bases in 37727 reads from 74 samples will be used for learning the error rates.errR <- learnErrors(filtRs)## 9809020 total bases in 37727 reads from 74 samples will be used for learning the error rates.#multithread=TRUELe nombre minimum de séquences à utiliser pour l’apprentissage des taux d’erreur peut être précisé avec le paramètre nreads.
Visualisation du taux d’erreur
plotErrors(errF, nominalQ=TRUE)
plotErrors(errR, nominalQ=TRUE)
Les taux d’erreur pour chaque transition (A->C, A->G,…) sont affichés. Chaque point est un taux d’erreur observé pour chaque score de qualité consensuel. La ligne noire montre l’erreur après convergence. La ligne rouge montre l’erreur sous la définition nominale de la valeur Q (pour la technologie Illumina).
Déreplication
Dans cette étape, toutes les séquences identiques vont être regroupées en séquences uniques auxquelles sont attribuées des abondances. Cela va diminuer les temps de calcul subséquants en éliminant des comparaisons redondantes. Les séquences dérépliquées prennent le nom des échantillons d’où elles proviennent.
derepFs <- derepFastq(filtFs)
names(derepFs) <- sample.names
derepRs <- derepFastq(filtRs)
names(derepRs) <- sample.namesL’avantage de DADA2 réside dans la conservation d’un résumé des informations de qualité associées à chaque séquence unique. Le profil de qualité consensuel d’une séquence unique est la moyenne des qualités de position des lectures dérépliquées. Ces profils de qualité informent le modèle d’erreur de l’étape suivante d’inférence d’échantillon, ce qui augmente considérablement la précision de DADA2.
Unique sequences
Inférence des échantillons
dadaFs <- dada(derepFs,
err = errF,
#multithread=TRUE,
pool=TRUE)## 74 samples were pooled: 37727 reads in 20364 unique sequences.dadaRs <- dada(derepRs,
err=errR,
#multithread=TRUE,
pool=TRUE)## 74 samples were pooled: 37727 reads in 18908 unique sequences.dadaFs[[1]]## dada-class: object describing DADA2 denoising results
## 50 sequence variants were inferred from 238 input unique sequences.
## Key parameters: OMEGA_A = 1e-40, OMEGA_C = 1e-40, BAND_SIZE = 16dadaRs[[1]]## dada-class: object describing DADA2 denoising results
## 57 sequence variants were inferred from 273 input unique sequences.
## Key parameters: OMEGA_A = 1e-40, OMEGA_C = 1e-40, BAND_SIZE = 16#save(dadaRs, file="data/dadaRs.rdata")
#save(dadaFs, file="data/dadaFs.rdata")Fusion des séquences pairées
FUSION !
Tout l’intérêt du séquençage en paire réside dans le but de fusionner les deux brins afin d’accroitre notre confiance en leur fiabilité. La fusion permet également d’obtenir des amplicons plus long.
La fonction mergePairs nécessite qu’on lui fournisse les objets calculés dans les deux étapes précédentes (derep et dada). Puis, les paramètres que nous pouvons modifiés librement sont:
- minOverlap : définit la taille minimale du chevauchement des brins sens et anti-sens pour que leur fusion soit accéptée. Les séquences qui ne fusionnent pas sont éliminées.
- maxMismatch : définit le nombre maximal d’incompatibilité nucléotidique dans le chevauchement.
D’autres paramètres peuvent également être modifiés, ils sont accessibles à la page d’aide de la fonction: ?mergePairs. Par exemple, si returnRejects = TRUE, les paires qui ont été rejetées en raison de discordances dans la région de chevauchement sont conservés dans la sortie.
mergers <- mergePairs(dadaFs, derepFs, dadaRs, derepRs,
minOverlap = 12,
maxMismatch = 0)Examinons les résultats!
head(mergers[[1]])## sequence
## 1 ATGCGATACTTGGTGTGAATTGCAGAATCCCGTGAATCATCGAGTCTTTGAACGCAAGTTGCGCCCGACGCCATTCGGCCGAGGGCACGTCTGCCTGGGTGTCACGCATCGTTGCCCCCAAACGCCTCCCGCCTCGCAAGGGGCGCGGGATCTCGTTTGGCGGCGGAAGTTGGCCTCCCGTGGGCCTGTGCTCGCGGTTAGCCTAAAAAGGAGTCCTCGGCGACGAGCGCCACGACAATCGGTGGTTGATGAGACCTCGGTCCCCGTCGTGCGTGTCTGGTCGCCACACGGTGTGACTCGTCGACCCTAACGCGTCGTACCCACGTCGCTCCCGACGCGACCCCAGGTCAGGCGGGACTACCCGCTGAGTTTAA
## 2 ATGCGATACTTGGTGTGAATTGCAGAATCCCGTGAATCATCGAGTCTTTGAACGCAAGTTGCGCCCGACGCCATTCGGCCGAGGGCACGTCTGCCTGGGTGTCACGCATCGTTGCCCCCAAACGCCTCCCGCCTCGCAAGGGGCGCGGGATCTCGTTTGGCGGCGGAAGTTGGCCTCCCGTGGGCCTGTGCTCGCGGTTAGCCTAAAAAGGAGTCCTCGGCGACGAGCGCCACGACAATCGGTGGTTGATGAGACCTCTGTCCCCGTCGTGCGTGTCTGGTCGCCACACGGTGTGACTCGTCGACCCTAACGCGTCGTACCCACGTCGCTCCCGACGCGACCCCAGGTCAGGCGGGACTACCCGCTGAGTTTAA
## 3 ATGCGATACTTGGTGTGAATTGCAGAATCCCGTGAATCATCGAGTCTTTGAACGCAAGTTGCGCCCGACGCCATTCGGCCGAGGGCACGTCTGCCTGGGTGTCACGCATCGTTGCCCCCAAACGCCCCCGCCTCGCAAGCGGCGCGGGATCTCGTTTGGTGGCGGAAGTTGGCCTCCCGTGGGCCTGTGCTCGCGGTTAGCCTAAAAAGGAGTCCTCGGCGACGAGCGCCACGACAATCGGTGGTTGATGAGACCTCGGTCCCCGTCGTGCGTGTCTGGTCGCCACACGGTGTGACTCGTCGACCCTAACGCGTCGTACCCACGTCGTTCCCAACGCGACCCCAGGTCAGGCGGGACTACCCGCTGAGTTTAA
## 4 ATGCGATACTTGGTGTGAATTGCAGAATCCCGTGAATCATCGAGTCTTTGAACGCAAGTTGCGCCCGACGCCATTCGGCCGAGGGCACGTCTGCCTGGGTGTCACGCATCGTTGCCCCCAAACGCCCCCGCCTCGCAAGGGGCGAGGGATCTCGTTTGGTGGCGGAAGTTGGCCTCCCGTGGGCCTGTGCTCGCGGTTAGCCTAAAAAGGAGTCCTCGGCGACGAGCGCCACGACAATCGGTGGTTGATGAGACCTCGGTCCCCGTCGTGCGTGTCTGGTCGCCACACGGTGTGACTCGTCGACCCTAACGCGTCGTACCCACGTCGCTCCCGACGCGACCCCAGGTCAGGCGAGACTACCCGCTGAGTTTAA
## 5 ATGCGATACTTGGTGTGAATTGCAGAATCCCGTGAATCATCGAGTCTTTGAACGCAAGTTGCGCCCGACGCCATTCGGCCGAGGGCACGTCTGCCTGGGTGTCACGCATCGTTGCCCCCAAACGCCCCCGCCTCGCAAGGGGCGCGGGATCTCGTTTGGTGGCGGAAGTTGGCCTCCCGTGGGCCTGTGCTCGCGGTTAGCCTAAAAAAGAGTCCTCGGCGACGAGCGCCACGACAATCGGTGGTTGATGAGACCTCGGTCCCCGTCGTGCGTGTCTGGTCGCCACACGGTGTGACTCGTCGACCCTAACGCGTCGTACCCACGTCGCTCCCGACGCGACCCCAGGTCAGGCGGGACTACCCGCTGAGTTTAA
## 6 ATGCGATACTTGGTGTGAATTGCAGAATCCCGTGAATCATCGAGTCTTTGAACGCAAGTTGCGCCCGACGCCATTCGGCCGAGGGCACGTCTGCCTGGGTGTCACGCATCGTTGCCCCCAAACGCCCCCGCCTCGCAAGGGGCGCGGGATCTCGTTTGGTGGCGGAAGTTGGCCTCCCGTGGGCCTGTGCTCGCGGTTAGCCTAAAAAGGAGTCCTCGGCGACGAGCGCCACGACAATCGGTGGTTGATGAGACCTCGGTCCCCGTCGTGCGTGTCTGGTCGCCACACGGTGTGACTCGTCGACCCTAACGCGTCGTACCCACGTCGCTCCCGACGCGACCCCAGGTCAGGCGGGACTACCCGCTGAGTTTAA
## abundance forward reverse nmatch nmismatch nindel prefer accept
## 1 63 1 3 148 0 0 1 TRUE
## 2 33 41 39 148 0 0 1 TRUE
## 3 16 37 16 149 0 0 1 TRUE
## 4 16 43 42 149 0 0 1 TRUE
## 5 14 6 10 149 0 0 1 TRUE
## 6 13 5 36 149 0 0 1 TRUEmax(mergers[[1]]$nmatch) # Taille du plus grand chevauchement (overlap)## [1] 258min(mergers[[1]]$nmatch) # Taille du plus petit chevauchement ## [1] 115Tableau des ASVs
Nous avons enfin nos ASVs que nous allons stoquer dans l’objet seqtab.
seqtab <- makeSequenceTable(mergers)
dim(seqtab)## [1] 74 551seqtab[,1]## S10 S11 S12 S13 S14 S15 S16 S17 S18 S1 S20 S21 S22 S23 S24 S25 S26 S27
## 0 0 0 0 0 0 46 207 127 0 0 1 0 0 1 0 0 532
## S28 S29 S2 S30 S31 S33 S34 S35 S36 S37 S39 S3 S41 S42 S43 S44 S45 S46
## 0 0 0 81 7 0 0 0 0 0 0 0 0 0 1 0 0 0
## S47 S48 S49 S4 S50 S51 S52 S53 S54 S55 S56 S58 S59 S5 S60 S61 S62 S63
## 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## S64 S65 S66 S67 S68 S69 S6 S72 S73 S74 S75 S76 S77 S78 S79 S7 S80 S81
## 0 0 0 0 0 0 1 0 0 0 0 0 6 0 0 0 0 0
## S8 S9
## 0 0Nous obtenons 551 ASVs à partir des 74000 séquences brutes que nous avions au début.
seqtab[,1] renseigne le nombre de fois qu’on retrouve la séquence du premier ASV dans chaque échantillon.
Inspectons la longueur de nos ASVs
hist(nchar(getSequences(seqtab)),xlab="Size", ylab="Frequency", main = "ASVs length", xlim=c(250,450), ylim=c(0,250)) 
Suppression des chimères
Cette étape vise à éliminer toutes les séquences non biologiques, les échantillons du tableau des ASV sont regroupés (method=“pooled”) pour l’identification. D’autres méthodes peuvent être utilisées comme la méthode consensus où chaque échantillon est vérifié individuellement pour identifier les bimères.
seqtab.nochim <- removeBimeraDenovo(seqtab,
method = "pooled",
#multithread = TRUE,
verbose = TRUE) ## Identified 5 bimeras out of 551 input sequences.#save(seqtab.nochim, file="data/seqtab.nochim.rdata")Examinons les résultats !
round((sum(seqtab.nochim)/sum(seqtab)*100),2) # Pourcentage du nombre total de séquences non-chimériques / nombre total de séquences## [1] 99.4hist(nchar(getSequences(seqtab.nochim)),xlab="Size", ylab="Frequency", main = "Non-chimeric ASVs length", xlim=c(250,450), ylim=c(0,250)) # Longueur de séquences non-chimériques
Maintenant nous pouvons transformer les occurences des ASVs en présence / absence. Ceci permet de quantifier le nombre d’ASVs par échantillons.
seqtab.nochim.bin <- ifelse(seqtab.nochim>0,1,0) Tableau de suivi
Le tableau qui suit permet de voir combien de séquences ont été éliminées à chaque étape. Une chute trop grande du nombre de séquences peut indiquer un problème. Par exemple, si peu de séquences passent l’étape de la suppresion des bimères cela peu indiquer qu’il reste des bouts d’amorces avec de nucléotides ambigües.
getN <- function(x) sum(getUniques(x))
track <- data.frame(Input=as.numeric(out[,1]), # Séquences brutes
Filtered=as.numeric(out[,2]), # Séquences filtrées
"Filt//In"=as.numeric(round(((out[,2]/out[,1])*100),2)),# % (Filtrées / Brutes)
Merge = as.numeric(sapply(mergers, getN)), # Mergés
"Mer//In"=as.numeric(round(((sapply(mergers, getN)/out[,1])*100),2)),# % (Mergés / Brutes)
Nonchim = as.numeric(rowSums(seqtab.nochim)),# Non-chimériques
"Nonchim//In"=as.numeric(round(((rowSums(seqtab.nochim)/out[,1])*100),2)),# % (Non-chimériques / Brutes)
ASV = as.numeric(rowSums(seqtab.nochim.bin))) # Nombre d'ASVs par échantillons
rownames(track) <- sample.names # Noms des lignes
head(track)## Input Filtered Filt..In Merge Mer..In Nonchim Nonchim..In ASV
## S10 1000 366 36.6 273 27.3 269 26.9 43
## S11 1000 432 43.2 342 34.2 342 34.2 58
## S12 1000 546 54.6 489 48.9 484 48.4 55
## S13 1000 417 41.7 340 34.0 338 33.8 71
## S14 1000 229 22.9 130 13.0 109 10.9 28
## S15 1000 512 51.2 460 46.0 459 45.9 68Une image vaut mille mots !
library(ggplot2)
library(reshape2)
gtrack<- track[,c(1,2,4,6)]
gtrack$ID <- rownames(gtrack)
lgtrack <- melt(gtrack, id.vars="ID")
bar_track <- ggplot(lgtrack ,aes(x=ID, y=as.numeric(value), fill=variable)) +
geom_bar(stat="identity", position = "identity") +
theme_classic() + # Thème
theme(axis.ticks.length=unit(0.3,"cm")) + # Taille de taquets
theme(axis.text.x = element_text(angle=45) , legend.title = element_blank())+ # Change l'orientation des x labels & supprime le titre de la légende
scale_x_discrete(name ="Sample ID", limits=rownames(track))+ # Change le titre de l'axe x et le pas de l'axe
scale_y_continuous(name="Abundance", breaks=seq(from = 0, to = 1000, by = 100))+ # Change le titre de l'axe y et le pas de l'axe.
ggtitle("Track")# Titre principal
bar_track 
Assignation taxonomique
Nous voyons enfin le bout du pipeline avec cette importante étape d’assignation taxonomique. Grâce à l’implémentation de la méthode de classification naïve bayésienne, la fonction assignTaxonomy prend en entrée l’ensemble de séquences à classer ainsi qu’un ensemble de séquences de référence avec une taxonomie connue. Les assignations taxonomiques sont données avec une confiance minimale de bootstrap spécifié avec le paramètre minBoot. La base de données de références (UNITE) est accessible sur ce lien https://unite.ut.ee/repository.php. D’autres bases de données sont également disponibles.
taxotab <- assignTaxonomy(seqtab.nochim,
refFasta = "reference_database/sh_general_release_dynamic_01.12.2017.fasta",
minBoot = 50, # Par défaut = 50. Bootsrap minimum (représente le niveau de confiance pour l'assignation à un rang taxonomique).
multithread=TRUE)## UNITE fungal taxonomic reference detected.save(taxotab, file = "data/taxotab.rdata")
#load("data/taxotab.rdata")Examinons les résultats !
# view(taxotab) # Tableau taxonomique au complet
write.table(taxotab[1:5,], row.names = FALSE) # Taxonomy des 5 premiers ASV sans la séquence (row.nqmes=FALSE)## "Kingdom" "Phylum" "Class" "Order" "Family" "Genus" "Species"
## "k__Fungi" "p__Basidiomycota" "c__Agaricomycetes" "o__Agaricales" "f__Tricholomataceae" "g__Mycena" NA
## "k__Fungi" "p__Mucoromycota" "c__Umbelopsidomycetes" "o__Umbelopsidales" "f__Umbelopsidaceae" "g__Umbelopsis" "s__dimorpha"
## "k__Plantae" NA NA NA NA NA NA
## "k__Fungi" "p__Ascomycota" "c__Dothideomycetes" "o__Mytilinidales" "f__Gloniaceae" "g__Cenococcum" "s__geophilum"
## "k__Fungi" "p__Basidiomycota" "c__Agaricomycetes" "o__Agaricales" "f__Hygrophoraceae" "g__Hygrocybe" "s__flavescens"unique(unname(taxotab[,7])) # Nombre d'espèces unique. ## [1] NA "s__dimorpha" "s__geophilum"
## [4] "s__flavescens" "s__algarvense" "s__terricola"
## [7] "s__punicea" "s__elongatum" "s__humilis"
## [10] "s__sylvestris" "s__opacum" "s__subvinosa"
## [13] "s__lucidum" "s__abramsii" "s__rexiana"
## [16] "s__ericae" "s__terminalis" "s__variata"
## [19] "s__conica" "s__zollingeri" "s__chlamydosporicum"
## [22] "s__asperellum" "s__album" "s__kathiae"
## [25] "s__deciduus" "s__saponaceum" "s__rufescens"
## [28] "s__populi" "s__auratus" "s__reidii"
## [31] "s__podzolica" "s__brunneoviolacea" "s__simile"
## [34] "s__heterochroma" "s__subsulphurea" "s__pseudozygospora"
## [37] "s__mutabilis" "s__microspora" "s__humicola"
## [40] "s__camphoratus" "s__acicola" "s__verrucosa"
## [43] "s__fuckelii" "s__var._bulbopilosa" "s__miniata"
## [46] "s__maius" "s__splendens" "s__paczoskii"
## [49] "s__pseudoartocreas" "s__echinulatum" "s__chlorophana"
## [52] "s__australis" "s__citrina" "s__lignicola"
## [55] "s__lubrica" "s__phyllophila" "s__cantharellus"
## [58] "s__trabinellum" "s__pygmaeum" "s__isabellina"
## [61] "s__changbaiensis" "s__fragilis" "s__carneum"
## [64] "s__fuscella" "s__bicolor" "s__nigrella"
## [67] "s__bulbillosa" "s__difforme" "s__spirale"
## [70] "s__finlandica" "s__myriocarpa" "s__reginae"
## [73] "s__fellea" "s__vrijmoediae" "s__lacmus"
## [76] "s__spurius" "s__arachnoidea" "s__dioscoreae"
## [79] "s__laetior" "s__coccinea" "s__spadicea"
## [82] "s__pyriforme" "s__macrocystis" "s__risigallina"
## [85] "s__pura" "s__cinereus" "s__foliicola"
## [88] "s__rebaudengoi" "s__fusiformis" "s__metachroides"
## [91] "s__diversispora" "s__fumosa" "s__hymenocystis"
## [94] "s__sublilacina" "s__rufum" "s__atrovirens"
## [97] "s__mors-panacis" "s__acerinum" "s__skinneri"
## [100] "s__glacialis" "s__cygneicollum" "s__crocea"
## [103] "s__pullulans" "s__delicatulum" "s__oblongisporum"
## [106] "s__globulifera" "s__silvestris" "s__variabilis"
## [109] "s__griseoviride" "s__nitrata" "s__renispora"
## [112] "s__sphaeroides" "s__eucalyptorum" "s__sindonia"
## [115] "s__grovesii" "s__piceae" "s__stuposa"
## [118] "s__moravica" "s__anomalovelatus" "s__miyagiana"
## [121] "s__pilosella" "s__flavidum" "s__glutinosum"
## [124] "s__fallax" "s__falcata" "s__rhododendri"
## [127] "s__fusispora" "s__scaurus" "s__schulzeri"
## [130] "s__lauri" "s__serotinus" "s__alliacea"
## [133] "s__cylichnium" "s__miyabei" "s__rosea"
## [136] "s__alpina" "s__albicastaneus" "s__sepiacea"
## [139] "s__fumosibrunneus" "s__spinosum" "s__hyalocuspica"
## [142] "s__boeremae" "s__epicalamia" "s__terrestris"
## [145] "s__rimosissimus" "s__bombacina" "s__physodes"
## [148] "s__verzuoliana" "s__acerina" "s__entomopaga"
## [151] "s__xylopini" "s__umbrosum" "s__physaroides"Nous obtenons 546 ASV différents dont 152 identifiés à l’espèce.
CHALLENGE
Pouvez-vous trouvez le nombre de famille différentes?
Nota Bene
Dada2 ne supprime pas les séquences singleton. Cependant, le pipeline n’est pas supposé infèrer des variants de séquence biologique qui ne sont supportés que par une seule séquence - les singletons sont supposés être trop difficiles à différencier des erreurs.
DADA2 mesure la diversité de façon cohérente à travers les différents paramètres de filtrage et de taux d’erreur. Les méthodes OTU ne le font pas.
Les ASV sans assignation d’espèces n’ont pas été assignié à la même espèce dans plus de 50% des affectations basées sur le kmer bootstrap (voir Wang et al., 2007 pour plus d’informations sur la méthode de classification bayesienne naïve).