Análisis ecológico de comunidades de peces

Golfo de California — Metabarcoding 12S rDNA

Este reporte presenta el flujo completo de análisis ecológico de comunidades de peces en el Golfo de California a partir de datos de metabarcoding del gen 12S rDNA. Partiendo de los archivos .RData generados por el pipeline DADA2, construimos el objeto phyloseq, evaluamos la calidad de los datos y aplicamos análisis de diversidad alfa, composición taxonómica y diversidad beta.

1. Preparación

Advertencia

Bug conocido: la columna sample_id en metadata.csv puede estar vacía (NA). Si es tu caso, restaura el archivo original desde GitHub antes de abrir R:

cd ~/metabarcoding-code
git checkout origin/main -- metadata/metadata.csv

Carga de paquetes

Abre R desde ~/metabarcoding-code/. Cargamos todos los paquetes necesarios al inicio del documento. Usamos suppressPackageStartupMessages() para mantener el output limpio. Los paquetes están organizados por función: manipulación de datos, análisis ecológico, visualización y tablas.

# Evitar que kable intente abrir navegador en el servidor
options(browser = "false")

suppressPackageStartupMessages({
  # Manipulación de datos
  library(dplyr)
  library(tidyr)
  library(tibble)

  # Bioinformática y ecología de comunidades
  library(phyloseq)
  library(Biostrings)
  library(vegan)

  # Visualización
  library(ggplot2)
  library(ggforce)
  library(gridExtra)
  library(viridis)
  library(RColorBrewer)
  library(pheatmap)
  library(colorspace)

  # Tablas
  library(knitr)
  library(kableExtra)
})

Definición de rutas

setwd("~/metabarcoding-code")

ruta_rdata       <- "final_data/rdata_output"
ruta_metadata    <- "metadata/Fish_speciesGC.csv"
ruta_csv_output  <- "final_data/csv_output/taxon_table.csv"
ruta_figs        <- "final_data/reporte_ecologico"

Carga de datos

Cargamos el archivo .RData más reciente generado por el pipeline DADA2, que contiene las matrices de abundancias, taxonomía y metadata.

# Buscar todos los .RData en el directorio de salida del pipeline
rdata_files <- list.files(ruta_rdata, pattern = "\\.RData$", full.names = TRUE)
stopifnot("No se encontró ningún .RData en rdata_output" = length(rdata_files) > 0)

# Usar el más reciente
rdata_file <- rdata_files[which.max(file.info(rdata_files)$mtime)]
message("Cargando (más reciente): ", basename(rdata_file))
load(rdata_file)

# --- Alinear sample_meta con rownames de asv_mat ---
# El RData puede traer sample_meta con rownames NA (bug en metadata.csv).
# Siempre recargamos desde el CSV y usamos rename_to como fuente confiable del ID.
sample_meta <- readr::read_csv(
  file.path("metadata", "metadata.csv"), show_col_types = FALSE
)
sample_meta <- as.data.frame(sample_meta)

# rename_to tiene el nombre completo del FASTQ: "12S-GC-001-d1_1_S1_L001"
# asv_mat usa una versión simplificada:          "12S-GC-001-d1_S1"
# Transformamos rename_to para que coincida con asv_mat
if ("rename_to" %in% colnames(sample_meta)) {
  meta_names <- sub("_L001$", "", sample_meta$rename_to)
  meta_names <- sub("(d1)_1_(S)", "\\1_\\2", meta_names)
  rownames(sample_meta) <- meta_names
} else {
  stop("metadata.csv no tiene columna 'rename_to'. Restaura desde GitHub:\n",
       "  cd ~/metabarcoding-code && git checkout origin/main -- metadata/metadata.csv")
}

sample_meta <- sample_meta[rownames(asv_mat), , drop = FALSE]

message("Objetos cargados: ", paste(ls(), collapse = ", "))
message("Muestras en asv_mat: ", nrow(asv_mat))
message("Muestras en sample_meta: ", nrow(sample_meta))
message("NAs en sample_meta: ", sum(is.na(sample_meta[,1])))

---

2. Construcción del objeto phyloseq

phyloseq es un paquete de R diseñado para el análisis de datos de secuenciación de amplicones (16S, 18S, ITS, 12S, etc.). Su función principal es integrar en un único objeto todas las capas de información de un estudio de microbiomas o comunidades: las abundancias de secuencias, la taxonomía y los metadatos de las muestras.

¿Qué es el objeto phyloseq?

El objeto phyloseq es un contenedor que agrupa tres tablas que deben estar siempre sincronizadas entre sí:

Componente	Función	Dimensiones
otu_table	Abundancias de cada ASV en cada muestra	ASVs × muestras
tax_table	Clasificación taxonómica de cada ASV	ASVs × niveles taxonómicos
sample_data	Metadatos de las muestras (sitio, fecha, tratamiento, etc.)	muestras × variables

Opcionalmente puede incluir un cuarto componente:

Componente	Función
refseq	Secuencia nucleotídica de cada ASV

La ventaja de este contenedor es que cualquier operación de filtrado o subsetting afecta a todas las tablas simultáneamente y de forma consistente. Por ejemplo, al filtrar para conservar solo muestras de un sitio, la tabla de abundancias, la taxonomía y los metadatos se reducen juntas — sin riesgo de desincronización.

---

Matrices de datos

Construimos las tres matrices que necesita phyloseq: abundancias (ASVs × muestras), taxonomía y metadata de muestras.

# Matriz de abundancias: ASVs como filas, muestras como columnas
otu_raw <- if (!is.null(asv_mat)) asv_mat else seqtab_clean
stopifnot(!is.null(otu_raw))
ASV_mat <- t(as.matrix(otu_raw))

# Tabla taxonómica
if (!is.null(tax_mat)) {
  taxa_df <- as.data.frame(tax_mat)
  for (rk in c("Kingdom","Phylum","Class","Order","Family","Genus","Species")) {
    if (!rk %in% colnames(taxa_df)) taxa_df[[rk]] <- NA_character_
  }
} else {
  stopifnot(!is.null(taxonomy_df), "Sequence" %in% names(taxonomy_df))
  keep <- intersect(
    c("Sequence","Kingdom","Phylum","Class","Order","Family","Genus","Species"),
    names(taxonomy_df)
  )
  taxa_df <- taxonomy_df[, keep, drop = FALSE] %>%
    dplyr::distinct(Sequence, .keep_all = TRUE) %>%
    tibble::column_to_rownames("Sequence")
}

# Alinear filas con ASV_mat y filtrar clases no-peces
taxa_df <- taxa_df[rownames(ASV_mat), , drop = FALSE]
taxa_df <- taxa_df[!(taxa_df$Class %in% c("Mammalia", "Aves") | is.na(taxa_df$Class)), ]

# Metadata de muestras (sample_meta ya fue alineada con asv_mat en cargar-datos)
md <- as.data.frame(sample_meta)

# Verificar que la metadata se cargó correctamente
if (all(is.na(md[, 1]))) {
  stop("sample_meta tiene puros NAs. Revisa que el chunk 'cargar-datos' se haya ejecutado correctamente.\n",
       "Si el problema persiste, restaura metadata.csv desde GitHub:\n",
       "  cd ~/metabarcoding-code && git checkout origin/main -- metadata/metadata.csv")
}

