Bioconductor

Paquetes esenciales para análisis ómico en R

r
bioconductor
rna-seq
microbiome
sequences
differential-expression
Referencia comentada de los paquetes Bioconductor que estructuran el flujo de análisis genómico en R: estructura de datos, manipulación de secuencias, expresión diferencial y microbioma.

Sobre Bioconductor

Bioconductor es el ecosistema de paquetes R orientado al análisis de datos biomédicos. Comparte infraestructura con CRAN pero opera con un release cycle propio (dos releases anuales, en abril y octubre, sincronizados con la versión de R). Esto garantiza coherencia entre paquetes pero obliga a fijar la versión correcta cuando reproduces análisis ajenos.

Instalación canónica:

if (!require("BiocManager", quietly = TRUE))
  install.packages("BiocManager")

BiocManager::install("DESeq2")          # un paquete
BiocManager::install(version = "3.20")  # fijar versión del release
BiocManager::valid()                    # verificar coherencia de versiones

Tres principios del ecosistema que conviene tener interiorizados:

  • Clases S4 compartidas. La mayoría de paquetes consume y devuelve SummarizedExperiment o derivados (DESeqDataSet, SingleCellExperiment, TreeSummarizedExperiment). Esto permite encadenar herramientas sin pegar datos a mano.
  • Vignettes como documentación principal. vignette("DESeq2") abre el manual canónico de cada paquete. En la práctica casi siempre es preferible a ?función.
  • Integración con GenomicRanges. Las coordenadas genómicas se manejan como objetos GRanges, no como data.frame con chr / start / end.

Esta página cataloga seis paquetes que estructuran la mayor parte de los flujos ómicos en R. El orden refleja jerarquía conceptual: primero la infraestructura (estructura de datos y secuencias), después los modelos de expresión diferencial, y finalmente la aplicación al dominio del microbioma.

SummarizedExperiment

SummarizedExperiment es la clase S4 que estandariza la representación de un experimento en Bioconductor: una matriz de assay (típicamente features × muestras), metadatos de filas (rowData) y de columnas (colData), todo en un objeto coherente.

Es la base sobre la que se construyen casi todos los workflows modernos. DESeqDataSet, SingleCellExperiment, TreeSummarizedExperiment y MultiAssayExperiment heredan o extienden esta clase. Internalizar su API ahorra horas a lo largo de cualquier proyecto serio de Bioconductor.

Cuándo usarlo

Siempre que tengas una matriz de medidas (counts, intensidades, abundancias) asociada a metadatos de muestra y/o de feature. Si tu output natural es un data.frame largo, casi siempre conviene reestructurar a SummarizedExperiment antes de seguir: es el formato que el resto del ecosistema espera.

Conceptos clave

  • assay(se) accede a la matriz principal. assays(se) si hay varias capas en paralelo (counts, normalized, logcounts…).
  • colData(se) son los metadatos de las muestras (DataFrame S4, no data.frame base).
  • rowData(se) son los metadatos de las features (genes, taxones, peaks).
  • Garantía de alineación: el subsetting (se[1:100, 1:5]) mantiene la coherencia entre assay y metadatos automáticamente, este es el valor real frente a hacerlo con matrices sueltas.
  • RangedSummarizedExperiment extiende la clase añadiendo coordenadas genómicas en rowRanges.

Patrón mínimo

library(SummarizedExperiment)

counts <- matrix(rpois(200, 10), nrow = 20)
rownames(counts) <- paste0("gene", 1:20)
colnames(counts) <- paste0("sample", 1:10)

col_meta <- DataFrame(
  condition = rep(c("ctrl", "treat"), each = 5),
  batch     = rep(c("A", "B"), 5)
)

se <- SummarizedExperiment(
  assays  = list(counts = counts),
  colData = col_meta
)

# Subset preservando la coherencia
se_treated <- se[, se$condition == "treat"]

