Fin de la segunda parte de la serie "De FASTQ a artículo"

Con este episodio, se completa la fase de cuantificación de la segunda parte de la serie "De FASTQ a artículo". El gráfico de volcanes de DESeq2 en el próximo episodio (S19) será el clímax de esta serie. Gracias por llegar hasta aquí.

Hoy en día, con el alineamiento similar de Salmon/kallisto como estándar, ¿por qué todavía se utilizan las herramientas de conteo tradicionales (RSEM, featureCounts y HTSeq)? Hay tres razones.

• Reproducibilidad de artículos: Cuando es necesario reproducir exactamente los flujos de trabajo de estudios anteriores.
• Alineamiento ya existente: Si ya se dispone de archivos BAM alineados con STAR/HISAT2, sería una pérdida descartar ese alineamiento.
• Requisitos especiales: Estimación por máxima verosimilitud (EM) a nivel de isoformas (RSEM) o conteo estricto (HTSeq).

En este episodio, compararemos estas tres herramientas y estableceremos reglas prácticas para la toma de decisiones.

featureCounts: el rey de la velocidad

Parte del paquete Subread. Desarrollado por Yang Liao y Wei Shi.

Principio

Es el enfoque más sencillo.

1. Se escanea cada lectura en el archivo BAM alineado.
2. Si una lectura se superpone con un gen específico (la unión de exones del archivo GTF), se suma 1 al conteo de ese gen.
3. Si una lectura se superpone con varios genes, se distribuye o descarta según la opción elegida.

Está implementado mediante árboles de intervalos, lo que lo hace extremadamente rápido.

Comando de ejecución

featureCounts \
    -a gencode.v44.annotation.gtf.gz \
    -o counts.tsv \
    -p --countReadPairs \
    -s 2 \
    -T 16 \
    -Q 30 \
    sample1.bam sample2.bam sample3.bam

• -a: Archivo de anotación GTF.
• -p --countReadPairs: Lecturas pareadas; cada par corresponde a un conteo.
• -s 2: Orientación de las lecturas. 0=sin orientación, 1=sentido, 2=antisentido.
• -Q 30: MAPQ mínimo.
• -T: Número de hilos.

Cada fila del archivo TSV de resultados corresponde a un gen y cada columna a una muestra.

Orientación de las lecturas — Opción que suele causar errores

El protocolo Illumina TruSeq Stranded lee las lecturas en la dirección antisentido del ARNm. featureCounts -s 2 refleja esta inversión. Si se configura incorrectamente, el conteo de genes se reduce a la mitad y pueden obtenerse resultados anómalos, donde solo los genes antisentido (por ejemplo, XIST en ratón) se identifican como genes.

Método para determinar la orientación de las lecturas:

# infer_experiment.py (paquete RSeQC)
infer_experiment.py -i sample.bam -r hg38_gencode.bed

Si este resultado es "1++,1--,2+-,2-+", entonces se invierte la orientación (=-s 2).

Ventajas

• Velocidad: De 5 a 10 veces más rápido que otras herramientas.
• Procesamiento simultáneo de múltiples muestras: Permite procesar varios archivos BAM en una sola ejecución.
• Ligereza: Al estar escrito en C, tiene un bajo consumo de memoria.

Desventajas

• No admite isoformas: Solo funciona a nivel de gen.
• Manejo ambiguo de lecturas con múltiples mapeos: Ofrece muchas opciones, por lo que se requiere precaución.

HTSeq: conteo estricto

Desarrollado por Simon Anders y Wolfgang Huber. El mismo equipo que los autores de DESeq2.

Principio

Similar a featureCounts, pero con opciones más estrictas.

• Modo "union" (predeterminado): Se cuenta una lectura solo si cubre un único gen; si cubre varios, se descarta.
• Modo "intersection-strict": Todas las partes de la lectura deben estar dentro de un único gen para ser contadas.
• Modo "intersection-nonempty": Un punto intermedio.

El modo "union" predeterminado suele ser adecuado para la mayoría de los casos prácticos.

Comando de ejecución

htseq-count \
    --format=bam \
    --order=name \
    --stranded=reverse \
    --mode=union \
    sample.name_sorted.bam \
    gencode.v44.annotation.gtf.gz > sample.counts.tsv

