11 Introducción al análisis por 16S rRNA usando DADA2 vs. Mothur

En esta sección vamos a explorar un set de datos de 98 muestras de piel de ballenas, jorobadas, azules y otras. Las muestras fueron obtenidas al amplificar y secuenciar un fragmento de la región variable 4 (V4) del gen 16S rRNA.

Las secuencias están disponibles en el directorio de trabajo del curso /MCV502_2019/data/whale_pe/ y también en la base de datos SRA de NCBI bajo el BioProject PRJNA428495.

11.1 DADA2

Vamos a analizar nuestros datos siguiendo la pipeline DADA2 y luego comparar esos resultados con el análisis de los mismos datos pero según la pipeline mothur. DADA2 ofrece ventajas con respecto a la estrategia de formar clusters (OTUs) en varios aspectos que incluyen mayor resolución, nombres consistentes entre diferentes análisis, mejor estimación de abundancias relativas, etc. Los siguientes artículos ofrecen una discusión más detallada al respecto ( artículo y artículo).

• Lo primero que vamos a hacer es configurar el directorio donde están las reads:

set.seed(100) 
# Configuramos un número de semilla para que los resultados 
# sean comparables entre estudiantes

# Simplemente cambiamos este path al que corresponda según donde están las reads
miseq_path <- "/Users/katterinne/Dropbox/CastroLab/MCV502_BCM634_BIO625/MCV502_2019/data/whale_pe/" 
list.files(miseq_path)

• Ahora necesitamos cargar los paquetes que vamos a utilizar. Recuerda que no necesitamos instalarlos porque ya lo hicimos en la sección ??.

# Generamos 3 objetos con la lista de paquetes provenientes de los repositorios CRAN, Bioconductor, y GitHub
cran_packages <- c("knitr", "qtl", "bookdown", "ggplot2", "grid", "gridExtra", "tidyverse", "xtable", "devtools", "kableExtra", "remotes")
bioc_packages <- c("phyloseq", "dada2", "DECIPHER", "phangorn", "ggpubr", "DESeq2", "genefilter", "philr")
git_packages <- c("btools", "fantaxtic")

# Cargar paquetes en la actual sesión de R
sapply(c(cran_packages, bioc_packages, git_packages), require, character.only = TRUE)

• En DADA2 las reads se trabajan inicialmente por separado, es decir, la copia forward o R1 separada de la copia reverse o R2. Por esto, nos tenemos que asegurar que ambos archivos estén ordenados. Luego, manipulamos el nombre de los archivos para generar automáticamente el nombre de las muestras en nuestro análisis:

fnFs <- sort(list.files(miseq_path, pattern="_1.fastq"))
fnRs <- sort(list.files(miseq_path, pattern="_2.fastq"))

# Extracción del nombre de las muestras
sampleNames <- sapply(strsplit(fnFs, "_"), `[`, 1)
# Especificamos el path completo para evitar errores de ambigüedad
fnFs <- file.path(miseq_path, fnFs)
fnRs <- file.path(miseq_path, fnRs)

• El paquete DADA2 tiene muchas funciones útiles para el preprocesamiento de los datos. Acá lo que hacemos es visualizar el perfil de calidad de las muestras para cada par de 5’ a 3’:

plotQualityProfile(fnFs[1:2])
plotQualityProfile(fnRs[1:2])

• Ahora vamos a proceder con el filtrado y corte de las reads de acuerdo a lo observado en los gráficos de calidad. Primero creamos un directorio donde vamos a poner las reads una vez realizado el control de calidad, y luego realizamos el dicho control. En específico usamos los argumentos trunLen =, maxN=, maxEE= y rm.phix= para indicar el corte promedio de cada read, el número máximo de bases indeterminadas, el número máximo de errores, y si es que queremos remover secuencias pertenecientes al control interno de Illumina.

# Creamos un directorio para poner las reads "limpias"
filt_path <- file.path(miseq_path, "filtered") 
if(!file_test("-d", filt_path)) dir.create(filt_path)
filtFs <- file.path(filt_path, paste0(sampleNames, "_F_filt.fastq.gz"))
filtRs <- file.path(filt_path, paste0(sampleNames, "_R_filt.fastq.gz"))

# Y finalmente procedemos con el control de calidad
out <- filterAndTrim(fnFs, filtFs, fnRs, filtRs, truncLen=c(250,200),
              maxN=0, maxEE=c(2,2), truncQ=2, rm.phix=TRUE,
              compress=TRUE, multithread=TRUE) 
# Si usan Windows, configuren multithread=FALSE

• DADA2 genera un modelo probabilístico de errores con el cual puede filtrar reads erróneas y así usar las restantes directamente para la etapa de clasificación taxonómica. Esta parte del método es la que nos permite tener una mayor resolución en comparación a los análisis basados en OTUs. Como las muestras probablemente tienen reads idénticas, no es eficiente usar cada una de ellas para los pasos río abajo. Por esto, DADA2 recomienda de-replicar las muestras para así disminuir la redundancia y avanzar eficientemente.

