Metagenómica por Amplicones

### --- Definimos los paquetes a instalar ——————###
## Desde el repositorio CRAN ##
cran_packages <- c("bookdown", "knitr", "tidyverse", "plyr", "grid", "gridExtra", "kableExtra", "xtable", "ggpubr")

## Desde el repositorio de Bioconductor ##
bioc_packages <- c("phyloseq", "dada2", "DECIPHER", "phangorn", "ggpubr", "BiocManager", "DESeq2", "microbiome", "philr")

## Desde el repositorio GitHub ##
git_source <- c("twbattaglia/btools", "gmteunisse/Fantaxtic", "MadsAlbertsen/ampvis2", "opisthokonta/tsnemicrobiota")
git_packages <- c("btools", "fantaxtic", "ampvis2", "tsnemicrobiota") 

### --- Instalación de paquetes ——————###
## Instalación de paquetes desde CRAN ##
.inst_cran <- cran_packages %in% installed.packages()
if(any(!.inst_cran)) {
  install.packages(cran_packages[!.inst_cran])
}

## Instalación de paquetes desde Bioconductor ##
if (!requireNamespace("BiocManager", quietly = TRUE)) {
  install.packages("BiocManager")
}
.inst_bioc <- bioc_packages %in% installed.packages()
if(any(!.inst_bioc)) {
  BiocManager::install(bioc_packages[!.inst_bioc])
}

## Instalación de paquetes desde GitHub ##
if (!requireNamespace("devtools", quietly = TRUE)) {
  install.packages("devtools")
}
.inst_git <- git_packages %in% installed.packages()
if(any(!.inst_git)) {
  devtools::install_github(git_source[!.inst_git])
}

stringi (1.8.3 -> 1.8.4) [CRAN]
xfun    (0.43  -> 0.44 ) [CRAN]
package 'stringi' successfully unpacked and MD5 sums checked
package 'xfun' successfully unpacked and MD5 sums checked

The downloaded binary packages are in
    C:\Users\Mario\AppData\Local\Temp\Rtmp0iCL0s\downloaded_packages
── R CMD build ─────────────────────────────────────────────────────────────────
* checking for file 'C:\Users\Mario\AppData\Local\Temp\Rtmp0iCL0s\remotes23aec71af5673\twbattaglia-btools-fa21d4c/DESCRIPTION' ... OK
* preparing 'btools':
* checking DESCRIPTION meta-information ... OK
* checking for LF line-endings in source and make files and shell scripts
* checking for empty or unneeded directories
Omitted 'LazyData' from DESCRIPTION
* building 'btools_0.0.1.tar.gz'

### --- Cargar los paquetes ——————###
sapply(c(cran_packages, bioc_packages, git_packages), require, character.only = TRUE)

      bookdown          knitr      tidyverse           plyr           grid 
          TRUE           TRUE           TRUE           TRUE           TRUE 
     gridExtra     kableExtra         xtable         ggpubr       phyloseq 
          TRUE           TRUE           TRUE           TRUE           TRUE 
         dada2       DECIPHER       phangorn         ggpubr    BiocManager 
          TRUE           TRUE           TRUE           TRUE           TRUE 
        DESeq2     microbiome          philr         btools      fantaxtic 
          TRUE           TRUE           TRUE          FALSE           TRUE 
       ampvis2 tsnemicrobiota 
          TRUE           TRUE 

library(reshape2)

1 Introducción

La metagenómica por amplicones es una técnica poderosa y ampliamente utilizada para estudiar la diversidad microbiana en diferentes ambientes. Este enfoque se basa en la secuenciación de regiones específicas del ADN ribosómico, como el gen 16S rRNA en bacterias, para identificar y clasificar organismos presentes en una muestra. La capacidad de obtener perfiles taxonómicos detallados y precisos permite investigar las comunidades microbianas en una variedad de contextos, desde la microbiota humana hasta los ecosistemas ambientales.

En este estudio, se ha llevado a cabo un análisis metagenómico de muestras de piel de ballenas recolectadas en diferentes ubicaciones geográficas de Chile en 2017. Utilizando la plataforma de secuenciación Illumina y el pipeline DADA2, se han procesado y analizado los datos de secuenciación para evaluar la diversidad microbiana en estas muestras. Los objetivos principales incluyen la evaluación de la calidad de las lecturas de secuenciación, la identificación de variantes de secuencia de amplicones (ASVs), la asignación taxonómica de las secuencias, y el análisis de la diversidad alfa y beta.

2 Directorio de trabajo

Se establece el directorio de trabajo donde se encuentra el conjunto de datos de secuenciación. Luego, se define la carpeta específica que contiene estos datos. Se listan los ficheros presentes en dicha carpeta para verificar su contenido. A continuación, se cargan en objetos separados (fnFs y fnRs) las rutas de los ficheros correspondientes a las lecturas forward (R1) y reverse (R2), respectivamente. Finalmente, se extraen y almacenan los nombres de las muestras en un objeto (sample.names), utilizando la primera parte del nombre de los ficheros como identificadores.

### --- Directorio de trabajo y datos ——————###
## Directorio de trabajo ##
setwd("C:/Users/Mario/Documents/Proyectos R/TecModAlgpractica8")

## Definimos la carpeta donde se encuentran los datos de secuenciación ##
path <- "data_set"

## Listamos los ficheros ##
list.files(path)

 [1] "filtered"           "SRR6442721_1.fastq" "SRR6442721_2.fastq"
 [4] "SRR6442722_1.fastq" "SRR6442722_2.fastq" "SRR6442723_1.fastq"
 [7] "SRR6442723_2.fastq" "SRR6442724_1.fastq" "SRR6442724_2.fastq"
[10] "SRR6442733_1.fastq" "SRR6442733_2.fastq" "SRR6442734_1.fastq"
[13] "SRR6442734_2.fastq" "SRR6442735_1.fastq" "SRR6442735_2.fastq"
[16] "SRR6442736_1.fastq" "SRR6442736_2.fastq"