Trampas habituales

  • No es data.frame. se$condition funciona como atajo para colData(se)$condition, pero el colData subyacente es DataFrame (S4): operaciones de dplyr no actúan directamente sobre él.
  • Múltiples assays. Si necesitas conviven counts y logcounts, usa assays(se) <- list(counts = ..., logcounts = ...) en lugar de duplicar objetos.
  • assayNames(se) es lo primero que se inspecciona cuando trabajas con código de terceros. Sin ese hábito, te encontrarás multiplicando ramas if.

Enlaces

Biostrings

Biostrings es la librería que representa eficientemente secuencias biológicas (ADN, ARN, proteína) en R. Define las clases DNAString, RNAString, AAString y sus colecciones (DNAStringSet, RNAStringSet, AAStringSet), con búsqueda de patrones, alineamiento básico y operaciones idiomáticas sobre secuencias.

Es la pieza más baja del stack de secuencias en Bioconductor. ShortRead (FASTQ), Rsamtools (BAM) y GenomicAlignments se apoyan en sus clases.

Cuándo usarlo

Manipulación directa de secuencias en R: leer FASTA, búsqueda de motivos, complemento reverso, traducción, k-mers. Cualquier preprocesamiento que no quieras delegar a herramientas externas como seqkit, bedtools o herramientas de alineamiento.

Cuándo NO usarlo

Análisis de lecturas a gran escala (decenas de millones de reads). Ahí entran herramientas fuera de R (BWA, minimap2, STAR) o paquetes que iteran sobre el FASTQ por bloques (ShortRead). En memoria de R, mantén las colecciones por debajo de unos pocos millones de elementos.

Conceptos clave

  • DNAStringSet es la colección eficiente. DNAString individual rara vez se usa solo.
  • API básica: width(seqs), subseq(seqs, ...), reverseComplement(seqs), translate(seqs).
  • Búsqueda: matchPattern() para una secuencia diana, vmatchPattern() para vectorizar contra una colección.
  • Soporta códigos IUPAC de ambigüedad y operaciones case-sensitive opcionales.
  • letterFrequency() y alphabetFrequency() para composición y conteos de k-meros.

Patrón mínimo

library(Biostrings)

# Leer FASTA
seqs <- readDNAStringSet("genome.fasta")

# Estadísticos básicos
width(seqs)
alphabetFrequency(seqs, baseOnly = TRUE)

# Búsqueda de motivo en toda la colección
hits <- vmatchPattern("GAATTC", seqs)        # diana EcoRI

# Complemento reverso y traducción a aminoácidos (frame 1)
rc       <- reverseComplement(seqs)
proteins <- translate(seqs, if.fuzzy.codon = "solve")

Trampas habituales

  • as.character() es caro. Si vas a iterar sobre una colección, evita convertir a cadena R. Opera con la API S4.
  • matchPattern vs vmatchPattern. El primero acepta una sola secuencia diana. El segundo, una colección. Confundirlos da errores opacos.
  • Marcos de lectura. translate() asume frame 1 por defecto. Para los 6 marcos hay que combinar subseq() + reverseComplement() manualmente.

Enlaces

DESeq2

DESeq2 es el estándar de facto para análisis de expresión diferencial en datos de RNA-seq bulk. Modela los conteos crudos con una distribución binomial negativa, estima la dispersión por gen con shrinkage hacia una curva tendencial, y aplica MAP shrinkage al log-fold-change para estabilizar las estimaciones en genes de baja expresión.

Desarrollado por Michael Love y Wolfgang Huber. Pareja natural con tximport para importar resultados de Salmon, Kallisto o RSEM a conteos a nivel de gen.

Cuándo usarlo

Bulk RNA-seq con réplicas (mínimo 3 por condición). Diseños factoriales con interacciones. Datasets de hasta ~100 muestras. Más allá, limma-voom es notablemente más rápido y suele ser preferible.

Cuándo NO usarlo

  • Single-cell. No está diseñado para zero-inflation extremo. Usa glmGamPoi, MAST o el patrón de pseudobulk + DESeq2 sobre el agregado.
  • Datos ya normalizados. DESeq2 exige conteos crudos enteros. Pasarle TPM, RPKM o FPKM rompe los supuestos.
  • Una sola réplica por condición. El modelo necesita estimar dispersión. N = 1 no es viable.

