Cargar y filtrar counts

r
bioconductor
rnaseq
Tres formas de cargar datos RNA-seq en R: matriz cruda desde CSV, importación moderna con tximeta desde Salmon y descarga de datasets públicos. Cómo construir el colData fiable y el filtrado de genes de baja expresión que de verdad importa.

Tres orígenes habituales

Tres situaciones de partida que cubren el 95 % de los casos:

  1. Matriz de counts en CSV/TSV: te dieron una tabla ya cuantificada.
  2. Output de Salmon (un directorio por muestra), el pipeline moderno estándar.
  3. Dataset público: airway, pasilla, recount3, o uno de GEO/SRA.

Cada uno tiene su patrón. Conociéndolos, casi cualquier dataset se carga en 10 líneas.

Opción 1: matriz CSV + metadata CSV

El caso “te paso esto por email”:

library(DESeq2)
library(readr)

counts <- read.csv("counts.csv", row.names = 1, check.names = FALSE) |>
    as.matrix()

meta <- read.csv("metadata.csv", row.names = 1)

Detalles importantes:

  • row.names = 1: la primera columna del CSV son los IDs de gen (filas) o de muestra. Sin esto, queda como columna normal.
  • check.names = FALSE: por defecto R sustituye caracteres “raros” en los nombres (espacios, guiones). En IDs de muestra puede romper la correspondencia con el metadata.
  • as.matrix(): DESeq2 espera matriz, no data.frame.

Verificar alineación, siempre:

stopifnot(all(colnames(counts) == rownames(meta)))

Si falla, reordena el meta:

meta <- meta[colnames(counts), ]      # mismo orden que las columnas del counts

Construir el DESeqDataSet:

dds <- DESeqDataSetFromMatrix(
    countData = counts,
    colData   = meta,
    design    = ~ condicion
)

Opción 2: Salmon con tximeta

Es el patrón moderno. Cada muestra tiene su carpeta de Salmon (con quant.sf dentro). Necesitas un coldata que apunte a esos archivos:

library(tximeta)

coldata <- data.frame(
    names      = c("m1", "m2", "m3", "m4", "m5", "m6"),
    files      = file.path("salmon", c("m1", "m2", "m3", "m4", "m5", "m6"), "quant.sf"),
    condicion  = c("ctrl", "ctrl", "ctrl", "tto", "tto", "tto")
)

stopifnot(all(file.exists(coldata$files)))

se  <- tximeta(coldata)              # carga a nivel de transcrito
gse <- summarizeToGene(se)           # agrega a nivel de gen
dds <- DESeqDataSet(gse, design = ~ condicion)

Qué hace tximeta:

  • Detecta el transcriptoma que usó Salmon (mediante un hash en los outputs).
  • Descarga la anotación correspondiente (gene IDs, símbolos) la primera vez y la cachea.
  • Devuelve un SummarizedExperiment con rowRanges ya anotados.

La diferencia frente a tximport (la versión anterior): tximport requiere que tú aportes el mapeo transcrito → gen. tximeta lo hace solo. Si tu Salmon usó un transcriptoma estándar (Ensembl, GENCODE, RefSeq), tximeta es la opción.

Opción 3: dataset de paquete

Para aprender, probar pipelines, o tener un caso de referencia:

library(airway)
data(airway)
se <- airway

se
# class: RangedSummarizedExperiment
# dim: 63677 8
# assays(1): counts
# colData names: SampleName cell dex albut Run avgLength Experiment Sample BioSample

8 muestras de 4 líneas celulares pulmonares (cell) tratadas con dexametasona (dex = "trt") o no (dex = "untrt"). Es el dataset canónico para aprender DESeq2, lo usaremos como ejemplo en varios tutoriales.

Conversión a DESeqDataSet:

library(DESeq2)
dds <- DESeqDataSet(se, design = ~ cell + dex)

Datasets parecidos en otros paquetes: pasilla (Drosophila), parathyroidSE, geneLevelData de recount3. Para datos públicos serios, recount3 descarga datasets uniformemente procesados de GEO/SRA.

El colData: cómo no equivocarse

El metadata es lo que más bugs introduce. Tres reglas:

1. Tipos de variables explícitos