## Cargamos en sendos objetos el path de los ficheros de la R1 y de la R2 ##
fnFs <- sort(list.files(path, pattern="_1.fastq", full.names = TRUE))
fnRs <- sort(list.files(path, pattern="_2.fastq", full.names = TRUE))

## Extraemos en un objeto los nombres de las muestras ##
sample.names <- sapply(strsplit(basename(fnFs), "_"), `[`, 1)

2.1 Preguntas

1. ¿Qué muestran los objetos fnFs y fnRs?

fnFs

[1] "data_set/SRR6442721_1.fastq" "data_set/SRR6442722_1.fastq"
[3] "data_set/SRR6442723_1.fastq" "data_set/SRR6442724_1.fastq"
[5] "data_set/SRR6442733_1.fastq" "data_set/SRR6442734_1.fastq"
[7] "data_set/SRR6442735_1.fastq" "data_set/SRR6442736_1.fastq"

fnRs

[1] "data_set/SRR6442721_2.fastq" "data_set/SRR6442722_2.fastq"
[3] "data_set/SRR6442723_2.fastq" "data_set/SRR6442724_2.fastq"
[5] "data_set/SRR6442733_2.fastq" "data_set/SRR6442734_2.fastq"
[7] "data_set/SRR6442735_2.fastq" "data_set/SRR6442736_2.fastq"

Se muestran los objetos fnFs y fnRs, los cuales contienen las rutas completas de los archivos de secuenciación en formato FASTQ para las lecturas forward (R1) y reverse (R2), respectivamente. Cada entrada en estos objetos representa un archivo correspondiente a una muestra específica, ordenados alfabéticamente.

2. ¿Qué muestra el objeto sample.names?

sample.names

[1] "SRR6442721" "SRR6442722" "SRR6442723" "SRR6442724" "SRR6442733"
[6] "SRR6442734" "SRR6442735" "SRR6442736"

Se muestra el objeto sample.names, que contiene los nombres de las muestras extraídos de los archivos de secuenciación. Estos nombres corresponden a la primera parte de los nombres de los archivos antes del primer guion bajo, y se utilizan para identificar de manera única cada muestra en los análisis subsecuentes.

3 Calidad

3.1 Inspección de Calidad de los Datos

3.1.1 Gráfica 1.1: Calidad de SRR6442721_1.fastq y SRR6442722_1.fastq

plotQualityProfile(fnFs[1:2])

Se muestra la calidad de las lecturas de secuenciación para los archivos SRR6442721_1.fastq y SRR6442722_1.fastq. En ambos gráficos, el eje Y representa la puntuación de calidad, mientras que el eje X representa el ciclo de secuenciación (posición a lo largo de la lectura). La línea roja horizontal parece indicar el umbral de calidad mínima.

En general, se observa que la calidad de las lecturas es alta en las primeras posiciones (con puntuaciones de calidad cercanas a 40) y tiende a disminuir ligeramente hacia el final de las lecturas, lo cual, como se ha visto a lo largo de las prácticas de la asignautra, es un patrón común en datos de secuenciación. La dispersión de los puntos muestra la variabilidad de la calidad a lo largo de los ciclos. Para SRR6442721_1.fastq, se registraron 51868 lecturas, mientras que para SRR6442722_1.fastq se registraron 70315 lecturas. Ambas muestras presentan una buena calidad general de las lecturas.

3.1.2 Gráfica 1.2: Calidad de SRR6442721_2.fastq y SRR6442722_2.fastq

plotQualityProfile(fnRs[1:2])

Se muestra la calidad de las lecturas de secuenciación para los archivos SRR6442721_2.fastq y SRR6442722_2.fastq. Parecido al fichero anterior, se observa que la calidad de las lecturas es alta en las primeras posiciones (con puntuaciones de calidad cercanas a 40) y tiende a disminuir hacia el final de las lecturas.

Para SRR6442721_2.fastq, se registraron 51868 lecturas, mientras que para SRR6442722_2.fastq se registraron 70315 lecturas. Ambas muestras presentan una buena calidad general de las lecturas en las primeras posiciones, aunque se observa una mayor variabilidad y una disminución significativa de la calidad hacia el final de las lecturas, especialmente en SRR6442722_2.fastq, donde hay una caída abrupta de la calidad después del ciclo 200. Esto indica que puede ser necesario aplicar un filtrado adicional para remover las lecturas de menor calidad antes de proceder con los análisis subsecuentes.

3.1.3 Tabla 1: Resumen de Calidad

Archivo	Número de Lecturas	Ciclo (Rango)	Calidad Promedio (Inicial)	Calidad Promedio (Final)
SRR6442721_1.fastq	51868	0-250	38	30
SRR6442722_1.fastq	70315	0-250	38	25
SRR6442721_2.fastq	51868	0-250	38	28
SRR6442722_2.fastq	70315	0-250	38	20

• Archivo: Nombre del archivo de secuenciación.
• Número de Lecturas: Número total de lecturas en cada archivo.
• Ciclo (Rango): Rango de ciclos de secuenciación (posiciones) representados en el gráfico de calidad.
• Calidad Promedio (Inicial): Calidad promedio de las lecturas al inicio de la secuenciación.
• Calidad Promedio (Final): Calidad promedio de las lecturas al final de la secuenciación.

3.2 Filtrado

Se realiza un control de calidad de las lecturas de secuenciación creando primero un directorio para almacenar las lecturas filtradas (“limpias”). Para cada muestra, se generan rutas para los archivos filtrados (filtFs y filtRs) y se asignan los nombres de las muestras a estos archivos. A continuación, se aplica el control de calidad utilizando la función filterAndTrim. Esta función filtra y recorta las lecturas basándose en varios parámetros:

• maxN = 0: No se permiten bases ‘N’ (desconocidas).
• truncQ = 2: Se recortan las lecturas en el primer lugar donde la calidad cae por debajo de 2.
• maxEE = c(2, 2): Se permite un máximo de 2 errores esperados por lectura (forward y reverse).
• rm.phix = TRUE: Se eliminan las secuencias PhiX, que se utilizan como control en muchas corridas de secuenciación.
• compress = TRUE: Se comprimen las lecturas filtradas.
• multithread = FALSE: Se desactiva el uso de múltiples hilos para el procesamiento paralelo.

El resultado del filtrado y recorte se almacena en el objeto out, que contiene información sobre el número de lecturas que pasaron el control de calidad para cada muestra.

### --- Control de calidad ——————###
## Creamos un directorio para poner las reads "limpias" ##
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

## Aplicamos el control de calidad ##
out <- filterAndTrim(fnFs, filtFs, fnRs, filtRs,
                     maxN = 0, truncQ = 2, maxEE = c(2, 2),
                     rm.phix = TRUE, compress = TRUE,
                     multithread = FALSE)
out

                   reads.in reads.out
SRR6442721_1.fastq    51868     44378
SRR6442722_1.fastq    70315     57126
SRR6442723_1.fastq    94944     81018
SRR6442724_1.fastq      118        60
SRR6442733_1.fastq     3406      2003
SRR6442734_1.fastq     2022      1533
SRR6442735_1.fastq    35029     29382
SRR6442736_1.fastq    78863     67349

Se realiza un control de calidad de las lecturas de secuenciación creando primero un directorio para almacenar las lecturas filtradas (“limpias”). Para cada muestra, se generan rutas para los archivos filtrados (filtFs y filtRs) y se asignan los nombres de las muestras a estos archivos. A continuación, se aplica el control de calidad utilizando la función filterAndTrim. Esta función filtra y recorta las lecturas basándose en varios parámetros:

• maxN = 0: No se permiten bases ‘N’ (desconocidas).
• truncQ = 2: Se recortan las lecturas en el primer lugar donde la calidad cae por debajo de 2.
• maxEE = c(2, 2): Se permite un máximo de 2 errores esperados por lectura (forward y reverse).
• rm.phix = TRUE: Se eliminan las secuencias PhiX, que se utilizan como control en muchas corridas de secuenciación.
• compress = TRUE: Se comprimen las lecturas filtradas.
• multithread = FALSE: Se desactiva el uso de múltiples hilos para el procesamiento paralelo.

3.2.1 Gráfica 2: Porcentaje de Lecturas que pasaron el Filtro de Calidad

# EXTRA: Mostrar los resultados en una gráfica. 
# Guardamos los datos manualmente en un data.frame
results <- data.frame(
  sample = c("SRR6442721_1.fastq", "SRR6442722_1.fastq", "SRR6442723_1.fastq", "SRR6442724_1.fastq",
             "SRR6442733_1.fastq", "SRR6442734_1.fastq", "SRR6442735_1.fastq", "SRR6442736_1.fastq"),
  reads.in = c(51868, 70315, 94944, 118, 3406, 2022, 35029, 78863),
  reads.out = c(44378, 57126, 81018, 60, 2003, 1533, 29382, 67349)
)

# Calcular el porcentaje de lecturas que pasaron el filtro y mostrar
results$percentage <- (results$reads.out / results$reads.in) * 100
ggplot(results, aes(x = sample, y = percentage)) +
  geom_bar(stat = "identity", fill = "steelblue") +
  geom_text(aes(label = sprintf("%.1f%%", percentage)), vjust = -0.3, size = 3.5) +
  labs(title = "Porcentaje de Lecturas que Pasaron el Filtro de Calidad",
       x = "Muestra",
       y = "Porcentaje de Lecturas Filtradas") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Los resultados de filterAndTrim muestran el número de lecturas que entraron (reads.in) y el número de lecturas que salieron (reads.out) después del proceso de filtrado y recorte para cada archivo de secuenciación. A continuación se presenta un resumen de los resultados:

3.2.2 Tabla 2: Lecturas Originales y Filtradas por Archivo

Archivo	Lecturas Originales	Lecturas Filtradas
SRR6442721_1.fastq	51868	44378
SRR6442722_1.fastq	70315	57126
SRR6442723_1.fastq	94944	81018
SRR6442724_1.fastq	118	60
SRR6442733_1.fastq	3406	2003
SRR6442734_1.fastq	2022	1533
SRR6442735_1.fastq	35029	29382
SRR6442736_1.fastq	78863	67349

Se observa que una proporción significativa de lecturas pasa el filtro en la mayoría de los archivos, indicando que la calidad de las lecturas es generalmente alta. Sin embargo, para SRR6442724_1.fastq, solo 60 de 118 lecturas pasaron el filtro, lo que sugiere una calidad baja para esta muestra.

3.3 Identificación de Errores

### --- Identificación de errores ——————###
## Aprender a identificar los errores ##
errF <- learnErrors(filtFs, multithread = TRUE)

70993893 total bases in 282849 reads from 8 samples will be used for learning the error rates.

errR <- learnErrors(filtRs, multithread = TRUE)

70914388 total bases in 282849 reads from 8 samples will be used for learning the error rates.

## Aplicamos el aprendizaje a las lecturas ##
dadaFs <- dada(filtFs, err = errF, multithread = TRUE)

Sample 1 - 44378 reads in 7265 unique sequences.
Sample 2 - 57126 reads in 9626 unique sequences.
Sample 3 - 81018 reads in 15490 unique sequences.
Sample 4 - 60 reads in 47 unique sequences.
Sample 5 - 2003 reads in 799 unique sequences.
Sample 6 - 1533 reads in 462 unique sequences.
Sample 7 - 29382 reads in 5938 unique sequences.
Sample 8 - 67349 reads in 9183 unique sequences.

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

