Resumen de la serie "De FASTQ a artículo" — ¿En qué punto nos encontramos?

• En S17, realizamos un alineamiento aproximado con Salmon sobre los archivos FASTQ sin procesar para obtener recuentos por gen.
• En S18, comparamos los métodos de cuantificación featureCounts y RSEM.
• En este episodio (S19), transformaremos esos recuentos en una figura para un artículo: un diagrama de dispersión (volcano plot).
• En el próximo S20, aprenderemos las alternativas edgeR/limma-voom y, en S21, realizaremos el análisis de vías GSEA.

En otras palabras, al completar el ejercicio de este episodio, dominarás la parte más importante del proceso estadístico que transforma los archivos FASTQ sin procesar en figuras listas para su publicación. Esta es una dificultad inevitable que todo investigador de laboratorio húmedo (wet-lab) encuentra al escribir un artículo.

¿Por qué una distribución binomial negativa y no una normal?

Primero, entendamos los datos. Una matriz de recuentos de RNA-seq tiene el siguiente aspecto:

Lo primero que notarás: son recuentos. Enteros. Sin valores negativos. Lo segundo: la varianza es proporcional a la media. Si ACTB tiene una media de 8000, su varianza también es grande; si TP53 tiene una media de 90, su varianza es pequeña. Lo tercero: la relación entre la media y la varianza es mucho más amplia de lo que predice la distribución normal (sobredispersión).

Enfoque ingenuo: Poisson

Al ser datos de recuentos, naturalmente podrías pensar en la distribución de Poisson. Poisson tiene Var(X) = E(X) = μ. Sin embargo, en los datos reales de RNA-seq, la varianza entre réplicas es abrumadoramente mayor que la media. Esto es especialmente cierto para los genes de baja expresión. Si usas Poisson directamente, las estimaciones estadísticas colapsarán.

Pasando a la distribución binomial negativa

La distribución binomial negativa (Negative Binomial, NB) es Poisson a la que se le añade un parámetro extra de dispersión α (dispersion).

Var(X) = μ + α · μ²

Si α es 0, converge a Poisson; si α es grande, la dispersión es mucho mayor. Al graficar un diagrama de dispersión (μ, Var) con datos reales de RNA-seq, esta relación encaja sorprendentemente bien. DESeq2 comienza desde aquí.

Reducción de la dispersión — La magia de DESeq2

Aquí es donde nos encontramos con las limitaciones prácticas. Los experimentos en los artículos suelen tener solo de 3 a 6 repeticiones. Al intentar estimar α para cada uno de los 20.000 genes con un número reducido de muestras, las estimaciones fluctúan drásticamente. Para algunos genes, α = 0,01; para otros, α = 5,0; en la mayoría de los casos, se trata simplemente de un error de muestreo.

Mike Love, el autor de DESeq2, recurre aquí al ajuste Bayesiano empírico. La idea es la siguiente:

1. Estimar el α_g sin ajustar para cada gen (presenta mucha fluctuación).
2. Representar el α de todos los genes en un gráfico de media-varianza y trazar una línea de ajuste suave (esta es la "tendencia").
3. El α_g* final para cada gen se determina como la media ponderada entre la estimación sin ajustar y el valor de la tendencia. Cuanto menor sea el número de muestras, más fuerte será la atracción hacia la tendencia.

¿Por qué es correcto? Porque introduce una probabilidad a priori biológica: "los genes con niveles de expresión similares tienden a tener una dispersión similar". Con este ajuste, los falsos positivos se reducen drásticamente y la reproducibilidad de los experimentos con pocas repeticiones aumenta notablemente.

Aclarar conceptos con StatQuest

Esta sección presenta conceptos complejos, por lo que es esencial contar con ayuda visual. Se recomienda familiarizarse con el tema a través de las conferencias de StatQuest con subtítulos en coreano.

• StatQuest: una introducción sencilla a RNA-seq y DESeq2 (Joshua Starmer, 26 min, subtítulos en coreano totalmente compatibles). Al activar los subtítulos en coreano en la configuración de YouTube, se obtiene una traducción precisa. La analogía del juego de emparejar tarjetas es muy memorable.

Prueba de Wald: ¿cómo determinar la significancia?