Conceptos clave

  • DESeqDataSet hereda de SummarizedExperiment. El design (fórmula) y el assay "counts" son obligatorios.
  • DESeq() ejecuta tres pasos en uno: estimateSizeFactors(), estimateDispersions(), nbinomWaldTest().
  • results() extrae la tabla. El contraste se especifica con contrast = c("variable", "nivel_num", "nivel_denom") o con name = "...".
  • lfcShrink() (apeglm, ashr o normal) corrige el LFC para ranking y visualización. Imprescindible para MA-plots, heatmaps y gene set enrichment downstream.
  • Independent filtering. results() filtra automáticamente genes de bajo conteo para mejorar el poder estadístico. Ese filtro aparece reportado en summary(res).

Patrón mínimo

library(DESeq2)

# Suponiendo `counts` (matriz integer genes × samples) y `coldata` (DataFrame)
coldata$condition <- factor(coldata$condition, levels = c("ctrl", "treat"))

dds <- DESeqDataSetFromMatrix(
  countData = counts,
  colData   = coldata,
  design    = ~ batch + condition   # variable de interés al final
)

# Prefiltrado opcional (acelera; no es obligatorio)
keep <- rowSums(counts(dds) >= 10) >= 3
dds  <- dds[keep, ]

# Modelo
dds <- DESeq(dds)

# Resultados con shrinkage del LFC
res <- lfcShrink(dds, coef = "condition_treat_vs_ctrl", type = "apeglm")
summary(res)

Trampas habituales

  • Orden de niveles del factor. El nivel de referencia (típicamente ctrl) debe ser el primero. Fijarlo explícitamente con factor(..., levels = ...) evita sorpresas en el signo del LFC.
  • Variable de interés al final del design. Convención de DESeq2: cuando no se pasa contrast ni name, results() devuelve el contraste del último término. Coloca siempre la variable principal al final.
  • Filtrar antes vs después. Filtrar genes con counts muy bajos antes de DESeq() es opcional. El independent filtering dentro de results() lo hace después, basado en baseMean. No filtres con criterios sesgados (p. ej. varianza alta), invalidan la corrección de p-values.
  • Batch effects en design, no en los datos. Incluirlos como covariable (~ batch + condition) es correcto. Aplicar removeBatchEffect() y pasar esos counts al modelo no lo es.

Enlaces

Relacionados en esta página

edgeR

edgeR resuelve el mismo problema que DESeq2, expresión diferencial sobre conteos, pero con una filosofía empírica Bayes ligeramente distinta: moderación de la dispersión por gen mediante estimación trended / tagwise, y test cuasi-likelihood F (QL F-test) en la versión moderna.

Desarrollado por el grupo de Gordon Smyth en WEHI. Suele dar resultados muy similares a DESeq2 en términos prácticos. Las diferencias importan en diseños complejos y en datasets de tamaño medio-alto donde el test QL F rinde mejor que el Wald de DESeq2.

Cuándo usarlo

  • Mismo nicho que DESeq2 (bulk RNA-seq con réplicas).
  • Diseños factoriales complejos con interacciones y contrastes múltiples, el QL F-test es generalmente más conservador y robusto.
  • Cuando ya estás integrado en el ecosistema limma para visualización o pasos posteriores.

Conceptos clave

  • DGEList es la estructura propia (no hereda de SummarizedExperiment, aunque hay conversión bidireccional).
  • Pipeline moderno recomendado: calcNormFactors() (TMM) → estimateDisp() → glmQLFit() → glmQLFTest().
  • La función legacy exactTest() solo aplica a comparaciones de dos grupos. Evítala en diseños reales con más de una covariable.
  • TMM (Trimmed Mean of M-values) ajusta por composición. Es diferente a la mediana de ratios de DESeq2 pero filosóficamente equivalente, ambos intentan corregir el efecto de unos pocos genes muy expresados sobre el tamaño efectivo de la librería.

Patrón mínimo

library(edgeR)