Sample 1 - 44378 reads in 11286 unique sequences.
Sample 2 - 57126 reads in 17187 unique sequences.
Sample 3 - 81018 reads in 21192 unique sequences.
Sample 4 - 60 reads in 49 unique sequences.
Sample 5 - 2003 reads in 1074 unique sequences.
Sample 6 - 1533 reads in 659 unique sequences.
Sample 7 - 29382 reads in 8729 unique sequences.
Sample 8 - 67349 reads in 14713 unique sequences.

## Mostramos los resultados de la primera muestra ##
dadaFs[[1]]

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

dadaRs[[1]]

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

Se procede a identificar y aprender los modelos de error de las lecturas filtradas mediante la función learnErrors de DADA2, la cual estima las tasas de error para cada conjunto de datos de secuencias. Esta función se aplica a las lecturas forward (filtFs) y reverse (filtRs), utilizando múltiples hilos (multithread = TRUE) para acelerar el procesamiento.

Una vez que se han aprendido los modelos de error (errF para las lecturas forward y errR para las lecturas reverse), se aplican estos modelos a las lecturas utilizando la función dada. Esta función corrige errores y realiza inferencias sobre las secuencias de las muestras. Los resultados de este proceso se almacenan en los objetos dadaFs y dadaRs para las lecturas forward y reverse, respectivamente. Finalmente, se muestran los resultados de la primera muestra en cada uno de estos objetos para verificar el éxito del proceso de corrección de errores y la inferencia de secuencias.

3.3.1 Tabla 3: Secuencias Únicas y Variantes

Muestra	Lecturas	Secuencias Únicas (Forward)	Secuencias Únicas (Reverse)	Variantes Inferidas (Forward)	Variantes Inferidas (Reverse)
1	44378	7265	11286	30	-
2	57126	9626	17187	27	-
3	81018	15490	21192	-	-
4	60	47	49	-	-
5	2003	799	1074	-	-
6	1533	462	659	-	-
7	29382	5938	8729	-	-
8	67349	9183	14713	-	-

3.3.2 Gráfia 3: Número de Secuencias Únicas por Muestra y Tipo de Lectura

# EXTRA: Visualizar los resultados con una gráfica
# Rellenamos un data.frame manualmente
results <- data.frame(
  sample = rep(c("Sample 1", "Sample 2", "Sample 3", "Sample 4", "Sample 5", "Sample 6", "Sample 7", "Sample 8"), 2),
  reads = c(44378, 57126, 81018, 60, 2003, 1533, 29382, 67349, 44378, 57126, 81018, 60, 2003, 1533, 29382, 67349),
  unique_sequences = c(7265, 9626, 15490, 47, 799, 462, 5938, 9183, 11286, 17187, 21192, 49, 1074, 659, 8729, 14713),
  type = rep(c("Forward", "Reverse"), each = 8)
)

# Generamos la visualización
ggplot(results, aes(x = sample, y = unique_sequences, fill = type)) +
  geom_bar(stat = "identity", position = "dodge") +
  geom_text(aes(label = unique_sequences), vjust = -0.3, position = position_dodge(0.9), size = 3.5) +
  labs(title = "Número de Secuencias Únicas por Muestra y Tipo de Lectura",
       x = "Muestra",
       y = "Número de Secuencias Únicas") +
  scale_fill_manual(values = c("steelblue", "coral")) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Los resultados muestran una variabilidad considerable en el número de secuencias únicas y variantes inferidas entre las muestras.

En términos de secuencias únicas, las muestras Sample 3 y Sample 8 presentan el mayor número de secuencias únicas tanto en las lecturas forward como en las reverse, con Sample 3 teniendo 15490 y 21192 secuencias únicas, y Sample 8 teniendo 9183 y 14713 secuencias únicas, respectivamente. Esto indica una alta diversidad de secuencias en estas muestras.

Por otro lado, la muestra Sample 4 presenta un número significativamente menor de lecturas y secuencias únicas, con solo 60 lecturas en total y 47-49 secuencias únicas, lo que puede sugerir problemas de calidad o cantidad insuficiente de muestra.

Adicionalmente, se identificaron variantes de secuencia en las muestras Sample 1 y Sample 2 durante el proceso de denoising. En Sample 1, se infirieron 30 variantes de 7265 secuencias únicas forward y 27 variantes de 11286 secuencias únicas reverse. Esto muestra la capacidad del algoritmo DADA2 para corregir errores y determinar variantes precisas en las muestras analizadas.

3.4 Reconstrucción de ASV

Se procede a reconstruir las Amplicon Sequence Variants (ASVs) solapando las lecturas forward (R1) y reverse (R2) mediante la función mergePairs de DADA2. Las ASVs son secuencias únicas de amplicones que representan variantes de secuencia específicas detectadas en las muestras. Este proceso permite ensamblar lecturas emparejadas en secuencias contiguas, mejorando la precisión y la fidelidad de las secuencias resultantes.

### --- Reconstrucción de los ASVs ——————###
## Reconstruimos los ASVs solapando las lecturas R1 y R2 ##
mergers <- mergePairs(dadaFs, filtFs, dadaRs, filtRs, verbose = TRUE)

## Mostramos las primeras líneas del objeto mergers ##
# No lo haremos directamente ya que las secuencias son demasiado largas, 
# pero se incluye una tabla con el contenido resumido
# head(mergers[[1]])

## Generamos una tabla de secuencias con los ASVs reconstruidos ##
seqtab <- makeSequenceTable(mergers)

El objeto mergers es un vector de data.frames. Cada data.frame contiene las siguientes columnas:

• sequence: Las secuencias resultantes después del emparejamiento, no se muestran aquí debido a su longitud.
• abundance: La abundancia de cada secuencia en la muestra.
• forward: El número de la lectura forward que se emparejó.
• reverse: El número de la lectura reverse que se emparejó.
• nmatch: El número de posiciones que coinciden perfectamente entre las lecturas forward y reverse.
• nmismatch: El número de posiciones que no coinciden entre las lecturas forward y reverse.
• nindel: El número de inserciones o deleciones entre las lecturas forward y reverse.
• prefer: La preferencia del emparejamiento basado en el puntaje de calidad.
• accept: Un valor lógico que indica si el emparejamiento fue aceptado o no.

