SummarizedExperiment y DESeqDataSet

r
bioconductor
rnaseq
El contenedor estándar de Bioconductor para datos ómicos. Cómo se estructura, por qué counts + colData + rowData en un solo objeto evita errores, y cómo DESeqDataSet extiende SummarizedExperiment con el diseño experimental.

La idea

En análisis ómico tienes tres piezas que tienen que ir juntas:

  • Una matriz de medidas: filas = genes, columnas = muestras (counts, intensidades, lo que sea).
  • Metadata de columnas (colData): información de cada muestra, condición, batch, sexo del donante.
  • Metadata de filas (rowData): información de cada gen, símbolo, cromosoma, longitud.

Si tienes las tres en data.frames sueltos, te equivocarás antes o después: filtrar la matriz pero olvidar filtrar el colData, reordenar muestras y desalinear, perder un gen al pasar de una versión de anotación a otra…

SummarizedExperiment resuelve esto: un único objeto que mantiene las tres sincronizadas. Filtras o reordenas, y todo se mueve a la vez. Es lo que comparten DESeq2, edgeR, limma, single-cell, ChIP-seq y casi todo Bioconductor moderno.

Estructura interna

SummarizedExperiment
├── assays         (lista de matrices: "counts", "tpm", "logcounts"...)
│     filas = genes
│     columnas = muestras
├── colData         (DataFrame: una fila por muestra, columnas = covariables)
├── rowData         (DataFrame: una fila por gen, columnas = anotaciones)
└── metadata        (lista libre: cualquier cosa adicional)

rowRanges es una variante de rowData para cuando las filas son rangos genómicos (no solo metadata): genes con coordenadas, picos de ChIP-seq, etc. En RNA-seq estándar lo verás cuando uses tximeta.

Construir uno desde cero

library(SummarizedExperiment)

counts <- matrix(rpois(60, 50), nrow = 10, ncol = 6)
rownames(counts) <- paste0("gen_", 1:10)
colnames(counts) <- paste0("m", 1:6)

coldata <- DataFrame(
    condicion = c("ctrl", "ctrl", "ctrl", "tto", "tto", "tto"),
    batch     = c("A", "A", "B", "A", "B", "B"),
    row.names = colnames(counts)
)

rowdata <- DataFrame(
    simbolo = paste0("GEN_", 1:10),
    cromosoma = sample(1:22, 10, replace = TRUE),
    row.names = rownames(counts)
)

se <- SummarizedExperiment(
    assays   = list(counts = counts),
    colData  = coldata,
    rowData  = rowdata
)

se
# class: SummarizedExperiment
# dim: 10 6
# assays(1): counts
# rownames(10): gen_1 gen_2 ... gen_9 gen_10
# colnames(6): m1 m2 ... m5 m6

Tres bloques (matriz, colData, rowData) en un solo objeto.

Acceso al contenido

assay(se)                  # primera matriz (counts)
assay(se, "counts")        # por nombre
assays(se)                 # lista con todas las matrices
colData(se)                # metadata de muestras (DataFrame)
rowData(se)                # metadata de genes (DataFrame)

dim(se)                    # 10 6
rownames(se)               # ids de genes
colnames(se)               # ids de muestras

Importante: colData(se) devuelve un DataFrame (con D mayúscula), clase de Bioconductor (S4Vectors::DataFrame), no el data.frame clásico. Comportamiento muy parecido, pero soporta columnas de objetos complejos (lists, Rle…) que un data.frame normal no.

Para convertirlo a data.frame clásico cuando sea útil:

as.data.frame(colData(se))

Subsetting: filas y columnas

La sintaxis es como un data.frame con dos índices: [filas, columnas]. Pero filas = genes, columnas = muestras, no es la convención de tidyverse (filas = observaciones).

se[1:3, ]                                       # primeros 3 genes, todas las muestras
se[, se$condicion == "tto"]                     # todas las filas, solo muestras tto
se[, se$batch == "A"]                           # solo batch A
se[rowSums(assay(se)) > 100, ]                  # genes con suma > 100 counts
se[rowData(se)$cromosoma == 1, ]                # genes del cromosoma 1

El acceso se$columna funciona como atajo para acceder a una columna del colData. Es lo que permite escribir filtros legibles.