y <- DGEList(counts = counts, group = coldata$condition)

# Filtrado recomendado por el propio paquete
keep <- filterByExpr(y, group = coldata$condition)
y    <- y[keep, , keep.lib.sizes = FALSE]

# Normalización por composición
y <- calcNormFactors(y, method = "TMM")

# Modelo
design <- model.matrix(~ batch + condition, data = coldata)
y      <- estimateDisp(y, design)
fit    <- glmQLFit(y, design)
qlf    <- glmQLFTest(fit, coef = "conditiontreat")

topTags(qlf)

Trampas habituales

  • exactTest() no escala a diseños complejos. Si tu modelo tiene más de una covariable, usa glmQLFit + glmQLFTest.
  • filterByExpr() ya considera el design. Diseñado para filtrar de forma estadísticamente válida. No filtres a mano sin entender qué hace.
  • TMM no produce conteos normalizados que se pasen a otra herramienta, produce factores de normalización que el modelo aplica internamente. Si necesitas exportar conteos normalizados para visualización, usa cpm(y, normalized.lib.sizes = TRUE).

Enlaces

Relacionados en esta página

  • DESeq2, alternativa con MAP shrinkage del LFC.
  • limma, del mismo grupo (WEHI), preferible para n grande.

limma

limma es la librería original de modelos lineales aplicados a expresión génica, primero microarrays, después adaptada a RNA-seq vía la transformación voom (variance modeling at the observational level). El combo limma + voom es la opción natural cuando tienes muchas muestras o diseños complejos: rinde más rápido que DESeq2 / edgeR y aguanta contrastes sofisticados (makeContrasts, eBayes) con elegancia.

Es además la pieza con la que se hacen los contrastes “raros”: correlaciones within-subject (duplicateCorrelation), ANOVA factorial con interacciones, ajustes por confounders, modelos con efectos aleatorios moderados.

Cuándo usarlo

  • n alto (>20 muestras) o muchas condiciones. La transformación voom convierte counts en log-CPM con pesos de precisión y permite aplicar la maquinaria lineal clásica.
  • Microarrays (caso histórico, aún común en datasets antiguos o plataformas Affymetrix HTA).
  • Diseños donde necesitas makeContrasts con múltiples comparaciones nombradas y simétricas.

Conceptos clave

  • voom(counts, design) produce un EList con log-CPM y matriz de pesos. A partir de ahí, lmFit + contrasts.fit + eBayes es el pipeline estándar.
  • eBayes() modera el error estándar de los t-stats mediante empirical Bayes, comparte información entre genes para estabilizar los tests.
  • topTable() extrae la tabla de resultados (equivalente a results() en DESeq2 o topTags() en edgeR).
  • Para correlaciones dentro de sujeto (datos pareados), usa duplicateCorrelation() antes de lmFit.

Patrón mínimo

library(limma)
library(edgeR)   # para DGEList y filterByExpr

y    <- DGEList(counts = counts)
keep <- filterByExpr(y, group = coldata$condition)
y    <- y[keep, , keep.lib.sizes = FALSE]
y    <- calcNormFactors(y)

design <- model.matrix(~ 0 + condition + batch, data = coldata)
colnames(design) <- gsub("condition", "", colnames(design))

v   <- voom(y, design, plot = FALSE)
fit <- lmFit(v, design)

contrasts <- makeContrasts(treat_vs_ctrl = treat - ctrl, levels = design)
fit2 <- contrasts.fit(fit, contrasts)
fit2 <- eBayes(fit2)

topTable(fit2, coef = "treat_vs_ctrl", number = 20)

Trampas habituales

  • No mezcles voom con datos ya normalizados. La entrada debe ser counts crudos integrados en DGEList + calcNormFactors. Pasarle TPM o log-counts rompe el modelo de pesos.
  • ~ condition (con intercept) vs ~ 0 + condition (sin). El segundo facilita makeContrasts con nombres claros (treat - ctrl). El primero contrasta automáticamente contra el nivel de referencia. Para diseños factoriales, el patrón sin intercept es más legible.
  • eBayes(trend = TRUE) activa el modo limma-trend, alternativa a voom cuando la varianza no depende fuertemente del nivel de expresión. Útil con counts ya transformados (p. ej. CPM logarítmicos) cuando no quieres pesos.