3.4.1 Tabla 4: Reconstrucción de ASV. Primer Elemento.

A continuación, se incluye el contenido de head(mergers[[1]]).

sequence	abundance	forward	reverse	nmatch	nmismatch	nindel	prefer	accept
TACGGAGG…	19806	1	1	249	0	0	1	TRUE
TACGGAGG…	8025	2	2	249	0	0	1	TRUE
TACGGAGG…	7663	3	3	249	0	0	1	TRUE
TACGGAGG…	3674	4	4	246	0	0	1	TRUE
TACGGAGG…	2231	5	5	249	0	0	1	TRUE
TACGGAGG…	1248	6	6	249	0	0	1	TRUE

3.5 Identificación de Qumieras

### --- Identificación de quimeras ——————###
## Eliminamos las secuencias quiméricas de los ASVs reconstruidos ##
seqtab.nochim <- removeBimeraDenovo(seqtab, method = "consensus", multithread = TRUE, verbose = TRUE)

Identified 78 bimeras out of 399 input sequences.

## Relación entre resultados sin quimeras y resultados brutos ##
proportion_non_chimeric <- sum(seqtab.nochim) / sum(seqtab)
proportion_non_chimeric

[1] 0.9879414

## Recopilamos la información del proceso de filtrado a un fichero ##
getN <- function(x) sum(getUniques(x))
track <- cbind(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

## Guardamos la tabla de seguimiento en un archivo
write.table(track, file = "summary_ASV.txt", quote = FALSE, sep = "\t")

## Mostramos la tabla de seguimiento
track

           input filtered denoisedF denoisedR merged nonchim
SRR6442721 51868    44378     44335     44322  44218   44025
SRR6442722 70315    57126     56701     56623  56053   56043
SRR6442723 94944    81018     80579     80683  79377   76516
SRR6442724   118       60        21        25     21      21
SRR6442733  3406     2003      1991      1593   1580    1580
SRR6442734  2022     1533      1409      1394   1345    1335
SRR6442735 35029    29382     29320     28494  28368   28366
SRR6442736 78863    67349     67201     67184  66847   66573

Se procede a identificar y eliminar las secuencias quiméricas de los ASVs reconstruidos utilizando la función removeBimeraDenovo con el método de consenso. Las secuencias quiméricas son artefactos de la PCR que pueden interferir con la precisión de los análisis de diversidad. La eliminación de estas secuencias es crucial para asegurar la integridad de los datos.

Para recopilar y documentar la información del proceso de filtrado, se crea una función getN que suma las secuencias únicas en cada etapa del proceso. Se construye una tabla track que resume el número de lecturas en cada etapa: lecturas originales (input), lecturas filtradas (filtered), lecturas desnoisadas forward (denoisedF) y reverse (denoisedR), lecturas emparejadas (merged), y secuencias no quiméricas (nonchim). Esta tabla se guarda en un archivo summary_ASV.txt y se imprime para su inspección.

Se han identificado y eliminado 78 secuencias quiméricas de un total de 399 secuencias de entrada utilizando la función removeBimeraDenovo. Esto representa una proporción de secuencias quiméricas muy baja, lo cual es un indicador positivo de la calidad de los datos de secuenciación. Después de este paso, se calcula que aproximadamente el 98.79% de las secuencias originales son no quiméricas y por tanto, utilizables para análisis posteriores.

El seguimiento del proceso de filtrado y corrección de errores se resume en la siguiente tabla, que muestra el número de lecturas en cada etapa del procesamiento para cada muestra. Esta información es crucial para evaluar la eficiencia de cada paso y la calidad final de los datos.

3.5.1 Tabla 5.1: Proporción de Secuencias No Quiméricas

Medida	Valor
Secuencias de Entrada	399
Secuencias Quiméricas Identificadas	78
Proporción de Secuencias No Quiméricas	98.79%

3.5.2 Tabla 5.2: Resumen del Proceso de Filtrado y Corrección de Errores

Muestra	input	filtered	denoisedF	denoisedR	merged	nonchim
SRR6442721	51868	44378	44335	44322	44218	44025
SRR6442722	70315	57126	56701	56623	56053	56043
SRR6442723	94944	81018	80579	80683	79377	76516
SRR6442724	118	60	21	25	21	21
SRR6442733	3406	2003	1991	1593	1580	1580
SRR6442734	2022	1533	1409	1394	1345	1335
SRR6442735	35029	29382	29320	28494	28368	28366
SRR6442736	78863	67349	67201	67184	66847	66573

3.5.3 Gráfica 5.1: Número de Lecturas por Etapa para cada Muestra. Barras

# Crear un data frame con los resultados
results <- data.frame(
  sample = rownames(track),
  input = track[, "input"],
  filtered = track[, "filtered"],
  denoisedF = track[, "denoisedF"],
  denoisedR = track[, "denoisedR"],
  merged = track[, "merged"],
  nonchim = track[, "nonchim"]
)

# Transformar los datos para el gráfico de barras
results_long <- melt(results, id.vars = "sample")

# Gráfico de barras
ggplot(results_long, aes(x = sample, y = value, fill = variable)) +
  geom_bar(stat = "identity", position = "dodge") +
  labs(title = "Proceso de Filtrado y Corrección de Errores",
       x = "Muestra",
       y = "Número de Lecturas",
       fill = "Etapa") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

3.6 Gráfica 5.2: Número de Lecturas por Etapa para cada Muestra. Mapa de Calor

# Crear un data frame para el mapa de calor
results_heatmap <- results_long
results_heatmap$value <- log10(results_heatmap$value + 1)  # Usar log10 para mejor visualización

# Mapa de calor
ggplot(results_heatmap, aes(x = sample, y = variable, fill = value)) +
  geom_tile() +
  scale_fill_gradient(low = "white", high = "steelblue") +
  labs(title = "Mapa de Calor del Proceso de Filtrado y Corrección de Errores",
       x = "Muestra",
       y = "Etapa",
       fill = "Log10(Número de Lecturas)") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

4 Asignación Taxonómica de los ASV’s

### --- Realización de la asignación taxonómica de los ASVs ——————###
## Cargamos en un objeto la base de datos taxonómica ##
fastaRef <- "DB/silva_nr_v132_train_set.fa"

## Realizamos la asignación taxonómica e incorporamos las especies ##
taxTab <- assignTaxonomy(seqtab.nochim, refFasta = fastaRef, multithread = TRUE)

## En el caso de querer agregar el rango taxonómico de especies, simplemente usamos una base de datos extra que contiene esta información ##
taxTabExtra <- addSpecies(taxTab, "DB/silva_species_assignment_v132.fa", verbose = TRUE)

33 out of 321 were assigned to the species level.
Of which 32 had genera consistent with the input table.

Se procede a la asignación taxonómica de las Amplicon Sequence Variants (ASVs) reconstruidas y filtradas (sin quimeras). Primero, se carga la base de datos taxonómica de referencia SILVA (versión 132) en el objeto fastaRef. Esta base de datos contiene las secuencias de referencia necesarias para la asignación taxonómica.

Utilizando la función assignTaxonomy, se asignan los rangos taxonómicos a cada ASV en el objeto seqtab.nochim. Esta función utiliza la base de datos de referencia para clasificar cada secuencia en taxones específicos (desde el dominio hasta el género) y se ejecuta en múltiples hilos (multithread = TRUE) para acelerar el proceso.

Para agregar el rango taxonómico de especies, se utiliza la función addSpecies con una base de datos adicional (silva_species_assignment_v132.fa) que contiene información específica de las especies. Esta función incorpora información de especies en la tabla de taxonomía taxTab, proporcionando una clasificación taxonómica más detallada para cada ASV.

4.1 FASTA de los ASV

### --- Extracción del fasta de los ASVs ——————###
## Creamos un objeto con las dimensiones del fichero fasta ##
asv_seqs <- colnames(seqtab.nochim)
asv_headers <- vector("character", length = dim(seqtab.nochim)[2])

## En un bloque añadimos el nombre a las secuencias fasta de los ASVs (>Seq_1, >Seq_2, …) ##
for (i in 1:dim(seqtab.nochim)[2]) {
  asv_headers[i] <- paste(">Seq", i, sep = "_")
}

## Escribimos el fichero de secuencias ##
asv_fasta <- c(rbind(asv_headers, asv_seqs))
write(asv_fasta, "ASVs.fa")

Se procede a la extracción de las secuencias de los Amplicon Sequence Variants (ASVs) y a la creación de un archivo fasta con estas secuencias. Primero, se crea un objeto asv_seqs que contiene los nombres de las columnas de seqtab.nochim, las cuales corresponden a las secuencias ASV.

A continuación, se inicializa un vector asv_headers que almacenará los encabezados de las secuencias fasta. Utilizando un bucle, se asignan nombres a cada secuencia en el formato “>Seq_1”, “>Seq_2”, y así sucesivamente, asegurando que cada ASV tenga un identificador único.

Finalmente, se combinan los encabezados y las secuencias en el objeto asv_fasta y se escribe este objeto en un archivo llamado “ASVs.fa”. Este archivo fasta contiene todas las secuencias ASV con sus correspondientes identificadores, y puede ser utilizado para análisis posteriores o como referencia en otros estudios.

4.2 Creación del Objeto PHYLOSEQ

### --- Creación del objeto phyloseq para la inspección ——————###
## Cargamos los metadatos (diseño experimental) ##
data_sample <- read.table("metadata.txt", header = TRUE, sep = "\t", row.names = 1)
knitr::kable(data_sample)

	sample_ID	bioproject_accession	study	biosample_accession	experiment	run	SRA_Sample	geo_loc_name	collection_date	sample_type	species	common_name	AvgSpotLen
SRR6442724	CHiO1	PRJNA428495	SRP128093	SAMN08292265	SRX3533958	SRR6442724	SRS2809233	Chile: Chiloe	2017	skin	Balaenoptera musculus	blue whale	501
SRR6442723	CHiO2	PRJNA428495	SRP128093	SAMN08292266	SRX3533959	SRR6442723	SRS2809232	Chile: Chiloe	2017	skin	Balaenoptera musculus	blue whale	501
SRR6442722	CHiO3	PRJNA428495	SRP128093	SAMN08292267	SRX3533960	SRR6442722	SRS2809234	Chile: Chiloe	2017	skin	Balaenoptera musculus	blue whale	501
SRR6442721	CHiO4	PRJNA428495	SRP128093	SAMN08292268	SRX3533961	SRR6442721	SRS2809235	Chile: Chiloe	2017	skin	Balaenoptera musculus	blue whale	501
SRR6442736	RNPH11	PRJNA428495	SRP128093	SAMN08292233	SRX3533946	SRR6442736	SRS2809220	Chile: Reserva_Nacional_Pinguino_de_Humboldt	2017	skin	Balaenoptera physalus	fin whale	501
SRR6442735	RNPH12	PRJNA428495	SRP128093	SAMN08292234	SRX3533947	SRR6442735	SRS2809221	Chile: Reserva_Nacional_Pinguino_de_Humboldt	2017	skin	Balaenoptera physalus	fin whale	501
SRR6442734	RNPH13	PRJNA428495	SRP128093	SAMN08292235	SRX3533948	SRR6442734	SRS2809222	Chile: Reserva_Nacional_Pinguino_de_Humboldt	2017	skin	Balaenoptera physalus	fin whale	501
SRR6442733	RNPH14	PRJNA428495	SRP128093	SAMN08292236	SRX3533949	SRR6442733	SRS2809223	Chile: Reserva_Nacional_Pinguino_de_Humboldt	2017	skin	Balaenoptera physalus	fin whale	501

## Creamos el objeto phyloseq ##
ps_orig <- phyloseq(
  otu_table(as.matrix(seqtab.nochim), taxa_are_rows = FALSE),
  sample_data(data_sample),
  tax_table(as.matrix(taxTabExtra))
)

## Mostramos el objeto phyloseq creado ##
ps_orig

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 321 taxa and 8 samples ]
sample_data() Sample Data:       [ 8 samples by 13 sample variables ]
tax_table()   Taxonomy Table:    [ 321 taxa by 7 taxonomic ranks ]