# Corrección de site_id donde sea necesario (dinámica, sin hardcoding)
md <- md %>%
  dplyr::mutate(
    site_id = dplyr::case_when(
      name == "Isla San Marcos (North)" ~ "MAR",
      TRUE ~ site_id
    )
  )

# Mover controles (MOCK/CTR) al final para tablas
is_ctrl <- grepl("MOCK|CTR", md$site_id, ignore.case = TRUE)
is_ctrl[is.na(is_ctrl)] <- FALSE
md <- rbind(md[!is_ctrl, ], md[is_ctrl, ])

stopifnot(nrow(md) > 0)

Curación taxonómica y geográfica

La base de datos de referencia contiene especies de todo el mundo. DADA2 asigna la especie con mayor similitud de secuencia, pero esa especie puede no distribuirse en el Golfo de California. La curación se hace en dos pasos:

Paso 1 — Correcciones de Order: Algunas familias necesitan reclasificación a nivel de orden siguiendo la taxonomía actualizada (Betancur-R et al. 2017), ya que la base de datos de referencia puede tener clasificaciones desactualizadas.

Paso 2 — Curación geográfica automática: Contrastamos cada asignación de especie contra una lista de especies documentadas para el Golfo de California (Fish_speciesGC.csv):

• Especie en la lista del GC → se mantiene tal cual (ej. Mycteroperca rosacea)
• Especie NO en GC, género SÍ, y el género tiene una sola especie en la lista → se asigna esa especie (ej. Elops solo tiene E. affinis en el GC)
• Especie NO en GC, género SÍ, pero con varias especies → se deja como Género sp. (no podemos resolver a especie)
• Ni especie ni género en la lista del GC → se limpia hasta familia (la asignación a género/especie no es confiable para esta región)

# Curación geográfica con lista de especies del Golfo de California
fish_sp      <- read.csv(ruta_metadata, row.names = NULL)
fish_species <- fish_sp$Species2

taxa2 <- as.data.frame(taxa_df)

# Correcciones de Order según taxonomía actualizada (Betancur-R et al. 2017)
taxa2$Order[!is.na(taxa2$Family) & taxa2$Family == "Sphyraenidae"]  <- "Carangiformes"
taxa2$Order[!is.na(taxa2$Family) & taxa2$Family == "Sciaenidae"]    <- "Perciformes"
taxa2$Order[!is.na(taxa2$Family) & taxa2$Family == "Malacanthidae"] <- "Perciformes"
taxa2$Order[!is.na(taxa2$Family) & taxa2$Family == "Pomacanthidae"] <- "Acanthuriformes"

# Filtrar artefactos (conservar Species=NA — la curación geográfica los resolverá)
taxa2 <- taxa2 %>%
  dplyr::filter(
    is.na(Species) |
      (Species != "Eukaryotic synthetic" & !grepl("^NA sp$", Species))
  )

# --- Curación geográfica ---
# Para cada género en la lista del GC, contar cuántas especies tiene
gc_by_genus <- split(fish_species, sub(" .*", "", fish_species))
fish_genera <- names(gc_by_genus)

# Géneros con una sola especie en el GC → se puede resolver la especie directamente
gc_genus_unico <- Filter(function(x) length(x) == 1, gc_by_genus)
gc_genus_unico_vec <- setNames(sapply(gc_genus_unico, `[[`, 1), names(gc_genus_unico))

# Para cada ASV, decidir qué nivel taxonómico mantener:
# 1. Especie en la lista del GC → mantener tal cual
# 2. Especie NO en GC, género SÍ con una sola especie → asignar esa especie
# 3. Especie NO en GC, género SÍ con varias especies → dejar como "Género sp."
# 4. Ni especie ni género en GC → limpiar hasta familia
taxa2 <- taxa2 %>%
  dplyr::mutate(
    in_gc_species  = Species %in% fish_species,
    in_gc_genus    = Genus %in% fish_genera,
    genus_unico    = Genus %in% names(gc_genus_unico_vec),
    Species = dplyr::case_when(
      in_gc_species               ~ Species,
      in_gc_genus & genus_unico   ~ gc_genus_unico_vec[Genus],
      in_gc_genus                 ~ paste0(Genus, " sp."),
      TRUE                        ~ NA_character_
    ),
    Genus = dplyr::case_when(
      is.na(Genus)   ~ NA_character_,
      in_gc_genus    ~ Genus,
      TRUE           ~ NA_character_
    )
  ) %>%
  dplyr::select(-in_gc_species, -in_gc_genus, -genus_unico)

# Resumen de la curación
n_sp_gc    <- sum(!is.na(taxa2$Species) & !grepl(" sp\\.$", taxa2$Species))
n_genus_sp <- sum(grepl(" sp\\.$", taxa2$Species), na.rm = TRUE)
n_familia  <- sum(is.na(taxa2$Species) & !is.na(taxa2$Family))
message("ASVs con especie confirmada en GC: ", n_sp_gc)
message("ASVs reducidas a Género sp.: ", n_genus_sp)
message("ASVs reducidas a familia: ", n_familia)

taxa_mat <- as.matrix(taxa2)

Construcción del objeto phyloseq

phyloseq integra en un solo objeto la tabla de abundancias, la taxonomía y la metadata de muestras. Una vez construido, todas las operaciones — filtrado, transformación, visualización — se realizan de manera consistente y reproducible. También guardamos las secuencias originales de cada ASV en el slot refseq y renombramos los ASVs como ASV1, ASV2, etc., para mayor legibilidad.

ASV_ps <- phyloseq::otu_table(ASV_mat, taxa_are_rows = TRUE)
TAX_ps <- phyloseq::tax_table(taxa_mat)
SMD_ps <- phyloseq::sample_data(md)

ps <- phyloseq::phyloseq(ASV_ps, TAX_ps, SMD_ps)

# Guardar secuencias originales y renombrar ASVs
dna <- Biostrings::DNAStringSet(phyloseq::taxa_names(ps))
names(dna) <- phyloseq::taxa_names(ps)
ps <- phyloseq::merge_phyloseq(ps, dna)
phyloseq::taxa_names(ps) <- paste0("ASV", seq_len(phyloseq::ntaxa(ps)))

message("Objeto phyloseq construido: ",
        phyloseq::nsamples(ps), " muestras, ",
        phyloseq::ntaxa(ps), " ASVs")
print(ps)

Guardar phyloseq y FASTA de ASVs (para VSEARCH)

Antes de continuar con los análisis ecológicos, guardamos el objeto phyloseq con las ASVs individuales (sin colapsar por especie) y exportamos las secuencias a FASTA. Estos archivos se usarán en el Día 4 para el agrupamiento en OTUs con VSEARCH.

dir.create("final_data/otus_vsearch", recursive = TRUE, showWarnings = FALSE)

# Guardar objeto phyloseq
save(ps, file = "final_data/otus_vsearch/phyloseq_pre_collapse.RData")

# Exportar secuencias ASV a FASTA
asv_seqs <- as.character(phyloseq::refseq(ps))
asv_headers <- phyloseq::taxa_names(ps)
asv_fasta <- c(rbind(paste0(">", asv_headers), asv_seqs))
writeLines(asv_fasta, "final_data/otus_vsearch/fish_ASV.fasta")

message("Phyloseq y FASTA guardados en final_data/otus_vsearch/")
message("Total ASVs exportadas: ", length(asv_seqs))