• --order=name: Se requiere que los archivos BAM estén ordenados alfabéticamente (preparados con samtools sort -n).
• --stranded: sí/no/inverso.
• --mode: unión/intersección estricta/intersección no vacía.

Ventajas

• Rigurosidad: Cumple estrictamente con los supuestos estadísticos, ya que fue desarrollado por los autores de DESeq2.
• Implementación de referencia estándar: Es el estándar de referencia comparado con los resultados de featureCounts.

Desventajas

• Velocidad: Es hasta cinco veces más lento que otras herramientas.
• Requisito de ordenación alfabética: Reduce la flexibilidad del flujo de trabajo.

RSEM: EM a nivel de isoforma

Desarrollado por Bo Li en Wisconsin. Es el estándar para la cuantificación de isoformas antes de la era del alineamiento aproximado.

Principio

Similar al algoritmo EM de Salmon, pero realiza el alineamiento primero. Diferencias con Salmon:

• RSEM: Alinea las lecturas con los transcritos usando STAR o Bowtie2, y luego distribuye las isoformas mediante EM.
• Salmon: Realiza un hash de k-mers sin alineamiento previo, y luego distribuye las isoformas mediante EM.

La etapa de alineamiento lo hace considerablemente más lento.

Comando de ejecución

# Índice de referencia
rsem-prepare-reference \
    --gtf gencode.v44.annotation.gtf.gz \
    --star \
    --star-path $(which STAR) \
    -p 16 \
    hg38.fa hg38_rsem

# Cuantificación
rsem-calculate-expression \
    --paired-end --star --star-path $(which STAR) \
    -p 16 \
    --output-genome-bam \
    clean_R1.fq.gz clean_R2.fq.gz \
    hg38_rsem \
    SAMPLE01

Archivos a nivel de gen (.genes.results) y a nivel de isoforma (.isoforms.results).

Ventajas

• Cuantificación a nivel de isoforma: Calcula el IsoPct (porcentaje de isoforma) para cada transcrito.
• Implementación de referencia: Muchos artículos citan los resultados de RSEM.

Desventajas

• Velocidad: Muy lento. Entre 20 y 50 veces más lento que Salmon.
• Requisito de alineamiento: Requiere una infraestructura de alineamiento.

Matriz de uso práctico de las tres herramientas

Conclusión clave: Para proyectos nuevos, usar Salmon (S17); si ya existe el alineamiento, usar featureCounts; para la reproducción de referencias, usar la herramienta original.

Comparación de resultados: Salmon vs. featureCounts

Al ejecutar ambas herramientas y comparar los resultados, los conteos a nivel de gen coinciden aproximadamente en un 90-95%. Las principales causas de las diferencias son:

• Manejo de lecturas con múltiples mapeos: Salmon distribuye mediante EM, mientras que featureCounts las descarta.
• Sesgo de isoforma: Salmon es más preciso cuando una isoforma específica es predominante.

Aunque afecta los resultados del análisis DGEA, la mayoría de las señales robustas son consistentes entre los resultados de ambas herramientas.

Cálculo manual: qué hace realmente el modo Union

Escenario hipotético. La lectura R se superpone con dos genes, g1 y g2, en el archivo GTF (exones superpuestos, etc.).

• HTSeq, modo union: Descarta R. Conteo: 0.
• featureCounts, configuración predeterminada: Descarta R (multi-overlap).
• featureCounts --allowMultiOverlap: Suma +1 tanto a g1 como a g2 (o distribución fraccionaria según la opción).

La convención predeterminada es descartar. Nota práctica: puede haber pérdida de conteos cuando un gen está dentro de otro (por ejemplo, tRNA mitocondrial dentro de rRNA). Estos genes requieren un tratamiento especial.

Creación de entradas para DESeq2 en R

Las tres herramientas pueden transferirse a DESeq2 desde R.

featureCounts

library(DESeq2)

counts <- read.table("counts.tsv", header=TRUE, skip=1, row.names=1)
counts <- counts[, 6:ncol(counts)]  # Eliminar columna meta
colnames(counts) <- gsub(".bam", "", colnames(counts))

coldata <- data.frame(row.names=colnames(counts),
                      condition=c("normal","normal","tumor","tumor"))

