Comprender completamente Salmon y kallisto: ¿por qué es tan rápido el pseudo-mapeo?

Salmon y kallisto: principios y práctica del pseudoalineamiento que omite el alineamiento

¿De dónde proviene la idea de que es posible realizar cuantificación sin alineamiento? Desde el principio de pseudoalineamiento que mapea lecturas a transcritos mediante hash de k-mer, pasando por la estadística para estimar la abundancia a nivel de transcrito con EM, hasta la integración con DESeq2 mediante tximport.

Intermedio

20min

Verificado (2026-07-22)

RNA-seq quantificationpseudo-alignmentSalmon

Ingreso a la segunda parte de la serie FASTQ to Paper

Dejamos atrás la tubería de ADN (S06～S11, GATK4) para adentrarnos en la cuantificación de RNA-seq. El objetivo de esta serie es el gráfico de volcanes DESeq2 de S19. Este episodio (S17) es el primer paso de ese viaje.

El objetivo del RNA-seq es la cantidad de expresión por gen. ¿Cómo se calcula esa expresión? Tradicionalmente, se alineaban las lecturas al genoma de referencia con STAR y se contaban las lecturas por gen usando featureCounts o RSEM. Sin embargo, alrededor de 2016, un enfoque completamente diferente se convirtió en el estándar: pseudoalineamiento.

¿Por qué pseudoalineamiento?

Observemos el desperdicio del alineamiento tradicional. El resultado final del RNA-seq es "¿cuántas lecturas se adhirieron al gen X?". Sin embargo, STAR calcula exactamente en qué coordenadas de la referencia se adhiere cada lectura.

Las coordenadas exactas son necesarias para la llamada de variantes. Para la cuantificación de RNA-seq, esa precisión no es necesaria; solo necesitamos saber a qué gen o transcrito pertenece.

Esta intuición es el punto de partida del pseudoalineamiento. Si en lugar de coordenadas exactas solo determinamos a qué transcrito pertenece una lectura, la cuantificación es mucho más rápida y requiere mucha menos memoria.

Principio — Hashing de k-mer y clases de compatibilidad

Aprendamos a intuir la estructura de datos central de kallisto.

1. Indexación de transcritos

Descomponer cada transcrito (~200,000 en GENCODE) en k-mers (k=31 es el estándar). Almacenar en un mapa hash a qué transcrito pertenece cada k-mer.

\text{index}[\text{kmer}] = \{t_1, t_5, t_{209}\}

Un k-mer puede pertenecer a múltiples transcritos (por ejemplo, exones compartidos entre isoformas).

2. Búsqueda de hash de k-mer en las lecturas

Extraer varios k-mers de una lectura y buscar cada uno en el hash. Determinar a qué conjunto de transcritos pertenece cada k-mer.

3. Clase de compatibilidad

La intersección de los conjuntos de transcritos obtenidos de los k-mers de una lectura es la clase de compatibilidad (clase de equivalencia) de esa lectura.

Ejemplo: si los k-mers de la lectura R devuelven {t1, t5}, {t1, t5, t9} y {t1, t5}, entonces la clase de compatibilidad de R es {t1, t5}. Esto significa que R proviene de t1 o t5.

No se calculan las coordenadas exactas.

Algoritmo EM — Estimación de abundancia de transcritos

Aún queda un problema con solo las clases de compatibilidad: ¿de dónde vino exactamente R, de t1 o de t5? Con una sola lectura es imposible determinarlo.

Sin embargo, con millones de lecturas, la estadística funciona. Se utiliza el algoritmo EM (Maximización de Expectación).

Paso E

Calcula la probabilidad de que cada lectura provenga de cada transcrito utilizando las estimaciones actuales de abundancia.

