¿Por qué es necesario este capítulo?

En el capítulo M14 analizamos la BWT. Dijimos que era una estructura de datos con una doble función: compresión y búsqueda. Sin embargo, aún no habíamos explicado en detalle por qué la búsqueda es tan rápida.

La respuesta es el índice FM, publicado en 2000 por Paolo Ferragina y Giovanni Manzini. Si añadimos a la BWT una matriz de rango y una matriz de conteo, podemos encontrar las posiciones de coincidencia del patrón P en tiempo O(|P|). Para buscar un patrón de 100 pb en el genoma humano de 3.000 millones de pb, solo se necesitan 100 consultas a la matriz.

Una vez que este capítulo esté completo, quedará claro por qué BWA es tan rápido y por qué el índice del genoma humano se almacena en un tamaño similar al original.

Los tres componentes del índice FM

El índice FM se compone de tres elementos:

1. BWT: La transformación de Burrows-Wheeler de la cadena original.
2. Matriz C: Para cada carácter, el número total de caracteres que son lexicográficamente menores que él en la cadena ordenada. Es decir, la posición en la que comienza ese carácter en la primera columna.
3. Matriz Occ (Rango): Occ(c, i) = número de veces que aparece el carácter c entre los primeros i caracteres de la BWT.

Veamos un ejemplo de T = BANANA$.

• BWT = ANNB$AA
• Primera columna ordenada = $AAABNN

Matriz C:

C['`'] = 0  # ` es el primero en orden lexicográfico
C['A'] = 1  # Antes de A: un $
C['B'] = 4  # Antes de B: $, A, A, A (4)
C['N'] = 5  # Antes de N: $, A, A, A, B (5)

Orden de ocurrencia (BWT = ANNB$AA):

i:      0 1 2 3 4 5 6 7
BWT:    A N N B $ A A
Occ(A): 1 1 1 1 1 2 3 3
Occ(N): 0 1 2 2 2 2 2 2
Occ(B): 0 0 0 1 1 1 1 1
Occ($): 0 0 0 0 1 1 1 1

Búsqueda hacia atrás: búsqueda del patrón desde el final

La búsqueda en el índice FM avanza desde el final del patrón. Esto se denomina búsqueda hacia atrás.

Busquemos el patrón P = ANA dentro de BANANA.

Inicialización: último carácter A

Se encuentra el rango de la primera columna en el BWT donde aparece A. En el arreglo C, A comienza en 1 y, dado que el recuento total de A es 3, el rango es [1, 4).

• start = C['A'] = 1
• end = C['A'] + count('A') = 1 + 3 = 4

Este rango apunta a las rotaciones ordenadas que comienzan con A.

1: A$BANAN
2: ANA$BAN
3: ANANA$B

Siguiente: NA (el penúltimo carácter desde el final: N)

Ahora, conservamos solo las rotaciones dentro del rango que tienen a N inmediatamente antes. Utilizamos el mapeo LF.

Nuevo rango: [6, 7). Es decir, solo queda una rotación ordenada, la número 6.

6: NA$BANA

Última: ANA (el tercer carácter contando desde el final es A)

Nuevo rango [3, 4). Una rotación ordenada, la tercera.

3: ANANA$B

Dado que SA[3] = 1, la cadena original comienza en la posición 1 con ANA.

Implementación completa en Python

from collections import Counter

def build_fm_index(text: str):
    text = text + "$"
    n = len(text)
    sa = sorted(range(n), key=lambda i: text[i:])
    bwt = "".join(text[(sa[i] - 1) % n] for i in range(n))

    # Arreglo C
    sorted_chars = sorted(text)
    C = {}
    for i, c in enumerate(sorted_chars):
        if c not in C:
            C[c] = i

    # Arreglo Occ
    alphabet = sorted(set(text))
    Occ = {c: [0] * (n + 1) for c in alphabet}
    for i, c in enumerate(bwt):
        for ch in alphabet:
            Occ[ch][i+1] = Occ[ch][i]
        Occ[c][i+1] += 1

    return {"bwt": bwt, "C": C, "Occ": Occ, "sa": sa, "n": n}

def fm_search(index, pattern: str) -> list[int]:
    C, Occ, sa, n = index["C"], index["Occ"], index["sa"], index["n"]
    start, end = 0, n
    for c in reversed(pattern):
        if c not in C:
            return []
        start = C[c] + Occ[c][start]
        end   = C[c] + Occ[c][end]
        if start >= end:
            return []
    return sorted(sa[start:end])

# Uso
index = build_fm_index("BANANA")
print(fm_search(index, "ANA"))  # [1, 3]

Se repite tantas veces como la longitud del patrón (3 veces), y en cada iteración solo se realizan unas pocas consultas al arreglo para obtener la posición de coincidencia. Tiempo total O(|P|), independiente del tamaño de referencia.

Optimizaciones reales de BWA

