Análisis de microbiota

Author

OMC Lab

Published

November 1, 2023

Introducción

Existen varias opciones para analizar los datos de secuenciación de microbiota. En este ejercicio utilizaremos herramientas basadas en R y Python para realizar el análisis de secuenciación de amplicones del gen ribosomal 16S. Estas herramientas generan resultados robustos, con tiempos de análisis cortas, y sin mucha necesidad de poder de cómputo.

Prerequisitos

R

DADA2

El proceso de secuenciación de amplicones introduce errores en los datos de secuenciación, lo cual complican la interpretación de los resultados. DADA2 implementa un algoritmo que modela los errores introducidos durante la secuenciación de amplicones y utiliza ese modelo de error para inferir la verdadera composición de la muestra. DADA2 reemplaza el paso tradicional de “selección de OTU” en los flujos de trabajo de secuenciación de amplicones, produciendo en su lugar tablas de mayor resolución de variantes de secuencia de amplicones (ASV).

Como se ve en el articulo de DADA2 y en evaluaciones comparativas adicionales, el método DADA2 es más sensible y específico que los métodos OTU tradicionales: DADA2 detecta variaciones biológicas reales que los métodos OTU pasan por alto y genera menos secuencias espurias. Otro artículo describe cómo reemplazar OTU con ASV mejora la precisión, exhaustividad y reproducibilidad del análisis de datos de genes marcadores.

Paquetes

BioConductor

dada2 esta en el repositorio de bioinformática BioConductor. Para instalar paquetes de este repositorio se necesita tener instalado BiocManager.

Code

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

El comando install.packages("BiocManager") es el necesario para instalarlo, pero la sección if (!requireNamespace("BiocManager", quietly = TRUE)) revisa se está instalado antes de descargar una nueva instalación.

Nota

Siempre es recomendable mantener actualizados los paquetes!

El siguiente comando actualiza la version de BioConductor.

Code

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install()

Bioconductor version 3.18 (BiocManager 1.30.22), R 4.3.1 (2023-06-16)

Old packages: 'ALL', 'Biostrings', 'evaluate', 'GenomicRanges',
  'RcppArmadillo', 'S4Vectors', 'xfun'

dada2

El siguiente comando instala la version 3.18 de dada2.

Warning

Si la consola marca algún error, lo mas probable es que esté relacionado con la version. El mensaje de salida les sugiere que versión es la actual, y en el comando se debe reemplazar 3.18 con la versión actual.

Code

BiocManager::install("dada2", version = "3.18")

Bioconductor version 3.18 (BiocManager 1.30.22), R 4.3.1 (2023-06-16)

Warning: package(s) not installed when version(s) same as or greater than current; use
  `force = TRUE` to re-install: 'dada2'

Old packages: 'ALL', 'Biostrings', 'evaluate', 'GenomicRanges',
  'RcppArmadillo', 'S4Vectors', 'xfun'

phyloseq

El siguiente comando instala phyloseq.

Code

BiocManager::install("phyloseq")

Bioconductor version 3.18 (BiocManager 1.30.22), R 4.3.1 (2023-06-16)

Warning: package(s) not installed when version(s) same as or greater than current; use
  `force = TRUE` to re-install: 'phyloseq'

Old packages: 'ALL', 'Biostrings', 'evaluate', 'GenomicRanges',
  'RcppArmadillo', 'S4Vectors', 'xfun'

Proceso general

El punto de partida del análisis con dada2 es un conjunto de archivos fastq demultiplexados correspondientes a las muestras de su estudio de secuenciación de amplicones. Es decir, dada2 espera que haya un archivo fastq individual para cada muestra (o dos archivos fastq, uno forward y otro reverse, para cada muestra).

Desmultiplexación

La desmultiplexación a menudo se realiza en el sitio de secuenciación, pero si eso no se ha hecho, hay una variedad de herramientas que lo hacen, incluidos los populares scripts de Python QIIME split_libraries_fastq.py seguido de split_sequence_file_on_sample_ids.py, y la utilidad idemp, entre otras.

