DADA2 pipeline - Español
Alexis Carteron & Simon Morvan
2021-09-13 16:21:59
Introducción
DADA2 es un pipeline bioinformático creado por Callahan et al., 2016. Consiste es una serie de pasos que filtran las secuencias crudas obtenidas con la secuenciación de Illumina. El último paso es obtener la taxonomía de las secuencias que han sido filtradas para estudiar la comunidad microbiana.
DADA2 tiene dos características principales que lo distinguen de otros pipelines comúnmente utilizados. Por un lado, procede a la modelización del error de secuenciación que se supone permite distinguir las secuencias mutantes de las erróneas. Por otro lado, a diferencia de otros pipelines como QIIME o Mothur, DADA2 no agrupa el 97% de las secuencias similares en unidades taxonómicas operativas (Operational Taxonomy Units, OTU). Sus variantes de secuencia de ampliación (Amplicon Sequence Variants, ASVs) no se agrupan si las secuencias no son 100% idénticas. Véase la figura anterior.
Construido originalmente para secuencias del gen marcador 16S (Bacterias), lo utilizaremos con secuencias del gen marcador ITS (Hongos) procedentes de la secuenciación Illumina MiSEQ 2x300 bp paired-end. Para acelerar la ejecución de cada paso, submuestreamos aleatoriamente un conjunto de datos para conservar sólo 1000 secuencias por muestra. Por último, Redde Caesari quae sunt Caesaris : este tutorial se ha inspirado en gran medida en el original tutorial DADA2.
En general, antes de iniciar este pipeline, debemos tomar algunas precauciones:
- Las muestras deben ser demultiplexadas: divididas en archivos fastq individuales por muestra.
- Si las secuencias están emparejadas, las secuencias directa e inversa deben estar en archivos fastq distintos pero deben contener lecturas en orden coincidente.
- Si las secuencias están emparejadas, las secuencias directa e inversa deben estar en archivos fastq distintos pero deben contener lecturas en orden coincidente.
- Los nucleótidos que no forman parte del amplicón (cebadores y adaptadores) deben ser eliminados. También se pueden eliminar en el paso de filtrado.
- La mayoría de las funciones tienen una opción de multihilo que permite un tiempo de cálculo más rápido al acceder a múltiples procesadores. Sólo hay que especificar multithread = TRUE para activarla. Atención, esta opción no funciona bajo Windows.
Esta figura tomada de Hugerth y Andersson, 2017 ilustra la diferencia teórica entre OTUs y ASV. Cada color representa un clado. Las estrellas amarillas indican mutaciones, las rojas indican errores de amplificación o secuenciación. El tamaño del espacio entre las secuencias indica su agrupación.
(A) OTUs agrupadas con un 100% de identidad
La más mínima variación en las secuencias provoca la creación de un nuevo grupo. Las secuencias mutantes y erróneas se tratan de la misma manera.
(B) OTUs agrupadas con un 97% de identidad
Una agrupación más amplia permite descartar las secuencias erróneas, aunque las secuencias mutantes también se agruparán en el grupo de consenso.
(C) ASVs
En teoría, el aprendizaje de la tasa de error permite agrupar las secuencias erróneas con las secuencias de consenso. Sin embargo, las secuencias mutantes se consideran por derecho propio.
¡Comencemos!
Primero cargaremos la biblioteca DADA2. Debería tener la última versión: packageVersion('dada2').
A continuación, crearemos una variable path que indique la ruta de acceso a los objetos que necesitaremos.
library(dada2); packageVersion("dada2")## Loading required package: Rcpp## [1] '1.18.0'path <- "data/ITS_sub/"Comprobemos qué hay al final del camino…
list.files(path)## [1] "cutadapt" "filtered_pairedend" "filtN"
## [4] "S1_R1.fastq.gz" "S1_R2.fastq.gz" "S10_R1.fastq.gz"
## [7] "S10_R2.fastq.gz" "S11_R1.fastq.gz" "S11_R2.fastq.gz"
## [10] "S12_R1.fastq.gz" "S12_R2.fastq.gz" "S13_R1.fastq.gz"
## [13] "S13_R2.fastq.gz" "S14_R1.fastq.gz" "S14_R2.fastq.gz"
## [16] "S15_R1.fastq.gz" "S15_R2.fastq.gz" "S16_R1.fastq.gz"
## [19] "S16_R2.fastq.gz" "S17_R1.fastq.gz" "S17_R2.fastq.gz"
## [22] "S18_R1.fastq.gz" "S18_R2.fastq.gz" "S2_R1.fastq.gz"
## [25] "S2_R2.fastq.gz" "S20_R1.fastq.gz" "S20_R2.fastq.gz"
## [28] "S21_R1.fastq.gz" "S21_R2.fastq.gz" "S22_R1.fastq.gz"
## [31] "S22_R2.fastq.gz" "S23_R1.fastq.gz" "S23_R2.fastq.gz"
## [34] "S24_R1.fastq.gz" "S24_R2.fastq.gz" "S25_R1.fastq.gz"
## [37] "S25_R2.fastq.gz" "S26_R1.fastq.gz" "S26_R2.fastq.gz"
## [40] "S27_R1.fastq.gz" "S27_R2.fastq.gz" "S28_R1.fastq.gz"
## [43] "S28_R2.fastq.gz" "S29_R1.fastq.gz" "S29_R2.fastq.gz"
## [46] "S3_R1.fastq.gz" "S3_R2.fastq.gz" "S30_R1.fastq.gz"
## [49] "S30_R2.fastq.gz" "S31_R1.fastq.gz" "S31_R2.fastq.gz"
## [52] "S33_R1.fastq.gz" "S33_R2.fastq.gz" "S34_R1.fastq.gz"
## [55] "S34_R2.fastq.gz" "S35_R1.fastq.gz" "S35_R2.fastq.gz"
## [58] "S36_R1.fastq.gz" "S36_R2.fastq.gz" "S37_R1.fastq.gz"
## [61] "S37_R2.fastq.gz" "S39_R1.fastq.gz" "S39_R2.fastq.gz"
## [64] "S4_R1.fastq.gz" "S4_R2.fastq.gz" "S41_R1.fastq.gz"
## [67] "S41_R2.fastq.gz" "S42_R1.fastq.gz" "S42_R2.fastq.gz"
## [70] "S43_R1.fastq.gz" "S43_R2.fastq.gz" "S44_R1.fastq.gz"
## [73] "S44_R2.fastq.gz" "S45_R1.fastq.gz" "S45_R2.fastq.gz"
## [76] "S46_R1.fastq.gz" "S46_R2.fastq.gz" "S47_R1.fastq.gz"
## [79] "S47_R2.fastq.gz" "S48_R1.fastq.gz" "S48_R2.fastq.gz"
## [82] "S49_R1.fastq.gz" "S49_R2.fastq.gz" "S5_R1.fastq.gz"
## [85] "S5_R2.fastq.gz" "S50_R1.fastq.gz" "S50_R2.fastq.gz"
## [88] "S51_R1.fastq.gz" "S51_R2.fastq.gz" "S52_R1.fastq.gz"
## [91] "S52_R2.fastq.gz" "S53_R1.fastq.gz" "S53_R2.fastq.gz"
## [94] "S54_R1.fastq.gz" "S54_R2.fastq.gz" "S55_R1.fastq.gz"
## [97] "S55_R2.fastq.gz" "S56_R1.fastq.gz" "S56_R2.fastq.gz"
## [100] "S58_R1.fastq.gz" "S58_R2.fastq.gz" "S59_R1.fastq.gz"
## [103] "S59_R2.fastq.gz" "S6_R1.fastq.gz" "S6_R2.fastq.gz"
## [106] "S60_R1.fastq.gz" "S60_R2.fastq.gz" "S61_R1.fastq.gz"
## [109] "S61_R2.fastq.gz" "S62_R1.fastq.gz" "S62_R2.fastq.gz"
## [112] "S63_R1.fastq.gz" "S63_R2.fastq.gz" "S64_R1.fastq.gz"
## [115] "S64_R2.fastq.gz" "S65_R1.fastq.gz" "S65_R2.fastq.gz"
## [118] "S66_R1.fastq.gz" "S66_R2.fastq.gz" "S67_R1.fastq.gz"
## [121] "S67_R2.fastq.gz" "S68_R1.fastq.gz" "S68_R2.fastq.gz"
## [124] "S69_R1.fastq.gz" "S69_R2.fastq.gz" "S7_R1.fastq.gz"
## [127] "S7_R2.fastq.gz" "S72_R1.fastq.gz" "S72_R2.fastq.gz"
## [130] "S73_R1.fastq.gz" "S73_R2.fastq.gz" "S74_R1.fastq.gz"
## [133] "S74_R2.fastq.gz" "S75_R1.fastq.gz" "S75_R2.fastq.gz"
## [136] "S76_R1.fastq.gz" "S76_R2.fastq.gz" "S77_R1.fastq.gz"
## [139] "S77_R2.fastq.gz" "S78_R1.fastq.gz" "S78_R2.fastq.gz"
## [142] "S79_R1.fastq.gz" "S79_R2.fastq.gz" "S8_R1.fastq.gz"
## [145] "S8_R2.fastq.gz" "S80_R1.fastq.gz" "S80_R2.fastq.gz"
## [148] "S81_R1.fastq.gz" "S81_R2.fastq.gz" "S9_R1.fastq.gz"
## [151] "S9_R2.fastq.gz"Debería ver los nombres de los archivos fastq.
Ahora leeremos los nombres de los archivos fastq, y manipularemos su cadena variable para obtener una lista de fastqs sentido y antisentido. La función ordenar dará el mismo orden entre fastq sentido y antisentido.
fnFs <- sort(list.files(path, pattern="_R1.fastq"))
fnRs <- sort(list.files(path, pattern="_R2.fastq"))Como los pares de archivos fastq sentido y antisentido pertenecen a la misma muestra, crearemos una variable que extraiga el nombre de esa muestra. En este caso, suponemos que los nombres de los archivos fastq tienen un formato: SAMPLENAME_XXX.fastq.
sample.names <- sapply(strsplit(fnFs, "_"), `[`, 1)
sample.names ## [1] "S1" "S10" "S11" "S12" "S13" "S14" "S15" "S16" "S17" "S18" "S2" "S20"
## [13] "S21" "S22" "S23" "S24" "S25" "S26" "S27" "S28" "S29" "S3" "S30" "S31"
## [25] "S33" "S34" "S35" "S36" "S37" "S39" "S4" "S41" "S42" "S43" "S44" "S45"
## [37] "S46" "S47" "S48" "S49" "S5" "S50" "S51" "S52" "S53" "S54" "S55" "S56"
## [49] "S58" "S59" "S6" "S60" "S61" "S62" "S63" "S64" "S65" "S66" "S67" "S68"
## [61] "S69" "S7" "S72" "S73" "S74" "S75" "S76" "S77" "S78" "S79" "S8" "S80"
## [73] "S81" "S9"Ahora especificaremos la ruta de acceso a los objetos fnFs y fnRS.
fnFs <- file.path(path, fnFs)
fnRs <- file.path(path, fnRs)Perfil de calidad
Este primer paso nos permite visualizar la calidad de las secuencias gracias a la puntuación Q asociada a cada nucleótido.
plotQualityProfile(fnFs[1]) # 1er echantillon sens
plotQualityProfile(fnRs[1]) # 1er echantillon anti-sens
En estas figuras, la mediana aparece en verde, y los cuartiles en naranja punteado. Aquí hemos elegido visualizar la primera dirección de la muestra (fnFs[1]) y el antisentido (fnRs[1]), pero es posible visualizar varios gráficos al mismo tiempo (fnFs[x:y]) o agregarlos de la siguiente manera.
plotQualityProfile(fnFs, aggregate = TRUE)
plotQualityProfile(fnRs, aggregate = TRUE)
El análisis de estos gráficos nos permite elegir los parámetros de filtrado y recorte para el siguiente paso. De hecho, el índice de puntuación Q nos da información sobre la precisión de la secuenciación (véase la tabla siguiente).
| Resultado Q | Precisión |
|---|---|
| 10 | 90 % |
| 20 | 99 % |
| 30 | 99.9 % |
| 40 | 99.99 % |
Otra herramienta más completa para evaluar la calidad de la secuencia es FastQC.
Filtrado y truncamiento
Primero crearemos una carpeta (filtered_pairedend) y objetos (filtFs y filtRs) para detener las secuencias filtradas.
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"))Proceda con la función filterAndTrim, su salida será almacenada en el objeto out.
out <- filterAndTrim(fnFs, filtFs, fnRs, filtRs,
truncQ=6,
truncLen = c(280,280),
trimLeft=c(18,20),
maxEE=c(2,2))## Creating output directory: data/ITS_sub//filtered_pairedend #multithread=TRUE)En primer lugar, la función necesita las secuencias no filtradas (fnFs y FnRs), así como los nombres de los objetos de las secuencias filtradas (filtFs y filtRs). A continuación, se pueden modificar varios parámetros según se desee:
- truncQ: define un índice de puntuación Q mínimo. En la primera instancia de una puntuación de calidad menor o igual a truncQ, la secuencia se trunca.
- truncLen: define a qué longitud se truncarán las secuencias. Se eliminan las secuencias más cortas que la longitud.
- trimLeft: define la longitud que se corta en el lado 5’ de las secuencias. Esto permite eliminar los cebadores si no se ha hecho antes.
(Cartillas:
ITS3_KYO2: GATGAAGAACGYAGYRAA = 18bp
ITS4: TCCTCCGCTTATTGATATGC = 20bp) - maxEE: define el número máximo de “errores esperados” permitidos en una lectura. Este filtro se basa en el índice de puntuación Q. Cuanto más alto sea el número, menos estricto será.
También se pueden modificar otros parámetros, que se pueden encontrar en la página de ayuda de la función: ?filterAndTrim.
Visualización del filtrado
pourc <- cbind((out[,2]/out[,1])*100) # Porcentaje de secuencias filtradas / no filtradas
pourc_disc <- cbind(out, pourc) # Combina el objeto out y pourc
pourc_disc ## reads.in reads.out
## S1_R1.fastq.gz 1000 595 59.5
## 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
## S2_R1.fastq.gz 1000 482 48.2
## 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
## S3_R1.fastq.gz 1000 556 55.6
## 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
## S4_R1.fastq.gz 1000 425 42.5
## 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
## S5_R1.fastq.gz 1000 387 38.7
## 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
## S6_R1.fastq.gz 1000 619 61.9
## 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
## S7_R1.fastq.gz 1000 521 52.1
## 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
## S8_R1.fastq.gz 1000 423 42.3
## S80_R1.fastq.gz 1000 457 45.7
## S81_R1.fastq.gz 1000 594 59.4
## S9_R1.fastq.gz 1000 492 49.2(mean(out[,2])/mean(out[,1]))*100 # Porcentaje medio## [1] 50.98243¡Casi la mitad de las secuencias no pasaron nuestros parámetros de filtrado!
CHALLENGE
Trace el perfil de calidad de la primera muestra filtrada y compárelo con su perfil de calidad sin filtrar. Utilice la función plotQualityProfile.
Existen otras herramientas de filtrado como el : Trimmomatic.
Aprender las tasas de error
Este paso consiste en estimar la tasa de error de secuenciación. El modelo de error se calcula alternando la estimación de las tasas de error y la inferencia de la composición de la muestra hasta que convergen en una solución coherente.
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=TRUEEl número mínimo de secuencias a utilizar para el entrenamiento de la tasa de error puede especificarse con el parámetro nreads.
Visualización de la tasa de error
plotErrors(errF, nominalQ=TRUE)
plotErrors(errR, nominalQ=TRUE)
Se muestran las tasas de error de cada transición (A->C, A->G,…). Cada punto es una tasa de error observada para cada puntuación de calidad consensuada. La línea negra muestra el error después de la convergencia. La línea roja muestra el error bajo la definición del valor Q nominal (para la tecnología Illumina).
De-replicación
En este paso, todas las secuencias idénticas se agruparán en secuencias únicas a las que se les asignan abundancias. Esto reducirá el tiempo de cálculo posterior al eliminar las comparaciones redundantes. Las secuencias desreplicadas reciben el nombre de las muestras de las que proceden.
derepFs <- derepFastq(filtFs)
names(derepFs) <- sample.names
derepRs <- derepFastq(filtRs)
names(derepRs) <- sample.namesLa ventaja de DADA2 es que mantiene un resumen de la información de calidad asociada a cada secuencia única. El perfil de calidad de consenso de una secuencia única es la media de las calidades posicionales de las lecturas desreplicadas. Estos perfiles de calidad informan al modelo de error del siguiente paso de inferencia de la muestra, lo que aumenta en gran medida la precisión de DADA2.
Unique sequences
Inferencia de muestras
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
## 70 sequence variants were inferred from 364 input unique sequences.
## Key parameters: OMEGA_A = 1e-40, OMEGA_C = 1e-40, BAND_SIZE = 16dadaRs[[1]]## dada-class: object describing DADA2 denoising results
## 73 sequence variants were inferred from 340 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")Fusión de secuencias emparejadas
¡FUSIÓN!
El objetivo de la secuenciación por pares es fusionar las dos cadenas para aumentar nuestra confianza en su fiabilidad. La fusión también permite obtener amplicones más largos.
La función mergePairs necesita que se le proporcionen los objetos calculados en los dos pasos anteriores (derep y dada). Entonces, los parámetros que podemos modificar libremente son:
- minOverlap: define el tamaño mínimo del solapamiento de las hebras sentido y antisentido para que su fusión sea exitosa. Las secuencias que no se fusionan se eliminan.
- maxMismatch: define el número máximo de incompatibilidades de nucleótidos en el solapamiento.
También se pueden cambiar otros parámetros, están accesibles en la página de ayuda de la función: ?mergePairs. Por ejemplo, si returnRejects = TRUE, los pares rechazados por falta de coincidencia en la región de solapamiento se mantienen en la salida.
mergers <- mergePairs(dadaFs, derepFs, dadaRs, derepRs,
minOverlap = 12,
maxMismatch = 0)¡Veamos los resultados!
head(mergers[[1]])## sequence
## 1 ATGCGATACGTAATGTGAATTGCAGAATTCAGTGAATCATCGAATCTTTGAACGCACCTTGCGCCCCTTGGTATTCCGAGGGGCACACCCGTTTGAGTGTCGTGAATACTCTCAACCTTCTTGGTTTCTTTGACCACGAAGGCTTGGACTTTGGAGGTTTTTCTTGCTGGCCTCTTTAGAAGCCAGCTCCTCCTAAATGAATGGGTGGGGTCCGCTTTGCTGATCCTCGACGTGATAAGCATCTCTTCTACGTCTCAGTGTCAGCTCGGAACCCCGCTTTCCAACCGTCTTTGGACAAAGACAATGTTCGAGTTGCGACTCGACCTTACAAACCTTGACCTCAAATCGGGTGAGACTACCCGCTGAACTTAA
## 2 ATGCGATAAGTAGTGTGAATTGCAGAATTCAGTGAATCATCGAATCTTTGAACGCACATTGCGCCCCTTGGTATCCCGAGGGGCATGCCTGTTCGAGCGTCATTTCACCACTCAAGCCTGGCTTGGTGTTGGGCGACGTCCCCTTTTGGGGACGCGTCTCGAAACGCTCGGCGGCGTGGCACCGGCTTTAAGCGTAGCAGAATCTTTCGCTTTGAAAGTCGGGGCCCCGTCTGCCGGAAGACCTACTCGCAAGGTTGACCTCGGATCAGGCAGGGATACCCGCTGAACTTAA
## 3 ATGCGATAAGTAATGTGAATTGCAGAATTCAGTGAATCATCGAATCTTTGAACGCACATTGCGCCCCTTGGTATTCCGAGGGGCATGCCTGTTCGAGCGTCATTATAACCACTCAAGCCCCGGCTTGGTCTTGGGGTTCGCGGTCCGCGGCCCTTAAACTCAGTGGCGGTGCCGTCTGGCTCTAAGCGCAGTAATTCTCTCGCTATAGTGTCTAGGTGGTTGCTTGCCATAATCCCCCAATTTTTTACGGTTGACCTCGGATCAGGTAGGGATACCCGCTGAACTTAA
## 4 ATGCGATACGTAATGTGAATTGCAGAATTCAGTGAATCATCGAATCTTTGAACGCACATTGCACTCCTTGGTATTCCGAGGAGTATGCCTGTTTCAGTATCATGAGCACTCTCACACCTAACCTTTGGGTTTATGGCGTGGAATTGGAATGCGCCGACTGTCATGGTTGGCCCTTCTAAAATGTAGTTCTTGGCTGTCACCTAATACAGCAGTTTGGCCTAATAGTTTTGGCATTCATTGTCAAATCTTTGGCTAACATTTGCTCCAGGAGTCAGTCTTGATAATACAGAAAACTCATTCAAATTTTGATCTGAAATCAGGTAGGGCTACCCGCTGAACTTAA
## 5 ATGCGATACGTAATGTGAATTGCAGAATTCCGTGAATCATTGAATCTTTGAACGCATCTTGCGCCTCTTGGTATTCCGAGGGGCATGCCTGTTTGAGTGTCATTAGAACTATCAAAAAAATAGATGATTTCAATCGTTAATTTTTTTGGAATTGGAGGTGGTGCTGGTCTTTTTCCATTAATGGCCCAAGCTCCTCCGAAATGCATTAGCGAATGCAGTGCACTTTTTCTCCTTGCTTTTTCTGGGCATTGATAGTTTACTCTCATGCCCTAAGCTGGTAGGGAGGAAGTCACAGAATGCTTCCCGCTCCTGAATGTAATACAAAACTTGACGATCAAACCCCTCAAATCAGGCAGGACTACCCGCTGAACTTAA
## 6 ATGCGATAAGTAATGCGAATTGCAGAATTCAGTGAGTCATCGAATCTTTGAACGCATATTGCGCCCTTTGGTATTCCGAAGGGCATGCCTGTTCGAGCGTCATGATCAACCATCAAGCCTGGCTTGTCGTTGGACCCTGTTGTCTCTGGGCGACAGGTCCGAAAGATAATGACGGTGTCATGGCAACCCCGAATGCAACGAGCTTTTTTATAGGCACGCATTTAGTGGTTGGCAAGGCCCCCTCGTGCGTTATTATTTTCTTACGGTTGACCTCGGATCAGGTAGGAATACCCGCTGAACTTAA
## abundance forward reverse nmatch nmismatch nindel prefer accept
## 1 172 14 11 150 0 0 2 TRUE
## 2 43 7 8 230 0 0 1 TRUE
## 3 37 21 17 234 0 0 2 TRUE
## 4 34 2 1 179 0 0 2 TRUE
## 5 32 3 2 147 0 0 1 TRUE
## 6 21 11 9 218 0 0 1 TRUEmax(mergers[[1]]$nmatch) # Taille du plus grand chevauchement (overlap)## [1] 257min(mergers[[1]]$nmatch) # Taille du plus petit chevauchement ## [1] 113Tabla de ASVs
Finalmente tenemos nuestros ASVs que pararemos en el objeto seqtab.
seqtab <- makeSequenceTable(mergers)
dim(seqtab)## [1] 74 551seqtab[,1]## S1 S10 S11 S12 S13 S14 S15 S16 S17 S18 S2 S20 S21 S22 S23 S24 S25 S26 S27 S28
## 0 0 0 0 0 0 0 46 207 127 0 0 1 0 0 1 0 0 532 0
## S29 S3 S30 S31 S33 S34 S35 S36 S37 S39 S4 S41 S42 S43 S44 S45 S46 S47 S48 S49
## 0 0 81 7 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0
## S5 S50 S51 S52 S53 S54 S55 S56 S58 S59 S6 S60 S61 S62 S63 S64 S65 S66 S67 S68
## 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0
## S69 S7 S72 S73 S74 S75 S76 S77 S78 S79 S8 S80 S81 S9
## 0 0 0 0 0 0 0 6 0 0 0 0 0 0Obtenemos 551 ASVs a partir de las 74000 secuencias en bruto que teníamos al principio.
seqtab[,1] da el número de veces que se encuentra la primera secuencia ASV en cada muestra.
Inspeccionar la longitud de nuestros ASV
hist(nchar(getSequences(seqtab)),xlab="Size", ylab="Frequency", main = "ASVs length", xlim=c(250,450), ylim=c(0,250)) 
Eliminación de las quimeras
Este paso tiene como objetivo eliminar todas las secuencias no biológicas, las muestras del array ASV se agrupan (method=“pooled”) para su identificación. Se pueden utilizar otros métodos, como el método de consenso, en el que cada muestra se comprueba individualmente para identificar los bímeros.
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")¡Veamos los resultados!
round((sum(seqtab.nochim)/sum(seqtab)*100),2) # Porcentaje del número total de secuencias no químicas / número total de secuencias## [1] 99.4hist(nchar(getSequences(seqtab.nochim)),xlab="Size", ylab="Frequency", main = "Non-chimeric ASVs length", xlim=c(250,450), ylim=c(0,250)) # Longitud de las secuencias no químicas
Ahora podemos transformar las apariciones de ASV en presencia/ausencia. Esto nos permite cuantificar el número de ASV por muestra.
seqtab.nochim.bin <- ifelse(seqtab.nochim>0,1,0) Cuadro de seguimiento
La siguiente tabla muestra cuántas secuencias se han eliminado en cada etapa. Si el número de secuencias desciende demasiado, puede indicar un problema. Por ejemplo, si pocas secuencias superan el paso de eliminación de bímeros, esto puede indicar que todavía hay trozos de cebadores con nucleótidos ambiguos.
getN <- function(x) sum(getUniques(x))
track <- data.frame(Input=as.numeric(out[,1]),
Filtered=as.numeric(out[,2]),
"Filt//In"=as.numeric(round(((out[,2]/out[,1])*100),2)),
Merge = as.numeric(sapply(mergers, getN)),
"Mer//In"=as.numeric(round(((sapply(mergers, getN)/out[,1])*100),2)),
Nonchim = as.numeric(rowSums(seqtab.nochim)),
"Nonchim//In"=as.numeric(round(((rowSums(seqtab.nochim)/out[,1])*100),2)),
ASV = as.numeric(rowSums(seqtab.nochim.bin)))
rownames(track) <- sample.names
head(track)## Input Filtered Filt..In Merge Mer..In Nonchim Nonchim..In ASV
## S1 1000 595 59.5 564 56.4 564 56.4 65
## 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 128 12.8 107 10.7 27Una imagen vale más que mil palabras.
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 
Asignación taxonómica
Por fin vemos el final de la tubería con este importante paso de asignación taxonómica. Gracias a la implementación del método de clasificación bayesiano ingenuo, la función assignTaxonomy toma como entrada el conjunto de secuencias a clasificar así como un conjunto de secuencias de referencia con una taxonomía conocida. Las asignaciones taxonómicas se dan con una confianza bootstrap mínima especificada con el parámetro minBoot. La base de datos de referencia (UNITE) es accesible en este enlace https://unite.ut.ee/repository.php. También hay otras bases de datos disponibles.
taxotab <- assignTaxonomy(seqtab.nochim,
refFasta = "reference_database/sh_general_release_dynamic_01.12.2017.fasta",
multithread=TRUE)## UNITE fungal taxonomic reference detected.save(taxotab, file = "data/taxotab.rdata")
#load("data/taxotab.rdata")¡Veamos los resultados!
# view(taxotab)
write.table(taxotab[1:5,], row.names = 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])) ## [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__microspora" "s__camphoratus" "s__acicola"
## [40] "s__verrucosa" "s__fuckelii" "s__var._bulbopilosa"
## [43] "s__miniata" "s__maius" "s__splendens"
## [46] "s__echinulatum" "s__chlorophana" "s__australis"
## [49] "s__citrina" "s__lignicola" "s__lubrica"
## [52] "s__phyllophila" "s__cantharellus" "s__trabinellum"
## [55] "s__pygmaeum" "s__isabellina" "s__changbaiensis"
## [58] "s__fragilis" "s__carneum" "s__fuscella"
## [61] "s__bicolor" "s__nigrella" "s__bulbillosa"
## [64] "s__difforme" "s__spirale" "s__finlandica"
## [67] "s__myriocarpa" "s__reginae" "s__fellea"
## [70] "s__vrijmoediae" "s__lacmus" "s__spurius"
## [73] "s__arachnoidea" "s__dioscoreae" "s__laetior"
## [76] "s__coccinea" "s__spadicea" "s__pyriforme"
## [79] "s__macrocystis" "s__risigallina" "s__pura"
## [82] "s__cinereus" "s__foliicola" "s__rebaudengoi"
## [85] "s__fusiformis" "s__metachroides" "s__diversispora"
## [88] "s__fumosa" "s__hymenocystis" "s__sublilacina"
## [91] "s__rufum" "s__atrovirens" "s__mors-panacis"
## [94] "s__acerinum" "s__skinneri" "s__glacialis"
## [97] "s__cygneicollum" "s__pullulans" "s__crocea"
## [100] "s__globulifera" "s__silvestris" "s__variabilis"
## [103] "s__griseoviride" "s__nitrata" "s__sphaeroides"
## [106] "s__renispora" "s__honrubiae" "s__sindonia"
## [109] "s__grovesii" "s__piceae" "s__stuposa"
## [112] "s__moravica" "s__anomalovelatus" "s__miyagiana"
## [115] "s__pilosella" "s__flavidum" "s__glutinosum"
## [118] "s__fallax" "s__falcata" "s__rhododendri"
## [121] "s__fusispora" "s__scaurus" "s__lauri"
## [124] "s__serotinus" "s__alliacea" "s__cylichnium"
## [127] "s__alpina" "s__miyabei" "s__rosea"
## [130] "s__albicastaneus" "s__sepiacea" "s__fumosibrunneus"
## [133] "s__spinosum" "s__hyalocuspica" "s__boeremae"
## [136] "s__epicalamia" "s__terrestris" "s__rimosissimus"
## [139] "s__bombacina" "s__physodes" "s__verzuoliana"
## [142] "s__acerina" "s__entomopaga" "s__xylopini"
## [145] "s__umbrosum" "s__physaroides"CHALLENGE
¿Puedes encontrar el número de familias diferentes?
Nota Bene
Dada2 no elimina las secuencias singleton. Sin embargo, se supone que la línea de producción no puede inferir variantes de secuencias biológicas que sólo se apoyan en una única secuencia, ya que se supone que las secuencias únicas son demasiado difíciles de diferenciar de los errores.
DADA2 mide la diversidad de forma consistente a través de los diferentes parámetros de filtrado y tasa de error. Los métodos OTU no lo hacen.
Las ASVs sin asignación de especie no fueron asignadas a la misma especie en más del 50% de las asignaciones kmer bootstrap (ver Wang et al., 2007 para más información sobre el método de clasificación bayesiano ingenuo).