---

OTUs: agrupamiento por similitud con VSEARCH

Las ASVs son secuencias únicas — entre dos ASVs puede haber solo 1 base de diferencia. Esto da gran resolución, pero nos lleva a preguntarnos: ¿existe 1 base de diferencia entre especies diferentes? ¿Cuál es el impacto en las estimaciones de diversidad?

Usamos vsearch (Rognes et al., 2016) para formar OTUs (Operational Taxonomic Units) a partir de las ASVs al 97% de similitud, reconstruimos un objeto phyloseq con la nueva información y comparamos la resolución taxonómica de ambas estrategias.

Agrupamiento de ASVs en OTUs

vsearch agrupa secuencias similares en clusters. El parámetro principal es el porcentaje de similitud (0.97 = 97%). A cada cluster se le asigna una secuencia representante (oturep).

En la terminal (no en R), ejecuta:

cd ~/metabarcoding-code
bash scripts/08b-vsearch.sh 0.97

Resultado: tres archivos generados:

Archivo	Contenido
_oturep.fasta	Secuencias representantes en formato FASTA
.uc	Relación de las agrupaciones en formato uclust
.log	Registro de la ejecución de vsearch

Procesar archivo .uc

El archivo .uc es una tabla separada por tabuladores con un renglón por cada secuencia analizada:

Columna	Descripción
1	Tipo: S (centroide/oturep), H (hit/miembro), C (registro de cluster)
2	Número de cluster (base 0)
3	Longitud de la secuencia
4	Porcentaje de identidad con el centroide (solo tipo H)
9	Etiqueta de la secuencia
10	Etiqueta del centroide (solo tipo H)

De vuelta en R, procesamos el archivo .uc:

uc_file <- "final_data/otus_vsearch/otus_0.97.uc"

uc_names <- c("Type", "clus_ID", "len", "id_pct", "strand",
              "NA1", "NA2", "aln", "ASV_ID", "centroid")

uc.df <- readr::read_tsv(uc_file, col_names = uc_names, show_col_types = FALSE) |>
  dplyr::select(-NA1, -NA2) |>
  dplyr::filter(!Type %in% "C") |>
  dplyr::relocate(ASV_ID)

# Número de OTUs (otureps) vs miembros
table(uc.df$Type)

La tabla muestra dos tipos de registros:

• S (Seed/centroide): secuencias que son representantes de un cluster — cada una define un OTU. En este caso 981 OTUs.
• H (Hit/miembro): secuencias que fueron asignadas a un cluster existente porque tienen ≥97% de similitud con el centroide. Aquí 422 ASVs fueron absorbidas por un OTU existente.

Es decir, de las 1,403 ASVs originales (981 + 422), VSEARCH las agrupó en 981 OTUs. La reducción (~30%) refleja la variación intraespecífica (haplotipos, errores residuales) que se colapsa al agrupar al 97%.

# Identificar otureps y su cluster
oturep.df <- uc.df |> dplyr::filter(Type %in% "S") |> dplyr::select(clus_ID, ASV_ID)
oturepID <- oturep.df |> dplyr::pull(ASV_ID)
message("Total OTUs (otureps): ", length(oturepID))

Reconstruir objeto phyloseq con OTUs

Extraemos los componentes del phyloseq original correspondientes a los otureps y colapsamos los conteos por cluster:

ps.asv <- ps

# Subset: solo los otureps
ps.sub <- phyloseq::subset_taxa(ps.asv, phyloseq::taxa_names(ps.asv) %in% oturepID)

oturep.seq <- phyloseq::refseq(ps.sub)
oturep.tax <- phyloseq::tax_table(ps.sub)
oturep.md  <- phyloseq::sample_data(ps.sub)

# Tabla de conteos original: asegurar formato muestras × ASVs (filas = muestras)
otu_mat <- as.matrix(phyloseq::otu_table(ps.asv))
if (phyloseq::taxa_are_rows(ps.asv)) otu_mat <- t(otu_mat)
conteos.ini <- data.frame(otu_mat, stringsAsFactors = FALSE, check.names = FALSE)

conteos.ini.long <- conteos.ini |>
  tibble::rownames_to_column(var = "SampleID") |>
  tidyr::pivot_longer(-SampleID, names_to = "ASV_ID", values_to = "conteo")

conteos.ini.long <- dplyr::left_join(conteos.ini.long, uc.df, by = "ASV_ID")

# Sumar conteos por muestra y cluster
conteos.sum.long <- conteos.ini.long |>
  dplyr::group_by(SampleID, clus_ID) |>
  dplyr::summarise(conteo_cluster = sum(conteo), .groups = "drop")

# Reemplazar cluster_ID por ASV_ID del oturep
conteos.sum.long <- dplyr::left_join(conteos.sum.long, oturep.df, by = "clus_ID")

# Convertir a matriz muestras × OTUs
conteos.sum <- conteos.sum.long |>
  dplyr::select(-clus_ID) |>
  tidyr::pivot_wider(names_from = ASV_ID, values_from = conteo_cluster, values_fill = 0) |>
  tibble::column_to_rownames(var = "SampleID")
conteos.sum <- as.matrix(conteos.sum)

message("Dimensiones OTU table: ", nrow(conteos.sum), " muestras × ", ncol(conteos.sum), " OTUs")

# Construir nuevo objeto phyloseq
ASV_otu  <- phyloseq::otu_table(conteos.sum, taxa_are_rows = FALSE)
TAXA_otu <- phyloseq::tax_table(oturep.tax)
SEQ_otu  <- phyloseq::refseq(oturep.seq)

ps.otu <- phyloseq::phyloseq(ASV_otu, TAXA_otu, SEQ_otu, oturep.md)
ps.otu

Comparar resolución taxonómica: ASVs vs OTUs

asv.tax <- data.frame(phyloseq::tax_table(ps.asv), stringsAsFactors = FALSE)
otu.tax <- data.frame(phyloseq::tax_table(ps.otu), stringsAsFactors = FALSE)

# --- Panel 1: Total ASVs vs OTUs ---
n_asv <- phyloseq::ntaxa(ps.asv)
n_otu <- phyloseq::ntaxa(ps.otu)
pct_red <- round((1 - n_otu / n_asv) * 100, 1)

df_total <- data.frame(
  Estrategia = factor(c("ASVs", "OTUs (97%)"), levels = c("ASVs", "OTUs (97%)")),
  Total = c(n_asv, n_otu)
)

p1 <- ggplot(df_total, aes(x = Estrategia, y = Total, fill = Estrategia)) +
  geom_col(width = 0.6, alpha = 0.85) +
  geom_text(aes(label = scales::comma(Total)), vjust = -0.5, fontface = "bold", size = 4.5) +
  scale_fill_manual(values = c("ASVs" = "#1b9e77", "OTUs (97%)" = "#d95f02")) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.15))) +
  labs(
    title = paste0("Reducción: ", pct_red, "%"),
    x = NULL, y = "Total de secuencias"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(face = "bold", hjust = 0.5),
    axis.title = element_text(face = "bold"),
    legend.position = "none"
  )

# --- Panel 2: Taxones únicos por nivel ---
df_comp <- data.frame(
  ASVs = apply(asv.tax, 2, function(x) length(unique(na.omit(x)))),
  OTUs = apply(otu.tax, 2, function(x) length(unique(na.omit(x))))
)