La implementación anterior es con fines didácticos, por lo que el arreglo Occ ocupa varias veces el tamaño original. BWA real utiliza las siguientes optimizaciones:

• Árbol de Wavelet: Almacena Occ en forma comprimida (especializado para el alfabeto de 4 letras del ADN).
• Muestreo del arreglo SA: No almacena todo el arreglo de sufijos, sino solo muestras a intervalos regulares. La posición de coincidencia se reconstruye moviéndose a las muestras adyacentes mediante el mapeo LF.
• Puntos de control: Solo se guardan puntos de control del arreglo Occ a intervalos regulares.

Resultado: El índice BWA para el genoma humano de 3 GB ocupa aproximadamente 3 GB (un tamaño similar al original). Una compresión muy impresionante.

Aplicaciones prácticas: alineamiento de lecturas cortas

BWA (2009), Bowtie (2009) y Bowtie2 (2012) utilizan todos el índice FM. El tiempo de alineamiento de una sola lectura de Illumina de 100 pb está en el orden de decenas de miles por segundo. Alinea los resultados de secuenciación con una cobertura de 30× del genoma humano (aproximadamente 300 millones de lecturas) en cuestión de horas.

• BWA-backtrack (primera generación de BWA): Para lecturas cortas (<70 pb).
• BWA-SW: Para longitudes intermedias (70 pb a 10 kb). Combina la extensión con Smith-Waterman.
• BWA-MEM (estándar práctico): Desde lecturas cortas hasta lecturas largas intermedias. Busca semillas con el índice FM y extiende los huecos con una función de penalización afín.

El episodio S03 aborda el uso práctico de BWA-MEM.

Problemas comunes en la práctica

• Archivos del índice: BWA guarda el índice en varios archivos (.amb, .ann, .bwt, .pac, .sa). Todos deben estar presentes y, por convención, colocarse en el mismo directorio que el archivo del genoma de referencia.
• Reutilización del índice: Mientras el genoma de referencia no se actualice, el índice se reutiliza continuamente. El comando bwa index solo debe ejecutarse una vez.
• Compatibilidad de versiones: Los índices de BWA suelen ser compatibles, pero si cambia el método de compresión BWT, podría ser necesario reconstruir el índice.
• Casos límite: Los genomas con muchas secuencias repetitivas pueden tener demasiadas coincidencias de semillas, lo que ralentiza la búsqueda. BWA-MEM limita el número máximo de candidatos mediante la opción -c.

Mapeo con ciencias de la computación

• Estructuras de datos Rank/Select: La base de la teoría de la información y la teoría de estructuras de datos. Permite acceder rápidamente a caracteres y conteos en posiciones arbitrarias dentro de una representación comprimida.
• Árbol de Wavelet: Extensión de Rank/Select para cadenas con alfabetos grandes. Presentado por Grossi, Gupta y Vitter en 2003.
• Índices comprimidos: Lograr que la suma del tamaño del original y del índice sea menor o igual al doble del original es un logro teóricamente interesante desde el punto de vista de la teoría de la información.
• Búsqueda inversa: Aplicación iterada del mapeo LF. La elegancia de la estructura de datos es clave en el diseño del algoritmo.

Ramificaciones hacia los siguientes capítulos

• Capítulo siguiente (M16): minimap2 · minimizador — Diseño de índices para lecturas largas. Un enfoque diferente al FM-index.
• Dos capítulos después (M17): DIAMOND — Optimización de índices para la búsqueda de proteínas.
• Seis capítulos después (S03): Práctica con BWA-MEM — Aplicación práctica de las estructuras de datos aprendidas.
• Diez capítulos después (S07～S12): GATK4 — Detección de variantes a partir de los resultados alineados con BWA.

Para profundizar

El texto principal es una narración reconstruida internamente por BPD. Para obtener más información, consulte lo siguiente:

• CMU 02-510 — El profesor Ben Langmead (autor de Bowtie) imparte FM-index and Read Alignment (con subtítulos completos y una excelente traducción automática). Un curso fundamental impartido por el propio autor de Bowtie.
• Artículo original: Ferragina, P. & Manzini, G. (2000), Opportunistic data structures with applications, FOCS. El origen del FM-index.
• Artículo original: Li, H. & Durbin, R. (2009), Fast and accurate short read alignment with Burrows-Wheeler transform, Bioinformatics 25, 1754–1760. El origen de BWA.
• Artículo original: Langmead, B. et al. (2009), Ultrafast and memory-efficient alignment of short DNA sequences to the human genome, Genome Biology 10, R25. El origen de Bowtie.
• Texto de referencia: Mäkinen et al. Genome-Scale Algorithms Capítulo 8. Un libro de texto de referencia.

Resuelva los problemas relacionados con la búsqueda de cadenas en Rosalind. Así comprenderá de forma natural por qué minimap2, en el siguiente capítulo M16, no utiliza el FM-index.