Bioconductor: el ecosistema de R para bioinformática y análisis genómico

Bioconductor es el repositorio de paquetes R especializado en bioinformática y biología computacional. Surgió en 2001 como respuesta a la necesidad de herramientas estadísticas robustas para analizar los datos del Proyecto Genoma Humano, y desde entonces ha crecido hasta superar los 2.200 paquetes en 2026. Su estructura es diferente a CRAN: los paquetes tienen ciclos de release semestrales sincronizados con versiones de R, lo que garantiza compatibilidad entre ellos.

Bioconductor 3.19, publicado en 2024, es la versión más reciente con soporte activo. La instalación de paquetes no se hace con install.packages() sino con BiocManager:

# Instalar el gestor de Bioconductor
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

# Instalar paquetes de Bioconductor
BiocManager::install("DESeq2")
BiocManager::install("limma")
BiocManager::install("GenomicRanges")
BiocManager::install("Biostrings")
BiocManager::install("edgeR")

# Ver versión de Bioconductor instalada
BiocManager::version()

# Comprobar paquetes desactualizados
BiocManager::valid()

Las estructuras de datos de Bioconductor

Bioconductor define sus propias clases S4 para representar datos genómicos. La más usada es SummarizedExperiment, que combina la matriz de conteos, los metadatos de las muestras y la anotación de los genes en un solo objeto:

library(SummarizedExperiment)

# Crear un SummarizedExperiment de ejemplo
conteos <- matrix(
  sample(1:1000, 500, replace = TRUE),
  nrow  = 100,
  ncol  = 5,
  dimnames = list(
    paste0("gen_", 1:100),
    paste0("muestra_", 1:5)
  )
)

metadatos_muestras <- data.frame(
  muestra  = paste0("muestra_", 1:5),
  condicion = c("control", "control", "tratado", "tratado", "tratado"),
  batch     = c(1, 2, 1, 2, 1),
  row.names = paste0("muestra_", 1:5)
)

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

# Acceso a los componentes
assay(se, "counts")[1:5, ]    # matriz de conteos
colData(se)                    # metadatos de muestras
rownames(se)[1:10]             # nombres de genes
dim(se)                        # dimensiones: 100 genes x 5 muestras

DESeq2: expresión diferencial en RNA-seq

DESeq2 es el paquete estándar para detectar genes con expresión diferencial entre condiciones a partir de datos de RNA-seq. Su modelo estadístico usa distribución binomial negativa y normaliza automáticamente por profundidad de secuenciación:

library(DESeq2)

# Crear objeto DESeqDataSet
dds <- DESeqDataSet(se, design = ~ condicion)

# Filtrar genes con pocos conteos (opcional pero recomendado)
keep <- rowSums(counts(dds) >= 10) >= 3
dds  <- dds[keep, ]

# Análisis completo: normalización + estimación de dispersión + test
dds <- DESeq(dds)

# Resultados de expresión diferencial
res <- results(dds, contrast = c("condicion", "tratado", "control"))
summary(res)

# Genes significativos
genes_sig <- subset(res, padj < 0.05 & abs(log2FoldChange) > 1)
nrow(genes_sig)

# MA plot (fold change vs expresión media)
plotMA(res, ylim = c(-5, 5))

# Volcano plot con ggplot2
library(ggplot2)
as.data.frame(res) |>
  dplyr::mutate(
    significativo = !is.na(padj) & padj < 0.05,
    regulacion    = dplyr::case_when(
      log2FoldChange >  1 & padj < 0.05 ~ "subido",
      log2FoldChange < -1 & padj < 0.05 ~ "bajado",
      TRUE ~ "sin cambio"
    )
  ) |>
  ggplot(aes(x = log2FoldChange, y = -log10(pvalue), color = regulacion)) +
    geom_point(alpha = 0.6, size = 1) +
    scale_color_manual(values = c("subido" = "red", "bajado" = "blue",
                                   "sin cambio" = "grey50")) +
    theme_minimal() +
    labs(title = "Volcano plot RNA-seq", x = "Log2 Fold Change")

GenomicRanges: intervalos genómicos

GenomicRanges es la clase fundamental para representar regiones del genoma: genes, exones, picos de ChIP-seq, variantes. Opera de forma vectorizada sobre millones de intervalos genómicos:

library(GenomicRanges)

# Crear GRanges (regiones genómicas)
genes <- GRanges(
  seqnames = c("chr1", "chr1", "chr2", "chr3"),
  ranges   = IRanges(
    start = c(1000, 5000, 2000, 8000),
    end   = c(2500, 7000, 4000, 9500)
  ),
  strand   = c("+", "-", "+", "+"),
  nombre   = c("BRCA1", "TP53", "EGFR", "MYC")
)

# Operaciones sobre rangos
width(genes)                   # longitud de cada región
genes[seqnames(genes) == "chr1"]  # filtrar por cromosoma

# Encontrar solapamientos entre dos conjuntos de regiones
picos_chip <- GRanges(
  seqnames = c("chr1", "chr2"),
  ranges   = IRanges(start = c(1200, 3000), end = c(1800, 3500))
)
findOverlaps(genes, picos_chip)
subsetByOverlaps(genes, picos_chip)

limma: microarrays y RNA-seq

limma (Linear Models for Microarray data) es anterior a DESeq2 y sigue siendo muy usado, especialmente para datos de microarrays y para RNA-seq cuando el tamaño de muestra es pequeño. Su modelo lineal es más flexible que DESeq2 para diseños experimentales complejos:

library(limma)

# Matriz de expresión (genes x muestras), ya normalizada y en log2
expr <- matrix(rnorm(500, mean = 6), nrow = 100, ncol = 5)
rownames(expr) <- paste0("gen_", 1:100)

# Matriz de diseño
grupo    <- factor(c("control", "control", "tratado", "tratado", "tratado"))
design   <- model.matrix(~ grupo)

# Ajuste del modelo lineal
fit   <- lmFit(expr, design)
fit2  <- eBayes(fit)

# Tabla de resultados con p-valores ajustados
topTable(fit2, coef = 2, n = 20, adjust.method = "BH")

Bioconductor es el ecosistema que hace que R sea insustituible en ciencias de la vida. No hay equivalente en Python con la misma integración, calidad estadística y trayectoria académica. Si te interesa cómo R interacciona con Python para combinar ambos ecosistemas, el artículo sobre reticulate cubre esa integración. Y para entender la base estadística que subyace a DESeq2 y limma, el artículo sobre modelos lineales en R es el punto de partida.

Imagen: Pexels / Google DeepMind