Ensamblaje de secuencias: unir fragmentos cortos mediante un grafo k-mer

Al finalizar este tema

Podrán construir una herramienta de ensamblaje que una cientos de lecturas de secuenciación cortas en un único contig largo, combinando la lista de adyacencia de grafos y los recorridos DFS/BFS aprendidos en los libros de texto. Comprenderán mediante código el principio interno del ensamblaje —el grafo de De Bruijn— que realizan ensambladores reales como SPAdes o Velvet, los cuales parecen operar por arte de magia.

Este texto es un ejemplo educativo general. El ensamblaje de secuencias real implica flujos de trabajo mucho más complejos, pero aborda con precisión los conceptos algorítmicos fundamentales en su núcleo.

"¿Por qué un genoma de 3 Gb no se obtiene de lecturas de 300 pb?" — La trampa del ensamblaje

Imaginemos que secuencian el genoma humano. Los secuenciadores modernos (Illumina) generan cientos de millones de lecturas cortas de aproximadamente 300 pb de una sola vez. Deben volver a unir estas piezas en una única secuencia continua para reconstruir el genoma humano de 3.000 millones de pb.

El enfoque más ingenuo consiste en encontrar las regiones superpuestas entre las lecturas y unirlas progresivamente.

def naive_assemble(reads: list[str]) -> str:
    result = reads[0]
    for read in reads[1:]:
        overlap = find_overlap(result, read)
        result += read[overlap:]
    return result

Este enfoque presenta dos problemas críticos.

Problema 1: El orden de las lecturas no tiene sentido. Las lecturas obtenidas del secuenciador provienen de posiciones aleatorias en el genoma, por lo que es muy poco probable que reads[0] y reads[1] sean adyacentes en el genoma real.

Problema 2: El cálculo de las superposiciones se vuelve explosivo. Para determinar qué pares de lecturas, entre 100 millones, se superponen realmente, en el peor de los casos, se necesitarían (100 millones)² = 10 mil billones de comparaciones. Esto no terminaría antes del fin del universo.

El enfoque real utiliza un grafo. Las lecturas se dividen en fragmentos cortos (k-mers), creando un grafo de De Bruijn donde cada k-mer es un nodo y las relaciones de conexión son las aristas. Luego, al encontrar un camino en este grafo, se obtiene la secuencia ensamblada.

De la caja negra a los componentes: explorando el grafo de De Bruijn

Para comprender este marco de trabajo, se requieren tres componentes.

Componente 1: k-mer — Fragmentación de cadenas

Un k-mer es una secuencia corta de longitud k. Si fragmentamos ATGCAT con k=3 (3-mer), obtenemos cuatro fragmentos: ATG, TGC, GCA y CAT.

def get_kmers(sequence: str, k: int) -> list[str]:
    return [sequence[i:i+k] for i in range(len(sequence) - k + 1)]

La característica principal es que los k-mers consecutivos comparten un prefijo/sufijo de (k-1)-mer. El sufijo de TG ATG = el prefijo de TGC. Esta coincidencia constituye la arista del grafo.

Componente 2: Lista de adyacencia del grafo

Los nodos del grafo de De Bruijn son (k-1)-mers. Las aristas son k-mers. La existencia de una arista de un nodo a otro significa que esos dos (k-1)-mers están conectados mediante una relación de prefijo y sufijo de un k-mer.

Este grafo se almacena como una lista de adyacencia (diccionario).

from collections import defaultdict

def build_de_bruijn(reads: list[str], k: int) -> dict[str, list[str]]:
    graph = defaultdict(list)
    for read in reads:
        for kmer in get_kmers(read, k):
            prefix = kmer[:-1]
            suffix = kmer[1:]
            graph[prefix].append(suffix)
    return graph

Al consultar graph["ATG"], se obtiene la lista de los nodos de destino de las aristas que salen del nodo ATG.

Componente 3: Recorrido DFS/BFS

Una vez construido el grafo, se debe encontrar un camino que recorra cada arista exactamente una vez (camino euleriano). Al concatenar este camino, se reconstruye la secuencia original.

Utilicemos DFS. Desde el nodo inicial, visitamos recursivamente los nodos vecinos, pero eliminamos la arista del grafo cada vez que se consume.

