¿Por qué es necesario este episodio?

Ya hemos aprendido BLAST; ahora busquemos en todo el proteoma. Supongamos que tenemos 1 millón de secuencias de proteínas predichas a partir del ensamblaje de un metagenoma microbiano y una base de datos NCBI nr con 500 millones de secuencias. Ejecutar BLASTP directamente tomaría varios días. En la práctica, estos procesos por lotes deben ejecutarse una vez al día. Por eso, en 2015, surgió Diamond, una herramienta que comenzó a producir los mismos resultados que BLAST de 100 a 10.000 veces más rápido.

En este episodio, analizaremos tres ideas de rediseño que permitieron a Diamond superar las limitaciones de BLAST a nivel de principios. Comprendamos cómo reducir únicamente el factor constante, manteniendo la precisión del algoritmo, puede ser un cambio radical en la práctica profesional.

La barrera que enfrenta BLASTP

Repasemos brevemente la estructura fundamental de BLAST que aprendimos en M10.

1. La secuencia de consulta se divide en semillas k-mer (para proteínas, normalmente k=3).
2. Se buscan las posiciones donde coinciden estas semillas en el índice de la base de datos.
3. A partir de las posiciones coincidentes, se expande en ambas direcciones para crear HSP (High-scoring Segment Pair).
4. Si la puntuación de expansión supera un umbral, se mantiene como candidato.

Este flujo es elegante en sí mismo. El problema radica en el fallo de caché del paso 2. A medida que la base de datos crece, el número de posiciones donde coincide una sola semilla aumenta exponencialmente, y cada vez se leen secuencias de la base de datos mediante accesos aleatorios. Más de la mitad del pipeline de la CPU se desperdicia en estado de espera.

En la era de los metagenomas, esta barrera es aún más crítica. Las consultas son grandes (más de 1 millón de proteínas predichas) y la base de datos también lo es (nr con más de 500 millones). El resultado es que un solo proceso por lotes tarda días. La primera herramienta en reducir esta brecha fue Diamond. Benjamin Buchfink et al. lo publicaron inicialmente en Nature Methods en 2015, y en 2021 lograron una aceleración de hasta 10.000 veces con la v2.

Tres ideas de rediseño de Diamond

Idea 1 — Doble indexación

BLAST solo indexa la consulta. La base de datos se escanea cada vez. Diamond indexa tanto la consulta como la base de datos, extrayendo las coincidencias de semillas a partir de la intersección de ambos índices.

¿Por qué esto es más rápido? El índice de la base de datos solo necesita crearse una vez. Además, el índice de la consulta se crea por adelantado para cada lote. La coincidencia de semillas se completa con un solo merge join (una estructura similar al sort-merge join de SQL) entre dos listas ordenadas. La localidad de caché mejora drásticamente y los accesos aleatorios se convierten en accesos secuenciales.

Veámoslo brevemente en una tabla. Veamos solo algunas semillas de ejemplo con k=3, para 3 secuencias de consulta y 4 secuencias de la base de datos.

El método BLAST recorre toda la base de datos D para cada semilla de consulta. El método Diamond ordena ambas listas por semilla para extraer las coincidencias mediante una única operación de combinación.

Lista ordenada de semillas de consulta:       (AGT, Q1) (AGC, Q3) (CAG, Q3) (GTC, Q1) (GTC, Q2) (TCA, Q2)
Semillas ordenadas de la base de datos:       (AGT, D1) (AGC, D3) (CAG, D2) (GTA, D1) (GCT, D3) (TCA, D2)

Emparejar semillas idénticas (AGT: Q1-D1, AGC: Q3-D3, CAG: Q3-D2, TCA: Q2-D2). Dado que ambas listas ya están ordenadas, basta con desplazar dos punteros en paralelo y recorrerlas una sola vez.

Idea 2 — Alfabeto reducido y semilla espaciada (spaced seed)

El alfabeto de proteínas tiene 20 caracteres. Incluso con k=3, existen 20³ = 8.000 combinaciones. La especificidad de la semilla es baja, lo que genera muchos falsos positivos. Diamond utiliza un alfabeto reducido agrupado por similitud química (de unos 11 caracteres). Aunque la semilla sigue teniendo un ancho de 4 caracteres, emplea una semilla espaciada que no requiere coincidencia en las dos posiciones centrales. El patrón de la semilla es **·**.

