Objetivos

Paquetes

Instalar los paquetes que no se encuentran instalados en el sistema.

library(ggplot2)
library(phyloseq)
library(ranacapa)
library(vegan)
library(ggsignif)

1. Creación del objeto phyloseq

Situarse en el directorio que contiene los objetos obtenidos del análisis con DADA2.

setwd("ruta/a/objetos")

Cargamos la tabla de secuencias (ASVs) después de eliminar las quimeras y la tabla con la asignación taxonómica de las ASVs.

load("seqtab.nochim")
load("taxa")

Cargar la tabla de metadats

sample.metadata <-read.csv("sampleData_alimentos.csv")
row.names(sample.metadata) <- sample.metadata$Run

Ahora construimos un objeto phyloseq directamente desde las salidas DADA2.

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

Función phyloseq()

La función phyloseq() del paquete phyloseq se utiliza para crear un objeto integrado para el análisis, que combina en una sola estructura la información de:

  • otu_table(): matriz de abundancias de ASVs u OTUs por muestra.
  • sample_data(): metadatos de las muestras (ej. sitio de muestreo, tratamiento).
  • tax_table(): tabla taxonómica con la clasificación jerárquica (Reino, Filo, Clase, Orden, Familia, Género, Especie).

El objeto resultante (physeq) es un contenedor que facilita el manejo conjunto de abundancias, taxonomía y metadatos, permitiendo realizar análisis estadísticos y gráficos complejos de forma más sencilla.

2. Curva de rarefacción

p1 <- ggrare(ps_subset,
             step  = 50,
             label = NULL,    
             color = "Sampling.position",
             se    = FALSE)

p1 <- p1 +
        geom_line(size = 1, linetype = "solid") +
        guides(
                color = guide_legend(override.aes = list(shape = NA)),  
                label = "none"                                          
        ) +
        theme_minimal() +
        labs(
                x     = "Number of sampled reads",
                y     = "Number of observed TAXA",
                title = ""
        )
p1

Función ggrare()

La función ggrare(), del paquete ranacapa, genera curvas de rarefacción utilizando ggplot2 para visualizar cómo la riqueza observada (número de ASVs/OTUs) cambia al submuestrear lecturas de una muestra.

La función ggrare()utiliza el mismo algoritmo que vegan::rarecurve().La diferencia es que ggrare() devuelve los resultados en un formato compatible con ggplot2, lo que permite personalizar y mejorar la visualización.

Los conceptos de diveridad de las especies fueron introducidos por Whittaker (1972).

3. Diversidad alfa


alpha <- plot_richness(
        ps_subset,
        x = "Sampling.position",
        measures = c("Shannon", "Simpson", "Observed"),
        color = "Sampling.position"
) +
        geom_point(size =2 , aes(color = Sampling.position)) +
        theme_minimal() +
        theme(
                axis.title.x = element_blank(),
                axis.text.x = element_text(size = 10, angle = 75, hjust = 1, vjust = 1),
                axis.line.x = element_line(color = "black", size = 0.5),
                axis.line.y = element_line(color = "black", size = 0.5),
                panel.background = element_blank(),
                panel.grid.major = element_blank(),
                panel.grid.minor = element_blank()
        ) +
        geom_signif(
                comparisons = list(c("Sticking", "Anal swab")),
                test = "wilcox.test",           
                test.args = list(exact = FALSE), 
                map_signif_level = TRUE,          
                step_increase = 0.1,
                textsize = 5,
                tip_length = 0,
                color = "black"                
        ) + ggtitle("Alpha Diversity")

alpha

Función plot_richness()

La función plot_richness() del paquete phyloseq genera gráficos de diversidad alfa (riqueza dentro de cada muestra), utilizando distintos índices ecológicos.

Parámetros principales:

physeq: objeto de clase phyloseq. x: variable de metadatos que se usará en el eje X (en este caso, “Sampling.position”).
measures: vector con los índices de diversidad a calcular. Los más comunes son:
- “Shannon”: combina riqueza y equidad de especies.
- “Simpson”: probabilidad de que dos lecturas elegidas al azar correspondan a la misma especie.
- “Observed”: número de ASVs/OTUs observadas (riqueza pura).
color: variable de metadatos para colorear los puntos o cajas del gráfico (aquí “Sampling.position”).

3. Diversidad beta

#Transformar a abundancias relativas

ps.prop <- transform_sample_counts(ps_subset, function(x) x / sum(x))

Ordinar con NMDS (Bray-Curtis) y graficar
set.seed(123)  # reproducibilidad del NMDS
ord.nmds.bray <- ordinate(ps.prop, method = "NMDS", distance = "bray")

ordplot <- plot_ordination(
  ps.prop, ord.nmds.bray,
  title = "",
  color = "Sampling.position",
  shape = "Sampling.position"
) +
  geom_point(size = 5) +
  theme_bw() +
  theme(text = element_text(size = 16))

#Agregar elipses por grupo (en este caso, Sampling.position)
p_ord <- ordplot +
  labs(colour = "Sampling.position", fill = "Sampling.position", shape = "Sampling.position") +
  ggforce::geom_mark_ellipse(aes(label = Sampling.position, fill = Sampling.position)) +
  theme(legend.position = "none")

#PERMANOVA (Bray-Curtis) con metadatos
metadata <- as(sample_data(ps_subset), "data.frame")

metadata <- metadata[sample_names(ps_subset), , drop = FALSE]

dist.bc <- phyloseq::distance(ps_subset, method = "bray")
permanova <- vegan::adonis2(dist.bc ~ Sampling.position, data = metadata)
pval <- permanova$`Pr(>F)`[1]

#Función para agregar asteriscos
p_stars <- function(p) {
  if (is.na(p)) return("")
  if (p < 0.001) "***" else if (p < 0.01) "**" else if (p < 0.05) "*" else "ns"
}

#Etiqueta con p-valor
lab <- sprintf("PERMANOVA (Bray–Curtis): p = %.3g  %s", pval, p_stars(pval))

p_ord_annot <- p_ord +
  annotate(
    "label",
    x = Inf, y = Inf, label = lab,
    hjust = 1.02, vjust = 1.5, size = 3.5, label.size = 0.3
  ) +
  coord_cartesian(clip = "off") +
  theme(plot.margin = margin(10, 40, 10, 10))

p_ord_annot

Función ordinate()

La función ordinate() del paquete phyloseq se utiliza para realizar análisis de ordenación multivariada a partir de un objeto phyloseq.

La ordenación es una técnica que permite resumir la complejidad de las comunidades microbianas en un espacio de pocas dimensiones, facilitando la visualización y comparación entre muestras.

Función adonis2()

La función adonis2() del paquete vegan se utiliza para realizar un PERMANOVA (Permutational Multivariate Analysis of Variance), una prueba estadística no paramétrica que evalúa si las comunidades microbianas difieren significativamente entre grupos definidos por una variable categórica (por ejemplo, sitio, tratamiento, condición experimental).