derepFs <- derepFastq(filtFs, verbose=TRUE)
derepRs <- derepFastq(filtRs, verbose=TRUE)

# Acá simplemente agregamos los nombres de las muestras al objeto de-replicado
names(derepFs) <- sampleNames
names(derepRs) <- sampleNames

• Una vez con las muestras de-replicadas procedemos a generar un modelo de errores. Para mayor detalle sobre este crucial paso revisen el artículo original aquí.

errF <- learnErrors(filtFs, multithread=TRUE)
errR <- learnErrors(filtRs, multithread=TRUE)

# Graficamos los errores para cada par
plotErrors(errF)
plotErrors(errR)

Cada panel en el gráfico nos indica la frecuencia de error de nuecleótido a nucleótido para todas las combinaciones. Naturalmente, las bases con mayor puntaje de calidad exhiben una menor frecuencia de error.

• Con esta información entonces procedemos al paso más importante de la pipeline, i.e., la inferencia de las Amplicon Sequence Variants (ASVs).

dadaFs <- dada(derepFs, err=errF, multithread=TRUE, pool = TRUE)
dadaRs <- dada(derepRs, err=errR, multithread=TRUE, pool = TRUE)

• Ahora nos queda hacer un poco de aseo: primero unimos las reads R1 y R2, luego generamos una tabla de secuencias que vamos a utilizar más tarde, y finalmente removemos secuencias artefactuales producto de la amplificación por PCR.

# Unimos las reads
mergers <- mergePairs(dadaFs, derepFs, dadaRs, derepRs)
# Generamos una tabla de secuencias
seqtabAll <- makeSequenceTable(mergers[!grepl("Mock", names(mergers))])
table(nchar(getSequences(seqtabAll)))
# Removemos las secuencias quiméricas
seqtabNoC <- removeBimeraDenovo(seqtabAll)

La tabla de secuencias sin las quimeras, es la tabla que usamos para realizar la clasificación taxonómica. En principio, podríamos usar cualquiera de las tres bases de datos más populares para clasificación de secuencias del 16S rRNA, i.e., GreenGenes, RDP o SILVA.

• En nuestro práctico vamos a utilizar SILVA, en particular, la versión 132.

fastaRef <- "/Users/katterinne/Dropbox/CastroLab/MCV502_BCM634_BIO625/MCV502_2019/data/silva_nr_v132_train_set.fa" # Simplemente cambiamos este path al que corresponda según donde está archivo "silva_nr_v132_train_set.fa"
taxTab <- assignTaxonomy(seqtabNoC, refFasta = fastaRef, multithread=TRUE)

# En el caso de querer agregar el rango taxonómico de especies, simplemente usamos una base de datos extra, la cual contiene esta información. 
taxTabExtra <- addSpecies(taxTab, "/Users/katterinne/Dropbox/CastroLab/MCV502_BCM634_BIO625/MCV502_2019/data/silva_species_assignment_v132.fa", verbose=TRUE) # Simplemente cambiamos este path al que corresponda según donde está archivo "silva_species_assignment_v132.fa"
unname(head(taxTab)) -> tabla
colnames(tabla) <- c("Kingdom", "Phylum", "Order", "Class", "Family", "Genus")

Ya que nuestro perfil taxonómico se construye con secuencias homólogas del 16S rRNA gene, podemos usar estas secuencias para inferir un árbol filogenético. Existen muchos paquetes de r que pueden hacer esto y aquí escogemos phangorn para hacer una inferencia basada en Maximum Likelihood.

*Primero alineamos las secuencias:

seqs <- getSequences(seqtabNoC)
# Este paso propaga los nombres de las secuencias al árbol
names(seqs) <- seqs
alignment <- AlignSeqs(DNAStringSet(seqs), anchor=NA,verbose=FALSE)

• Y con ese alineamiento inferimos un árbol de partida o starting tree para inicializar la búsqueda por ML. También ajustamos un modelo de sustitución nucleotídica para parametrizar la tasa de cambio de un nucleótido a otro y asó poder inferir correctamente el largo de las ramas y la topología del árbol.

phangAlign <- phyDat(as(alignment, "matrix"), type="DNA")
dm <- dist.ml(phangAlign)
treeNJ <- NJ(dm) 
fit = pml(treeNJ, data=phangAlign)
fitGTR <- update(fit, k=4, inv=0.2)
fitGTR <- optim.pml(fitGTR, model="GTR", optInv=TRUE, optGamma=TRUE,
        rearrangement = "stochastic", control = pml.control(trace = 0))

detach("package:phangorn", unload=TRUE)