df_comp <- df_comp |>
  tibble::rownames_to_column(var = "Nivel") |>
  tidyr::pivot_longer(-Nivel, names_to = "Estrategia", values_to = "Conteo") |>
  dplyr::mutate(Nivel = factor(Nivel,
    levels = c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species"),
    ordered = TRUE))

p2 <- ggplot(df_comp, aes(x = Nivel, y = Conteo, fill = Estrategia)) +
  geom_col(position = "dodge", width = 0.7, alpha = 0.85) +
  geom_text(aes(label = Conteo, group = Estrategia),
            position = position_dodge(width = 0.7), vjust = -0.5, size = 3) +
  scale_fill_manual(values = c("ASVs" = "#1b9e77", "OTUs" = "#d95f02")) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.15))) +
  labs(
    title = "Taxones únicos por nivel",
    x = "Nivel taxonómico", y = "Número de taxones"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(face = "bold", hjust = 0.5),
    axis.title = element_text(face = "bold"),
    axis.text.x = element_text(angle = 45, hjust = 1),
    legend.position = "bottom"
  )

# Combinar paneles
gridExtra::grid.arrange(p1, p2, ncol = 2, widths = c(1, 2.5),
                        top = grid::textGrob("Comparación ASVs vs OTUs (97% similitud)",
                                             gp = grid::gpar(fontface = "bold", fontsize = 14)))

ggsave("final_data/reporte_ecologico/comparacion_asvs_vs_otus.png",
       gridExtra::arrangeGrob(p1, p2, ncol = 2, widths = c(1, 2.5)),
       width = 12, height = 5, dpi = 300)

Nota

Interpretación: De 1,403 ASVs, VSEARCH las agrupó en 981 OTUs (reducción del ~30%), pero la resolución taxonómica prácticamente no cambió: de Kingdom a Family los números son idénticos, y a nivel de Species solo se pierden 2. Esto indica que las ASVs colapsadas eran variación intraespecífica (haplotipos de la misma especie), no especies distintas. Al colapsar por especie con tax_glom(), ambas estrategias producirán resultados muy similares.

Tip

¿Cuándo sí importa la elección ASVs vs OTUs?

En este caso los resultados fueron casi idénticos porque 12S es un marcador conservado, los peces del Golfo de California están bien representados en la base de datos, y la variación intraespecífica es baja. Sin embargo, los resultados pueden diferir significativamente cuando:

• El marcador tiene mayor variación intraespecífica — COI (barcoding animal) o ITS (hongos) generan más ASVs por especie. Mayor variación intraespecífica = mayor diferencia entre ASVs y OTUs.
• La comunidad es más diversa — en invertebrados, hongos o bacterias (16S), con cientos de géneros cercanos, el umbral de 97% puede colapsar especies distintas en el mismo OTU, perdiendo resolución taxonómica real.
• La base de datos es incompleta — si muchas ASVs quedan sin asignación a especie, no se pueden colapsar con tax_glom(). Los OTUs retienen esa diversidad “sin nombre” mejor que las ASVs sueltas.
• La divergencia interespecífica es baja — en grupos donde especies distintas difieren <1% (ej. Vibrio en 16S), el 97% colapsa especies reales en un solo OTU.

En resumen: para análisis ecológicos de comunidades a nivel de especie, ambas estrategias suelen dar resultados equivalentes. La diferencia importa cuando se trabaja sin colapsar por especie, o cuando el marcador y el grupo taxonómico lo demandan.

Filtrado de muestras y transformación

Separamos controles y muestras problemáticas antes de los análisis ecológicos. No aplicamos rarefacción — usamos abundancias relativas para composición y β-diversidad, y Chao1/ACE para estimar riqueza.

# Separar controles (MOCK, CTR) para análisis independiente
ps_ctrl <- phyloseq::subset_samples(ps, grepl("MOCK|CTR", site_id, ignore.case = TRUE))
ps_ctrl <- phyloseq::prune_taxa(phyloseq::taxa_sums(ps_ctrl) > 0, ps_ctrl)
message("Muestras control separadas: ", phyloseq::nsamples(ps_ctrl))

# Excluir controles y sitios SLG, LIB, BSN
ps_filt <- phyloseq::subset_samples(
  ps,
  !site_id %in% c("SLG", "LIB", "BSN") & !grepl("MOCK|CTR", site_id, ignore.case = TRUE)
)
ps_filt <- phyloseq::prune_taxa(phyloseq::taxa_sums(ps_filt) > 0, ps_filt)
message("Muestras para análisis: ", phyloseq::nsamples(ps_filt))

# Transformar a abundancias relativas (proporciones)
ps_rel <- phyloseq::transform_sample_counts(ps_filt, function(x) x / sum(x))

---

3. Diversidad alfa

La diversidad alfa (α) describe la riqueza y equitatividad de especies dentro de cada sitio. Calculamos los índices directamente sobre ps_filt (conteos crudos). Chao1 y ACE corrigen el sesgo por diferencias de profundidad de secuenciación sin necesidad de rarefactar.

Curva de acumulación de especies

La curva de acumulación muestra cómo aumenta el número de especies únicas detectadas conforme se incorporan más muestras. Si alcanza una asíntota, el esfuerzo de muestreo fue suficiente para capturar la mayor parte de la riqueza.

ps_sp <- phyloseq::tax_glom(ps_filt, "Species")

comm_sp <- as.data.frame(as.matrix(phyloseq::otu_table(ps_sp)))
if (phyloseq::taxa_are_rows(ps_sp)) comm_sp <- as.data.frame(t(as.matrix(comm_sp)))
comm_pa <- (comm_sp > 0) * 1

set.seed(123)
acc <- vegan::specaccum(comm_pa, method = "random", permutations = 1000)

acc_df <- data.frame(
  muestras = acc$sites,
  riqueza  = acc$richness,
  lower    = acc$richness - acc$sd,
  upper    = acc$richness + acc$sd
)

chao_est <- vegan::specpool(comm_pa)$chao

ggplot(acc_df, aes(x = muestras, y = riqueza)) +
  geom_ribbon(aes(ymin = lower, ymax = upper), fill = "#1b9e77", alpha = 0.25) +
  geom_line(color = "#1b9e77", linewidth = 1.1) +
  geom_point(color = "#1b9e77", size = 2) +
  geom_hline(yintercept = chao_est, linetype = "dashed", color = "#d95f02", linewidth = 0.9) +
  annotate("text", x = max(acc_df$muestras), y = chao_est,
           label = paste0("Chao = ", round(chao_est, 1)),
           hjust = 1, vjust = -0.7, size = 3.5, color = "#d95f02") +
  scale_x_continuous(expand = expansion(mult = c(0.01, 0.02))) +
  scale_y_continuous(expand = expansion(mult = c(0.02, 0.08))) +
  labs(title = "Curva de acumulacion de especies",
       subtitle = "Banda: +/- 1 SD entre 1,000 permutaciones aleatorias",
       x = "Numero de muestras", y = "Especies acumuladas") +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold"),
        plot.subtitle = element_text(color = "gray40"),
        axis.title = element_text(face = "bold"),
        panel.grid.minor = element_blank())

ggsave("final_data/reporte_ecologico/curva_acumulacion.png",
       width = 8, height = 5, dpi = 300)

Tip

