Entrega Final

Introducción

En este notebook presenta el procedimiento llevado a cabo para el análisis de datos de metabarcoding, a partir el artículo “The sources and transmission routes of microbial populations throughout a meat processing facility”. Este artículo busca analizar e identificar fuentes potenciales de contaminación microbiana en la carne de cerdos, durante el proceso de faena en el matadero.

Objetivo

El objetivo de nuestro estudio se reduce a analizar cómo cambia la composición microbiana en la carne de los cerdos de los distintos granjeros, e identificar si alguna composición presenta un peligro mayor para el matadero.

1. Preprocesamiento de los datos

Para llevar a cabo el análisis, primero se deben cargar, filtrar y preprocesar los datos tomados en el matadero, los cuales fueron descargados del European Nucleotide Archive.

1.1 Setup del ambiente de trabajo

Se hará uso de paquetes de R, y librerías opensource de Bioconductor para el procesamiento de datos.

library(ampvis2)
library(mia)
library(miaViz)
library(ggplot2)
library(Biostrings)
library(phyloseq)
library(ANCOMBC)
library(tidyr)
library(patchwork)
library(ggsignif)
library(scater)
library(ggpubr)
library(reshape2)
library(RColorBrewer)
library(dada2)
library(vegan)

Además, se debe setear el directorio que contenga los archivos (.fastq) con los datos descargados:

# Reemplazar por directorio utilizado
setwd("~/UCU/alimentos1/")
path <- getwd()

1.2 Registrar metainformación de las muestras

Diferentes muestras fueron tomadas durante el proceso de faena, en diferentes cerdos, lugares y herramientas. Las lecturas realizadas en cada muestra fueron almacenadas en archivos .fastq, los cuales contienen las secuencias de ADN y calidad de sus respectivas lecturas.

Debido a que los archivos contienen demasiada información, nunca se cargan en memoria. Sino que las funciones utilizadas se encargan de buscar los datos en cada archivo cuando sea necesario, a partir de su ruta. Entonces, se guarda el listado de rutas hacia los archivos de todas las lecturas forward (f) y reverse (r).

# Archivos que contengan '_1.fastq' representan lecturas f
fnFs <- sort(list.files(path, pattern="_1.fastq", full.names = TRUE)) 
# Archivos que contengan '_2.fastq' representan lecturas r
fnRs <- sort(list.files(path, pattern="_2.fastq", full.names = TRUE))

Luego, se extraen los nombres de cada muestra a partir del nombre de los archivos:

# Formato: SAMPLENAME_XXX.fastq
sample.names <- sapply(strsplit(basename(fnFs), "_"), `[`, 1)

Finalmente se importa la metadata de cada muestra:

sample.metadata <- read.csv("sampleData_alimentos.csv")

1.3 Chequeo y filtrado de calidad

Para filtrar lecturas de mala calidad, se grafica la calidad de las mismas para definir donde truncarlas.

Es característico de una lectura Illumina que al final de la secuencia, el score de calidad decaiga, y que la lectura reverse tenga menor calidad que la forward. Entonces se trunca las regiones de mala calidad.

En este caso, se utilizarán las regiones donde el score sea mayor a 30 aproximadamente. Por lo tanto se truncaran las bases posteriores a 290 en lecturas forward, y 200 para reverse.

Para realizar el filtrado y truncado es útil la función filterAndTrim, la cual permite:

• Cortar registros de mala calidad

• Filtrar lecturas con datos nulos

• Cortar primers

Utilizando los parámetros:

• truncLen=c(290, 200) ⟶ Lecturas con tamaño menor a 290 (forward) o 200 (reverse) son descartadas

• maxEE=(2,2) ⟶ Lecturas con mas de 2 errores esperados son descartadas

• maxN=0 ⟶ Lecturas con datos nulos son descartadas

Se debe tener cuidado de no truncar en exceso, debido a que dada2 solo tiene en cuenta la secuencia solapada entre forward y reverse, entonces al truncar de más habrá poco solapamiento, y menos datos serán procesados por dada2.

# Archivos donde se almacenarán las muestras filtradas y truncadas
filtFs <- file.path(path, "filtered", paste0(sample.names, "_F_filt.fastq.gz")) 
filtRs <- file.path(path, "filtered", paste0(sample.names, "_R_filt.fastq.gz"))