Además, dada2 espera que no haya bases no biológicas en los datos de secuenciación. Esto puede requerir un procesamiento previo con software externo si, como es relativamente común, los primers de PCR se incluyeron en la región del amplicón secuenciada.

Una vez que los archivos fastq demultiplexados sin nucleótidos no biológicos están disponibles, el proceso dada2 procede de la siguiente manera:

Filtrar y recortar: filterAndTrim()
Desreplicar: derepFastq()
Tasas de error: learnErrors()
Inferir la composición de la muestra: dada()
Fusionar lecturas emparejadas: mergePairs()
Crear tabla de secuencia: makeSequenceTable()
Eliminar quimeras: removeBimeraDenovo()
Asignar taxonomia: assignTaxonomy()

dada2 genera una tabla de características de variantes de secuencia de amplicones (una tabla ASV): una matriz, en la que cada fila corresponde a una muestra procesada, y cada columna corresponde a una secuencia inferida, no quimérica, en la que el valor de cada entrada es el número de veces que se observó ASV en ese muestra. Esta tabla es análoga a la tabla OTU tradicional, excepto que tiene una resolución más alta, es decir, variantes de secuencia de amplicones exactas en lugar de grupos (generalmente 97%) de lecturas de secuenciación.

Ejercicio

Partiendo de un grupo de secuencias de microbiota 16S 2x250 de Illumina se ejecutará la secuencia dada2-phyloseq.

bibliotecas

El primer paso siempre es cargar las bibliotecas que se utilizarán. El comando packageVersion("paquete") les indica la versión del paquete que se está utilizando, y es opcional incluirlo.

Code

library("dada2"); packageVersion("dada2")

Loading required package: Rcpp

[1] '1.30.0'

Code

library(phyloseq)
library(Biostrings)

Loading required package: BiocGenerics


Attaching package: 'BiocGenerics'

The following objects are masked from 'package:stats':

    IQR, mad, sd, var, xtabs

The following objects are masked from 'package:base':

    anyDuplicated, aperm, append, as.data.frame, basename, cbind,
    colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
    get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
    match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
    Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort,
    table, tapply, union, unique, unsplit, which.max, which.min

Loading required package: S4Vectors

Loading required package: stats4


Attaching package: 'S4Vectors'

The following object is masked from 'package:utils':

    findMatches

The following objects are masked from 'package:base':

    expand.grid, I, unname

Loading required package: IRanges


Attaching package: 'IRanges'

The following object is masked from 'package:phyloseq':

    distance

Loading required package: XVector

Loading required package: GenomeInfoDb


Attaching package: 'Biostrings'

The following object is masked from 'package:base':

    strsplit

Code

library(ggplot2)

archivos de secuenciación

Una vez descargados los archivos de secuenciación hay que especificar la dirección de la carpeta donde se encuentran las secuencias. list.files permite ver el contenido de la carpeta.

Code

path <- "/Users/oscar/Library/CloudStorage/GoogleDrive-o.medinac@gmail.com/My Drive/Clase microbiota" # CAMBIAR al directorio que contiene los archivos fastq descomprimidos.
list.files(path)

 [1] "F3D143_S209_L001_R1_001.fastq"          
 [2] "F3D143_S209_L001_R2_001.fastq"          
 [3] "F3D144_S210_L001_R1_001.fastq"          
 [4] "F3D144_S210_L001_R2_001.fastq"          
 [5] "F3D145_S211_L001_R1_001.fastq"          
 [6] "F3D145_S211_L001_R2_001.fastq"          
 [7] "F3D146_S212_L001_R1_001.fastq"          
 [8] "F3D146_S212_L001_R2_001.fastq"          
 [9] "F3D147_S213_L001_R1_001.fastq"          