Se procede a la creación de un objeto phyloseq para la inspección y análisis de los datos de secuenciación. El paquete phyloseq en R facilita la integración y el análisis de datos de microbiomas.

Primero, se cargan los metadatos experimentales desde un archivo llamado metadata.txt. Este archivo contiene información adicional sobre las muestras, como las condiciones experimentales, y se carga en el objeto data_sample.

A continuación, se crea el objeto phyloseq (ps_orig) utilizando tres componentes principales:

• otu_table: Se crea una tabla de Unidades Taxonómicas Operacionales (OTUs) a partir de seqtab.nochim, con las secuencias como columnas y las muestras como filas.
• sample_data: Se incorpora la información de los metadatos de las muestras cargada en data_sample.
• tax_table: Se añade la tabla de taxonomía taxTabExtra que contiene la clasificación taxonómica de las secuencias.

Los datos proporcionados son metadatos experimentales para un conjunto de muestras de secuenciación de microbiomas de piel de ballenas recolectadas en diferentes ubicaciones de Chile en 2017. Cada muestra tiene un identificador único de secuenciación (por ejemplo, SRR6442724) y está asociada con diversos atributos como el proyecto biológico, el estudio, el tipo de muestra, la localización geográfica y la especie del huésped.

Las muestras se dividen en dos grupos principales basados en su localización geográfica: Chiloe y Reserva Nacional Pingüino de Humboldt. Además, cada grupo representa diferentes especies de ballenas, específicamente la ballena azul (Balaenoptera musculus) y la ballena de aleta (Balaenoptera physalus). Estos metadatos son esenciales para realizar análisis que relacionen la diversidad microbiana con factores como el hábitat y la especie del huésped.