Interpretación: Si la curva observada (verde) se aproxima al estimador Chao (naranja discontinuo), el muestreo capturó la mayor parte de la riqueza. Una brecha grande indica que podrían detectarse más especies con mayor esfuerzo.

Índices de diversidad

Calculamos seis índices complementarios:

Índice	Qué mide	¿Usa abundancia?	Rango
Observed	Número de ASVs detectadas (riqueza observada)	No — solo presencia/ausencia	0 – ∞
Chao1	Riqueza estimada incluyendo especies raras no detectadas, basado en singletons y doubletons	Parcial — usa singletons/doubletons	≥ Observed
ACE	Estimador de cobertura de abundancia; usa especies con <10 lecturas para estimar riqueza total	Parcial — distingue raras (<10) vs abundantes	≥ Observed
Shannon	Combina riqueza y equitatividad; sensible a especies raras	Sí — usa proporciones de cada especie	0 – ln(S)
Simpson	Probabilidad de que dos lecturas al azar sean de la misma especie (dominancia)	Sí — pesa más las especies abundantes	0 – 1 (1 = alta diversidad)
InvSimpson	Inverso de Simpson; se interpreta como el número efectivo de especies dominantes	Sí — derivado de Simpson	1 – S

Cálculo de índices y boxplots por región

Calculamos todos los índices con phyloseq::estimate_richness() y graficamos la distribución de valores por región del Golfo de California. La prueba de Kruskal-Wallis evalúa si hay diferencias significativas entre regiones.

alpha_df <- phyloseq::estimate_richness(
  ps_filt,
  measures = c("Observed", "Chao1", "ACE", "Shannon", "Simpson", "InvSimpson")
)
alpha_df$site_id   <- phyloseq::sample_data(ps_filt)$site_id
alpha_df$region_gc <- phyloseq::sample_data(ps_filt)$region_gc

# Reshape a formato largo para ggplot
alpha_long <- alpha_df %>%
  dplyr::select(site_id, region_gc, Observed, Chao1, ACE, Shannon, Simpson, InvSimpson) %>%
  tidyr::pivot_longer(
    cols      = -c(site_id, region_gc),
    names_to  = "index",
    values_to = "value"
  ) %>%
  dplyr::mutate(
    index = factor(index,
                   levels = c("Observed","Chao1","ACE","Shannon","Simpson","InvSimpson"))
  )

# Paleta de regiones (colorblind-friendly: RColorBrewer Set2)
regions     <- unique(alpha_long$region_gc)
pal_regions <- RColorBrewer::brewer.pal(max(3, length(regions)), "Set2")[seq_along(regions)]
names(pal_regions) <- regions

# Kruskal-Wallis por índice
kw_results <- alpha_long %>%
  dplyr::group_by(index) %>%
  dplyr::summarise(
    p_value = tryCatch(
      kruskal.test(value ~ region_gc)$p.value,
      error = function(e) NA_real_
    ),
    .groups = "drop"
  ) %>%
  dplyr::mutate(
    label = dplyr::case_when(
      p_value < 0.001 ~ "p < 0.001 ***",
      p_value < 0.01  ~ paste0("p = ", round(p_value, 3), " **"),
      p_value < 0.05  ~ paste0("p = ", round(p_value, 3), " *"),
      !is.na(p_value) ~ paste0("p = ", round(p_value, 3), " ns"),
      TRUE            ~ "p = NA"
    )
  )

# Unir resultados al dataframe largo para anotar
alpha_long <- alpha_long %>%
  dplyr::left_join(kw_results, by = "index")

# Calcular posición y del texto de anotación (máximo + 5%)
y_pos <- alpha_long %>%
  dplyr::group_by(index) %>%
  dplyr::summarise(y_ann = max(value, na.rm = TRUE) * 1.08, .groups = "drop")

alpha_long <- alpha_long %>% dplyr::left_join(y_pos, by = "index")

ggplot(alpha_long, aes(x = region_gc, y = value, fill = region_gc)) +
  geom_boxplot(alpha = 0.75, outlier.shape = 21, outlier.size = 2) +
  geom_jitter(aes(color = region_gc), width = 0.15, size = 2, alpha = 0.7) +
  geom_text(
    data = dplyr::distinct(alpha_long, index, label, y_ann),
    aes(x = 1.5, y = y_ann, label = label),
    inherit.aes = FALSE,
    size = 3, fontface = "italic", color = "gray30"
  ) +
  scale_fill_manual(values  = pal_regions) +
  scale_color_manual(values = pal_regions) +
  facet_wrap(~index, scales = "free_y", ncol = 3) +
  labs(
    title    = "Diversidad alfa por región del Golfo de California",
    subtitle = "Conteos crudos - Kruskal-Wallis entre regiones",
    x        = "Región",
    y        = "Valor del índice"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title      = element_text(face = "bold"),
    axis.title      = element_text(face = "bold"),
    strip.text      = element_text(face = "bold"),
    legend.position = "none",
    panel.grid.minor = element_blank()
  )

ggsave("final_data/reporte_ecologico/diversidad_alfa_regiones.png",
       width = 12, height = 6, dpi = 300)

Tip

Interpretación: Un índice de Shannon entre 1 y 3 es típico para comunidades de peces. Valores de Simpson cercanos a 1 indican alta diversidad (baja dominancia). Si el p-valor de Kruskal-Wallis es < 0.05, hay evidencia de diferencias significativas en diversidad alfa entre las regiones del Golfo.

# Tabla resumen por región (ASVs)
alpha_summary <- alpha_df %>%
  dplyr::group_by(region_gc) %>%
  dplyr::summarise(
    dplyr::across(
      c(Observed, Chao1, Shannon, InvSimpson),
      list(media = ~round(mean(.x, na.rm=TRUE), 2),
           sd    = ~round(sd(.x, na.rm=TRUE), 2)),
      .names = "{.col}_{.fn}"
    ),
    .groups = "drop"
  )

write.csv(alpha_summary, "final_data/reporte_ecologico/diversidad_alfa_asvs_resumen.csv", row.names = FALSE)
alpha_summary

Diversidad alfa a nivel de especie

El análisis anterior mide diversidad a nivel de ASVs — cada variante de secuencia única cuenta como una unidad. Pero múltiples ASVs pueden pertenecer a la misma especie (variantes intraespecíficas, alelos diferentes, haplotipos). Para obtener una medida más interpretable ecológicamente, colapsamos las ASVs por especie con tax_glom().

A nivel de especie usamos solo Observed, Shannon e InvSimpson — Chao1 y ACE no son aplicables aquí porque al colapsar se suman los conteos de todas las ASVs de cada especie, eliminando los singletons que estos estimadores necesitan para funcionar.

# Colapsar ASVs por especie (ps_sp ya existe de la sección anterior)
message("ASVs antes de colapsar: ", phyloseq::ntaxa(ps_filt))
message("Especies despues de colapsar: ", phyloseq::ntaxa(ps_sp))

# Solo indices aplicables tras tax_glom (sin Chao1/ACE)
alpha_sp <- suppressWarnings(phyloseq::estimate_richness(
  ps_sp,
  measures = c("Observed", "Shannon", "Simpson", "InvSimpson")
))
alpha_sp$site_id   <- phyloseq::sample_data(ps_sp)$site_id
alpha_sp$region_gc <- phyloseq::sample_data(ps_sp)$region_gc