[10] "F3D147_S213_L001_R2_001.fastq"          
[11] "F3D148_S214_L001_R1_001.fastq"          
[12] "F3D148_S214_L001_R2_001.fastq"          
[13] "F3D149_S215_L001_R1_001.fastq"          
[14] "F3D149_S215_L001_R2_001.fastq"          
[15] "F3D150_S216_L001_R1_001.fastq"          
[16] "F3D150_S216_L001_R2_001.fastq"          
[17] "filtered"                               
[18] "HMP_MOCK.v35.fasta"                     
[19] "Mock_S280_L001_R1_001.fastq"            
[20] "Mock_S280_L001_R2_001.fastq"            
[21] "silva_nr99_v138.1_train_set.fa"         
[22] "silva_nr99_v138.1_wSpecies_train_set.fa"
[23] "silva_species_assignment_v138.1.fa"

Para generar la lista de archivos a analizar se crea un archivo para las secuencias forward y uno para las reverse. Para que todas las secuencias incluidas hay que indicar el patrón del nombre para ambos casos.

Code

# Los nombres de los archivos fastq forward y reverse tienen el formato: SAMPLENAME_R1_001.fastq y SAMPLENAME_R2_001.fastq
fnFs <- sort(list.files(path, pattern="_R1_001.fastq", full.names = TRUE))
fnRs <- sort(list.files(path, pattern="_R2_001.fastq", full.names = TRUE))
# Extraer los nombres de las muestras, asumiendo que los archivos tienen el formato: SAMPLENAME_XXX.fastq
sample.names <- sapply(strsplit(basename(fnFs), "_"), `[`, 1)

calidad de la secuenciación

Visualizar la calidad de lectura de las secuencias contenidas en los archivos generados previamente.

Code

plotQualityProfile(fnFs[1:2])

Code

plotQualityProfile(fnRs[1:2])

La gráfica en escala de grises resultante es un mapa de calor de la frecuencia de cada puntuación de calidad en cada base. La linea verde representa la calidad promedio, y las lineas naranja representan los cuartiles.

En este ejemplo se puede observar que las secuencias forward tiene buena calidad, mientras que las reverse tienen peor calidad, específicamente al final (esto es común para secuencias de Illumina). En estas secuencias utilizaremos los primeros 240nt y 160nt.

Note

Siempre es recomendable recortar los últimos 10 nucleótidos para evitar errores de cualquier naturaleza debida a la secuenciación.

Warning

Después de truncar las secuencias aún deben poder traslaparse para realizar la fusión (paso 5).

filtrar

Asignar los nombres de archivo para los archivos fastq.gz filtrados.

Code

filtFs <- file.path(path, "filtered", paste0(sample.names, "_F_filt.fastq.gz"))
filtRs <- file.path(path, "filtered", paste0(sample.names, "_R_filt.fastq.gz"))
names(filtFs) <- sample.names
names(filtRs) <- sample.names

Los parámetros estandar para filtrar son: maxN=0 (DADA2 no requiere Ns), truncQ=2, rm.phix=TRUE y maxEE=2. El parámetro maxEE establece el número máximo de “errores esperados” permitidos en la lectura.

Code

out <- filterAndTrim(fnFs, filtFs, fnRs, filtRs, truncLen=c(240,160),
              maxN=0, maxEE=c(2,2), truncQ=2, rm.phix=TRUE,
              compress=TRUE, multithread=TRUE) # En Windows dejar multithread=FALSE
head(out)

                              reads.in reads.out
F3D143_S209_L001_R1_001.fastq     3178      2941
F3D144_S210_L001_R1_001.fastq     4827      4312
F3D145_S211_L001_R1_001.fastq     7377      6741
F3D146_S212_L001_R1_001.fastq     5021      4560
F3D147_S213_L001_R1_001.fastq    17070     15637
F3D148_S214_L001_R1_001.fastq    12405     11413