P(t | R) = \frac{a_t \cdot \mathbb{1}[t \in \text{class}(R)]}{\sum_{t'} a_{t'} \cdot \mathbb{1}[t' \in \text{class}(R)]}

Paso M

La nueva abundancia de cada transcrito = número de lecturas esperadas que provienen de ese transcrito / longitud del transcrito.

a_t^{\text{new}} = \frac{\sum_R P(t|R)}{L_t}

Repite estos dos pasos → convergencia → abundancia final.

Este es el algoritmo de cuantificación común a kallisto y Salmon. RSEM también aplica el mismo EM después del alineamiento de lecturas.

Diferencias entre Salmon y kallisto

Aunque las dos herramientas son similares en principio, existen diferencias en su implementación.

Característica	Salmon	kallisto
Creador	Rob Patro (Stony Brook)	Nicolas Bray, Lior Pachter (UC Berkeley)
Hash de k-mer	Mapeo cuasi	Pseudoalineamiento
Corrección de sesgo	Integración de sesgo de secuencia, sesgo GC y sesgo posicional	Por separado
Bootstrap	Soportado	Soportado
Velocidad	Un poco más lento pero preciso	Un poco más rápido pero con corrección de sesgo limitada

En la práctica, se utiliza principalmente Salmon. Dado que integra la corrección de sesgos, es estable en bibliotecas con alto sesgo GC.

Comando de ejecución de Salmon

Generación del índice

bash

salmon index \
    -t gencode.v44.transcripts.fa.gz \
    -i gencode_v44_salmon \
    -k 31 \
    --gencode

-t: Transcritos en formato FASTA.
-i: Salida del índice.
--gencode: Limpieza de los encabezados de GENCODE (solo IDs de transcritos).
El tamaño aproximado del índice de transcritos humanos es de 2 GB.

Cuantificación

bash

salmon quant \
    -i gencode_v44_salmon \
    -l A \
    -1 clean_R1.fq.gz -2 clean_R2.fq.gz \
    --validateMappings \
    --gcBias --seqBias --posBias \
    -p 16 \
    -o SAMPLE01_salmon

-l A: Detección automática del tipo de biblioteca.
--validateMappings: Mejora de la precisión del quasi-mapeo.
--gcBias --seqBias --posBias: Corrección de tres tipos de sesgo.

En la carpeta de resultados, se encuentra quant.sf — abundancia por transcrito (TPM · longitud efectiva · número de lecturas).

Comando de ejecución de kallisto

bash

# Índice
kallisto index -i gencode_v44_kallisto gencode.v44.transcripts.fa.gz

# Cuantificación
kallisto quant \
    -i gencode_v44_kallisto \
    -o SAMPLE01_kallisto \
    -b 100 \
    -t 16 \
    clean_R1.fq.gz clean_R2.fq.gz

-b 100: bootstrap 100 veces (estimación de la varianza).

TPM · CPM · counts — ¿qué es qué?

Indicadores de cuantificación de RNA-seq.

counts (lecturas sin procesar): número de lecturas asignadas a un transcrito. Sesgado por el tamaño de la biblioteca y la longitud del transcrito.
CPM (Counts Per Million): normalización por el tamaño de la biblioteca. count / total_reads * 1M.
RPKM / FPKM: normalización por el tamaño de la biblioteca y la longitud del transcrito.
TPM (Transcripts Per Million): normalización primero por la longitud del transcrito y luego por el tamaño de la biblioteca. Versión mejorada de RPKM. Adecuado para comparaciones entre muestras.

Para DGEA (S19 DESeq2) se utilizan los counts sin procesar. Los valores correspondientes son los generados por Salmon/kallisto, es decir, est_counts o NumReads.

tximport — Convertir resultados de Salmon/kallisto a R

Para ejecutar DESeq2 en R, es necesario convertir los resultados de Salmon/kallisto a un data.frame de R. La herramienta estándar para ello es tximport.

library(tximport)

# Archivos quant.sf de Salmon
files <- file.path("data", samples$id, "quant.sf")
names(files) <- samples$id

# Correspondencia transcrito-gen
tx2gene <- read.table("tx2gene.tsv", header=TRUE)

# Integrar como recuentos a nivel de gen
txi <- tximport(files, type = "salmon", tx2gene = tx2gene)

head(txi$counts)

Se inserta este objeto txi directamente en DESeq2 (continuación desde S19).

Crear tx2gene

Extraer el mapeo transcrito-gen del GTF de GENCODE.

bash

zcat gencode.v44.annotation.gtf.gz | awk -F'\t' '$3=="transcript"' | \
    grep -oP 'transcript_id "[^"]+"|gene_id "[^"]+"|gene_name "[^"]+"' | \
    ... crear tx2gene.tsv mediante distintos métodos

En la práctica, se maneja mediante funciones del paquete R GenomicFeatures:

library(GenomicFeatures)
txdb <- makeTxDbFromGFF("gencode.v44.annotation.gtf.gz")
k <- keys(txdb, keytype = "TXNAME")
tx2gene <- select(txdb, k, "GENEID", "TXNAME")

Cálculo manual — Lo que realmente hace la iteración de EM

Un escenario sencillo. Los transcritos t1 y t2 comparten un exón. La lectura R pertenece a la clase de compatibilidad {t1, t2} derivada del k-mer de ese exón compartido.

Abundancia inicial: a_t1 = a_t2 = 0.5.
Paso E: Probabilidad de que R provenga de t1 = 0.5, probabilidad de que provenga de t2 = 0.5. Se distribuye por igual.
Paso M: Recálculo de la abundancia de cada t.

A medida que múltiples lecturas se distribuyen de esta manera, los valores convergen hacia las abundancias reales de cada t. En particular, si t1 está mucho más expresado, las lecturas únicas de t1 impulsan al alza la abundancia de t1, y la EM asigna también las lecturas del exón compartido predominantemente a t1.

La belleza de este algoritmo radica en que permite la cuantificación a nivel de transcrito sin necesidad de coordenadas exactas.

Práctica en Colab

Salmon es fácil de ejecutar en Colab.

bash

!pip install pyroe -q
!apt-get install -y salmon >/dev/null

!wget -q https://sra-pub-src-1.s3.amazonaws.com/SRR1039508/SRR1039508_1.fastq.gz -O r1.fq.gz
!wget -q https://sra-pub-src-1.s3.amazonaws.com/SRR1039508/SRR1039508_2.fastq.gz -O r2.fq.gz
!wget -q https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_44/gencode.v44.transcripts.fa.gz

!salmon index -t gencode.v44.transcripts.fa.gz -i gencode_idx -k 31 --gencode
!salmon quant -i gencode_idx -l A -1 r1.fq.gz -2 r2.fq.gz --validateMappings -p 2 -o SAMPLE01
!head SAMPLE01/quant.sf

Treinta minutos bastan para familiarizarse con el procedimiento.

Mapeo CS

Correspondencia aproximada mediante hash: una consulta hash de k-mer es O(1), frente al coste dependiente de la longitud de lectura y del genoma de un alineamiento genómico.
Algoritmo EM: consulte el episodio de DryBench «Fundamentos del algoritmo EM» (em-algorithm-basics).
Clase de compatibilidad: combina mediante operaciones de conjuntos los resultados de la búsqueda hash.

Conclusión

La idea de cuantificar sin realizar un alineamiento completo es sorprendente. Este enfoque redujo el análisis de RNA-seq de horas a minutos y facilitó la migración de muchos pipelines a la nube. En el siguiente episodio (S18) veremos la alternativa al pseudoalineamiento: realizar un alineamiento convencional y después contar con featureCounts, RSEM o HTSeq. Compararemos cuándo conviene elegir cada enfoque.

Para profundizar

Patro R. et al. (2017), Salmon provides fast and bias-aware quantification of transcript expression. Nature Methods 14:417. Artículo original de Salmon.
Bray N.L. et al. (2016), Near-optimal probabilistic RNA-seq quantification. Nature Biotechnology 34:525. Artículo original de kallisto.
Soneson C. et al. (2015), Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences. F1000 Research 4:1521. Artículo original de tximport.
Clase Xiaole Shirley Liu Harvard STAT115 W4 — RNA-seq quantification, con subtítulos completos y una buena introducción intuitiva a Salmon y kallisto.
Blog de Pachter Lab (https://liorpachter.wordpress.com/): comentarios metodológicos del autor de kallisto.

Abra un archivo quant.sf de Salmon y examine la columna TPM para familiarizarse con la cuantificación. En el siguiente episodio pasaremos a las herramientas de recuento posteriores al alineamiento convencional.

BioPlayground

Salmon y kallisto: principios y práctica del pseudoalineamiento que omite el alineamiento

Ingreso a la segunda parte de la serie FASTQ to Paper

¿Por qué pseudoalineamiento?

Principio — Hashing de k-mer y clases de compatibilidad

1. Indexación de transcritos

2. Búsqueda de hash de k-mer en las lecturas

3. Clase de compatibilidad

Algoritmo EM — Estimación de abundancia de transcritos

Paso E

Paso M

Diferencias entre Salmon y kallisto

Comando de ejecución de Salmon

Generación del índice

Cuantificación

Comando de ejecución de kallisto

TPM · CPM · counts — ¿qué es qué?

tximport — Convertir resultados de Salmon/kallisto a R

Crear tx2gene

Cálculo manual — Lo que realmente hace la iteración de EM

Práctica en Colab

Mapeo CS

Conclusión

Para profundizar

💬 Preguntas y comentarios