El alfabeto reducido reduce el número de coincidencias de semillas sin perder la homología remota. La semilla espaciada disminuye la tasa de falsos positivos aleatorios y mantiene las señales de alineamiento reales. La combinación de estos dos trucos reduce drásticamente el número de semillas candidatas en comparación con BLAST, aligerando la carga de la etapa de extensión posterior.

Idea 3 — Disposición por bloques y vectorización SIMD

BLAST realiza la extensión inmediatamente al encontrar un hit. Diamond divide la base de datos en bloques de tamaño fijo (2 GB por defecto) para su procesamiento. Dentro de cada bloque, recopila todas las coincidencias de semillas y luego procesa la extensión en lotes (batch).

¿Por qué es beneficioso el procesamiento por lotes? La extensión Smith-Waterman solo requiere tres celdas vecinas (arriba, izquierda y diagonal) para calcular una celda de la matriz. Este flujo de datos se adapta perfectamente a las instrucciones SIMD (SSE·AVX2). Al procesar la extensión de múltiples pares de secuencias en paralelo, se pueden completar entre 16 y 32 celdas simultáneamente con una sola instrucción. Esto multiplica teóricamente el rendimiento de un solo núcleo de CPU.

A esto se suma la paralelización por hilos, lo que permite una aceleración de hasta cuatro órdenes de magnitud en comparación con BLAST en un servidor de 8 núcleos. En la v2, la brecha se amplió aún más con la adición de la extensión mediante GPU.

Compromisos (trade-offs) del modo de sensibilidad

El interruptor más importante al usar Diamond en entornos prácticos es el modo de sensibilidad.

Si desea reproducir exactamente los resultados de BLAST, utilice --ultra-sensitive. Aun así, sigue siendo de 2 a 5 veces más rápido. Por el contrario, si solo necesita un cribado de candidatos, la configuración predeterminada ofrece una aceleración de dos órdenes de magnitud de forma gratuita.

En lugar de memorizar esta tabla, recordemos el principio. La velocidad y la sensibilidad se ajustan mediante la especificidad de la semilla. Si la semilla es densa, hay menos candidatos y la velocidad es mayor, pero se pierden las homologías distantes. Si la semilla es dispersa, ocurre lo contrario. En la práctica, la elección de este modo influye más en la calidad del resultado que la selección de la herramienta.

Práctica: ejecute su primera búsqueda con Diamond CLI

Funciona directamente en instancias de CPU gratuitas de Colab o SageMaker Studio Lab.

# 1) Instalación (también funciona sin conda)
wget https://github.com/bbuchfink/diamond/releases/download/v2.1.9/diamond-linux64.tar.gz
tar xzf diamond-linux64.tar.gz
./diamond version

# 2) Crear una base de datos pequeña (por ejemplo, 5000 secuencias de SwissProt)
wget -O sprot.fasta.gz https://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/complete/uniprot_sprot.fasta.gz
gunzip sprot.fasta.gz
./diamond makedb --in sprot.fasta -d sprot

# 3) Preparar la consulta (por ejemplo, insulina humana P01308)
cat > insulin.fasta <<'EOF'
>P01308 Human Insulin
MALWMRLLPLLALLALWGPDPAAAFVNQHLCGSHLVEALYLVCGERGFFYTPKTRREAED
LQVGQVELGGGPGAGSLQPLALEGSLQKRGIVEQCCTSICSLYQLENYCN
EOF

# 4) Ejecutar la búsqueda — modo default
./diamond blastp -q insulin.fasta -d sprot -o hits.tsv \
  --threads 4 --outfmt 6 qseqid sseqid pident length evalue bitscore

# 5) Comparar modos de sensibilidad — ultra-sensitive
./diamond blastp -q insulin.fasta -d sprot -o hits_us.tsv \
  --ultra-sensitive --threads 4 --outfmt 6 qseqid sseqid pident length evalue bitscore

# 6) Comparar el número de líneas de resultados
wc -l hits.tsv hits_us.tsv

Las proteínas altamente conservadas, como la insulina, generan decenas de coincidencias entre especies incluso en el modo predeterminado. Al cambiar a --ultra-sensitive, se observa que la detección de homólogos distantes aumenta varias veces.

Formato de los resultados

La salida predeterminada de Diamond hereda directamente el formato outfmt 6 de BLAST. Se puede especificar explícitamente el orden de las columnas, lo cual es esencial al escribir scripts de posprocesamiento.

# Filtro sencillo con Python
import csv

MIN_PIDENT = 30.0
MIN_LENGTH = 50
MAX_EVALUE = 1e-10