Enlaces

Relacionados en esta página

phyloseq

phyloseq es el framework dominante para análisis integrado de microbioma en R: combina tabla de abundancias (OTU/ASV), taxonomía, metadatos de muestra y opcionalmente filogenia y secuencias representativas, en un objeto único (phyloseq-class). Cierra el ciclo desde la salida de DADA2 (o QIIME2) hasta los análisis estadísticos y la visualización.

Desarrollado por Paul J. McMurdie y Susan Holmes. Maduro y muy extendido en publicaciones, aunque su ritmo de desarrollo se ha enlentecido frente a alternativas modernas basadas en el stack de Bioconductor (ver más abajo).

Cuándo usarlo

  • Datos de amplicon (16S, 18S, ITS) o shotgun metagenómica representados como tabla de abundancia + taxonomía + metadatos.
  • Diversidad alfa (estimate_richness) y beta (distance + ordinate).
  • Visualización rápida con plot_bar, plot_ordination, plot_heatmap.

Cuándo NO usarlo (alternativas modernas)

  • Para un proyecto nuevo, considera mia / miaViz o MicrobiotaProcess, construidos sobre TreeSummarizedExperiment, son compatibles con el resto de Bioconductor moderno y mejor mantenidos.
  • Para downstream avanzado (modelos jerárquicos, longitudinal), microbiome, ANCOMBC y ALDEx2 ofrecen métodos más recientes que los disponibles dentro de phyloseq.

Conceptos clave

  • El objeto integra hasta cinco slots: otu_table, tax_table, sample_data, phy_tree, refseq. No todos son obligatorios.
  • taxa_are_rows(physeq) indica la orientación de la matriz. Compruébalo siempre antes de operar sobre filas/columnas, es la fuente número uno de errores opacos.
  • prune_samples() / prune_taxa() hacen subsetting preservando la coherencia entre todos los slots.
  • Rarefacción: rarefy_even_depth() está disponible, pero el propio McMurdie publicó en 2014 una crítica argumentando que es estadísticamente subóptima. Para tests formales, prefiere modelos sobre counts crudos (DESeq2, edgeR aplicados al microbioma, ANCOMBC, ALDEx2).

Patrón mínimo

library(phyloseq)

# Asumiendo otu_mat (matriz numérica), tax_mat (matriz de taxonomía con
# rangos como columnas) y sample_df (data.frame con metadatos).
OTU <- otu_table(otu_mat, taxa_are_rows = TRUE)
TAX <- tax_table(tax_mat)
SAM <- sample_data(sample_df)

ps <- phyloseq(OTU, TAX, SAM)

# Diversidad alfa
alpha <- estimate_richness(ps, measures = c("Shannon", "Simpson", "Observed"))

# Ordenación beta (Bray-Curtis + PCoA)
ord <- ordinate(ps, method = "PCoA", distance = "bray")
plot_ordination(ps, ord, color = "condition") + ggplot2::theme_bw()

Trampas habituales

  • taxa_are_rows. Confundir la orientación de la matriz lleva a errores opacos. Convención típica de DADA2: ASVs en filas. Convención típica de QIIME (BIOM): muestras en filas. Compruébalo siempre tras un import.
  • Importar desde QIIME2. qiime2R (fuera de Bioconductor) facilita el puente. El import directo del .qza es frágil y dependiente de la versión.
  • Rarefacción por ritual. Muy común en literatura pero criticable. Considera tests sobre counts crudos con normalización dentro del modelo (DESeq2, edgeR) o pruebas adaptadas a composición (ALDEx2, ANCOMBC).

Enlaces

Relacionados en esta página

  • SummarizedExperiment, base de las alternativas modernas tipo TreeSummarizedExperiment.
  • DESeq2, útil también para tests de abundancia diferencial en datos de microbioma con counts crudos.