Lo crucial: al hacer subset, se mueve todo a la vez. Si filtras a 3 muestras, el assay se queda con 3 columnas Y el colData con 3 filas Y el rowData no se toca. Si filtras a 100 genes, los assays se reducen a 100 filas Y el rowData también. Sin sincronizar a mano.

Múltiples assays

Una de las ventajas del contenedor: puedes tener varias matrices sobre el mismo objeto:

assays(se)$logcounts <- log2(assay(se) + 1)

assayNames(se)
# [1] "counts" "logcounts"

assay(se, "logcounts")[1:3, 1:3]

Útil para guardar “counts crudos” + “counts normalizados” + “vst” en un solo objeto. Lo verás al normalizar con DESeq2.

DESeqDataSet: la extensión para DESeq2

DESeq2 no usa SummarizedExperiment tal cual, usa DESeqDataSet, una subclase que añade:

  • Un diseño (design): la fórmula que define qué quieres modelar.
  • Validación: los counts deben ser enteros no negativos.
  • Slots adicionales para size factors, dispersions, etc.

Construcción desde matriz cruda:

library(DESeq2)

dds <- DESeqDataSetFromMatrix(
    countData = counts,
    colData   = coldata,
    design    = ~ batch + condicion       # modelar batch como covariable
)

O desde un SummarizedExperiment ya construido:

dds <- DESeqDataSet(se, design = ~ batch + condicion)

Y desde Salmon/tximeta (la forma idiomática moderna):

dds <- DESeqDataSet(gse, design = ~ condicion)

Hereda todo de SummarizedExperiment: assay(), colData(), rowData(), subsetting con [ , ]. Más los métodos propios de DESeq2: sizeFactors(), dispersions(), results(), etc.

La fórmula design

La fórmula define qué se modela. Sintaxis estándar de R:

~ condicion                  # solo condición
~ batch + condicion          # condición ajustando por batch
~ sex + age + condicion      # tres covariables, condición la última
~ condicion * tiempo         # interacción

Convención DESeq2: la variable de interés va al final de la fórmula. Las anteriores se tratan como covariables a ajustar. Lo veremos a fondo en el tutorial 7 (diseño y contrastes).

El patrón idiomático

El flujo típico al recibir datos:

# 1. Cargar matriz + metadata
counts <- read.csv("counts.csv", row.names = 1) |> as.matrix()
meta   <- read.csv("meta.csv", row.names = 1)

# 2. Verificar alineación
stopifnot(all(colnames(counts) == rownames(meta)))

# 3. Construir DESeqDataSet directamente
dds <- DESeqDataSetFromMatrix(
    countData = counts,
    colData   = meta,
    design    = ~ condicion
)

stopifnot(all(colnames(counts) == rownames(meta))) es la línea más importante de cualquier análisis RNA-seq. Si las columnas del counts no están en el mismo orden que las filas del metadata, todo el análisis está mal. Pero no falla, produce resultados con sentido aparente pero totalmente equivocados. Verifica siempre.

Trampas habituales

  • Counts decimales (de Salmon sin tximport). DESeq2 exige enteros, si te quejas pasa lo que pasa. round() en su sitio o, mejor, usa tximeta que ya lo hace.
  • colData con el mismo orden de filas que columnas del counts pero no los mismos nombres. Pasa cuando se procesan en orden distinto. Siempre verifica con all(colnames(counts) == rownames(coldata)).
  • Confundir rowData con rowRanges. Si tu objeto viene de tximeta o summarizeOverlaps, suele venir con rowRanges (con coordenadas genómicas) además de rowData. Para análisis estándar es transparente, pero al exportar y reimportar hay que conservar.
  • Olvidar el design al construir un DESeqDataSet. Lo permite (design = ~ 1) pero luego no podrás llamar a DESeq(). Mejor definirlo desde el principio aunque después lo cambies.
  • Asumir que filas = observaciones. En tidyverse sí. En SummarizedExperiment filas = genes, columnas = muestras. Cuando llevas el objeto a tidyverse (con broom::tidy() o similar), recuerda transponer mentalmente.

En la siguiente entrega

Ya entiendes el contenedor. Ahora toca cargar datos reales: importar de Salmon con tximeta, leer counts de un CSV, ensamblar el colData y filtrar genes poco expresados antes del análisis. Lo siguiente.