Comenzando desde el muro de OLC

En M22, confirmamos que OLC (Overlap-Layout-Consensus) se topa con la pared NP-completa del problema del camino hamiltoniano. Con millones de lecturas, los métodos exactos no son viables.

Lo que supera esta pared es el Grafo de De Bruijn. Al cambiar solo una definición de nodo, el problema pasa a resolverse en tiempo polinomial. Esta es la estética del diseño algorítmico. En esta sección derivamos la esencia de esa transformación.

Definición de k-mer

La primera redefinición de De Bruijn se centra en los nodos. Los nodos son k-mers (subsecuencias de longitud k), no lecturas.

Si fragmentamos la lectura ACGTACG con k=4, obtenemos los siguientes 4-mers:

ACGT · CGTA · GTAC · TACG

De un solo read (L - k + 1) k-mers se generan. Si los reads son de millones, los k-mers serían aún más numerosos, pero al eliminar las duplicaciones se comprimen al número real de k-mers distintos que existen en el genoma. Esta es la ventaja clave.

Definición del Grafo de De Bruijn

• Nodos: (k-1)-mers distintos
• Aristas: cada k-mer uno (una arista dirigida desde los primeros k-1 caracteres del k-mer hacia los últimos k-1 caracteres)

Por ejemplo, el 4-mer ACGT se convierte en una arista que conecta el nodo de 3-mer ACG → CGT.

Si insertamos los cuatro 4-mers del read ACGTACG en un Grafo de De Bruijn.

ACGT: ACG → CGT
CGTA: CGT → GTA
GTAC: GTA → TAC
TACG: TAC → ACG

4 nodos (ACG, CGT, GTA, TAC), 4 aristas. La secuencia original se representa como un camino que recorre cada arista del grafo exactamente una vez. Esto es un camino euleriano.

Camino euleriano — la vía rápida de tiempo polinomial

Teorema del camino euleriano (Euler, siglo XVIII, problema de los puentes de Königsberg):

La condición necesaria y suficiente para que exista un camino euleriano en un grafo dirigido es: (1) el grafo sea conexo, y (2) el nodo inicial tenga outdeg - indeg = 1, el nodo final tenga indeg - outdeg = 1, y los demás nodos tengan outdeg = indeg.

Si se cumplen estas condiciones, se puede encontrar un camino euleriano en tiempo O(número de aristas) mediante el algoritmo de Hierholzer.

La idea de Hierholzer es sencilla:

1. Comenzar en un nodo arbitrario y seguir las aristas disponibles hasta donde sea posible.
2. Cuando no se pueda avanzar más, se ha formado un ciclo. Desde ese ciclo, volver a un nodo que aún tenga aristas sin usar y formar otro ciclo.
3. Fusionar los dos ciclos.
4. Repetir el proceso hasta haber utilizado todas las aristas.

Complejidad O(N + E), lineal respecto a la suma del número de nodos y aristas. Al pasar del problema NP-completo del camino hamiltoniano al camino euleriano, que se resuelve en tiempo polinomial, es posible procesar millones de lecturas en cuestión de minutos.

Ensamblaje manual paso a paso

Supongamos que tenemos 4 lecturas observadas.

R1: ACGTAC
R2: CGTACG
R3: GTACGT
R4: TACGTA

k=4 para fragmentos de 4-meros. Tres por lectura, un total de 12 (con duplicados).

R1: ACGT · CGTA · GTAC R2: CGTA · GTAC · TACG R3: GTAC · TACG · ACGT R4: TACG · ACGT · CGTA

Tras eliminar duplicados, hay 4 4-meros distintos: ACGT, CGTA, GTAC, TACG.

Grafo de De Bruijn con nodos de 3-meros:

ACG → CGT (from ACGT)
CGT → GTA (from CGTA)
GTA → TAC (from GTAC)
TAC → ACG (from TACG)

4 nodos y 4 aristas. La topología del grafo es un ciclo completo. Existe un circuito de Euler.

Ruta: ACG → CGT → GTA → TAC → ACG.

Al concatenar los 3-meros en esta ruta (manteniendo los primeros 3 caracteres y conectando solo el último carácter de cada nodo).

ACG + T + A + C + ... = ACGTACG (genoma original)

Se reconstruye el genoma completo. El ensamblaje manual queda así.

Por qué aparece un tiempo polinómico