def find_eulerian_path(graph: dict[str, list[str]], start: str) -> list[str]:
    stack = [start]
    path = []
    graph = {k: list(v) for k, v in graph.items()}

    while stack:
        node = stack[-1]
        if graph.get(node):
            next_node = graph[node].pop()
            stack.append(next_node)
        else:
            path.append(stack.pop())

    return path[::-1]

Este es el algoritmo de Hierholzer. Su complejidad temporal es O(número de aristas).

Combinar las tres partes: pipeline de ensamblaje

Ahora, unamos las tres partes en un único pipeline.

from collections import defaultdict


def kmerize(sequence: str, k: int) -> list[str]:
    return [sequence[i:i+k] for i in range(len(sequence) - k + 1)]


def build_de_bruijn_graph(reads: list[str], k: int) -> dict[str, list[str]]:
    graph = defaultdict(list)
    for read in reads:
        for kmer in kmerize(read, k):
            prefix, suffix = kmer[:-1], kmer[1:]
            graph[prefix].append(suffix)
    return dict(graph)


def find_start_node(graph: dict[str, list[str]]) -> str:
    out_degree = {node: len(edges) for node, edges in graph.items()}
    in_degree: dict[str, int] = defaultdict(int)
    for edges in graph.values():
        for target in edges:
            in_degree[target] += 1

    for node in graph:
        if out_degree[node] > in_degree[node]:
            return node
    return next(iter(graph))


def eulerian_path(graph: dict[str, list[str]], start: str) -> list[str]:
    graph = {k: list(v) for k, v in graph.items()}
    stack, path = [start], []
    while stack:
        node = stack[-1]
        if node in graph and graph[node]:
            stack.append(graph[node].pop())
        else:
            path.append(stack.pop())
    return path[::-1]


def path_to_sequence(path: list[str]) -> str:
    if not path:
        return ""
    return path[0] + "".join(node[-1] for node in path[1:])


def assemble(reads: list[str], k: int = 5) -> str:
    graph = build_de_bruijn_graph(reads, k)
    start = find_start_node(graph)
    path = eulerian_path(graph, start)
    return path_to_sequence(path)

Veamos esto con un ejemplo sencillo.

original = "ATGCATGCATGATG"
reads = [original[i:i+8] for i in range(0, len(original) - 7)]
# ['ATGCATGC', 'TGCATGCA', 'GCATGCAT', 'CATGCATG', 'ATGCATGA', 'TGCATGAT', 'GCATGATG']

result = assemble(reads, k=5)
print(result)  # ATGCATGCATGATG (o una reconstrucción similar)

Fading: los tres espacios en blanco que debes completar

Ahora es tu turno. El código anterior contiene tres espacios en blanco que deben abordarse en la práctica. Solo te daré una pista para cada uno.

Espacio en blanco 1: Filtrado de errores de secuenciación

Las lecturas reales contienen caracteres erróneos mezclados. Los k-mers con errores crean nodos incorrectos de frecuencia extremadamente baja en el gráfico. Estos nodos contaminan el ensamblaje.

def filter_low_frequency_kmers(reads: list[str], k: int, threshold: int = 2) -> list[str]:
    from collections import Counter
    kmer_counts: Counter[str] = Counter()
    for read in reads:
        for kmer in kmerize(read, k):
            kmer_counts[kmer] += 1

    # TODO: filtrar lecturas que contengan k-mers con frecuencia por debajo del umbral
    filtered = []
    for read in reads:
        # Tu turno: verificar si la frecuencia de todos los k-mers en esta lectura es mayor o igual al umbral
        pass
    return filtered

Pista: all(kmer_counts[kmer] >= threshold for kmer in kmerize(read, k))

Espacio en blanco 2: Procesamiento de múltiples puntos de inicio

En los datos de secuenciación reales, es posible que el genoma no se conecte en un único camino euleriano. En la mayoría de los casos, se deben generar múltiples contigs.

def assemble_multi_contigs(reads: list[str], k: int) -> list[str]:
    graph = build_de_bruijn_graph(reads, k)
    contigs = []

    while graph:
        # TODO: seleccionar un nodo inicial en el grafo restante y ensamblar un contig
        # Eliminar las aristas utilizadas después del ensamblaje
        # Repetir mientras el grafo restante no esté vacío
        pass

    return contigs