tasas de error

DADA2 usa un modelo paramátrico de error, y cada amplicón tiene una tasa de error diferente. Este es un modelo de machine-learning, donde se intercala la estimación de la tasa de error y la inferencia de la comnposición de la muestra, hasta generar una solución consistente. Este es un proceso relativamente rápido, pero tomará mas tiempo a medida que el procesador sea menos reciente.

Code

errF <- learnErrors(filtFs, multithread=TRUE)

16072080 total bases in 66967 reads from 9 samples will be used for learning the error rates.

Code

errR <- learnErrors(filtRs, multithread=TRUE)

10714720 total bases in 66967 reads from 9 samples will be used for learning the error rates.

Al graficar los errores se muestran todas las posibles combinaciones (A→C, A→G, etc). Las lineas negras muestran las tasas de error estimado. Las lineas rojas muestran las tasas de error esperado. Se debe esperar un buen ajuste del error estimado (linea negra) con las tasas observadas (puntos rojos), y la tasa de error debe ser menor a medida la calidad sea mayor.

Code

plotErrors(errF, nominalQ=TRUE)

Warning: Transformation introduced infinite values in continuous y-axis
Transformation introduced infinite values in continuous y-axis

inferencia

El algoritmo de inferencia de las muestras se aplica en las secuencias filtradas y recortadas.

Code

dadaFs <- dada(filtFs, err=errF, multithread=TRUE)

Sample 1 - 2941 reads in 939 unique sequences.
Sample 2 - 4312 reads in 1267 unique sequences.
Sample 3 - 6741 reads in 1756 unique sequences.
Sample 4 - 4560 reads in 1438 unique sequences.
Sample 5 - 15637 reads in 3590 unique sequences.
Sample 6 - 11413 reads in 2762 unique sequences.
Sample 7 - 12017 reads in 3021 unique sequences.
Sample 8 - 5032 reads in 1566 unique sequences.
Sample 9 - 4314 reads in 897 unique sequences.

Code

dadaRs <- dada(filtRs, err=errR, multithread=TRUE)

Sample 1 - 2941 reads in 880 unique sequences.
Sample 2 - 4312 reads in 1286 unique sequences.
Sample 3 - 6741 reads in 1803 unique sequences.
Sample 4 - 4560 reads in 1265 unique sequences.
Sample 5 - 15637 reads in 3414 unique sequences.
Sample 6 - 11413 reads in 2522 unique sequences.
Sample 7 - 12017 reads in 2771 unique sequences.
Sample 8 - 5032 reads in 1415 unique sequences.
Sample 9 - 4314 reads in 732 unique sequences.

El resultado nos muestra el objeto dada-class. Este algoritmo infiere las variantes de secuencia verdaderas a partir de las secuencias únicas en cada muestra.

Code

dadaFs[[1]]

dada-class: object describing DADA2 denoising results
74 sequence variants were inferred from 939 input unique sequences.
Key parameters: OMEGA_A = 1e-40, OMEGA_C = 1e-40, BAND_SIZE = 16

Note

dada-class contiene mas información, como diagnósticos múltiples de la calidad. help("dada-class") muestra las diferentes opciones.

fusión

Se fusionan las lecturas forward y reverse para obtener secuencias completas sin ruido. Esto es el mismo principio que los “contig” usados en secuenciación de cromosomas. Los valores estandar necesitan un traslape de las secuencias de al menos 12 bases, con 100% de identidad en esta región (pueden modificarse estos parámetros en la función).

Code

mergers <- mergePairs(dadaFs, filtFs, dadaRs, filtRs, verbose=TRUE)

2549 paired-reads (in 59 unique pairings) successfully merged out of 2781 (in 118 pairings) input.

3624 paired-reads (in 53 unique pairings) successfully merged out of 4106 (in 156 pairings) input.

6080 paired-reads (in 81 unique pairings) successfully merged out of 6515 (in 197 pairings) input.