# Filtrado y truncado de datos
filterAndTrim_out <- filterAndTrim(fnFs, filtFs, fnRs, filtRs, truncLen=c(290,200), maxN=0, maxEE=c(2,2), truncQ=2, rm.phix=TRUE, compress=TRUE, multithread=TRUE)

# Almacenar resultado
save(filterAndTrim_out, file="filterAndTrim_out")

Al ser una operación costosa, se almacena el resultado en un archivo para evitar recalcular:

load("filterAndTrim_out")

El resultado de esta función es una matriz que contiene para cada muestra, la cantidad de secuencias que existían originalmente, y cuántas resultaron luego del filtrado.

Permite identificar cuantas lecturas fueron descartadas. Además, almacena en los archivos especificados (filtFs y filtRs) los resultados filtrados.

col_sum <- colSums(filterAndTrim_out)
filtered_reads <- col_sum[1] - col_sum[2]

message(sprintf("Fueron filtradas %s lecturas, un %%%.2f de las lecturas originales.", filtered_reads, filtered_reads * 100 / col_sum[1]))

Fueron filtradas 456939 lecturas, un %20.61 de las lecturas originales.

1.4 Aprendizaje de errores

Dada2 necesita conocer las probabilidades de errores para poder decidir si dos secuencias que se diferencian en 1 o 2 bases son la misma y haya ocurrido un error, o sean realmente dos secuencias diferentes.

Por lo tanto, se debe estimar los errores, y modelar las probabilidades de que hayan ocurrido transiciones (A ⟶ T, A ⟶ C, C ⟶ A, etc):

# forward
errF <- learnErrors(filtFs, multithread=TRUE)
save(errF, file="errF")

# reverse
errR <- learnErrors(filtRs, multithread=TRUE)
save(errR, file="errR")

Al ser una operación costosa, se almacena el resultado en un archivo para evitar recalcular:

load("errF")
load("errR")

Se puede observar que los errores son acordes a sus estimados, ya que presentan gráficos similares:

1.5 Derreplicacion

La derreplicacion consiste en contar las apariciones de cada secuencia en las muestras, lo cual es necesario para luego agrupar los ASVs encontrados:

# forward
derepFs <- derepFastq(filtFs, verbose=TRUE)
names(derepFs) <- sample.names
save(derepFs, file= "derepFs")

# reverse
derepRs <- derepFastq(filtRs, verbose=TRUE)
names(derepRs) <- sample.names
save(derepRs, file= "derepRs")

Para evitar el costoso cálculo:

load("derepFs")
load("derepRs")

Obteniéndose por muestra, una matriz con el conteo de cada secuencia, y sus scores promedio de calidad.

1.6 Identificar ASVs con DADA2

Ahora que se tienen los datos prontos para agrupar, se identifican los ASVs para las secuencias forward y reverse utilizando la función dada:

dadaFs <- dada(derepFs, err=errF, pool = TRUE , multithread=TRUE)
save(dadaFs, file="dadaFs")

#Reverse reads
dadaRs <- dada(derepRs, err=errR, pool = TRUE , multithread=TRUE)
save(dadaRs, file="dadaRs")

Para evitar el costoso cálculo:

1.7 Identificar cuantas veces aparece un ASV en cada muestra

Para esto, se debe unir los ASVs forward y reverse en una única secuencia:

mergers <- mergePairs(dadaFs, derepFs, dadaRs, derepRs, verbose=TRUE) 
save(mergers, file="mergers")

Para evitar el costoso cálculo:

El resultado de esta función es una lista de dataframes con información sobre la unión de las secuencias forward y reverse, por ejemplo: solapamiento, missmatch, etc.

Esta información es utilizada para contar las apariciones de cada ASV en cada muestra:

1.8 Filtrado de ASVs resultantes y quimeras

Las secuencias se esperan que tengan un largo dentro de un intervalo. Este intervalo se calcula a partir de los primers, largo de forward y reverse, porque se espera que las secuencias tengan un largo similar al solapamiento del forward y reverse.

En este caso el intervalo es entre 440 y 466, por lo tanto se filtran todos los ASV que no lo cumplan:

Además, es posible que existan quimeras en los datos. Las quimeras son secuencias malformadas y no existentes que pueden aparecer luego de la secuenciación, y deben ser filtradas:

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

Para evitar el costoso cálculo:

1.9 Resumen de información perdida durante el proceso

