De FASTQ a count matrix: el camino

Qué pasa entre el secuenciador y la matriz que cargas en R. QC con FastQC, pseudoalignment con Salmon vs alignment con STAR, cuantificación a nivel de gen y por qué cada elección afecta al análisis downstream.

Qué hay antes de R

Cuando recibes una “matriz de counts” para análisis, detrás hay un pipeline de varias horas que ya tomó decisiones por ti. Aunque no vayas a ejecutar tú esa parte, conviene entender qué se hizo. Las decisiones afectan al análisis downstream, saber qué entra evita interpretar mal.

El camino estándar:

FASTQ → QC → (trimming) → alineamiento o pseudoalineamiento → cuantificación → count matrix

Cada flecha es una herramienta, y cada herramienta tiene parámetros.

FASTQ: las lecturas crudas

El secuenciador (Illumina, MGI, PacBio, ONT…) produce archivos FASTQ, texto plano comprimido con cuatro líneas por lectura:

@SEQ_ID
ATCGATCGATCGATCG...
+
!''*((((***+))%%%...   ← calidades en codificación ASCII

Una muestra típica de RNA-seq con Illumina son 20-50 millones de lecturas por muestra, lecturas de 50-150 bp. Archivo: 1-10 GB comprimido.

Las lecturas pueden ser:

• Single-end (muestra.fastq.gz), una lectura por fragmento.
• Paired-end (muestra_R1.fastq.gz, muestra_R2.fastq.gz), los dos extremos del mismo fragmento. Más informativo, estándar actual.

QC: FastQC y MultiQC

Antes de cualquier cosa, control de calidad:

fastqc muestra_R1.fastq.gz muestra_R2.fastq.gz
multiqc .   # agrega FastQC + alineamiento + cuantificación en un HTML

Lo que se mira:

• Distribución de calidades por posición: calidad baja al final de la read es normal, baja en el medio es problema.
• Contenido GC: si se aparta mucho del esperado, hay contaminación.
• Sobre-representación de secuencias: adaptadores no recortados, contaminación de rRNA.
• Duplicación: alta en RNA-seq es habitual (genes muy expresados), no necesariamente patológico.

MultiQC es el agregador imprescindible. Combina reports de FastQC, Salmon, STAR, picard… en un único HTML navegable. Cuando recibes datos del centro de secuenciación, el MultiQC suele venir con ellos. Si no, pídelo.

Trimming: ¿hace falta?

El trimming de adaptadores y bases de baja calidad se hacía sistemáticamente hace años. Hoy hay matices:

• Para STAR con soft-clipping o Salmon/kallisto con pseudoalignment, el trimming no es necesario, los algoritmos lo manejan internamente.
• Para alineadores estrictos o protocolos especiales (UMI, smRNA), sí.

Herramientas: Trim Galore, fastp, cutadapt. Si vas a usar Salmon en el siguiente paso, normalmente puedes saltar trimming.

Alineamiento vs pseudoalineamiento

Hay dos filosofías para “asignar lecturas a genes”:

Alignment clásico (STAR, HISAT2)

Encuentra la posición exacta de cada lectura en el genoma. Output: archivo BAM con coordenadas.

STAR --runMode alignReads \
     --genomeDir indice_genoma/ \
     --readFilesIn muestra_R1.fq.gz muestra_R2.fq.gz \
     --outSAMtype BAM SortedByCoordinate \
     --quantMode GeneCounts          # opcional: counts a la vez

• ✓ Posición exacta, útil si quieres ver splicing, variantes, etc.
• ✓ Estándar histórico, todo el mundo entiende los BAMs.
• ✗ Lento (horas por muestra), necesita decenas de GB de RAM.
• ✗ Archivos BAM enormes (10-30 GB por muestra).

Pseudoalignment (Salmon, kallisto)

No alinea a posiciones, solo calcula la probabilidad de que cada lectura venga de cada transcrito. Output directo: cuentas por transcrito.

salmon quant -i indice_transcriptoma/ \
             -l A \
             -1 muestra_R1.fq.gz -2 muestra_R2.fq.gz \
             --gcBias --seqBias \
             -p 8 -o muestra_salmon/

• ✓ 10-100× más rápido que STAR. Una muestra en 5-15 minutos.
• ✓ No produce BAMs gigantes. Output ~ pocos MB.
• ✓ Modela sesgos (GC, posición) con --gcBias --seqBias.
• ✗ No te da posiciones, no apto si quieres ver splicing junctions.
• ✗ Requiere un buen transcriptoma anotado.

Para análisis diferencial estándar, Salmon es la elección actual. Es más rápido, igual o más preciso (en benchmarks recientes), y se integra muy bien con Bioconductor vía tximeta o tximport.

Cuantificación a nivel de gen

Importante: Salmon cuenta por transcrito, no por gen. Para análisis de expresión diferencial a nivel de gen (lo más común), hay que agregar transcritos al gen correspondiente.

Este paso lo hace tximport o tximeta dentro de R:

library(tximeta)

se <- tximeta(coldata)        # carga + anota automáticamente
gse <- summarizeToGene(se)    # agrega transcritos al gen

tximeta reconoce qué transcriptoma usó Salmon, descarga los metadatos, y produce un SummarizedExperiment listo. Es el patrón moderno, lo veremos a fondo en el tutorial 4.

Otra opción: featureCounts sobre BAM

Si tu pipeline produce BAMs (STAR o HISAT2), la cuantificación se hace a posteriori con featureCounts (paquete Rsubread):

library(Rsubread)
fc <- featureCounts(files = c("m1.bam", "m2.bam"),
                    annot.ext = "anotacion.gtf",
                    isPairedEnd = TRUE)
counts <- fc$counts

Más lento que Salmon, pero útil cuando ya tienes BAMs por otra razón.

¿Qué cuentas debe haber en la matriz?

La pregunta sutil: cuando la matriz dice gen_X = 1234, ¿qué unidad es?

• Raw counts (enteros): número de lecturas asignadas al gen. Esto es lo que necesita DESeq2 y casi cualquier modelo estadístico de RNA-seq. No uses counts normalizados como TPM o FPKM como entrada.
• TPM / FPKM / RPKM: counts normalizados por longitud del gen y profundidad. Útiles para visualizar y comparar genes dentro de una muestra, pero no para análisis diferencial.
• Estimated counts (decimales): Salmon devuelve estimaciones probabilísticas, no enteros (NumReads). tximport los redondea a enteros para downstream.

Regla: entrada a DESeq2 = raw integer counts. Si la matriz tiene decimales o se llama “TPM matrix”, está mal para análisis diferencial.

Metadata: la otra mitad

Una matriz de counts sin información de muestras no sirve. Necesitas un colData con al menos:

• Identificador de muestra (mismo que la columna del counts).
• Condición experimental (tratamiento, control, batch…).
• Idealmente: lote técnico, fecha de procesamiento, libreria, etc.

Sin metadata, no hay diseño experimental. Sin diseño, no hay análisis diferencial. Lo veremos en el tutorial 4.

Resumen del camino

Etapa	Herramienta típica	Output
Lecturas	secuenciador	FASTQ
QC	FastQC + MultiQC	HTML
(Trimming)	Trim Galore / fastp	FASTQ recortados
Cuantificación	Salmon	counts por transcrito
Importar a R	tximeta	SummarizedExperiment
Resumir a gen	summarizeToGene	counts por gen

A partir de “counts por gen”, entras a Bioconductor de lleno.

Trampas habituales

• Confundir counts con TPM. Si la matriz que recibes tiene decimales o totales que parecen normalizados (~10⁶ por columna en TPM), no es lo que necesita DESeq2.
• Datos de proveedor sin MultiQC. Si no hay HTML de QC, pídelo. Sin QC no se puede confiar en nada del downstream.
• Trimming agresivo cuando no hace falta. Quitar las bases de baja calidad del final con Trim Galore es razonable, pero si vas a usar Salmon, normalmente no aporta.
• Mezclar single-end y paired-end en un mismo experimento. Es una fuente de variabilidad técnica artificial. Si toca, modelar el tipo de lectura como covariable.
• Olvidar el strand-specificity. Si la librería es stranded, hay que decírselo al cuantificador (-l A en Salmon detecta automático. En featureCounts hay que indicarlo). Equivocarse aquí da resultados que parecen razonables pero están mal por la mitad.

En la siguiente entrega

Ya sabes de dónde vienen los counts. Bioconductor los almacena en un contenedor especial: SummarizedExperiment (y su primo DESeqDataSet). Entender bien este objeto es la clave para que el resto fluya, counts, metadata y anotación viven juntos en una sola estructura. Lo siguiente.