3954 paired-reads (in 91 unique pairings) successfully merged out of 4384 (in 190 pairings) input.

14238 paired-reads (in 143 unique pairings) successfully merged out of 15369 (in 351 pairings) input.

10528 paired-reads (in 120 unique pairings) successfully merged out of 11159 (in 277 pairings) input.

11155 paired-reads (in 137 unique pairings) successfully merged out of 11798 (in 299 pairings) input.

4350 paired-reads (in 85 unique pairings) successfully merged out of 4806 (in 180 pairings) input.

4269 paired-reads (in 20 unique pairings) successfully merged out of 4281 (in 28 pairings) input.

La visualización muestra una lista de las secuencias y su abundancia, entre otros valores. Cada muestra crea un data.frame independiente.

Code

head(mergers[[1]])

                                                                                                                                                                                                                                                      sequence
1 TACGGAGGATGCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGTAGGCGGCCTGCCAAGTCAGCGGTAAAATTGCGGGGCTCAACCCCGTACAGCCGTTGAAACTGCCGGGCTCGAGTGGGCGAGAAGTATGCGGAATGCGTGGTGTAGCGGTGAAATGCATAGATATCACGCAGAACCCCGATTGCGAAGGCAGCATACCGGCGCCCTACTGACGCTGAGGCACGAAAGTGCGGGGATCAAACAGG
2 TACGGAGGATGCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGTAGGCGGGCTTTTAAGTCAGCGGTAAAAATTCGGGGCTCAACCCCGTCCGGCCGTTGAAACTGGGGGCCTTGAGTGGGCGAGAAGAAGGCGGAATGCGTGGTGTAGCGGTGAAATGCATAGATATCACGCAGAACCCCGATTGCGAAGGCAGCCTTCCGGCGCCCTACTGACGCTGAGGCACGAAAGTGCGGGGATCGAACAGG
3 TACGGAGGATGCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGCAGGCGGAAGATCAAGTCAGCGGTAAAATTGAGAGGCTCAACCTCTTCGAGCCGTTGAAACTGGTTTTCTTGAGTGAGCGAGAAGTATGCGGAATGCGTGGTGTAGCGGTGAAATGCATAGATATCACGCAGAACTCCGATTGCGAAGGCAGCATACCGGCGCTCAACTGACGCTCATGCACGAAAGTGTGGGTATCGAACAGG
4 TACGGAGGATGCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGTAGGCGGGCTGTTAAGTCAGCGGTCAAATGTCGGGGCTCAACCCCGGCCTGCCGTTGAAACTGGCGGCCTCGAGTGGGCGAGAAGTATGCGGAATGCGTGGTGTAGCGGTGAAATGCATAGATATCACGCAGAACTCCGATTGCGAAGGCAGCATACCGGCGCCCGACTGACGCTGAGGCACGAAAGCGTGGGTATCGAACAGG
5 TACGGAGGATGCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGCAGGCGGACTCTCAAGTCAGCGGTCAAATCGCGGGGCTCAACCCCGTTCCGCCGTTGAAACTGGGAGCCTTGAGTGCGCGAGAAGTAGGCGGAATGCGTGGTGTAGCGGTGAAATGCATAGATATCACGCAGAACTCCGATTGCGAAGGCAGCCTACCGGCGCGCAACTGACGCTCATGCACGAAAGCGTGGGTATCGAACAGG
6 TACGGAGGATGCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGTAGGCGGGATGCCAAGTCAGCGGTAAAAAAGCGGTGCTCAACGCCGTCGAGCCGTTGAAACTGGCGTTCTTGAGTGGGCGAGAAGTATGCGGAATGCGTGGTGTAGCGGTGAAATGCATAGATATCACGCAGAACTCCGATTGCGAAGGCAGCATACCGGCGCCCTACTGACGCTGAGGCACGAAAGCGTGGGTATCGAACAGG
  abundance forward reverse nmatch nmismatch nindel prefer accept
