¿Por qué RNA-seq utiliza otros alineadores?

Las lecturas de ADN se alinean de forma continua al genoma de referencia. Sin embargo, el ARN ha perdido sus intrones durante la maduración del ARNm. Por lo tanto, para alinear las lecturas de ARN al genoma de referencia, es frecuente que una lectura deba dividirse en dos fragmentos que se unan a exones distintos. El segmento correspondiente al intrón está presente en la referencia pero ausente por completo en la lectura.

Este tipo de "alineación con huecos del orden de cientos a decenas de miles de pares de bases" no es manejado adecuadamente por BWA-MEM. La penalización de hueco afín de BWA-MEM aplica una penalidad logarítmica basada en la longitud del hueco, lo que resulta inadecuado para cubrir el rango típico de tamaños de intrones (de 200 pb a 100 kpb). Por ello, se requieren alineadores específicos para RNA-seq, como STAR e HISAT2.

¿Qué es una unión de empalme?

Se denomina unión de empalme al punto donde el extremo 3′ de un exón se une con el inicio 5′ del siguiente exón.

Genoma: ...[exón1==========]--------intrón (2 kb)--------[exón2==========]...
Read:                ACGTCC | la read atraviesa la unión | AATCGA

Cuando un read cruza un empalme, se divide de manera que 100 pb se asignan al extremo inicial (exón 1) y 50 pb al extremo final (exón 2), alineándose por separado. En la CIGAR resultante aparece una etiqueta N (región omitida) como 100M2000N50M. Esta N constituye la huella del empalme.

La detección de empalmes es difícil por tres razones:

1. No se conocen de antemano todas las posiciones de empalme: además de los empalmes registrados en GENCODE, existen empalmes nuevos (novel).
2. Los intrones tienen longitudes extremas: varían naturalmente desde 20 pb hasta 1 Mpb.
3. Las regiones cercanas al empalme son cortas: cuando solo unos pocos pb del read se solapan con el exón opuesto al empalme, la confianza de la alineación disminuye.

STAR y HISAT2 abordan estos tres problemas mediante enfoques distintos.

STAR — índice grande, alta precisión

STAR, publicado por Alex Dobin en 2013, emplea una estrategia basada en búsqueda con arreglo de sufijos (SA) + expansión del empalme. Los principios de su alineación se resumen en dos pasos:

1. Identificar el segmento continuo máximo del read que se une completamente al genoma (MMP, Maximum Mappable Prefix).
2. Localizar los fragmentos restantes del read en otras posiciones y confirmar el empalme si la distancia entre ellos corresponde a un candidato de intrón.

Este enfoque es ventajoso porque permite detectar empalmes de novo sin necesidad de cargar previamente información de empalmes desde archivos GTF.

Índice y requisitos de memoria

El costo de STAR es su gran índice. Para el genoma humano, es necesario cargar en memoria aproximadamente 30 GB del índice completo. Esta es la razón por la que STAR resulta costoso en entornos de nube; no cabe en la capa gratuita de Colab (12 GB).

Creación del índice:

STAR --runMode genomeGenerate \
     --genomeDir /data/star_index \
     --genomeFastaFiles hg38.fa \
     --sjdbGTFfile gencode.v44.annotation.gtf \
     --sjdbOverhang 100 \
     --runThreadN 16

Al insertar el archivo GTF de GENCODE mediante --sjdbGTFfile, se cargan previamente las informaciones de unión conocidas, lo que aumenta significativamente la precisión del primer alineamiento. Para --sjdbOverhang, es habitual establecerlo en longitud de lectura - 1 (por ejemplo, cualquier valor entre 100 y 149 para lecturas de 150 pb).

Ejecución básica + modo de dos pasadas

STAR --runMode alignReads \
     --genomeDir /data/star_index \
     --readFilesIn clean_R1.fq.gz clean_R2.fq.gz \
     --readFilesCommand zcat \
     --outSAMtype BAM SortedByCoordinate \
     --outSAMattrRGline ID:S1 SM:SAMPLE01 LB:lib1 PL:ILLUMINA \
     --twopassMode Basic \
     --runThreadN 16 \
     --outFileNamePrefix SAMPLE01_

--twopassMode Basic es la función destacada de STAR. El principio es el siguiente:

• Primera pasada: alineación solo con el genoma y los empalmes conocidos (GTF). Durante este proceso, se recopilan candidatos de empalme de novo.
• Segunda pasada: se añaden los empalmes de novo descubiertos en la primera pasada a los empalmes conocidos y se realiza una nueva alineación.

Al aprender nuevos empalmes incluso con una sola muestra y re-alinear, la sensibilidad para detectar isoformas nuevas aumenta notablemente. Es el estándar para análisis de RNA-seq novedosos.

HISAT2: una herramienta para la era de los grafos con índices pequeños

HISAT2 es un alineador basado en grafos desarrollado por Ben Langmead y Daehwan Kim. Aborda directamente la situación en la que el índice de 30 GB de STAR resulta oneroso.

¿Por qué es pequeño el índice?

HISAT2 utiliza un índice FM jerárquico (de ahí su nombre) junto con un GFM (Graph FM-index) que codifica los empalmes y variantes conocidos en un grafo. Gracias a esta estructura de datos, el índice del genoma humano más las variantes/empalmes conocidos se comprime a aproximadamente 8 GB.