Es posible obtener la cantidad de ASVs forward y reverse (denoised), cantidad de secuencias resultantes del merge, y cantidad de ASVs resultantes luego de eliminar largos no esperados y que contengan quimeras.

Estos cálculos permiten identificar la información recortada durante el preprocesamiento de los datos:

# Funcion auxiliar para contar los valores unicos de una colección
getN <- function(x) sum(getUniques(x))

track <- cbind(filterAndTrim_out, sapply(dadaFs, getN), sapply(dadaRs, getN), sapply(mergers, getN), rowSums(seqtab.nochim))
colnames(track) <- c("input", "filtered", "denoisedF", "denoisedR", "merged", "nonchim") 
rownames(track) <- sample.names

# Mostrar resultados
head(track)

           input filtered denoisedF denoisedR merged nonchim
ERR4144822 51293    44055     43140     43110  39762   27904
ERR4144823 37976    34005     33267     33215  30645   21949
ERR4144824 31514    27413     26851     26763  24647   17679
ERR4144825 31934    28665     28023     28071  26052   19236
ERR4144826 26243    21528     21016     21058  19637   13124
ERR4144827 37867    33094     32246     32539  30534   20689

total_nonchim <- sum(track[, 'nonchim'])
total_input <-sum(track[, 'input'])
message(sprintf("Se perdieron un total de %s muestras, un %%%.2f",  total_input - total_nonchim, sum(total_nonchim * 100 / total_input)))

Se perdieron un total de 1040460 muestras, un %53.06

1.10 Unificar datos con Phyloseq

Phyloseq es un framework que ofrece funcionalidades para analizar los datos, por ejemplo diversidad alfa, beta, abundancias diferenciales, etc, en base a toda la información de las muestras que se hayan recolectado.

Se debe cargar la base de datos que contiene la taxonomía (informacion biologica) de los ASVs:

taxa_gtdb <- assignTaxonomy(seqtab.nochim, "", multithread=TRUE)

Para evitar el costoso cálculo:

load("taxa_gtdb")

Luego, importar las muestras que se vayan a utilizar:

sample_data <- as.data.frame(read.table("sampleData_alimentos.csv", header=TRUE, sep=","))
rownames(sample_data) <- sample_data[,1]

Finalmente, cargar los datos a Phyloseq:

physeq <- phyloseq(otu_table(seqtab.nochim, taxa_are_rows=FALSE),sample_data(sample_data),tax_table(taxa_gtdb)) 
save(physeq, file="physeq")

2. Interpretación de resultados

2.1 Análisis de muestras

Para el objetivo del estudio no es necesario analizar todas las muestras tomadas, es posible tomar un subset útil para la investigación. En este caso, son las muestras tomadas inicialmente de la piel de los cerdos, y en el recto anal.

ps <- subset_samples(physeq, Sampling.position %in% c("Sticking", "Anal swab"))
ps.prop <- transform_sample_counts(ps, function(OTU) OTU/sum(OTU))

message("Cantidad de muestras: ", nsamples(ps), "\n",
        "Cantidad de ASVs: ", ntaxa(ps), "\n",
        "Cantidad de ASVs no presentes en las muestras: ", ntaxa(prune_taxa(taxa_sums(ps) == 0, ps)))

Cantidad de muestras: 8
Cantidad de ASVs: 15413
Cantidad de ASVs no presentes en las muestras: 7606

Se filtran los ASV que no estan presentes en ninguna muestra, porque no son utiles:

ps <- prune_taxa(taxa_sums(ps) > 0, ps)

La distribución de ASV resultantes por muestra se puede apreciar en el siguiente gráfico:

reads_per_ASV <- data.frame(nreads=sort(taxa_sums(ps)), TRUE, SVs= 1:ntaxa(ps))

ggplot(data=reads_per_ASV,aes(x = SVs, y = nreads)) +
  geom_area (stat="identity", fill = "#69BE28") +
  scale_x_log10()

Como del granjero A se muestrearon dos cerdos, existen cerca del doble de lecturas que los demás granjeros.

# Filter reads
filtered_reads <- data.frame(
    nreads=sort(sample_sums(ps)), 
    Farmer=sample_data(ps)$Farmer,
    Position=sample_data(ps)$Sampling.position, 
    TRUE, 
    samples=sample_data(ps)$Farmer, 
    type = "Samples"
)