with open("hits.tsv") as f:
    reader = csv.reader(f, delimiter="\t")
    for row in reader:
        qseqid, sseqid, pident, length, evalue, bitscore = row
        if (float(pident) >= MIN_PIDENT and
            int(length) >= MIN_LENGTH and
            float(evalue) <= MAX_EVALUE):
            print(qseqid, sseqid, pident, evalue)

En la práctica, en este punto se vuelve a insertar pandas, se realiza un join de la taxonomía con los resultados, se vuelve a realizar un join con los nombres de los genes y, finalmente, se organiza todo en una sola tabla. Lo importante aquí es el hábito de establecer y documentar de antemano el umbral de filtrado. Si no se definen conjuntamente el modo de sensibilidad y el umbral de postprocesamiento, la reproducibilidad no será posible.

Complejidad y percepción práctica

La complejidad temporal teórica es O((m + n) × número de candidatos de semilla), similar a BLAST, pero los factores constantes difieren en dos o tres órdenes de magnitud. El objetivo de Diamond es precisamente el rediseño de los factores constantes. Aunque la clase asintótica del algoritmo permanece igual, cuando el tiempo de ejecución real se optimiza según las características del hardware, la experiencia práctica cambia por completo.

Esto ilustra una lección común en CS: incluso con la misma complejidad asintótica, rediseñar los factores constantes cambia las reglas del juego. La localidad de datos, SIMD, procesamiento por lotes y paralelismo de hilos —al combinarse estas cuatro herramientas— han permitido que procesos que en la era de BLAST requerían días en un clúster de servidores, ahora se completen en pocas horas en una computadora portátil.

Mapeo en CS: la unión de hashing y procesamiento por lotes

Aquí se combinan el indexado basado en hash y el join por lotes tratados en DryBench.

• Indexado por hash: Mapeo de semilla → posición. Reducción del espacio de hash mediante un alfabeto abreviado.
• Merge join: Intersección de dos listas de semillas ordenadas (isomorfo al sort-merge join de SQL).
• Localidad de caché: Rediseño de acceso aleatorio a acceso secuencial.
• Vectorización SIMD: Cálculo paralelo de las celdas de la rejilla en la fase de expansión.
• Procesamiento por lotes: Agrupación de semillas por bloques para procesar las expansiones de forma masiva.

El hecho de que el diseñador de algoritmos sea consciente o no de que la arquitectura de la CPU ha evolucionado sobre tres ejes —pipeline, caché e instrucciones vectoriales— determina si herramientas como esta pueden existir o no.

Ramificaciones hacia los próximos capítulos

• Siguiente capítulo (M18): Fundamentos de HMM + algoritmo de Viterbi — La versión probabilística del alineamiento de secuencias. Trata el estado de cada posición mediante probabilidades.
• Dos capítulos después (M19): HMM Forward-Backward — Cálculo hacia adelante y hacia atrás de las probabilidades de estado.
• Extensión práctica: Clasificación taxonómica de metagenomas con la opción --taxon-classify y búsqueda de repeticiones similares a PSI-BLAST con la opción --iterate.

El legado de BLAST no ha desaparecido. Usar Diamond no es más que superponer un rediseño para la era del hardware sobre el esqueleto algorítmico establecido por BLAST. En el próximo capítulo, cambiaremos a un modelo probabilístico con una estructura fundamentalmente distinta.

Para profundizar más

El cuerpo del texto es una narrativa reconstruida por BPD. Si ya comprendes los principios, puedes profundizar mediante los artículos originales y las clases académicas correspondientes.

• Buchfink, Xie, Huson (2015), Fast and sensitive protein alignment using DIAMOND, Nature Methods — Artículo original. Conceptos de doble indexación y alfabeto abreviado.
• Buchfink, Reuter, Drost (2021), Sensitive protein alignments at tree-of-life scale using DIAMOND, Nature Methods — Artículo v2. Cinco modos de sensibilidad y detalles de SIMD.
• Documentación oficial de Diamond: https://github.com/bbuchfink/diamond/wiki — Todas las opciones de CLI · Benchmarks de los modos de sensibilidad.
• EMBL-EBI Training — Lista de reproducción de tutoriales prácticos. Trata Diamond en el contexto de pipelines metagenómicos.

Si tienes tiempo, compara directamente el número de hits y el tiempo de ejecución por cada modo de sensibilidad utilizando una base de datos de tamaño intermedio como UniRef50. El objetivo de este ejercicio es experimentar cómo y cuánto Diamond alivia los cuellos de botella en el pipeline.