meta$condicion <- factor(meta$condicion, levels = c("ctrl", "tto"))
meta$batch     <- factor(meta$batch)

DESeq2 usa el primer nivel del factor como referencia. Si dejas que R lo decida alfabéticamente, “ctrl” sale antes que “tto” (alfabético) y funciona. Pero "WT" vs "KO", alfabético da KO como referencia, que no es lo que quieres. Define los niveles explícitos.

2. Renombrar columnas problemáticas

Si una columna se llama Treatment Group con espacio, las fórmulas (~ Treatment Group) fallan. Renombra a treatment_group o treatment antes de empezar.

3. Documenta unidades y codificación

En el colData mismo, o en un documento aparte. ¿dex = "trt" es el grupo tratado o el control? Resulta obvio ahora, no lo será dentro de 6 meses.

Filtrar genes de baja expresión

La matriz de counts típica tiene 30.000-60.000 genes. Muchos están casi sin lecturas en ninguna muestra. Esos genes:

  • No aportan información (no se puede testear nada sin counts).
  • Inflan el coste de los ajustes de p-valor (multiple testing).
  • Hacen el análisis innecesariamente lento.

El filtrado idiomático:

keep <- rowSums(counts(dds) >= 10) >= 3
dds <- dds[keep, ]

Lectura: “al menos 3 muestras tienen ≥ 10 counts”. El criterio >= 10 y >= 3 (tamaño del grupo experimental más pequeño) es lo que recomienda la vignette de DESeq2, no algo dogmático, pero buen default.

Filtros menos cuidadosos que se ven a veces:

keep <- rowSums(counts(dds)) > 0          # ❌ deja basura: 1 count en una muestra
keep <- rowMeans(counts(dds)) >= 10       # ⚠️ ignora distribución por grupo

El primero es muy laxo. El segundo deja pasar genes que son altos en 1-2 muestras outlier.

No hace falta filtrar mucho: DESeq2 hace su propio “independent filtering” en results(). El pre-filtro de arriba es más por velocidad que por estadística. No te obsesiones con esto.

¿Qué NO filtrar antes?

  • No quitar muestras sin justificación clara (outlier real, falló QC). Cada muestra menos = menos potencia.
  • No quitar genes “que no esperas que sean DE”. La estadística necesita el conjunto completo para estimar la dispersión.
  • No transformar counts antes de DESeq2. Los modelos de DESeq2 esperan enteros crudos. La normalización la hace el paquete por dentro.

Inspección rápida

Después de cargar, dos comprobaciones de salud:

dim(dds)              # genes × muestras
colSums(counts(dds))  # library size por muestra
table(dds$condicion)  # tamaño de cada grupo

colSums(counts) debe estar en órdenes parecidos entre muestras, habitualmente de 5 a 50 millones. Si una muestra tiene 100k y el resto 30M, algo va mal (problema en la secuenciación, archivo corrupto, etc.). Investigar antes de seguir.

Trampas habituales

  • Counts y metadata en distinto orden sin verificar. Sintetiza todo bien hasta que los resultados son raros. La línea stopifnot(all(colnames(counts) == rownames(meta))) cuesta un segundo y salva análisis.
  • Decimales de Salmon sin tximport/tximeta. DESeqDataSetFromMatrix se queja. Convierte con round() o, mejor, usa el flujo tximeta.
  • Filtros agresivos antes de DESeq2. Hay tutoriales por ahí que sugieren rowSums > 100 u otros. No hace daño si los genes filtrados realmente no tienen señal, pero pierdes información de la distribución.
  • No definir niveles de factor. DESeq2 toma el primero alfabético como referencia, y a veces eso es contraintuitivo. Define con levels = explícito.
  • check.names = TRUE (default) que renombra columnas con guiones, espacios, etc. y rompe la correspondencia con metadata. Usa check.names = FALSE al leer counts.

En la siguiente entrega

Ya tienes un DESeqDataSet cargado y filtrado. La siguiente entrega es DESeq2 funcionando: la función DESeq(), qué hace por dentro (size factors, dispersiones, fit), cómo extraer resultados con results() y la interpretación básica. El núcleo del análisis. Lo siguiente.