# Plot
ggplot(filtered_reads, aes(x=samples, y=nreads)) +
  geom_bar(aes(color=Position, fill=Position), stat="identity") + 
  facet_wrap(~Farmer, scales="free_x")

2.2 Abundancias relativas

# Graficar abundancias relativas por filo
plot_bar(ps.prop, x="Sample", fill="Phylum") +
  facet_wrap(~Farmer, scales="free_x") +
  geom_bar(aes(color=Phylum, fill=Phylum), stat="identity", position="stack")

En el gráfico anterior se puede ver en todas las muestras una gran presencia de bacteroidetes, las mismas son comunes en el intestino de muchos animales. Estas bacterias pueden desempeñar un papel importante en la digestión y la absorción de nutrientes en los animales, así como en la salud general de su sistema digestivo.

Además, se puede ver presencia de fusobacteriotas en las muestras. Son conocidas por su capacidad para causar infecciones en los animales, incluyendo infecciones orales, abscesos y enfermedades respiratorias. Por ejemplo, en los cerdos se ha informado que las Fusobacterium pueden estar implicadas en casos de necrosis hepática y abscesos. La infección con Fusobacterium puede ser el resultado de una lesión en la piel o una infección previa, lo que permite que las bacterias entren en el cuerpo y causen una infección.

Otras bacterias presentes en todas las muestras son del filo Spirochaetota, las mismas son conocidas por ser patógenas pero es común en el tracto gastrointestinal de los cerdos, aunque no se consideran patógenas en este caso para ellos.

# Abundancias relativas por género
tse <- makeTreeSummarizedExperimentFromPhyloseq(ps) 
tse <- transformSamples(tse, method = "relabundance")
tse <- addTaxonomyTree(tse)
tse_gen <- agglomerateByRank(tse, rank = "Genus")
rownames(tse_gen)<-sub("Genus:","", rownames(tse_gen))

plotAbundanceDensity(tse_gen, layout = "jitter", abund_values="relabundance", n = 20, point_size=1, point_shape=19, point_alpha=0.5,                                               colour_by="Farmer") + 
  scale_x_log10(label=scales::percent)

En este otro gráfico, se puede ver los porcentajes de abundancia relativa entre cada género.

Como por ejemplo, campylobacter está presente en un más de 10% de abundancia relativa en una muestra del cerdo perteneciente al granjero B. Campylobacter puede ser un patógeno que causa enfermedades gastrointestinales en humanos, por lo que la presencia de Campylobacter en cerdos puede plantear preocupaciones en términos de seguridad alimentaria. La bacteria puede transmitirse a los humanos a través del consumo de alimentos contaminados con Campylobacter, como carne cruda o insuficientemente cocida, así como a través del contacto directo con animales infectados.

# Abundancia relativa por familia
tse_family <- agglomerateByRank(tse, rank = "Family")
rownames(tse_family)<-sub("Family:","", rownames(tse_family))

plotAbundanceDensity(tse_family, layout = "jitter", abund_values="relabundance", n = 20, point_size=1, point_shape=19, point_alpha=0.5,                                               colour_by="Farmer") + scale_x_log10(label=scales::percent)

Las bacterias pertenecientes a la familia Moraxellaceae se encuentran presente en todas las muestras pero mayormente en los chanchos del granjero A y C. Moraxellaceae pueden causar problemas en la calidad de la carne de cerdo, como el olor rancio, lo que puede ser un problema importante para la industria cárnica.

Las bacterias enterobacteriaceae son patógenas y pueden causar enfermedades como la salmonella. La misma está más presente en las muestras del chancho del granjero B.

2.3 Diversidad alfa

La diversidad alfa en metabarcoding se refiere a la diversidad de especies o tipos de organismos presentes en una muestra de ADN biológico y se puede medir utilizando varios índices, como el índice de Simpson.

rarefied<-rarefy_even_depth(ps,min(sample_sums(ps)),rngseed=1)
tse_rarefied <- makeTreeSummarizedExperimentFromPhyloseq(rarefied)
# Uniformidad
tse_rarefied <- estimateEvenness(tse_rarefied, 
                        abund_values = "counts",
                        index="simpson")

plots <- lapply(c("simpson"),
                plotColData,
                object = tse_rarefied,
                x = "Farmer",
                colour_by = "Farmer")

