Extraer las K lecturas de mayor calidad: mantener solo 100 de un millón usando un montículo

Al finalizar este tema

Podrán crear una herramienta que seleccione en tiempo real las K lecturas de mayor calidad de entre grandes volúmenes de lecturas de secuenciación, combinando el montículo (cola de prioridad) y el ordenamiento aprendidos en los libros de texto. Entenderán mediante código por qué el montículo es mucho más rápido que el ordenamiento completo y por qué es ideal para datos en tiempo real.

Este artículo es un ejemplo educativo general. El filtrado real de lecturas implica procesos más complejos, como el recorte de calidad y la eliminación de adaptadores, pero aborda con precisión el concepto central de top-K.

"Las 1.000 mejores de entre un millón": la trampa del enfoque ingenuo

Imaginen que tienen 1 millón de lecturas de secuenciación y han calculado la puntuación media de calidad para cada una. Quieren conservar solo las 1.000 mejores.

Enfoque ingenuo A: Ordenamiento completo

def topk_by_sort(reads: list[dict], k: int) -> list[dict]:
    return sorted(reads, key=lambda r: -r["avg_quality"])[:k]

Este enfoque es correcto, pero tiene una complejidad temporal de O(n log n). Para un millón de elementos, se requieren aproximadamente 2 × 10⁷ comparaciones. En Python, esto tarda unos segundos. Y, lo más importante, todos los datos deben estar en la memoria.

Enfoque ingenuo B: actualizar el valor mínimo en cada iteración

def topk_by_scan(reads: list[dict], k: int) -> list[dict]:
    result = []
    for read in reads:
        result.append(read)
        if len(result) > k:
            worst = min(range(len(result)), key=lambda i: result[i]["avg_quality"])
            result.pop(worst)
    return result

Esto es O(n × k). 1 millón × 1.000 = 10⁹. Es mucho más lento.

El enfoque óptimo utiliza un montículo. Se mantiene un min-heap de tamaño K. Cuando llega un nuevo dato, se compara con el valor mínimo del montículo y se reemplaza si es mejor. La complejidad temporal es O(n log K). 1 millón × log(1.000) ≈ 10⁷. Es varias veces más rápido que la ordenación y permite la transmisión en tiempo real.

De la abstracción a los componentes

Componente 1: Min-heap

Un montículo es un árbol binario en el que el padre siempre es menor o igual que sus hijos (min-heap). Python heapq simula un min-heap mediante una lista.

import heapq

nums = []
heapq.heappush(nums, 5)
heapq.heappush(nums, 3)
heapq.heappush(nums, 8)
heapq.heappush(nums, 1)

print(nums)  # [1, 3, 8, 5] — tiene forma de lista, pero cumple la propiedad de heap

smallest = heapq.heappop(nums)  # 1

Propiedad clave: heappush y heappop tienen una complejidad de O(log n). La consulta del valor mínimo es O(1) (nums[0]).

Componente 2: Los K elementos principales

Ahora, se obtienen los K elementos principales utilizando un min-heap de tamaño K.

def topk_streaming(scores: list[float], k: int) -> list[float]:
    heap: list[float] = []

    for score in scores:
        if len(heap) < k:
            heapq.heappush(heap, score)
        else:
            if score > heap[0]:
                heapq.heapreplace(heap, score)

    return sorted(heap, reverse=True)

¿Por qué funciona? El valor mínimo del montón (heap[0]) es el K-ésimo valor entre los K valores más grandes. Si un nuevo valor es mayor que este, debe reemplazar al K-ésimo valor; esto es exactamente lo que hace heapreplace.

Análisis de la complejidad temporal: Cada una de las n iteraciones tiene una complejidad de O(log K). En total, O(n log K). Si k = 1000 y n = 1.000.000, se realizan aproximadamente 10⁷ operaciones. Es ligeramente más rápido que ordenar, O(n log n) ≈ 2 × 10⁷, pero solo utiliza una cantidad de memoria proporcional a K.

Parte 3: Mantener la identidad mediante tuplas

read no es un único escalar, sino un diccionario. Al almacenarlo en el montón, se guarda como una tupla (score, read).

def topk_reads(reads_iter, k: int) -> list[dict]:
    heap: list[tuple[float, int, dict]] = []
    counter = 0

    for read in reads_iter:
        score = read["avg_quality"]
        counter += 1
        entry = (score, counter, read)

        if len(heap) < k:
            heapq.heappush(heap, entry)
        else:
            if score > heap[0][0]:
                heapq.heapreplace(heap, entry)

    return [entry[2] for entry in sorted(heap, reverse=True)]

¿Por qué necesitamos un contador? Cuando la calidad de dos lecturas es la misma, el diccionario read se compara con otros, lo que genera una excepción. Dado que las tuplas se comparan secuencialmente a partir del primer elemento, se utiliza counter como segundo elemento para desempatar en caso de que las calidades sean iguales.

Procesamiento en flujo continuo: la verdadera fortaleza del heap

¿Qué sucede si un millón de lecturas no caben en la memoria al mismo tiempo? ¿O si es necesario leerlas línea por línea desde un archivo? El acceso al heap sigue funcionando de la misma manera, y este es el punto clave que lo diferencia del enfoque de ordenamiento.