dds <- DESeqDataSetFromMatrix(countData=counts, colData=coldata, design=~condition)

HTSeq

library(DESeq2)

files <- list.files("counts_dir", pattern="*.tsv", full.names=TRUE)
sampleTable <- data.frame(
    sampleName = gsub(".counts.tsv", "", basename(files)),
    fileName = files,
    condition = c("normal","normal","tumor","tumor"))

dds <- DESeqDataSetFromHTSeqCount(sampleTable=sampleTable,
                                  directory=".",
                                  design=~condition)

RSEM

library(tximport)
files <- list.files("rsem", pattern="*.genes.results", full.names=TRUE)
txi <- tximport(files, type="rsem", txIn=FALSE, txOut=FALSE)

dds <- DESeqDataSetFromTximport(txi, colData=coldata, design=~condition)

Los tres enfoques convergen finalmente en el objeto DESeqDataSet de DESeq2.

Pequeño ejercicio práctico en Colab

!apt-get install -y subread >/dev/null

# Asumir que existe el BAM de alineación STAR (preparado en S04)
!featureCounts \
    -a gencode.v44.annotation.gtf.gz \
    -o counts.tsv \
    -p --countReadPairs \
    -s 2 \
    -T 2 \
    -Q 30 \
    SAMPLE01_Aligned.sortedByCoord.out.bam

!head -5 counts.tsv

Columna de resultados: Geneid, Chr, Start, End, Strand, Length, sample1.bam. El recuento por gen se encuentra en la última columna.

Mapeo CS

• Árbol de intervalos: Consulta de exones GTF en featureCounts. Consulte la sección "Fundamentos del árbol de intervalos" de DryBench (interval-tree-basics).
• Manejo de mapeos múltiples: Reglas de decisión cuando una lectura se asigna a varios intervalos.
• Reutilización de EM: RSEM aplica un algoritmo de maximización de esperanza (EM) después de la alineación, similar a Salmon.

Resumen del recorrido de la serie "FASTQ to Paper"

Las 13 entregas anteriores (S06～S18) han completado el flujo de trabajo de cuantificación y variantes de la serie "FASTQ to Paper".

• S06～S11 (6 entregas de GATK4): Flujo de trabajo estándar para variantes de ADN.
• S17～S18 (2 entregas de cuantificación de RNA-seq): Cálculo de niveles de expresión.
• Próximas S19～S21 (3 entregas de DGEA): Se completará el gráfico de volcanes real.

Ahora pueden comprender todo el recorrido, desde los archivos FASTQ sin procesar hasta las figuras del artículo, como un flujo continuo.

Conclusión

RSEM, featureCounts y HTSeq fueron los estándares previos al alineamiento y siguen siendo válidos en ciertas situaciones. Tengan en cuenta que la existencia de múltiples herramientas se debe a sus respectivas fortalezas. La siguiente entrega (S19) es el punto culminante de la serie: completar el gráfico de volcanes con DESeq2. Ya se ha publicado como un episodio piloto; continúen con los resultados de cuantificación aprendidos hoy para verificarlos.

Para profundizar más

• Liao Y. et al. (2014), featureCounts: un programa eficiente de propósito general para asignar lecturas de secuencia a características genómicas. Bioinformatics 30:923. Artículo original de featureCounts.
• Anders S. et al. (2015), HTSeq: un marco de trabajo de Python para trabajar con datos de secuenciación de alto rendimiento. Bioinformatics 31:166. Artículo original de HTSeq.
• Li B. & Dewey C.N. (2011), RSEM: cuantificación precisa de transcritos a partir de datos de RNA-Seq con o sin un genoma de referencia. BMC Bioinformatics 12:323. Artículo original de RSEM.
• Xiaole Shirley Liu Harvard STAT115 W4: Clase sobre cuantificación de RNA-seq. Con subtítulos completos. Caso práctico sobre las diferencias entre los resultados de featureCounts y Salmon.
• Tutorial de tximport de Bioconductor: Método estándar para integrar los resultados de las tres herramientas mediante DESeq2.

Ejecute featureCounts y Salmon en una sola muestra y compare los recuentos para obtener una mejor comprensión práctica. En el siguiente episodio, confirmaremos la conclusión de esta serie con un gráfico de volcanes de DESeq2 para el estudio piloto S19.