La clave está en la definición de arista = k-mer observado. En OLC, la superposición era una "relación posible" que podía existir o no de manera probabilística. En De Bruijn, cada arista es un "hecho observado". Encontrar un camino que utilice cada observación exactamente una vez se mapea con precisión al problema del camino euleriano.

Mapear con precisión la definición del problema a un problema de tiempo polinómico en teoría de grafos. Esa es la elegancia de De Bruijn. Solo al expresar la información de manera diferente, cambia la complejidad fundamental del algoritmo.

Compensación en la elección de k

El rendimiento práctico del Grafo de De Bruijn es extremadamente sensible al valor de k.

Si k es demasiado pequeño:

• Los mismos k-mers provenientes de diferentes posiciones se superponen, aumentando los ciclos repetitivos en el grafo
• En las regiones repetitivas, aparecen múltiples caminos como candidatos, impidiendo identificar la respuesta correcta
• El ensamblaje final se fragmenta en contigs más pequeños

Si k es demasiado grande:

• Los k-mers con errores quedan aislados en el grafo e inutilizables
• El número de k-mers generados por cada lectura disminuye, dejando el grafo disperso y fragmentado
• Se producen huecos donde la cobertura es insuficiente

En la práctica, herramientas como SPAdes crean grafos con múltiples valores de k (21, 33, 55, 77, 99, etc.) y los fusionan. Este enfoque de múltiples k fue la clave del éxito de SPAdes.

Burbujas y puntas — Recuperación de errores del grafo

En la práctica con el Grafo de De Bruijn, las lecturas con errores crean dos defectos topológicos en el grafo.

Burbuja: Si hay dos caminos paralelos entre dos nodos, uno es la respuesta correcta y el otro es un error. Los caminos paralelos cortos se determinan como errores y se eliminan.

Punta: Una rama que sobresale brevemente desde un extremo del grafo. Si la cobertura es baja, es probable que una rama creada por una lectura con error tenga una longitud inferior al umbral. Se eliminan las puntas de longitud inferior al umbral.

Estos dos trucos son responsables del 80% de la precisión de los ensambladores basados en el Grafo de De Bruijn.

Implementación mínima en Python (~40 líneas)

from collections import defaultdict

def build_de_bruijn(reads, k):
    """Construye un grafo de De Bruijn a partir de lecturas. Nodos = (k-1)-mers; aristas = k-mers."""
    graph = defaultdict(list)
    indeg = defaultdict(int)
    outdeg = defaultdict(int)
    for read in reads:
        for i in range(len(read) - k + 1):
            kmer = read[i:i+k]
            u, v = kmer[:-1], kmer[1:]
            graph[u].append(v)
            outdeg[u] += 1
            indeg[v] += 1
    return graph, indeg, outdeg

def find_euler_path(graph, indeg, outdeg):
    """Camino euleriano mediante el algoritmo de Hierholzer."""
    # Determinar el nodo inicial
    start = None
    for node in list(outdeg.keys()) + list(indeg.keys()):
        if outdeg[node] - indeg[node] == 1:
            start = node
            break
    if start is None:
        start = next(iter(graph))  # Inicio arbitrario si existe un circuito euleriano

    # Copiar el grafo para no modificar el original
    g = {u: list(vs) for u, vs in graph.items()}
    stack = [start]
    path = []
    while stack:
        node = stack[-1]
        if g.get(node):
            stack.append(g[node].pop())
        else:
            path.append(stack.pop())
    return path[::-1]

def assemble(path):
    """Reconstruye la cadena genómica a partir del camino euleriano."""
    if not path:
        return ""
    result = path[0]
    for node in path[1:]:
        result += node[-1]
    return result

reads = ["ACGTAC", "CGTACG", "GTACGT", "TACGTA"]
k = 4
graph, indeg, outdeg = build_de_bruijn(reads, k)
path = find_euler_path(graph, indeg, outdeg)
genome = assemble(path)
print("Path:", path)
print("Genome:", genome)

Al ejecutarlo, se obtiene un camino euleriano en forma de ciclo y el genoma original se reconstruye. Las 40 líneas constituyen la base completa del ensamblaje con De Bruijn. Herramientas prácticas como SPAdes son una ingeniería masiva que añade sobre esta base el manejo múltiple de k, la eliminación de burbujas/puntas y el procesamiento de repeticiones, junto con información de pares extremos (paired-end).

Complejidad