• STAR: 30 GB (cercano al genoma original)
• HISAT2: 8 GB (compresión por grafo)

En la capa gratuita de Colab, el único alineador de RNA-seq que puede ejecutarse prácticamente es HISAT2.

Ejecución básica

hisat2-build hg38.fa hg38_hisat2   # Crear el índice

hisat2 -x hg38_hisat2 \
       -1 clean_R1.fq.gz -2 clean_R2.fq.gz \
       --rg-id S1 --rg SM:SAMPLE01 --rg LB:lib1 --rg PL:ILLUMINA \
       --known-splicesite-infile gencode_ss.tsv \
       -p 16 \
       | samtools sort -@ 4 -o SAMPLE01.bam -

--known-splicesite-infile es la opción de entrada de unión conocida correspondiente al --sjdbGTFfile de STAR. Se extrae del archivo GTF mediante el siguiente comando:

hisat2_extract_splice_sites.py gencode.v44.annotation.gtf > gencode_ss.tsv

STAR vs. HISAT2: reglas de decisión para la práctica

Cálculo manual: lo que realmente hace MMP

Construyamos una intuición sobre el MMP de STAR. Supongamos que una lectura de 150 bp alinea perfectamente sus primeros 90 bp en chr7:100-190, pero los 60 bp restantes no alinean inmediatamente después en la referencia.

STAR procede de la siguiente manera:

1. Fija el MMP de los primeros 90 bp en chr7:100-190.
2. Busca de nuevo los 60 bp restantes en otras posiciones de la referencia.
3. Encuentra que esos 60 bp alinean en chr7:2190-2250.
4. Interpreta el intervalo de 2000 bp como un posible intrón, compatible con una distribución de longitudes de 100 bp a 1 Mbp.
5. Confirma la unión: chr7:190 → chr7:2190.
6. Registra 90M2000N60M en el CIGAR.

Esta lógica permite que STAR encuentre uniones de novo incluso sin información previa. El modo de dos pasadas reúne esas uniones y las reutiliza en un segundo alineamiento.

GTF y GENCODE: el diccionario de referencia de las uniones

La calidad del alineamiento de RNA-seq depende en gran medida de la calidad del archivo GTF.

• GENCODE (https://www.gencodegenes.org/): referencia habitual para humanos y ratón; se actualiza cada año.
• Ensembl (https://www.ensembl.org/): generalmente concordante con GENCODE y compatible con más especies.
• RefSeq (NCBI): utilizado con frecuencia en informes clínicos.

La gestión de versiones es importante: incluso entre GENCODE v44 y v46 cambia el catálogo de uniones. Fijar y documentar la versión del GTF al iniciar el proyecto facilita la reproducción posterior.

Práctica sencilla con HISAT2 en Colab

Como un índice del genoma humano completo sigue siendo pesado incluso con 8 GB, en el ejemplo crearemos un índice únicamente del cromosoma 22, de unos 50 MB.

!wget -q https://hgdownload.soe.ucsc.edu/goldenPath/hg38/chromosomes/chr22.fa.gz
!gunzip chr22.fa.gz
!hisat2-build chr22.fa chr22_hisat2

!wget -q https://sra-pub-src-1.s3.amazonaws.com/SRR1039508/SRR1039508_1.fastq.gz -O r1.fq.gz
!hisat2 -x chr22_hisat2 -U r1.fq.gz -p 2 -S out.sam
!samtools flagstat out.sam

El resultado de flagstat muestra la proporción de lecturas alineadas. Como solo se indexó el cromosoma 22, es normal que la mayoría no alinee.

Mapeo CS

• Consulta de array de sufijos (STAR): localiza posiciones de alineamiento ordenando los sufijos de la cadena. Consulte la unidad de DryBench sobre índices de grafos (graph-index-fundamentals).
• FM-index de grafo (HISAT2): codifica el genoma y las variantes conocidas como un grafo y busca sobre él; combina compresión de estructuras de datos y búsqueda en grafos.
• Aprendizaje en dos pasadas: mejora el segundo alineamiento utilizando parámetros aprendidos en el primero, una aplicación de la idea de autoentrenamiento al alineador.

Conclusión

BWA-MEM alinea el ADN de forma continua, mientras que STAR y HISAT2 pueden atravesar uniones de splicing. En la siguiente unidad (S05) utilizaremos samtools para ordenar, indexar y comprimir como CRAM los archivos BAM creados hasta ahora. La canalización de GATK4 comienza en S06 y la cuantificación de RNA-seq continúa en S17–S21.

Para profundizar

• Dobin A. et al. (2013), STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29:15. Artículo original de STAR, con detalles sobre MMP y las dos pasadas.
• Kim D. et al. (2019), Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nature Biotechnology 37:907. FM-index de grafo de HISAT2.
• Xiaole Shirley Liu, Harvard STAT115 W4: clase sobre alineamiento de RNA-seq con subtítulos y una interpretación útil de las proporciones de alineamientos únicos y múltiples.
• Broad Institute BroadE — RNA-seq alignment: serie de vídeos con casos prácticos de ajuste de opciones de STAR.
• Documentación de GENCODE (https://www.gencodegenes.org/human/): fuente principal del catálogo de uniones.

Ejecutar en Colab un alineamiento del cromosoma 22 con HISAT2 ayuda a adquirir intuición sobre las uniones de splicing. En la siguiente unidad trataremos este BAM con samtools como en un flujo de trabajo real.