Caso completo: análisis diferencial con airway

r
bioconductor
rnaseq
Un análisis RNA-seq de principio a fin sobre el dataset airway. Carga, filtrado, modelo con DESeq2 ajustando por línea celular, shrinkage, QC visual, volcano publicable, anotación y enriquecimiento, todo en el orden que aparece en un análisis real.

El dataset

Cerramos la ruta con el dataset airway, un estudio clásico de Himes et al. (2014) sobre el efecto de dexametasona en células del músculo liso de las vías aéreas humanas. 8 muestras: 4 líneas celulares × 2 tratamientos (untrt vs trt).

Es el dataset canónico de Bioconductor para aprender DESeq2, y como caso real, ilustra dos retos típicos: N pequeño (8 muestras) y batch técnico evidente (las líneas celulares).

1. Carga e inspección inicial

library(DESeq2)
library(airway)

data(airway)
se <- airway

se
# class: RangedSummarizedExperiment
# dim: 63677 8
# assays(1): counts
# colData names(9): SampleName cell ... Sample BioSample

Inspección del colData:

colData(se)[, c("cell", "dex")]
# DataFrame with 8 rows
#                    cell    dex
#            <character> <factor>
# SRR1039508      N61311   untrt
# SRR1039509      N61311     trt
# SRR1039512     N052611   untrt
# SRR1039513     N052611     trt
# SRR1039516     N080611   untrt
# SRR1039517     N080611     trt
# SRR1039520     N061011   untrt
# SRR1039521     N061011     trt

4 líneas celulares, cada una con tratado y no tratado. Diseño pareado, el más limpio posible para 8 muestras.

Comprobación de tamaño de librería:

colSums(assay(se)) / 1e6
# SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516 SRR1039517 SRR1039520 SRR1039521
#      20.6       18.8       25.3       15.2       24.4       30.8       19.1       21.2

Entre 15M y 31M reads. Variación habitual, nada preocupante.

2. Construcción del DESeqDataSet

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

El diseño ~ cell + dex significa:

  • Modela el efecto medio de cada línea celular como covariable.
  • Sobre lo que queda, mide el efecto de dex (tratamiento).
  • Por convención, la variable de interés (dex) va al final.

Niveles del factor dex:

levels(dds$dex)
# [1] "trt"   "untrt"

trt viene antes alfabéticamente, pero conceptualmente queremos untrt como referencia:

dds$dex <- relevel(dds$dex, ref = "untrt")
levels(dds$dex)
# [1] "untrt" "trt"

Ahora un LFC positivo significará “sube con dexametasona”.

3. Filtrado mínimo

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

dim(dds)
# [1] 18032     8

De 63.677 genes pasamos a ~18.000 con expresión razonable. Criterio >= 10 en >= 4 muestras (la mitad de las 8), más estricto que el default por el N pequeño.

4. Ejecución de DESeq2

dds <- DESeq(dds)

Tarda menos de 10 segundos en este dataset. Tres pasos internos: size factors, dispersión, Wald.

Diagnóstico de dispersión:

plotDispEsts(dds)

La curva debe bajar monotónicamente. Si tiene “joroba” o saltos, problema. Aquí se ve limpia.

5. Extracción de resultados con shrinkage

res <- lfcShrink(dds, coef = "dex_trt_vs_untrt", type = "apeglm")

summary(res)
# adjusted p-value < 0.1
# LFC > 0 (up)      : 1545, 8.6%
# LFC < 0 (down)    : 1379, 7.6%

~3000 genes con padj < 0.1. Magnitud razonable para un efecto biológico fuerte como dexametasona.

Versión más estricta:

res05 <- results(dds, alpha = 0.05, name = "dex_trt_vs_untrt")
summary(res05)

6. QC visual

PCA

vsd <- vst(dds, blind = FALSE)

library(ggplot2)
pca_data <- plotPCA(vsd, intgroup = c("dex", "cell"), returnData = TRUE)
pv <- round(100 * attr(pca_data, "percentVar"))

ggplot(pca_data, aes(PC1, PC2, color = dex, shape = cell)) +
    geom_point(size = 4, alpha = 0.85) +
    xlab(paste0("PC1: ", pv[1], "% varianza")) +
    ylab(paste0("PC2: ", pv[2], "% varianza")) +
    scale_color_manual(values = c(untrt = "#5F8575", trt = "#A24B4B")) +
    coord_fixed() +
    theme_minimal(base_size = 12)

Lo que se ve en este dataset: PC1 separa por línea celular (~50 % varianza), PC2 separa por tratamiento (~25 %). Es decir, la variación biológica entre individuos domina sobre la del tratamiento. Esto no es problema, el modelo ~ cell + dex ajusta exactamente eso. Sin ese ajuste, el efecto del tratamiento se “diluiría” en la variación entre líneas.

Heatmap de muestras

library(pheatmap)
library(RColorBrewer)

sample_dist <- dist(t(assay(vsd)))
mat <- as.matrix(sample_dist)
rownames(mat) <- paste(vsd$cell, vsd$dex, sep = " · ")
colnames(mat) <- NULL

pheatmap(mat,
         clustering_distance_rows = sample_dist,
         clustering_distance_cols = sample_dist,
         col = colorRampPalette(rev(brewer.pal(9, "Greens")))(255),
         border_color = NA)

Bloques 2×2 de cada línea celular, las dos muestras de N052611, las dos de N61311, etc. Confirma lo del PCA.

7. Anotación

library(AnnotationDbi)
library(org.Hs.eg.db)
library(dplyr)
library(tibble)