1       244       1       1    148         0      0      2   TRUE
2       231       2       2    148         0      0      2   TRUE
3       228       3       4    148         0      0      1   TRUE
4       204       4       5    148         0      0      1   TRUE
5       176       5       6    148         0      0      1   TRUE
6       152       6       3    148         0      0      2   TRUE

tabla de secuencias

La tabla de variantes de secuencia del amplicón (ASV) es una tabla de mayor resolución que las tablas de OTUs producidas por otros tipos de análisis.

Code

seqtab <- makeSequenceTable(mergers)
dim(seqtab)

[1]   9 224

Code

# Distribución de la longitud de las secuencias
table(nchar(getSequences(seqtab)))


252 253 254 255 
 73 143   6   2

Note

Las secuencias mas largas o cortas que lo esperado, resultado de inespecificidades, se pueden remover con el comando seqtab2 <- seqtab[,nchar(colnames(seqtab)) %in% 250:256]), donde se delimita la longitud de los fragmentos esperados.

quimeras

El algoritmo corrige sustituciones, inserciones y deleciones, pero no las quimeras. Después de quitar el ruido es muy simple identificar ASVs quiméricos.

Code

seqtab.nochim <- removeBimeraDenovo(seqtab, method="consensus", multithread=TRUE, verbose=TRUE)

Identified 54 bimeras out of 224 input sequences.

Code

dim(seqtab.nochim)

[1]   9 170

La frecuencia de secuencias no quiméricas varia en cada grupo de datos, y depende de factores como la complejidad de la muestra y el experimento.

Code

sum(seqtab.nochim)/sum(seqtab)

[1] 0.9512898

Note

La cantidad de quimeras debe ser muy baja. Si se eliminan muchas secuencias en este paso es necesario revisar los pasos anteriores. Lo mas probable es que se deba secuencias de primers con nucleotidos ambiguos.

resúmen de lecturas

Es posible obtener una tabla con el número de lecturas que se mantuvieron en cada uno de los pasos del análisis.

Code

getN <- function(x) sum(getUniques(x))
track <- cbind(out, sapply(dadaFs, getN), sapply(dadaRs, getN), sapply(mergers, getN), rowSums(seqtab.nochim))
# Si solo se esta procesando una muestra remover los comandos `sapply`: por ejemplo reemplazar sapply(dadaFs, getN) con getN(dadaFs)
colnames(track) <- c("input", "filtered", "denoisedF", "denoisedR", "merged", "nonchim")
rownames(track) <- sample.names
head(track)

       input filtered denoisedF denoisedR merged nonchim
F3D143  3178     2941      2822      2864   2549    2515
F3D144  4827     4312      4151      4224   3624    3485
F3D145  7377     6741      6593      6628   6080    5821
F3D146  5021     4560      4450      4467   3954    3865
F3D147 17070    15637     15446     15507  14238   13005
F3D148 12405    11413     11244     11269  10528    9936

Este es el paso final para revisar la integridad del análisis. No debe haber un paso donde haya una disminución súbita del número de lecturas.

asignar taxonomia

dada2 utiliza un metodo de clasificación Bayesiana para realizar la asignación taxonómica a las ASVs. La función assignTaxonomytoma la tabla de secuencias y un conjunto de secuencias de referencia.

El sitio web de dada2 contiene secuencias de referencia de las bases de datos RDP, Silva, y UNITE. Hay secuencias para asignación inicial a nivel de género o especie, así como asignación de especie secuencial (posterior a asignación de género) taxa <- addSpecies(taxa, "~/directorio/silva_species_assignment_.fa"). En este ejemplo en la carpera de secuencias se encuentran descargados los archivos silva_nr99_v138.1_train_set.

Warning

Es importante que los archivos silva esten en la misma carpeta que las secuencias.

Code