plots <- lapply(plots, "+", 
                theme(axis.text.x = element_blank(),
                      axis.title.x = element_blank(),
                      axis.ticks.x = element_blank()))

plots[[1]] + plot_layout(guides = "collect")

La gráfica de diversidad alfa de Simpson muestra la riqueza de especies y la equitabilidad en una comunidad. En esta gráfica, la diversidad alfa se representa en el eje de las ordenadas, mientras que el tamaño de la muestra se representa en el eje de abscisas.

La diversidad de Simpson se expresa como un índice de dominancia y se calcula restando la suma de los cuadrados de las proporciones de cada especie en la comunidad de uno. Cuanto más cercano sea el índice de Simpson a 1, menor será la diversidad alfa, lo que indica una mayor dominancia de una o unas pocas especies en la comunidad. Por otro lado, cuanto más cercano sea el índice de Simpson a 0, mayor será la diversidad alfa, lo que indica que las especies están distribuidas de manera uniforme en la comunidad.

message("Indice simpson promedio del granjero B: ", mean(subset(x, Farmer == "B")$simpson))

Indice simpson promedio del granjero B: 0.200961048738534

message("Indice simpson promedio del granjero A: ", mean(subset(x, Farmer == "A")$simpson))

Indice simpson promedio del granjero A: 0.402348623504202

En este caso, las muestras del granjero B presentan un menor índice simpson promedio, por lo tanto cuenta con una mayor diversidad alfa.

Por otro lado, las muestras del granjero A presentan un mayor índice simpson que las demás muestras, por lo tanto cuentan con una diversidad alfa menor. Esto puede significar que existe menor cantidad de bacterias diferentes, o una gran proporción de las bacterias presentes son similares.

2.4 Diversidad beta

En metabarcoding, la diversidad beta se refiere a la variación en la composición de especies entre diferentes muestras.

# Aglomeracion por gènero manteniendo los NA
tse_genus <- agglomerateByRank(tse, rank = "Genus")
tse_genus_notNA <- tse_genus[!is.na(rowData(tse_genus)$Genus), ]
tse_genus_notNA <- runMDS(tse_genus_notNA, FUN = vegan::vegdist, method = "bray", name = "PCoA_BC", exprs_values = "counts")

# Create ggplot object
p <- plotReducedDim(tse_genus_notNA, "PCoA_BC", colour_by = "Farmer")

# Add explained variance for each axis
e <- attr(reducedDim(tse_genus_notNA, "PCoA_BC"), "eig");
rel_eig <- e/sum(e[e>0])          
p + labs(x = paste("PCoA 1 (", round(100 * rel_eig[[1]],1), "%", ")", sep = ""),
              y = paste("PCoA 2 (", round(100 * rel_eig[[2]],1), "%", ")", sep = ""))

Los puntos en la gráfica representan las muestras y su posición se basa en su similitud biológica. Las muestras que están más cerca en la gráfica se consideran más similares en términos de su composición biológica.

Analizando la gráfica se puede concluir que las muestras coinciden con lo esperado ya que entre ellas hay cierta distancia lo que refiere a que las muestras no tienen composición biológica similar.

3. Conclusión

Luego de analizar las muestras bacterianas de los cerdos de los tres granjeros, se encontró que el cerdo del granjero B tenía la mayor cantidad de Campylobacter y Chloroflexi en su tracto gastrointestinal. Esto puede ser preocupante ya que ambas bacterias se sabe que pueden causar enfermedades en los seres humanos.

La presencia de Campylobacter en los cerdos es particularmente preocupante ya que se ha demostrado que esta bacteria puede ser un importante factor de riesgo para la transmisión de enfermedades alimentarias. Además, la resistencia a los antimicrobianos.

La presencia de Chloroflexi en los cerdos también puede ser preocupante ya que se sabe que algunas especies pueden causar enfermedades respiratorias en los seres humanos. Además, se sabe que algunas bacterias Chloroflexi son resistentes a los antimicrobianos, lo que puede complicar el tratamiento de las enfermedades causadas por estas bacterias.

Por otro lado, se puede ver que la carne de los granjeros A y C tienen una mayor cantidad de Moraxellaceae que la carne del granjero B. La misma presenta un riesgo en términos de desperdicio alimenticio pero no de riesgo humano.

Contestando la pregunta disparadora de la investigación se podría decir que la carne del granjero B es la que mayor peligrosidad tiene.