def topk_from_fastq(fastq_path: str, k: int) -> list[dict]:
    heap: list[tuple[float, int, dict]] = []
    counter = 0

    with open(fastq_path) as f:
        while True:
            header = f.readline().strip()
            if not header:
                break
            seq = f.readline().strip()
            plus = f.readline().strip()
            qual = f.readline().strip()

            avg_quality = sum(ord(c) - 33 for c in qual) / len(qual)
            counter += 1
            read = {"header": header, "seq": seq, "qual": qual, "avg_quality": avg_quality}
            entry = (avg_quality, counter, read)

            if len(heap) < k:
                heapq.heappush(heap, entry)
            elif avg_quality > heap[0][0]:
                heapq.heapreplace(heap, entry)

    return [entry[2] for entry in sorted(heap, reverse=True)]

Esta función calcula los K elementos más grandes utilizando una memoria equivalente a K lecturas, independientemente de si el tamaño del archivo es de 100 GB o 1 TB. También funciona directamente en entornos de transmisión de datos en tiempo real. Esta es una aplicación práctica importante de la estructura de datos de montón (heap).

Implementación concisa con heapq.nlargest

La biblioteca estándar de Python ya proporciona esta funcionalidad.

from heapq import nlargest

top_reads = nlargest(1000, reads_iter, key=lambda r: r["avg_quality"])

Internamente, se utiliza exactamente el mismo algoritmo. En la práctica, se aplica este método. Sin embargo, es fundamental comprender por qué esta sintaxis es tan eficiente: esto permite implementar directamente la función top-K en otros lenguajes cuando sea necesario y personalizarla cuando K sea grande o existan condiciones especiales.

Atenuación: los dos espacios en blanco que deben completarse

Espacio en blanco 1: top-K condicional

Selecciona los valores principales (top-K) entre las lecturas cuya calidad supera un umbral específico y cuya longitud también es superior a un valor mínimo determinado.

def topk_conditional(
    reads_iter,
    k: int,
    min_quality: float = 20.0,
    min_length: int = 100
) -> list[dict]:
    """
    Obtiene el top-K solo entre las reads que cumplen las condiciones.
    Las que no las cumplen no se insertan en el heap.
    """
    heap: list = []
    counter = 0

    for read in reads_iter:
        # TODO: comprobar las condiciones de quality/length; si no se cumplen, usar continue
        # Si se cumplen, insertar en el heap con push o heapreplace
        pass

    return [entry[2] for entry in sorted(heap, reverse=True)]

Sugerencia: Primero, if read["avg_quality"] < min_quality or len(read["seq"]) < min_length: continue.

Espacio en blanco 2: Ordenación por múltiples criterios

Si la calidad es la misma, priorizar las lecturas más cortas (o al revés). Esto se incluye como el segundo elemento de la tupla.

def topk_by_multiple(
    reads_iter,
    k: int
) -> list[dict]:
    """
    Orden primario: avg_quality descendente
    Orden secundario: longitud ascendente (con la misma quality, se prefiere la más corta)
    """
    heap: list = []
    counter = 0

    for read in reads_iter:
        score = read["avg_quality"]
        length_key = -len(read["seq"])
        counter += 1

        # TODO: construir la tupla que se guardará en el heap
        # Prioridad: (score descendente, length ascendente)
        # Como es un min-heap, conservar score e invertir el signo de length
        pass

    return [entry[3] for entry in sorted(heap, reverse=True)]

Pista: entry = (score, length_key, counter, read). El orden en que se extraen elementos de un min-heap es score ascendente; en caso de empate, length_key ascendente (= length descendente). Por lo tanto, para obtener el top, se debe ordenar de forma inversa.

Reflexión — Diferencias con un filtro de lectura (read filter) real

Refinamiento del cálculo de calidad: Su avg_quality es un promedio simple. En la práctica, se analiza la calidad de cada posición individualmente; dado que la degradación de la calidad al principio o al final es común, se requieren filtros por posición independientes.

Recorte de adaptadores (Adapter trimming): El primer paso de un filtro real es la eliminación de secuencias de adaptadores de secuenciación. fastp y trimmomatic son estándares. El filtro Top-K es una etapa posterior.

Procesamiento de paired-end: En la práctica, los dos reads deben procesarse juntos como un par (mate). Si uno de ellos es descartado por el filtro, su pareja también debe gestionarse de la misma manera.

Archivos de mapeo de memoria: Los archivos FASTQ de tamaño masivo se acceden mediante mmap o se procesan mediante parsing por streaming en estado comprimido (gzip). Si cambia su open a gzip.open, se añadirá soporte para gzip.

Aceleración por GPU: El procesamiento de reads masivos (miles de millones) utiliza toolchains de GPU como NVIDIA Parabricks. Aunque el heap en sí no es amigable con la GPU, el cálculo de calidad se puede paralelizar.

Proyecto de expansión

1. Soporte para paired-end: Filtro Top-K manteniendo los pares mientras se realiza el streaming de dos archivos FASTQ simultáneamente.

2. Soporte automático para gzip: Verificar la extensión del archivo y, si es .gz, usar gzip.open automáticamente.

3. Dashboard de distribución de calidad: Comparar la distribución de calidad antes y después del filtrado utilizando matplotlib.

4. Adición de bottom-K: Extraer reads de baja calidad y guardarlos por separado (para diagnóstico de problemas).

Mapa de componentes de esta unidad

• [F] Heap · Cola de prioridad: Mantiene el Top-K mediante un min-heap de tamaño K. Push/replace con O(log K).
• [F] Comparación con ordenamiento: sort O(n log n) vs heap O(n log K). Cuándo es ventajoso cada uno.
• [W] E/S de archivos: Parsing de FASTQ (se proporciona como script completo).

[F] = Implementación propia / [W] = Proporcionado como código completo.