4.3 Modificación del objeto PHYLOSEQ

### --- Incorporación del nombre de las secuencias a los ASVs reconstruidos ——————###
## Añadimos nombres a los ASVs ##
n_seqs <- seq(ntaxa(ps_orig))
taxa_names(ps_orig) <- paste("Seq", formatC(n_seqs, flag = "0"), sep = "_")
head(taxa_names(ps_orig)) # una vista rápida de los nuevos nombres

[1] "Seq_1" "Seq_2" "Seq_3" "Seq_4" "Seq_5" "Seq_6"

## Calculamos las distancias entre los taxones ##
library(ape)  # Asegúrate de tener el paquete ape instalado
random_tree <- rtree(ntaxa(ps_orig), rooted = TRUE, tip.label = taxa_names(ps_orig))
random_tree

Phylogenetic tree with 321 tips and 320 internal nodes.

Tip labels:
  Seq_103, Seq_104, Seq_74, Seq_94, Seq_314, Seq_276, ...

Rooted; includes branch lengths.

## Recalculamos el objeto phyloseq ##
otu <- otu_table(ps_orig)
tax <- tax_table(ps_orig)
data <- sample_data(ps_orig)
ps <- phyloseq(otu, tax, data, random_tree)
ps

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 321 taxa and 8 samples ]
sample_data() Sample Data:       [ 8 samples by 13 sample variables ]
tax_table()   Taxonomy Table:    [ 321 taxa by 7 taxonomic ranks ]
phy_tree()    Phylogenetic Tree: [ 321 tips and 320 internal nodes ]

Los resultados indican que se ha creado exitosamente un objeto phyloseq que integra varias capas de datos, incluyendo una tabla de OTUs, datos de muestras, una tabla de taxonomía y un árbol filogenético. Este objeto contiene 321 taxones y 8 muestras, con 13 variables de muestra. La tabla de taxonomía abarca 7 rangos taxonómicos, y el árbol filogenético tiene 321 tips y 320 nodos internos.

Este objeto phyloseq permite una inspección detallada y análisis de la composición del microbioma, las relaciones filogenéticas y las características ambientales de las muestras. La incorporación de nombres secuenciales a los ASVs (por ejemplo, “Seq_1”, “Seq_2”, etc.) facilita la identificación y el seguimiento de cada secuencia a lo largo de los análisis.

4.3.1 Tabla 6.1: Resumen del Objeto phyloseq

Elemento	Descripción
OTU Table	321 taxones y 8 muestras
Sample Data	8 muestras con 13 variables de muestra
Taxonomy Table	321 taxones con 7 rangos taxonómicos
Phylogenetic Tree	Árbol filogenético con 321 tips y 320 nodos internos

4.4 Cálculo de la Biodiversidad

### --- Cálculo y visualización de la diversidad alfa y beta ——————###
## Calculamos los índices de diversidad alfa ##
alpha_diversity <- estimate_richness(ps, measures = c("Observed", "Shannon", "Simpson"))
## Calculamos los índices de diversidad beta ##
du <- ordinate(ps, "PCoA", "unifrac", weighted = TRUE)

1. Diversidad Alfa:

• Se calculan los índices de diversidad alfa (Observed, Shannon, y Simpson) con la función estimate_richness.
• Se generan gráficos de diversidad alfa usando la función plot_richness:
  • El primer gráfico muestra la diversidad alfa por muestra (run), coloreada por el tratamiento (common_name).
  • El segundo gráfico muestra la diversidad alfa por tratamiento (common_name), utilizando un boxplot para representar la distribución del índice de Shannon.

2. Diversidad Beta:

• Se calcula la diversidad beta mediante un análisis de coordenadas principales (PCoA) utilizando la distancia Unifrac ponderada con la función ordinate.
• Se genera un gráfico de la ordenación utilizando la función plot_ordination, coloreado por el tratamiento (common_name) y ajustando el tamaño y la transparencia de los puntos para mejorar la visualización.

4.4.1 Gráfica 7.1: Diversidad Alfa por Tratamiento y Muestra

## Representamos la diversidad alfa por muestra coloreada por tratamiento ##
plot_richness(ps, x = "run", measures = c("Shannon", "Simpson"), color = "common_name") + 
  theme_bw()

Se observa que la diversidad alfa, tanto en términos de Shannon como de Simpson, varía entre las muestras. Las muestras de ballena azul tienden a mostrar una mayor diversidad alfa (mayores valores de los índices) en comparación con las muestras de ballena de aleta. Esto sugiere que las comunidades microbianas asociadas a las ballenas azules son más diversas que las asociadas a las ballenas de aleta. La diferencia en los índices de diversidad alfa puede reflejar variaciones en la composición y estructura de las comunidades microbianas entre las dos especies de ballenas.

4.4.2 Gráfica 7.2: Diversidad Alfa por Tratamiento

## Representamos la diversidad alfa por tratamiento ##
plot_richness(ps, x = "common_name", color = "common_name", measures = c("Shannon")) + 
  geom_boxplot()

Cada boxplot representa la distribución de los valores de diversidad alfa (índice de Shannon) para las muestras asociadas a cada especie.

• Para la ballena azul (blue whale), los valores de diversidad alfa muestran una mediana más alta y una mayor dispersión en comparación con la ballena de aleta (fin whale), indicando una mayor variabilidad en la diversidad de las comunidades microbianas asociadas a esta especie.
• Para la ballena de aleta (fin whale), la mediana de la diversidad alfa es ligeramente más baja y los valores están más concentrados, lo que sugiere una comunidad microbiana más homogénea en estas muestras.

4.4.3 Gráfica 7.3: Diversidad Beta

## Calculamos y representamos la diversidad beta ##
p <- plot_ordination(ps, du, color = "common_name") +
  geom_point(size = 7, alpha = 0.75)
p

Se observa una clara separación entre las muestras de ballena azul y las de ballena de aleta a lo largo del eje X, lo que sugiere que la composición de las comunidades microbianas es significativamente diferente entre estas dos especies de ballenas. La agrupación de puntos dentro de cada grupo de especies indica una mayor similitud en la composición microbiana entre las muestras de la misma especie.

Esta diferenciación en la diversidad beta entre las especies refuerza la idea de que diferentes especies de ballenas albergan comunidades microbianas distintas, probablemente influenciadas por factores específicos del huésped como la dieta, el comportamiento y el hábitat.

4.5 Filtrado por Abundancia

### --- Filtramos aquellos taxones menos abundantes ——————###
## Filtramos los taxones que aparecen menos de 10 veces en el 10% de las muestras ##
ps_filter <- filter_taxa(ps, function(x) sum(x > 10) > (0.1 * length(x)), TRUE)
ps_filter

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 238 taxa and 8 samples ]
sample_data() Sample Data:       [ 8 samples by 13 sample variables ]
tax_table()   Taxonomy Table:    [ 238 taxa by 7 taxonomic ranks ]
phy_tree()    Phylogenetic Tree: [ 238 tips and 237 internal nodes ]

## Calculamos la abundancia relativa ##
ps.relative <- transform_sample_counts(ps_filter, function(OTU) (OTU / sum(OTU)) * 100)
count_relative <- t(otu_table(ps_filter))
tax_filter <- tax_table(ps_filter)
abun_tax <- cbind(count_relative, tax_filter)

4.5.1 Gráfica 8: Taxones más abundantes

## Con un diagrama de barras representamos los taxones más abundantes ##
fantaxtic_bar(ps_filter, color_by = "Phylum", label_by = "Class", facet_by = NULL, grid_by = NULL, other_color = "Grey")

Se observa que la composición taxonómica varía significativamente entre las muestras. Algunas muestras, como SRR6442734 y SRR6442735, están dominadas por un pequeño número de taxones, mientras que otras, como SRR6442721 y SRR6442723, muestran una mayor diversidad de taxones con una distribución más equilibrada de abundancias.

Las muestras asociadas a las ballenas de diferentes especies y ubicaciones geográficas pueden presentar diferencias en la composición microbiana. Esta variabilidad podría estar influenciada por factores específicos del huésped y el entorno, como la dieta, el comportamiento, el hábitat y la interacción con otros organismos.

5 Conclusión

El análisis metagenómico realizado en este estudio ha permitido caracterizar la diversidad microbiana presente en las muestras de piel de ballenas recolectadas en Chile. Los resultados indican una alta calidad en las lecturas de secuenciación, con una proporción significativa de secuencias pasando los filtros de calidad establecidos. La identificación y eliminación de secuencias quiméricas ha mejorado la precisión de los datos, permitiendo una asignación taxonómica detallada utilizando la base de datos SILVA.

Los índices de diversidad alfa y beta han revelado diferencias significativas en la composición microbiana entre las especies de ballenas y sus respectivas ubicaciones geográficas. Las ballenas azules mostraron una mayor diversidad alfa en comparación con las ballenas de aleta, lo que sugiere una mayor heterogeneidad en las comunidades microbianas asociadas. Además, el análisis de la diversidad beta ha demostrado una clara separación en la composición microbiana entre las dos especies de ballenas, indicando que factores específicos del huésped y el entorno influyen en la estructura de las comunidades microbianas.

El filtrado de taxones menos abundantes y la transformación a abundancias relativas han permitido identificar los taxones dominantes en cada muestra, proporcionando una visión detallada de la composición taxonómica. Este enfoque ha destacado la variabilidad en la composición microbiana entre las muestras, reflejando la influencia de factores específicos del huésped y del ambiente.