• Construcción del grafo: O(número total de k-mers) = O(N × L). N es el número de lecturas.
• Camino euleriano: O(V + E). Lineal respecto a la suma del número de nodos y aristas.
• Total: Pocos minutos incluso con millones de lecturas.

Una ganancia decisiva frente al tiempo exponencial de OLC. No obstante, en la era de las lecturas largas (long-reads), las propias lecturas son más grandes por lo que se puede aumentar k; además, como se trató en M22, OLC resulta ser más estable.

SPAdes — El trono de la era de los short-reads

SPAdes (Bankevich et al. 2012) es el ensamblador De Bruijn emblemático. Es el estándar para el ensamblaje de genomas bacterianos y metagenomas. Características principales:

• k múltiple: Utiliza secuencialmente 21, 33, 55, 77 y 99 (con la opción --only-assembler).
• Modo híbrido: Combina De Bruijn para short-reads con long-reads (--pacbio, --nanopore) para rellenar huecos.
• Modo plásmido: Opción de ensamblaje específica para plásmidos.
• metaSPAdes: Modo especializado para microbiomas.

Práctica:

# Instalar SPAdes
wget https://github.com/ablab/spades/releases/download/v4.0.0/SPAdes-4.0.0-Linux.tar.gz
tar xzf SPAdes-4.0.0-Linux.tar.gz
./SPAdes-4.0.0-Linux/bin/spades.py --version

# Ensamblaje bacteriano pequeño (paired-end)
spades.py -1 reads_R1.fastq.gz -2 reads_R2.fastq.gz -o my_assembly \
          -k 21,33,55,77 -t 4

# Comprobar los resultados
grep -c "^>" my_assembly/scaffolds.fasta

Un genoma bacteriano (~5 Mb) se ensambla en decenas de minutos en un portátil. La escala del genoma humano requiere un clúster.

Mapeo de la informática — El poder de la redefinición de problemas

Este es un caso práctico de la transformación de problemas tratada en DryBench.

• Camino hamiltoniano (NP-completo) ← los nodos son lecturas
• Camino euleriano (tiempo polinómico) ← los nodos son k-mers y las aristas son lecturas

Al representar de otra manera el mismo problema físico —concatenar lecturas—, su complejidad computacional cambia de raíz. Una de las técnicas más poderosas del diseño de algoritmos es redefinir el problema. La clave es poder transformar una representación difícil en otra sencilla.

Esta lección aparece también en otros ámbitos. La FFT convierte la multiplicación de polinomios en una representación puntual y reduce O(n²) a O(n log n). El dual de la programación lineal traslada las partes difíciles del problema primal a una forma más sencilla. La representación es el algoritmo.

Obtener una puntuación en Rosalind

El problema DBG de Rosalind proporciona una lista de k-mers y pide construir el grafo de De Bruijn. Puede utilizarse directamente la función build_de_bruijn implementada en Python más arriba.

EULR de Rosalind: problema del camino euleriano. Puede utilizarse directamente la función find_euler_path.

Ramificaciones para el próximo artículo

• Próximo artículo (M24): hifiasm — ensamblaje diploide con lecturas PacBio HiFi, incluido el trio binning específico de los progenitores. La práctica moderna de OLC y culminación de la serie de ensamblaje del Micro Tier.
• Después de la serie de ensamblaje (M25): alineamiento múltiple de secuencias — alineamiento progresivo · MAFFT.
• Extensión: string graph (Myers) — teoría unificada de OLC y De Bruijn.

Si desea profundizar

• Pevzner, Tang, Waterman (2001), An Eulerian path approach to DNA fragment assembly, PNAS — primer artículo que aplicó De Bruijn al ensamblaje; la gran intuición de Pevzner.
• Bankevich et al. (2012), SPAdes: A New Genome Assembly Algorithm and Its Applications, J Comput Biol — artículo original de SPAdes y quintaesencia del enfoque multi-k.
• UC Berkeley CS176 — conferencia sobre ensamblaje del profesor Yun S. Song, con diagramas de la topología de De Bruijn desarrollados en la pizarra.
• Ben Langmead (JHU) — Lista de reproducción de YouTube sobre ensamblaje. Inducción manual de grafos.
• Compau & Pevzner Bioinformatics Algorithms Capítulo 3 — Tutorial de De Bruijn integrado con Rosalind.

Antes de pasar al siguiente capítulo, insertemos algunos errores en el código de Python anterior. Verificar manualmente cómo se forman las burbujas y las puntas, y cómo deben eliminarse, es la práctica real de este capítulo.