# Reshape para graficar
alpha_sp_long <- alpha_sp %>%
  dplyr::select(site_id, region_gc, Observed, Shannon, Simpson, InvSimpson) %>%
  tidyr::pivot_longer(
    cols      = -c(site_id, region_gc),
    names_to  = "index",
    values_to = "value"
  ) %>%
  dplyr::mutate(
    index = factor(index, levels = c("Observed", "Shannon", "Simpson", "InvSimpson"))
  )

# Kruskal-Wallis por índice
kw_sp <- alpha_sp_long %>%
  dplyr::group_by(index) %>%
  dplyr::summarise(
    p_value = tryCatch(kruskal.test(value ~ region_gc)$p.value, error = function(e) NA_real_),
    .groups = "drop"
  ) %>%
  dplyr::mutate(
    label = dplyr::case_when(
      p_value < 0.001 ~ "p < 0.001 ***",
      p_value < 0.01  ~ paste0("p = ", round(p_value, 3), " **"),
      p_value < 0.05  ~ paste0("p = ", round(p_value, 3), " *"),
      !is.na(p_value) ~ paste0("p = ", round(p_value, 3), " ns"),
      TRUE            ~ "p = NA"
    )
  )

alpha_sp_long <- alpha_sp_long %>% dplyr::left_join(kw_sp, by = "index")
y_pos_sp <- alpha_sp_long %>%
  dplyr::group_by(index) %>%
  dplyr::summarise(y_ann = max(value, na.rm = TRUE) * 1.08, .groups = "drop")
alpha_sp_long <- alpha_sp_long %>% dplyr::left_join(y_pos_sp, by = "index")

ggplot(alpha_sp_long, aes(x = region_gc, y = value, fill = region_gc)) +
  geom_boxplot(alpha = 0.75, outlier.shape = 21, outlier.size = 2) +
  geom_jitter(aes(color = region_gc), width = 0.15, size = 2, alpha = 0.7) +
  geom_text(
    data = dplyr::distinct(alpha_sp_long, index, label, y_ann),
    aes(x = 1.5, y = y_ann, label = label),
    inherit.aes = FALSE,
    size = 3, fontface = "italic", color = "gray30"
  ) +
  scale_fill_manual(values  = pal_regions) +
  scale_color_manual(values = pal_regions) +
  facet_wrap(~index, scales = "free_y", ncol = 3) +
  labs(
    title    = "Diversidad alfa por region - nivel de especie",
    subtitle = "Colapsado por especie - Kruskal-Wallis entre regiones",
    x        = "Región",
    y        = "Valor del índice"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title      = element_text(face = "bold"),
    axis.title      = element_text(face = "bold"),
    strip.text      = element_text(face = "bold"),
    legend.position = "none",
    panel.grid.minor = element_blank()
  )

ggsave("final_data/reporte_ecologico/diversidad_alfa_especies.png",
       width = 12, height = 6, dpi = 300)

# Tabla resumen por región (especies)
alpha_sp_summary <- alpha_sp %>%
  dplyr::group_by(region_gc) %>%
  dplyr::summarise(
    dplyr::across(
      c(Observed, Shannon, InvSimpson),
      list(media = ~round(mean(.x, na.rm=TRUE), 2),
           sd    = ~round(sd(.x, na.rm=TRUE), 2)),
      .names = "{.col}_{.fn}"
    ),
    .groups = "drop"
  )

write.csv(alpha_sp_summary, "final_data/reporte_ecologico/diversidad_alfa_especies_resumen.csv", row.names = FALSE)
alpha_sp_summary

Nota

¿Qué esperamos ver? El número de especies (“Observed”) será menor que el de ASVs, porque varias ASVs se colapsan en una sola especie. Sin embargo, los patrones de diversidad entre regiones deberían mantenerse — si la región Central es más diversa a nivel de ASVs, también debería serlo a nivel de especie. Si los patrones cambian, puede indicar que hay muchas ASVs sin asignar a nivel de especie que están inflando la diversidad en algunas regiones.

---

4. Composición taxonómica

Las barras apiladas muestran la abundancia relativa por especie en cada sitio. Son la visualización central del metabarcoding: qué detectamos y en qué proporción.

# Colapsar y calcular proporciones (ps_sp creado en curva de acumulacion)
ps_sp_prop <- phyloseq::transform_sample_counts(ps_sp, function(x) x / sum(x))

Tabla de abundancias por especie y muestra

La tabla de abundancias por especie y muestra es el producto central del análisis: registra cuántas lecturas corresponden a cada especie en cada sitio. Se exporta en CSV para su uso en otros programas.

ASV_df <- as.data.frame(phyloseq::otu_table(ps_sp)) %>%
  tibble::rownames_to_column(var = "seqID") %>%
  tidyr::pivot_longer(cols = -seqID, names_to = "Sample_name", values_to = "nReads")

TAXA_df <- as.data.frame(phyloseq::tax_table(ps_sp)) %>%
  tibble::rownames_to_column(var = "seqID")

# Mapeo de Sample_name a site_id
name_map <- data.frame(
  Sample_name = phyloseq::sample_names(ps_sp),
  site_id     = phyloseq::sample_data(ps_sp)$site_id,
  stringsAsFactors = FALSE
)

species.site.reads <- ASV_df %>%
  dplyr::left_join(TAXA_df, by = "seqID") %>%
  dplyr::left_join(name_map, by = "Sample_name") %>%
  dplyr::filter(!is.na(Species)) %>%
  dplyr::group_by(site_id, Species) %>%
  dplyr::summarise(nReads = sum(nReads), .groups = "drop") %>%
  tidyr::pivot_wider(names_from = site_id, values_from = nReads, values_fill = 0)

write.csv(species.site.reads, ruta_csv_output, row.names = FALSE)

species.site.reads[1:10, 1:3] %>%
  kable(format = "html", caption = "Primeros 10 renglones de la tabla de abundancias", align = "c") %>%
  kable_styling(full_width = TRUE, font_size = 10)

Barras apiladas por especie

Abundancia relativa de cada especie curada geográficamente en cada sitio del Golfo de California.

# Top N especies por abundancia total; el resto se agrupa en "Otros"
top_n <- 15

df_sp <- phyloseq::psmelt(ps_sp_prop) %>%
  dplyr::filter(site_id %in% md$site_id)

# Ranking de especies por abundancia media global
sp_rank <- df_sp %>%
  dplyr::group_by(Species) %>%
  dplyr::summarise(mean_abund = mean(Abundance, na.rm = TRUE), .groups = "drop") %>%
  dplyr::arrange(dplyr::desc(mean_abund))

top_species <- sp_rank$Species[1:min(top_n, nrow(sp_rank))]
n_otros <- sum(!df_sp$Species %in% top_species & !is.na(df_sp$Species))
message("Top ", top_n, " especies mostradas; ",
        length(unique(df_sp$Species)) - length(top_species),
        " especies agrupadas en 'Otros'")

# Reclasificar las menos abundantes como "Otros"
df_sp <- df_sp %>%
  dplyr::mutate(
    Species_plot = dplyr::if_else(Species %in% top_species, Species, "Otros"),
    Species_plot = factor(Species_plot, levels = c("Otros", rev(top_species)))
  ) %>%
  dplyr::group_by(site_id, Species_plot) %>%
  dplyr::summarise(Abundance = sum(Abundance), .groups = "drop") %>%
  dplyr::mutate(site_id = factor(site_id, levels = unique(md$site_id)))