Ahora, tenemos todos los ingredientes para formar un objeto de R que contenga todo lo que nos importa en un experimento metagenómico, i.e., una tabla de cuentas que indica el número de reads por seceuncia del 16S rRNA gene, una tabla con el linaje taxonómico de esas secuencias, un árbol filogenético que relaciona esas secuencias entre sí, y finalmente, una tabla con variables asociadas a nuestras muestras, también llamada metadata.

• Con estos elementos procedemos a generar un objeto phyloseq:

samdf <- read.csv("/Users/katterinne/Dropbox/CastroLab/MCV502_BCM634_BIO625/MCV502_2019/data/metadata.csv", header=TRUE, row.names = 1) # Simplemente cambiamos este path al que corresponda según donde está archivo "metadata.csv"
rownames(seqtabNoC) %in% rownames(samdf)
all(rownames(seqtabAll) %in% samdf$run)

ps <- phyloseq(otu_table(seqtabNoC, taxa_are_rows=FALSE), 
               sample_data(samdf), 
               tax_table(taxTab),
               phy_tree(fitGTR$tree))
ps <- prune_samples(sample_names(ps) != "Mock", ps) # Remover potenciales muestras sintéticas
ps

• Y el producto final del análisis por DADA2 es entonces un objeto phyloseq:

Slot	Descripción	Resultado
otu_table()	OTU Table:	[ 1476 taxa and 98 samples ]
sample_data()	Sample Data:	[ 98 samples by 13 sample variables ]
tax_table()	Taxonomy Table:	[ 1476 taxa by 6 taxonomic ranks ]
phy_tree()	Phylogenetic Tree:	[ 1476 tips and 1474 internal nodes ]

11.2 Mothur

El otro paradigma en análisis de 16S rRNA es agrupar las reads de acuerdo a su porcentaje de similitud a un umbral determinado. Típicamente, ese umbral es de 97% de identidad de secuencia, el cual es totalmente arbitrario. Cuando se comenzó a usar high-throughput sequencing para hacer perfiles taxonómicos de comunidades microbianas, fue necesario agrupar las secuencias usando este método para minimizar los efectos de errores en la secuenciación y clasificación de reads. A estas agrupaciones de reads se les llama Operational Taxonomic Units (OTUs). Es difícil asignar una taxonomía clara a cada OTU por al menos tres razones:

1. No existe una razón biológica de porqué especies de microorganismos deben agruparse en grupos discretos de secuencias al 97%.
2. Existen microorganismos de distintas especies que poseen la misma secuencia del gen 16S rRNA.
3. Al revés, existen microorgansimos de la misma especie que tienen diferente secuencia del gen 16S rRNA.

Por estos motivos una OTU no corresponde necesariamente a una especie sino que a un linaje fabricado pero operacional para poder investigar la diversidad y estructura de una comunidad microbiana. En este práctico, vamos a revisar los pasos necesarios en mothur para poder ir de reads hasta tabla de taxonomía lista para ser analizada en R usando los mismos datos que ya hemos usado en DADA2. Una diferencia clave entre DADA2 y mothur es que este último no está implementado en R sino que en la terminal directamente, por lo que para comenzar tenemos que abrir la terminal, cambiar directorio hasta donde están las reads.

• En mothur partimos uniendo los reads (a diferencia con DADA2) para formar contigs.

make.contigs(file=stability.files,processors=8)

• El comando summary.seqs en mothur es muy usado porque es el comando que te permite evaluar el efecto de otros comandos. Una vez que tenemos listos los contigs (stability.trim.contigs.fasta), aplicamos el comando:

summary.seqs(fasta=stability.trim.contigs.fasta)

• Un control de calidad básico es asegurarse que el largo de las reads sea el esperado y que éstas no contengan bases indeterminadas (Ns):

screen.seqs(fasta=stability.trim.contigs.fasta,group=stability.contigs.groups,maxambig=0,minlength=100,maxlength=300)

• En el paso anterior, nos aseguramos que las reads no tengan ambiguedades y que su largo vaya de 100 a 300 bases. Ahora le pedimos a mothur que extraiga las reads únicas, que las cuente, y que nos dé un resumen de esas operaciones:

unique.seqs(fasta=stability.trim.contigs.good.fasta)
count.seqs(name=stability.trim.contigs.good.names,group=stability.contigs.good.groups)
summary.seqs(count=stability.trim.contigs.good.count_table)

• Ahora que tenemos un conjunto de reads únicas podemos proceder al alineamiento en contra de una base de datos. Aquí otra vez usamos SILVA:

align.seqs(fasta=stability.trim.contigs.good.unique.fasta,reference=silva.nr_v123.v4.align.txt,flip=t)
summary.seqs(fasta=stability.trim.contigs.good.unique.align,count=stability.trim.contigs.good.count_table)