Pista: Llame a find_start_node repetidamente y, en cada iteración, ensamble uno con eulerian_path. Dado que las aristas utilizadas se eliminan (pop) tras cada ensamblaje, el grafo se reduce de forma natural.

Espacio vacío 3: Manejo de secuencias repetitivas

Los genomas contienen frecuentemente secuencias repetitivas. Cuando una misma secuencia corta aparece en múltiples lugares del genoma, el nodo correspondiente a esa secuencia en el grafo de De Bruijn se convierte en un punto de intersección de múltiples aristas, lo que genera ambigüedad en el ensamblaje.

La resolución completa de este problema requiere que los ensambladores reales utilicen diversas heurísticas sofisticadas (lecturas emparejadas, lecturas largas, etc.). El enfoque que deben intentar:

def detect_repeat_nodes(graph: dict[str, list[str]]) -> set[str]:
    """Los nodos con aristas que salen en múltiples direcciones son candidatos a secuencias repetidas"""
    # TODO: devolver los nodos con out-degree > 1 o in-degree > 1
    pass

Mostrar solo la ubicación de estos nodos repetidos permite emitir una advertencia al usuario: "El ensamblaje es ambiguo aquí".

Reflexión — ¿En qué se diferencia este código de un ensamblador de producción?

Aunque el ensamblador que han creado comparte raíces conceptuales con SPAdes o Velvet, los ensambladores de producción reales son mucho más sofisticados. A continuación, se resumen las principales diferencias.

Escalabilidad: Los ensambladores reales procesan cientos de millones de lecturas. Su implementación en Python no puede manejar esta escala debido a la sobrecarga de los diccionarios. Los ensambladores reales están escritos en C++ y comprimen la memoria codificando los k-meros como enteros (2 bits por base).

Lecturas emparejadas (paired-end reads): Los secuenciadores reales generan dos lecturas desde los extremos de un fragmento y conocen aproximadamente la distancia entre ellas. Esta información sobre el intervalo es crucial para resolver problemas de secuencias repetidas. Su implementación aún no incluye esto.

Lecturas largas (PacBio/Nanopore): Los secuenciadores recientes generan lecturas mucho más largas (10 kb+). En este caso, el enfoque cambia completamente (OLC — Overlap-Layout-Consensus), utilizando grafos de solapamiento de cadenas en lugar de De Bruijn.

Puntuaciones de calidad: Las lecturas reales incluyen una puntuación de calidad para cada base. Esta información permite realizar un filtrado de errores mucho más sofisticado.

Proyectos de ampliación

Es posible ampliar este ensamblador en varias direcciones.

1. Entrada de archivos FASTQ: Los datos de secuenciación reales llegan en formato FASTQ. Añada una función para analizar este archivo y leer las lecturas junto con sus puntuaciones de calidad.

2. Informe estadístico: Calcule y muestre el N50 (indicador de la distribución de longitudes de los contigs), la longitud total del ensamblaje, el número de contigs, etc., de los resultados del ensamblaje.

3. Visualización: Dibuje un pequeño grafo de De Bruijn con networkx + matplotlib. Verifique visualmente cómo los nodos repetidos complican el ensamblaje.

4. Prueba con datos reales: Descargue datos de secuenciación públicos de un genoma bacteriano (el de E. coli es de aproximadamente 4.6M bp) y ejecute su ensamblador, comparando los resultados con SPAdes.

Mapa de componentes de este capítulo

Se resume aquí qué conceptos individuales de informática se combinaron en esta aplicación práctica. Para una explicación detallada de cada concepto, consulte DryBench.

• [F] Lista de adyacencia de grafos: Representación mediante diccionario graph[node] = [neighbor1, neighbor2, ...]. Almacenamiento del grafo de De Bruijn.
• [F] DFS/BFS: Camino euleriano mediante el DFS basado en pilas de Hierholzer. Recorrido real del ensamblaje.
• [W] Tabla hash / diccionario: defaultdict(list) para búsqueda de nodos en O(1).
• [W] Segmentación k-mer: Slicing de cadenas y comprensión de listas.

[F] = Implementación propia / [W] = Concepto de herramienta proporcionado con el código completo.