Para cada gen, se estima el cambio logarítmico en la expresión β entre dos condiciones (por ejemplo, tejido de cáncer de mama frente a tejido normal) y se realiza una prueba de Wald para comprobar si β ≠ 0.

z = β̂ / SE(β̂)

Asumiendo que z sigue una distribución normal estándar, se obtiene el valor p bilateral. DESeq2 también aplica un ajuste al propio cambio logarítmico para estabilizar las estimaciones extremas de β en los genes con baja expresión. La función lfcShrink() cumple esta función.

Benjamini-Hochberg: el problema de realizar 20.000 pruebas

Aquí nos enfrentamos directamente al problema de las pruebas múltiples en el aprendizaje automático. Si se realizan pruebas simultáneas para 20.000 genes, incluso con un umbral p de 0,05, por pura casualidad aparecerán 1.000 genes como "significativos". Si se representan estos resultados directamente en un diagrama de volcanes, el artículo sería rechazado.

El ajuste de Benjamini-Hochberg (BH) funciona de la siguiente manera:

1. Ordenar los valores p de todos los genes de menor a mayor: p₍₁₎ ≤ p₍₂₎ ≤ ... ≤ p₍ₘ₎.
2. Definir el objetivo de la tasa de falsos descubrimientos (FDR) (por ejemplo, 0,05).
3. Encontrar el valor máximo de k que satisfaga la siguiente condición.

p₍ₖ₎ ≤ (k / m) × FDR

4. Se consideran significativos los genes hasta el valor k.

La columna padj que devuelve DESeq2 contiene los valores p ajustados mediante el método de Benjamini-Hochberg (BH). El umbral de significancia estándar para el gráfico de dispersión (volcano plot) es padj < 0.05.

Si aún no se comprende bien la razón del uso de BH, este breve video de StatQuest también tiene subtítulos excelentes en coreano.

• StatQuest — False Discovery Rates (FDR), clearly explained (Joshua Starmer, 18 minutos, subtítulos en coreano disponibles y de alta calidad). Utiliza una analogía con las cartas para mostrar de forma intuitiva por qué los falsos positivos aumentan al realizar 20.000 pruebas.

Creación de un gráfico de dispersión (volcano plot) de cáncer de mama de TCGA con 30 líneas de código en R

Ahora, pasemos a la práctica. Iniciemos el kernel de R en Colab.

Paso 1: Importar los datos de TCGA

# Instalación de Bioconductor (5 minutos en la primera ejecución de Colab)
if (!require("BiocManager")) install.packages("BiocManager")
BiocManager::install(c("DESeq2", "TCGAbiolinks", "EnhancedVolcano"))

library(DESeq2)
library(TCGAbiolinks)

# Datos de recuento crudo TCGA-BRCA (tejido tumoral vs tejido normal)
query <- GDCquery(
  project    = "TCGA-BRCA",
  data.category = "Transcriptome Profiling",
  data.type  = "Gene Expression Quantification",
  workflow.type = "STAR - Counts",
  sample.type = c("Primary Tumor", "Solid Tissue Normal")
)
GDCdownload(query)
data <- GDCprepare(query)

Se obtienen los recuentos iniciales de pacientes con cáncer de mama y tejidos normales del TCGA. La primera descarga puede tardar, pero el uso de cBioPortal o UCSC Xena agiliza el proceso.

Paso 2: Ejecutar DESeq2

dds <- DESeqDataSet(data, design = ~ sample_type)
dds <- DESeq(dds)  # Aquí ocurren la estimación de dispersión, la prueba de Wald y el ajuste de BH

res <- results(dds, contrast = c("sample_type", "Primary Tumor", "Solid Tissue Normal"))
res <- lfcShrink(dds, coef = "sample_type_Primary Tumor_vs_Solid Tissue Normal", type = "apeglm")
summary(res)

El resumen que genera summary(res) con las etiquetas "up-regulated / down-regulated / outlier / low-count" es una frase que se puede incluir directamente en la sección de métodos de un artículo científico.

Paso 3: Gráfico de dispersión de volcanes

library(EnhancedVolcano)
EnhancedVolcano(res,
    lab      = rownames(res),
    x        = "log2FoldChange",
    y        = "padj",
    pCutoff  = 0.05,
    FCcutoff = 1.0,
    title    = "TCGA-BRCA Tumor vs Normal"
)