# Paleta: colores distintos para top especies + gris para "Otros"
n_top <- length(top_species)
base_colors <- c(
  RColorBrewer::brewer.pal(12, "Paired"),
  RColorBrewer::brewer.pal(8, "Set2"),
  RColorBrewer::brewer.pal(8, "Dark2")
)
pal_sp <- c(base_colors[1:n_top], "Otros" = "#D3D3D3")
names(pal_sp)[1:n_top] <- top_species

gp_sp <- ggplot(df_sp, aes(x = site_id, y = Abundance, fill = Species_plot)) +
  geom_bar(stat = "identity") +
  scale_fill_manual(values = pal_sp) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  labs(
    title = paste0("Abundancia relativa por especie (Top ", top_n, " + Otros)"),
    x     = "Sitio",
    y     = "Abundancia relativa (%)",
    fill  = "Especie"
  ) +
  theme_linedraw(base_size = 11) +
  theme(
    axis.text.x      = element_text(angle = 90, vjust = 0.95, hjust = 1, size = 9),
    legend.position  = "bottom",
    legend.key.size  = unit(0.35, "cm"),
    legend.text      = element_text(face = "italic", size = 7),
    legend.title     = element_text(face = "bold"),
    plot.title       = element_text(face = "bold"),
    axis.title       = element_text(face = "bold")
  ) +
  guides(fill = guide_legend(nrow = ceiling((n_top + 1) / 4), byrow = TRUE, title.position = "top"))

print(gp_sp)
ggsave(
  file.path(ruta_figs, "abundancia_relativa_especie.png"),
  gp_sp, width = 12, height = 9, dpi = 300
)

Heatmap: Top 20 especies × sitios

El heatmap complementa los gráficos de barras mostrando simultáneamente la distribución de las 20 especies más abundantes en todos los sitios. Los valores están transformados con log(x+1) para reducir el efecto de species muy dominantes y hacer visibles las abundancias bajas. El clustering jerárquico (Ward D2) agrupa automáticamente sitios y especies con patrones similares.

# Top 20 especies por abundancia total
top20_names <- names(sort(phyloseq::taxa_sums(ps_sp), decreasing = TRUE))[1:20]
ps_top20    <- phyloseq::prune_taxa(top20_names, ps_sp)

# Construir matriz sites × species
mat_heat <- log1p(as.data.frame(as.matrix(phyloseq::otu_table(ps_top20))))
if (phyloseq::taxa_are_rows(ps_top20)) mat_heat <- as.data.frame(t(as.matrix(mat_heat)))

# Etiquetas: site_id como filas, Species como columnas
rownames(mat_heat) <- phyloseq::sample_data(ps_top20)$site_id
colnames(mat_heat) <- phyloseq::tax_table(ps_top20)[, "Species"]

png("final_data/reporte_ecologico/heatmap_top20_especies.png",
    width = 2400, height = 1800, res = 300)
pheatmap::pheatmap(
  t(mat_heat),
  scale              = "none",
  color              = viridis::viridis(50, option = "D"),
  clustering_method  = "ward.D2",
  fontsize_row       = 9,
  fontsize_col       = 9,
  main               = "Heatmap Top 20 especies — log(x+1)",
  border_color       = NA,
  cellwidth          = 18,
  cellheight         = 14
)
dev.off()

Tip

Interpretación: El clustering de columnas (sitios) muestra qué sitios tienen comunidades similares. El clustering de filas (especies) muestra qué especies co-ocurren. Bloques de color intenso en el mismo cluster de sitios indican asociaciones entre región geográfica y composición de especies.

---

5. Diversidad beta

La diversidad beta (β) mide las diferencias en composición de comunidades entre sitios. Usamos Bray-Curtis sobre abundancias relativas a nivel de especie — sensible a cambios en abundancia de especies dominantes.

Flujo: betadisper → PERMANOVA + ANOSIM → NMDS.

Preparación de matrices

# Abundancias relativas por especie (ps_sp creado en curva de acumulacion)
ps_sp_rel <- phyloseq::transform_sample_counts(ps_sp, function(x) x / sum(x))

comm_sp_rel <- as.data.frame(as.matrix(phyloseq::otu_table(ps_sp_rel)))
if (phyloseq::taxa_are_rows(ps_sp_rel)) comm_sp_rel <- as.data.frame(t(as.matrix(comm_sp_rel)))

meta_sp <- data.frame(phyloseq::sample_data(ps_sp_rel), stringsAsFactors = FALSE)
meta_sp$sample_id <- rownames(meta_sp)

samps_keep <- rownames(comm_sp_rel)
meta_filt  <- meta_sp[meta_sp$sample_id %in% samps_keep, ]
comm_filt  <- comm_sp_rel[samps_keep, , drop = FALSE]
rownames(comm_filt) <- meta_filt$site_id

dist_sp_bc <- vegan::vegdist(comm_filt, method = "bray")

message("Muestras en analisis beta: ", nrow(comm_filt))

Homogeneidad de dispersiones (betadisper) — Bray-Curtis

Antes de aplicar PERMANOVA es obligatorio verificar la homogeneidad de varianzas entre grupos. Si los grupos tienen dispersiones muy distintas, un PERMANOVA significativo puede reflejar esa diferencia de dispersión y no necesariamente diferencias en la posición de los centroides (composición media).

bd      <- vegan::betadisper(dist_sp_bc, meta_filt$region_gc)
bd_test <- vegan::permutest(bd, permutations = 999)

cat("=== Test de homogeneidad de dispersiones (betadisper — Bray-Curtis) ===\n")
print(bd_test)

# Visualizar dispersiones
bd_df <- data.frame(
  site_id   = rownames(bd$vectors),
  region_gc = meta_filt$region_gc,
  Distancia = bd$distances
)

ggplot(bd_df, aes(x = region_gc, y = Distancia, fill = region_gc)) +
  geom_boxplot(alpha = 0.7, outlier.shape = 21) +
  geom_jitter(aes(color = region_gc), width = 0.15, size = 2.5, alpha = 0.8) +
  scale_fill_manual(values  = pal_regions) +
  scale_color_manual(values = pal_regions) +
  annotate(
    "text",
    x     = 1.5,
    y     = max(bd_df$Distancia) * 1.05,
    label = paste0("betadisper p = ", round(bd_test$tab$`Pr(>F)`[1], 3)),
    size  = 3.8, fontface = "italic", color = "gray30"
  ) +
  labs(
    title    = "Homogeneidad de dispersiones por región (betadisper)",
    subtitle = "Distancia al centroide de cada grupo — Bray-Curtis",
    x        = "Región",
    y        = "Distancia al centroide"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title      = element_text(face = "bold"),
    axis.title      = element_text(face = "bold"),
    legend.position = "none"
  )

ggsave("final_data/reporte_ecologico/betadisper_regiones.png",
       width = 8, height = 5, dpi = 300)

Nota

Recuerda: Si el p-valor de betadisper es < 0.05, las dispersiones son heterogéneas entre grupos. En ese caso, interpreta el PERMANOVA con precaución: la significancia puede deberse a diferencias en dispersión, no solo en composición media.

PERMANOVA y ANOSIM