taxa <- assignTaxonomy(seqtab.nochim, "/Users/oscar/Library/CloudStorage/GoogleDrive-o.medinac@gmail.com/My Drive/Clase microbiota/silva_nr99_v138.1_train_set.fa", multithread=TRUE)

Las asignaciones se pueden visualizar en una tabla

Code

taxa.print <- taxa
# Se remueven los nombres de las secuencias solamente para visualizarlas
rownames(taxa.print) <- NULL
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" "Muribaculaceae"
[6,] "Bacteria" "Bacteroidota" "Bacteroidia" "Bacteroidales" "Muribaculaceae"
     Genus
[1,] NA   
[2,] NA   
[3,] NA   
[4,] NA   
[5,] NA   
[6,] NA

Note

Es importante que los archivos Silva esten en la misma carpeta que las secuencias.

Si la asignación no parece correcta (muchas secuencias estan asignadas como Eukaryota NA NA NA NA NA), la orientación de las lecturas podría esta en sentido opuesta a las de referencia. dada2 puede intentar realizar la asignación taxonómica con el reverso-complementario de las secuencias con el argumento tryRC=TRUE

phyloseq

A partir de esta tabla se genera un data.frame de las muestras. En las muestras del ejemplo el nombre incluye el sexo (G), número de ratón (X), y el dia en que fueron muestreadas (DY). Por ejemplo, F3D147_S213_L001_R2_001 es hembra, ratón 3, del dia 147.

Code

samples.out <- rownames(seqtab.nochim)
subject <- sapply(strsplit(samples.out, "D"), `[`, 1)
gender <- substr(subject,1,1)
subject <- substr(subject,2,999)
day <- as.integer(sapply(strsplit(samples.out, "D"), `[`, 2))
samdf <- data.frame(Subject=subject, Gender=gender, Day=day)
samdf$When <- "Early"
samdf$When[samdf$Day>200] <- "Late"
rownames(samdf) <- samples.out

Con estos datos de genera un objeto de phyloseq directamente de los resultados de dada2.

Code

ps <- phyloseq(otu_table(seqtab.nochim, taxa_are_rows=FALSE), 
               sample_data(samdf), 
               tax_table(taxa))
# Este comando elimina del análisis posterior las secuencias no deseadas.
ps <- prune_samples(sample_names(ps) != "Mock", ps)

Es más conveniente utilizar nombres cortos para las ASVs (por ejemplo ASV17) en lugar de la secuencia completa de DNA, con el fin de generar tablas y visualizaciones adecuadas, pero se debe mantener la secuencia para análisis en otras bases de datos. Estas secuencias completas se almacenan en la sección refseq del objeto phyloseq, y se pueden recuperar con el comando refseq(ps).

Code

# Renombrar los taxa a un nombre corto
dna <- Biostrings::DNAStringSet(taxa_names(ps))
names(dna) <- taxa_names(ps)
ps <- merge_phyloseq(ps, dna)
taxa_names(ps) <- paste0("ASV", seq(ntaxa(ps)))
ps

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 170 taxa and 8 samples ]
sample_data() Sample Data:       [ 8 samples by 4 sample variables ]
tax_table()   Taxonomy Table:    [ 170 taxa by 6 taxonomic ranks ]
refseq()      DNAStringSet:      [ 170 reference sequences ]

Finalmente se pueden visualizar los análisis usando las diferentes opciones de phyloseq.

diversidad alfa

Code

# plot_richness(physeq, x = "samples", color = NULL, shape = NULL, title = NULL, scales = "free_y", nrow = 1, shsi = NULL, measures = NULL, sortby = NULL)
# En nuestro ejemplo estan graficadas por dia, y separadas en temprano y tardio.
plot_richness(ps, x="Day", measures=c("Observed", "Chao1", "ACE", "Shannon", "Simpson", "InvSimpson", "Fisher"), color="When")