res_df <- as.data.frame(res) |>
    rownames_to_column("ensembl_id") |>
    mutate(
        simbolo = mapIds(org.Hs.eg.db, keys = ensembl_id,
                         column = "SYMBOL", keytype = "ENSEMBL",
                         multiVals = "first"),
        descripcion = mapIds(org.Hs.eg.db, keys = ensembl_id,
                             column = "GENENAME", keytype = "ENSEMBL",
                             multiVals = "first"),
        entrez_id = mapIds(org.Hs.eg.db, keys = ensembl_id,
                           column = "ENTREZID", keytype = "ENSEMBL",
                           multiVals = "first")
    ) |>
    arrange(padj)

Inspección rápida de los top 10:

res_df |>
    filter(!is.na(padj)) |>
    head(10) |>
    select(simbolo, log2FoldChange, padj, descripcion)

Los top genes incluyen reguladores conocidos de respuesta a glucocorticoides: FKBP5, KLF15, PER1. Validación biológica directa, sale lo que la literatura espera.

8. Volcano publicable

library(ggrepel)

df_plot <- res_df |>
    filter(!is.na(padj)) |>
    mutate(estado = case_when(
        padj < 0.05 & log2FoldChange >  1 ~ "Up",
        padj < 0.05 & log2FoldChange < -1 ~ "Down",
        TRUE                              ~ "NS"
    ))

etiquetas <- df_plot |>
    filter(estado != "NS") |>
    slice_min(padj, n = 12)

ggplot(df_plot, aes(log2FoldChange, -log10(padj), color = estado)) +
    geom_point(size = 1, alpha = 0.5) +
    geom_text_repel(data = etiquetas,
                    aes(label = simbolo),
                    size = 3.2, box.padding = 0.5,
                    max.overlaps = 20, seed = 42) +
    scale_color_manual(values = c(Up = "#A24B4B", Down = "#3D6FA8", NS = "grey70")) +
    geom_vline(xintercept = c(-1, 1), linetype = "dashed", color = "grey50") +
    geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "grey50") +
    labs(x = "log2 fold-change (trt vs untrt)",
         y = "-log10 p-valor ajustado",
         color = NULL) +
    theme_minimal(base_size = 12) +
    theme(panel.grid.minor = element_blank(),
          legend.position = "top")

Volcano simétrico con tendencia bilateral, efecto del tratamiento amplio, no sesgado. Las etiquetas destacan los genes conocidos de respuesta a glucocorticoides.

9. Enriquecimiento funcional

library(clusterProfiler)

genes_up <- res_df |>
    filter(padj < 0.05, log2FoldChange > 1) |>
    pull(entrez_id) |>
    na.omit() |>
    unique()

universo <- res_df |>
    pull(entrez_id) |>
    na.omit() |>
    unique()

go_bp <- enrichGO(gene          = genes_up,
                  universe      = universo,
                  OrgDb         = org.Hs.eg.db,
                  ont           = "BP",
                  pAdjustMethod = "BH",
                  pvalueCutoff  = 0.05)

dotplot(go_bp, showCategory = 12)

Términos esperables: “response to glucocorticoid”, “response to corticosteroid”, “regulation of inflammatory response”. Coherente con la biología del tratamiento.

10. Export del resultado

res_df |>
    select(ensembl_id, simbolo, log2FoldChange, lfcSE, padj, baseMean, descripcion) |>
    write.csv("airway_DE_results.csv", row.names = FALSE)

Una tabla limpia: ID, símbolo, LFC encogido, error, p-valor ajustado, expresión media, descripción. Es lo que acompaña a la figura en un informe.

11. Reporte de versiones

Esencial para reproducibilidad:

sessionInfo()

Captura las versiones de R, DESeq2, organismos, dependencias. Siempre al final del documento Quarto del análisis.

El flujo completo en una vista

Repasando, los pasos del análisis:

  1. Carga + inspección (data(airway), colData, colSums).
  2. Diseño con relevel (~ cell + dex, relevel(ref = "untrt")).
  3. Filtrado mínimo (rowSums(counts >= 10) >= 4).
  4. DESeq(): un solo comando, tres pasos internos.
  5. lfcShrink con apeglm: para visualización y ranking.
  6. PCA + heatmap: diagnóstico antes de interpretar.
  7. Anotación: ENSEMBL → símbolo + descripción.
  8. Volcano publicable: con etiquetas y umbrales claros.
  9. Enriquecimiento GO: sanity check biológico.
  10. Export + sessionInfo(): reproducibilidad.

Diez pasos. Esto es un análisis estándar moderno, lo que llena 60-100 líneas en un script o un Quarto.

Has terminado la ruta

Con esto cierras RNA-seq con Bioconductor. Has visto:

  • El ecosistema Bioconductor y su filosofía.
  • El camino FASTQ → counts (conceptual) y los contenedores SummarizedExperiment / DESeqDataSet.
  • El motor de DESeq2: size factors, dispersión, Wald, contrastes.
  • Transformaciones (vst, rlog) y por qué importan.
  • Diseño experimental: covariables, interacciones, contrastes explícitos.
  • Shrinkage de LFC: apeglm como estándar moderno.
  • Visualización: PCA, heatmap, volcano, MA.
  • Anotación y enriquecimiento funcional.
  • Un caso real completo.

A partir de aquí, las direcciones naturales: análisis multi-factor, batch effects con sva/RUV, single-cell (Seurat, Bioconductor scRNA-seq), tiempo (maSigPro), variantes (VariantAnnotation).

Cuando publiquemos el libro RNA-seq con DESeq2: del diseño al informe, cubriremos diseños complejos, gestión de batch effect, pipelines con targets y renv, e integración con Quarto para reports profesionales.