permanova_sp_bc <- vegan::adonis2(
  comm_filt ~ region_gc,
  data         = data.frame(meta_filt, stringsAsFactors = FALSE),
  permutations = 999,
  method       = "bray"
)
anosim_sp_bc <- vegan::anosim(dist_sp_bc, meta_filt$region_gc)

as.data.frame(permanova_sp_bc) %>%
  tibble::rownames_to_column("Fuente") %>%
  kable(format = "html", digits = 4, align = "c",
        caption = "PERMANOVA — Bray-Curtis, 999 permutaciones") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = TRUE, font_size = 10)

cat("ANOSIM R =", round(anosim_sp_bc$statistic, 3),
    "| p =", round(anosim_sp_bc$signif, 4), "\n")

NMDS — Bray-Curtis

El NMDS (Non-metric Multidimensional Scaling) representa las distancias entre comunidades en 2D preservando al máximo las relaciones originales. Un stress < 0.1 es excelente; < 0.2 es aceptable. Usamos Bray-Curtis como métrica principal.

set.seed(123)
nmds_sp_bc <- vegan::metaMDS(dist_sp_bc, k = 2, trymax = 500, trace = FALSE)

nmds_sp_bc_df <- as.data.frame(nmds_sp_bc$points) %>%
  tibble::rownames_to_column("site_id") %>%
  dplyr::left_join(
    meta_filt %>% dplyr::select(site_id, region_gc),
    by = "site_id"
  ) %>%
  dplyr::mutate(region_gc = as.factor(region_gc))

ggplot(nmds_sp_bc_df, aes(x = MDS1, y = MDS2, color = region_gc)) +
  geom_point(size = 3.5, alpha = 0.9) +
  geom_text(aes(label = site_id), vjust = -1, size = 3.2, fontface = "bold") +
  scale_color_manual(values = pal_regions) +
  labs(
    title    = "NMDS basado en especies (Bray-Curtis)",
    subtitle = "Métrica principal — sitios coloreados por región",
    x        = "NMDS1",
    y        = "NMDS2",
    color    = "Región"
  ) +
  annotate("text", x = -Inf, y = -Inf,
           label = paste0("Stress = ", round(nmds_sp_bc$stress, 3)),
           hjust = -0.1, vjust = -0.5, size = 3.8, color = "gray30") +
  annotate("text", x = -Inf, y = -Inf,
           label = paste0("PERMANOVA R² = ", round(permanova_sp_bc$R2[1], 2),
                          ", p = ", signif(permanova_sp_bc$`Pr(>F)`[1], 3)),
           hjust = -0.1, vjust = -2.2, size = 3.5, color = "#d95f02") +
  annotate("text", x = -Inf, y = -Inf,
           label = paste0("ANOSIM R = ", round(anosim_sp_bc$statistic, 2),
                          ", p = ", signif(anosim_sp_bc$signif, 3)),
           hjust = -0.1, vjust = -3.9, size = 3.5, color = "#1b9e77") +
  theme_minimal(base_size = 12) +
  theme(
    plot.title       = element_text(face = "bold"),
    plot.subtitle    = element_text(color = "gray40"),
    axis.title       = element_text(face = "bold"),
    panel.grid.minor = element_blank(),
    legend.position  = "bottom"
  )

ggsave("final_data/reporte_ecologico/nmds_especies_braycurtis.png",
       width = 9, height = 7, dpi = 300)

Tip

Interpretación: Sitios agrupados en el NMDS tienen composiciones de especies similares. Si el agrupamiento coincide con las regiones del Golfo, hay estructuración espacial de las comunidades.

---

6. Control positivo — MOCK

La muestra MOCK es un control positivo con ADN de especies conocidas. Permite validar el pipeline: ¿detectamos las especies esperadas? ¿aparecen especies inesperadas? Se aplica la misma curación geográfica que a las muestras ambientales, por lo que especies no documentadas en el GC pueden perder su asignación a especie.

# Detección dinámica del MOCK (usa ps_ctrl, creado al separar controles)
mock_meta <- data.frame(phyloseq::sample_data(ps_ctrl), stringsAsFactors = FALSE)
mock_ids <- mock_meta$site_id[grepl("MOCK", mock_meta$site_id, ignore.case = TRUE)]

if (length(mock_ids) > 0) {
  mock_sid <- mock_ids[1]
  message("Muestra MOCK (site_id): ", mock_sid)

  # Construir tabla de abundancias por especie desde ps_ctrl
  ps_mock <- phyloseq::subset_samples(ps_ctrl, site_id == mock_sid)
  ps_mock_sp <- phyloseq::tax_glom(ps_mock, "Species")
  mock_counts <- as.numeric(phyloseq::otu_table(ps_mock_sp))
  mock_taxa   <- as.data.frame(phyloseq::tax_table(ps_mock_sp), stringsAsFactors = FALSE)

  df_mock <- data.frame(
    Species = mock_taxa$Species,
    reads   = mock_counts,
    stringsAsFactors = FALSE
  ) %>%
    dplyr::filter(reads > 0) %>%
    dplyr::mutate(rel_abund = reads / sum(reads)) %>%
    dplyr::arrange(dplyr::desc(rel_abund))

  df_mock$Species <- factor(df_mock$Species, levels = df_mock$Species)
  n_mock <- nrow(df_mock)

  # Paleta
  if (n_mock <= 12) {
    pal_mock <- RColorBrewer::brewer.pal(max(3, n_mock), "Paired")[1:n_mock]
  } else {
    pal_mock <- c(RColorBrewer::brewer.pal(12, "Paired"),
                  viridis::viridis(n_mock - 12, option = "D"))
  }
  names(pal_mock) <- levels(df_mock$Species)

  # Gráfico
  print(
    ggplot(df_mock, aes(x = "MOCK", y = rel_abund, fill = Species)) +
      geom_bar(stat = "identity", width = 0.5) +
      scale_fill_manual(values = pal_mock) +
      scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
      labs(title = paste0("Control positivo: MOCK (", n_mock, " especies)"),
           x = NULL, y = "Abundancia relativa (%)") +
      theme_linedraw(base_size = 13) +
      theme(legend.position = "right", legend.key.size = unit(0.5, "cm"),
            legend.text = element_text(face = "italic", size = 9),
            legend.title = element_text(face = "bold"),
            plot.title = element_text(face = "bold"),
            axis.title = element_text(face = "bold")) +
      guides(fill = guide_legend(ncol = 2))
  )

  ggsave("final_data/reporte_ecologico/control_mock.png",
         width = 8, height = 6, dpi = 300)

  # Tabla
  df_mock %>%
    dplyr::mutate(rel_abund = scales::percent(rel_abund, accuracy = 0.1)) %>%
    dplyr::rename(Especie = Species, Lecturas = reads,
                  `Abundancia relativa` = rel_abund) %>%
    kable(format = "html", align = "c",
          caption = paste0("Especies detectadas en MOCK (", mock_col, ")")) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                  full_width = TRUE, font_size = 10)
} else {
  message("No se detectó ninguna muestra MOCK en este dataset.")
}

Nota

Un MOCK bien resuelto contiene únicamente las especies de la mezcla, en proporciones similares a las esperadas. Especies adicionales pueden indicar contaminación cruzada (cross-talk); distorsión de proporciones puede reflejar sesgos de amplificación.

---

Información de sesión

sessionInfo()