Warning in estimate_richness(physeq, split = TRUE, measures = measures): The data you have provided does not have
any singletons. This is highly suspicious. Results of richness
estimates (for example) are probably unreliable, or wrong, if you have already
trimmed low-abundance taxa from the data.

We recommended that you find the un-trimmed data and retry.

ordenadas

Code

# plot_ordination(physeq, ordination, type = "samples", axes = 1:2, color = NULL, shape = NULL, label = NULL, title = NULL, justDF = FALSE)
ps.prop <- transform_sample_counts(ps, function(otu) otu/sum(otu))
ord.nmds.bray <- ordinate(ps.prop, method="NMDS", distance="bray")

Run 0 stress 0.08020701 
Run 1 stress 0.08020702 
... Procrustes: rmse 1.84494e-05  max resid 2.467851e-05 
... Similar to previous best
Run 2 stress 0.08020701 
... Procrustes: rmse 2.59855e-06  max resid 4.364363e-06 
... Similar to previous best
Run 3 stress 0.08020701 
... Procrustes: rmse 5.509653e-07  max resid 7.300018e-07 
... Similar to previous best
Run 4 stress 0.08020701 
... Procrustes: rmse 2.528789e-06  max resid 4.241341e-06 
... Similar to previous best
Run 5 stress 0.08020701 
... Procrustes: rmse 1.938798e-06  max resid 3.222313e-06 
... Similar to previous best
Run 6 stress 0.09330431 
Run 7 stress 0.08020701 
... Procrustes: rmse 1.174677e-05  max resid 1.818458e-05 
... Similar to previous best
Run 8 stress 0.08020701 
... Procrustes: rmse 4.684691e-07  max resid 6.559921e-07 
... Similar to previous best
Run 9 stress 0.07193353 
... New best solution
... Procrustes: rmse 0.237632  max resid 0.3307452 
Run 10 stress 0.07193353 
... Procrustes: rmse 2.220985e-06  max resid 3.524983e-06 
... Similar to previous best
Run 11 stress 0.06847299 
... New best solution
... Procrustes: rmse 0.06776059  max resid 0.1286094 
Run 12 stress 0.08020701 
Run 13 stress 0.08020701 
Run 14 stress 0.07193353 
Run 15 stress 0.2520583 
Run 16 stress 0.07193353 
Run 17 stress 0.08020701 
Run 18 stress 0.06847297 
... New best solution
... Procrustes: rmse 3.417334e-05  max resid 4.45469e-05 
... Similar to previous best
Run 19 stress 0.07193353 
Run 20 stress 0.08020701 
*** Best solution repeated 1 times

Code

plot_ordination(ps.prop, ord.nmds.bray, color="When", title="Bray NMDS")

barras

Code

# plot_bar(physeq, x="Sample", y="Abundance", fill=NULL, title=NULL, facet_grid=NULL)
# Aqui se seleccionan las especies que se quieven visualizar. En este caso son las 20 mas abundantes solamente.
top20 <- names(sort(taxa_sums(ps), decreasing=TRUE))[1:20]
ps.top20 <- transform_sample_counts(ps, function(OTU) OTU/sum(OTU))
ps.top20 <- prune_taxa(top20, ps.top20)
# Estan graficadas por dia, agrupadas en temprano y tardío, y coloreadas por familia.
plot_bar(ps.top20, x="Day", fill="Family") + facet_wrap(~When, scales="free_x")

redes

Code

# plot_network(g, physeq=NULL, type="samples", color=NULL, shape=NULL, point_size=4, alpha=1, label="value", hjust = 1.35, line_weight=0.5, line_color=color, line_alpha=0.4, layout.method=layout.fruchterman.reingold, title=NULL)
data(enterotype)
ig <- make_network(enterotype, max.dist=0.3)
plot_network(ig, enterotype, color="Enterotype", line_weight=0.3, label=NULL)

Note

Siempre hay información extra en la documentación de los paquetes.