• El último comando nos muestra de dónde hasta dónde alinean nuestras secuencias en contra de la referencia de SILVA y su distribución. Con esta información, y para minimizar errores de alineamiento, podemos ajustar los siguientes comandos para cortar (start=13862,end=23444) y seguir usando solamente las reads que alinean de manera más congruente contra SILVA.

screen.seqs(fasta=stability.trim.contigs.good.unique.align,count=stability.trim.contigs.good.count_table,summary=stability.trim.contigs.good.unique.summary,start=13862,end=23444,maxhomop=8)
summary.seqs(fasta=current,count=current)
filter.seqs(fasta=stability.trim.contigs.good.unique.good.align,vertical=T,trump=.)
unique.seqs(fasta=stability.trim.contigs.good.unique.good.filter.fasta,count=stability.trim.contigs.good.good.count_table)
pre.cluster(fasta=stability.trim.contigs.good.unique.good.filter.unique.fasta,count=stability.trim.contigs.good.unique.good.filter.count_table,diffs=2)

Los últimos tres comandos filtran las secuecnias que alinean fuera del rango establecido, seleccionan nuevamente las reads únicas, y finalmente realizan un paso de pre agrupación o pre clustering que tolera hasta dos diferencias nucleotídicas en las secuencias.

• Al igual que con DADA2, necesitamos eliminar las secuencias que pueden ser quiméricas, es decir, producto de recombinación introducida durante el paso de PCR.

chimera.uchime(fasta=stability.trim.contigs.good.unique.good.filter.unique.precluster.fasta,count=stability.trim.contigs.good.unique.good.filter.unique.precluster.count_table,dereplicate=t)
remove.seqs(fasta=stability.trim.contigs.good.unique.good.filter.unique.precluster.fasta,accnos=stability.trim.contigs.good.unique.good.filter.unique.precluster.denovo.uchime.accnos)
summary.seqs(fasta=current,count=current)

Ahora tenemos un conjunto de secuencias filtradas de calidad suficiente para el paso crucial de clasificación.

• Usando el comando classify.seqs podemos alinear las secuencias en contra de la base de datos de SILVA a un cutoff determinado. Posteriormente, removemos alineamientos espurios que pueden corresponder a falsos positivos o, simplemente, a que los partidores universales de 16S eran también complementarios a ortólogos del 16S en plantas o eucariontes en general.

classify.seqs(fasta=stability.trim.contigs.good.unique.good.filter.unique.precluster.pick.fasta,count=stability.trim.contigs.good.unique.good.filter.unique.precluster.denovo.uchime.pick.count_table,reference=trainset9_032012.pds.fasta,taxonomy=trainset9_032012.pds.tax,cutoff=80)
remove.lineage(fasta=stability.trim.contigs.good.unique.good.filter.unique.precluster.pick.fasta,count=stability.trim.contigs.good.unique.good.filter.unique.precluster.denovo.uchime.pick.count_table,taxonomy=stability.trim.contigs.good.unique.good.filter.unique.precluster.pick.pds.wang.taxonomy,taxon=Chloroplast-Mitochondria-unknown-Eukaryota)

• Finalmente, calculamos una matriz de distancia entre las secuencias, las agrupamos y propagamos la clasificación taxonómica de tal forma de contar con una tabla de OTUs y una tabla de taxonomía.

dist.seqs(fasta=stability.trim.contigs.good.unique.good.filter.unique.precluster.pick.pick.fasta,cutoff=0.20)
cluster(column=stability.trim.contigs.good.unique.good.filter.unique.precluster.pick.pick.dist,count=stability.trim.contigs.good.unique.good.filter.unique.precluster.denovo.uchime.pick.pick.count_table,method=opti)
make.shared(list=stability.trim.contigs.good.unique.good.filter.unique.precluster.pick.pick.opti_mcc.list,count=stability.trim.contigs.good.unique.good.filter.unique.precluster.denovo.uchime.pick.pick.count_table,label=0.03)
classify.otu(list=stability.trim.contigs.good.unique.good.filter.unique.precluster.pick.pick.opti_mcc.list,count=stability.trim.contigs.good.unique.good.filter.unique.precluster.denovo.uchime.pick.pick.count_table,taxonomy=stability.trim.contigs.good.unique.good.filter.unique.precluster.pick.pds.wang.pick.taxonomy,label=0.03)

• Escribimos el output para su posterior uso en R en formato BIOM. El formato BIOM es el industry standard en el mundo metagenómico, reconocido por el Earth Microbiome Project y la organización Standards in Genomic Sciences.

make.biom(shared=stability.trim.contigs.good.unique.good.filter.unique.precluster.pick.pick.opti_mcc.shared,constaxonomy=stability.trim.contigs.good.unique.good.filter.unique.precluster.pick.pick.opti_mcc.0.03.cons.taxonomy)