Con estas tres líneas se obtiene un gráfico de volcanes adecuado para incluir en un artículo científico. Eje X = log2 fold change (los valores negativos indican mayor expresión en el tejido normal, los positivos en el tumor), Eje Y = −log10(padj). Los puntos que aparecen en las esquinas superior derecha e izquierda son candidatos a genes cuya expresión cambia drásticamente de manera específica del tumor.

Al observar esta figura en el laboratorio, probablemente se reconozcan nombres familiares: BRCA1, TP53, ERBB2, MKI67, entre otros. En la siguiente entrega (S21), al realizar un análisis de enriquecimiento de vías (GSEA) con esta lista de genes, aparecerán claramente las vías bien conocidas del ciclo celular y la respuesta al daño del ADN. En ese momento se completará el verdadero panorama del artículo científico.

¿Por qué DESeq2 en lugar de edgeR o limma?

• DESeq2: Aplica una contracción (shrinkage) fuerte en genes de baja expresión. Es robusto con un número pequeño de réplicas (3 a 6) y es el estándar en la comunidad académica.
• edgeR: Es similar a DESeq2. Su método de dispersión tag-wise difiere ligeramente. En datos bulk, su rendimiento es casi equivalente al de DESeq2.
• limma-voom: Convierte los datos a log CPM y asume una distribución normal. Es más potente cuando el número de muestras es grande (30 o más).

El taller práctico de BPD comienza con DESeq2 debido al tamaño reducido de la muestra, pero en la siguiente entrega (S20) se observará empíricamente cómo divergen los tres métodos.

Mapeo conceptual

• Prueba de hipótesis múltiple (Multiple hypothesis testing): Concepto tratado en DryBench. BH es el método de corrección más práctico.
• Bayes empírico (Empirical Bayes): En lugar de estimar la dispersión para cada uno de los 20.000 genes de forma independiente, este método comparte información entre ellos. Es exactamente la misma idea que la "regularización" en CS.
• Estimador de contracción (Shrinkage estimator): Es la versión biológica del estimador de James-Stein. Atrae las estimaciones individuales hacia un objetivo común.

Errores frecuentes

• Ingresar valores normalizados (TPM, RPKM) en lugar de conteos brutos → DESeq2 fallará. En DESeqDataSet se deben ingresar obligatoriamente los conteos brutos (raw counts).
• Juzgar la significancia con padj en lugar de pvalue → Provoca una explosión de falsos positivos.
• No eliminar genes de baja expresión (conteo < 10 en todas las muestras) → Los picos aleatorios se interpretan erróneamente como significativos. Utilicemos filterByExpr().
• No controlar el efecto de lote (batch effect) → Podríamos estar observando diferencias entre "lote de secuenciación A vs B" en lugar de tumor vs normal. Especifiquemos esto con design = ~ batch + sample_type.

Para profundizar más

El texto principal es una descripción reconstruida directamente por BPD. Para un estudio más profundo, consulten los siguientes materiales clásicos en orden.

• Harvard STAT115: Xiaole Shirley Liu, normalización y estabilización de la varianza con la biblioteca DESeq2 (52 minutos, con subtítulos completos). Se explica en detalle el cálculo del factor de tamaño y la derivación de la transformación VST.
• StatQuest: Dos videos de Joshua Starmer (ambos con subtítulos en coreano):
  • Introducción a RNA-seq y DESeq2 (26 minutos)
  • Tasas de descubrimiento falso (FDR), explicación clara (18 minutos)
• Libro web gratuito de referencia: Holmes & Huber, Modern Statistics for Modern Biology. Los principios estadísticos de este libro están disponibles gratuitamente en huber.embl.de/msmb. Se recomienda leer detenidamente el Capítulo 8, "High-throughput count data".
• Artículo original: Love, Huber, Anders (2014), Estimación moderada del cambio y la dispersión para datos de RNA-seq con DESeq2, Genome Biology 15:550. Es la fuente original de la derivación del método de contracción (shrinkage).

En el siguiente episodio S20, compararemos empíricamente edgeR y limma-voom, y en S21, utilizaremos estos resultados para ejecutar GSEA. Hemos recorrido exactamente la mitad del camino en el proceso "de FASTQ a publicación".