Búsqueda de secuencias de proteínas con FAISS: Capturando homólogos lejanos que BLAST pasa por alto mediante modelos fundacionales

BLAST ha sido el estándar para la búsqueda de homología de proteínas durante 30 años, pero su rendimiento se desploma drásticamente cuando la similitud de secuencia es inferior al 30 %. Sin embargo, existen innumerables homólogos lejanos (remote homologs) en la naturaleza que están evolutivamente relacionados pero tienen secuencias muy diferentes; pasar por alto estos casos impide la estimación de la función de nuevos genes. En este artículo, construiremos un pipeline práctico que resuelve este problema utilizando el espacio latente aprendido por el modelo fundacional de proteínas (ESM3) y añade una búsqueda a nivel de servicio en tiempo real sobre bases de datos masivas de secuencias mediante FAISS.

📚 Recomendación de lecturas previas (Altamente recomendada)

Este artículo es una guía avanzada y técnica de IA × Biología. Antes de comenzar, recomendamos encarecidamente revisar primero los siguientes artículos de DryBench:

• DryBench ai-native #3 Transformers y embeddings
• DryBench ai-native #5 Mecanismos de atención
• DryBench ai-native #13 HuggingFace y APIs comerciales

Si se adentra en este artículo sin haber leído las lecturas previas, el código práctico comenzará inmediatamente sin volver a explicar el significado de los vectores de embedding, la concentración de información de la atención o el principio de carga de modelos en HuggingFace, por lo que resultará difícil de seguir.

Lo que ya aprendimos en DryBench

En DryBench ai-native #3 aprendimos que los transformers aprenden las relaciones entre tokens y representan cada token como un vector denso (dense embedding) que refleja el contexto; en #5 aprendimos que la atención es un mecanismo que captura dependencias de largo alcance al aprender qué tan importantes son dos posiciones arbitrarias dentro de una secuencia entre sí; y en #13 aprendimos que HuggingFace ha estandarizado los pesos, tokenizadores y APIs de inferencia de estos modelos, permitiendo cargarlos con unas pocas líneas de código.

Pero, ¿qué sucede cuando aplicamos estos tres conceptos al dominio de las proteínas? A diferencia del lenguaje natural, las "oraciones" proteicas han sido moldeadas por presiones evolutivas durante 3 mil millones de años, y la información contextual (aminoácidos vecinos) determina la estructura de plegamiento tridimensional. La serie ESM de Meta AI (→ spin-off de EvolutionaryScale) es el primer modelo fundacional de proteínas a gran escala que aprendió esta información evolutiva mediante transformers, y múltiples estudios han confirmado que los embeddings aprendidos contienen potencialmente información sobre la estructura 3D, la función y la estabilidad de las proteínas [1][2][5]. Este artículo presenta el pipeline para convertir esos embeddings en una herramienta de búsqueda realmente útil.

Definición del problema complejo

Limitaciones fundamentales de BLAST

BLAST (Basic Local Alignment Search Tool, 1990) asigna puntuaciones a la similitud de secuencias mediante una matriz de sustitución y penalizaciones por huecos (gap penalty) para devolver coincidencias (hits) estadísticamente significativas. Este enfoque es extremadamente preciso para homólogos cercanos (close homologs) con una identidad de secuencia superior al 30–40 %, pero falla en las siguientes dos situaciones:

• Homólogos remotos (remote homologs): Cuando existe una relación evolutiva, pero la identidad de secuencia cae al 20–30 %. Aunque la estructura 3D y la función siguen conservadas, BLAST pierde significancia estadística. El rango del 20–35 %, conocido como "zona crepuscular" (twilight zone), es particularmente difícil, mientras que la "zona de medianoche" (midnight zone, < 20 %) es prácticamente indetectable.
• Emparejamiento de dominios específicos en proteínas multidominio: BLAST busca alineamientos locales de la secuencia completa; si la puntuación del alineamiento global es baja, puede pasar por alto dominios individuales que presentan un fuerte homólogo.

Las soluciones tradicionales para mitigar esto incluyen PSI-BLAST (perfil iterativo), HHblits (HMM-HMM) y Foldseek (basado en estructura), pero cada una presenta problemas de coste computacional y dependencia de datos.

Lo que resuelve el enfoque de embeddings

Los modelos fundacionales como ESM3 absorben los acoplamientos evolutivos (evolutionary coupling) en el espacio latente durante el proceso de entrenamiento. Como resultado, las embeddings finales se encuentran cerca en el espacio vectorial si la estructura y la función son similares, incluso si la similitud de secuencia es baja. Aprovechando esta propiedad:

• Es posible recuperar homólogos remotos que BLAST pasó por alto mediante búsqueda KNN con embeddings (según evidencia bibliográfica [2][3]).
• Al calcular una sola embedding por secuencia, todas las búsquedas posteriores solo requieren cálculos de similitud vectorial → ventaja en grandes volúmenes de datos frente al alineamiento por pares (pairwise alignment) de BLAST.
• La transferencia de función (por ejemplo, términos GO) puede predecirse inmediatamente mediante la votación mayoritaria de las etiquetas de los $k$ vecinos más cercanos encontrados.

Métricas objetivo de este proyecto

• Construcción de un índice de embeddings para un subconjunto de 500.000 secuencias de UniProt Swiss-Prot.
• Tasa de detección de homólogos remotos en 1000 consultas: Mejora del recall en +25 puntos porcentuales frente a BLAST (basado en el benchmark CAFA [4]).
• Latencia media de búsqueda: Menos de 100 ms por consulta (mediante ajuste de FAISS HNSW).
• F1 de transferencia de función GO: Superior a 0,65 (nivel de vanguardia en CAFA-3).

Stack de herramientas y requisitos de infraestructura

Requisitos de infraestructura:

• Fase de extracción de embeddings: GPU de consumo pequeña (RTX 4060 con 8 GB o más). ESM3 1.4B consume aprox. 3 GB de VRAM en fp16; se recomienda un batch de 4 a 8.
• Construcción y búsqueda del índice FAISS: Solo CPU es suficiente. RAM de 16 GB o más (500.000 secuencias × 1024 dimensiones float32 ≈ 2 GB).
• Disco: UniProt Swiss-Prot comprimido aprox. 200 MB, caché de embeddings aprox. 2 a 4 GB.

Coste estimado de reproducción para el estudiante: Coste de API: 0 (ejecución local). Tiempo de GPU: aprox. 2 a 4 horas (estimación para el embedding de 500.000 secuencias, basado en RTX 4060). Descarga de datos: aprox. 200 MB.

Implementación práctica del pipeline

Flujo completo:

Paso 1. Descarga y análisis de UniProt Swiss-Prot

Swiss-Prot es la sección de UniProt curada manualmente, en la que las anotaciones GO se incluyen como etiquetas fiables.

import gzip
from pathlib import Path
from typing import Iterator, NamedTuple

import requests
from Bio import SeqIO


SWISSPROT_URL = "https://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/complete/uniprot_sprot.fasta.gz"


class ProteinRecord(NamedTuple):
    """Registro de secuencia de proteínas."""
    accession: str
    sequence: str
    description: str
    go_terms: list[str]


def download_swissprot(dest: Path) -> Path:
    """Descarga de FASTA UniProt Swiss-Prot."""
    dest.parent.mkdir(parents=True, exist_ok=True)
    if dest.exists():
        return dest
    with requests.get(SWISSPROT_URL, stream=True, timeout=300) as resp:
        resp.raise_for_status()
        with open(dest, "wb") as f:
            for chunk in resp.iter_content(chunk_size=8192):
                f.write(chunk)
    return dest


def parse_swissprot(fasta_gz: Path, max_len: int = 1024) -> Iterator[ProteinRecord]:
    """Análisis en streaming de FASTA. Se omiten las secuencias con longitud >= max_len (límite de contexto de ESM)."""
    with gzip.open(fasta_gz, "rt") as f:
        for record in SeqIO.parse(f, "fasta"):
            seq = str(record.seq).upper()
            if len(seq) > max_len or len(seq) < 30:
                continue
            # Extraer el accession de la descripción FASTA de UniProt (ej: sp|P12345|GENE_ORG)
            parts = record.id.split("|")
            accession = parts[1] if len(parts) >= 2 else record.id
            yield ProteinRecord(
                accession=accession,
                sequence=seq,
                description=record.description,
                go_terms=[],  # GO se fusiona desde una fuente separada, ver paso a continuación
            )

La anotación GO se fusiona a partir de un archivo independiente (goa_uniprot_all.gaf.gz [6]). En la práctica, es más sencillo realizar la unión mediante SQLite o DuckDB.

import sqlite3
from collections import defaultdict


def build_go_index(gaf_path: Path, db_path: Path) -> None:
    """Indexar GAF (GO Annotation File) como mapeo accession → [IDs de GO]."""
    conn = sqlite3.connect(db_path)
    conn.execute("CREATE TABLE IF NOT EXISTS go (accession TEXT, go_id TEXT)")
    conn.execute("CREATE INDEX IF NOT EXISTS idx_acc ON go(accession)")
    with gzip.open(gaf_path, "rt") as f:
        batch = []
        for line in f:
            if line.startswith("!"):
                continue
            fields = line.split("\t")
            if len(fields) < 5:
                continue
            batch.append((fields[1], fields[4]))  # DB_Object_ID, GO_ID
            if len(batch) >= 10000:
                conn.executemany("INSERT INTO go VALUES (?, ?)", batch)
                batch.clear()
        if batch:
            conn.executemany("INSERT INTO go VALUES (?, ?)", batch)
    conn.commit()
    conn.close()


def lookup_go_terms(accession: str, db_path: Path) -> list[str]:
    conn = sqlite3.connect(db_path)
    cursor = conn.execute("SELECT go_id FROM go WHERE accession = ?", (accession,))
    terms = [row[0] for row in cursor]
    conn.close()
    return terms

Paso 2. Extracción de embeddings de ESM3

El modelo ESM3 de pesos abiertos de tamaño pequeño (1.4B, esm3-sm-open-v1) es un modelo de lenguaje con enmascaramiento de tres vías (secuencia-estructura-función). Incluso si solo se ingresa la secuencia, el estado oculto (hidden state) de la última capa de la vía de secuencia genera una representación significativa.

import torch
from esm.models.esm3 import ESM3
from esm.sdk.api import ESMProtein


class ESM3Embedder:
    """Envoltorio para extraer incrustaciones de secuencia ESM3."""

    def __init__(self, device: str = "cuda"):
        self.device = device
        self.model = ESM3.from_pretrained("esm3-sm-open-v1").to(device).eval()

    @torch.no_grad()
    def embed(self, sequence: str) -> torch.Tensor:
        """Promedio de pooling del último estado oculto del rastro de secuencia de una sola secuencia.

        Returns: vector (hidden_dim,). hidden_dim=1536 para el modelo de 1.4B.
        """
        protein = ESMProtein(sequence=sequence)
        encoded = self.model.encode(protein)
        # Extraer el último estado oculto de la secuencia codificada
        output = self.model.forward(sequence_tokens=encoded.sequence.unsqueeze(0).to(self.device))
        # (1, seq_len, hidden_dim) → promedio de pooling → (hidden_dim,)
        return output.embeddings.squeeze(0).mean(dim=0).cpu()

    def embed_batch(self, sequences: list[str], batch_size: int = 4) -> torch.Tensor:
        """Procesamiento por lotes de múltiples secuencias. Ajustar batch_size según los límites de VRAM."""
        embeddings = []
        for i in range(0, len(sequences), batch_size):
            chunk = sequences[i:i + batch_size]
            for seq in chunk:
                embeddings.append(self.embed(seq))
        return torch.stack(embeddings)

Alternativa a ESM2 (si el acceso a ESM3 es difícil, con licencia MIT de uso completamente libre):

from transformers import AutoTokenizer, AutoModel


class ESM2Embedder:
    """Línea base ESM2. hidden_dim=1280 (según t33_650M)."""

    MODEL_ID = "facebook/esm2_t33_650M_UR50D"

    def __init__(self, device: str = "cuda"):
        self.device = device
        self.tokenizer = AutoTokenizer.from_pretrained(self.MODEL_ID)
        self.model = AutoModel.from_pretrained(self.MODEL_ID).to(device).eval()

    @torch.no_grad()
    def embed(self, sequence: str) -> torch.Tensor:
        inputs = self.tokenizer(sequence, return_tensors="pt", truncation=True, max_length=1024).to(self.device)
        outputs = self.model(**inputs)
        # last_hidden_state: (1, seq_len, 1280) → promedio pooling → (1280,)
        return outputs.last_hidden_state.squeeze(0).mean(dim=0).cpu()

Paso 3. Construcción del índice FAISS: ajuste de Flat, IVF y HNSW

FAISS ofrece tres tipos principales de índices, cada uno con un equilibrio diferente entre precisión, velocidad y memoria.

• IndexFlatIP: búsqueda exacta. Precisión del 100%, velocidad O(N), memoria O(N × d).
• IndexIVFFlat: aproximación basada en clústeres. Ajuste de nlist (número de clústeres) y nprobe (número de clústeres visitados durante la búsqueda). Precisión del 95～99%, velocidad O(nprobe/nlist × N).
• IndexHNSWFlat: aproximación basada en grafos. Ajuste de M (número de vecinos del grafo), efConstruction y efSearch. Precisión superior al 98%, muy rápido, memoria algo elevada.

import faiss
import numpy as np


def build_flat_index(embeddings: np.ndarray) -> faiss.Index:
    """Precisión del 100%, para uso según estándares de benchmark."""
    d = embeddings.shape[1]
    index = faiss.IndexFlatIP(d)  # producto interno (equivalente a la similitud coseno después de normalizar las incrustaciones)
    faiss.normalize_L2(embeddings)
    index.add(embeddings)
    return index


def build_ivf_index(embeddings: np.ndarray, nlist: int = 4096) -> faiss.Index:
    """IVF: nlist cerca de sqrt(N) es óptimo empíricamente. Para 500k secuencias, se recomienda nlist~700."""
    d = embeddings.shape[1]
    quantizer = faiss.IndexFlatIP(d)
    index = faiss.IndexIVFFlat(quantizer, d, nlist, faiss.METRIC_INNER_PRODUCT)
    faiss.normalize_L2(embeddings)
    index.train(embeddings)
    index.add(embeddings)
    return index


def build_hnsw_index(embeddings: np.ndarray, M: int = 32) -> faiss.Index:
    """HNSW: M=16~64, mayor precisión y memoria con valores más altos. efConstruction=200 es estándar."""
    d = embeddings.shape[1]
    index = faiss.IndexHNSWFlat(d, M, faiss.METRIC_INNER_PRODUCT)
    index.hnsw.efConstruction = 200
    faiss.normalize_L2(embeddings)
    index.add(embeddings)
    return index

Paso 4. Búsqueda KNN de la secuencia de consulta

def search(
    index: faiss.Index,
    query_embedding: np.ndarray,
    k: int = 10,
    ef_search: int | None = None,
    nprobe: int | None = None,
) -> tuple[np.ndarray, np.ndarray]:
    """Búsqueda de vecinos top-k. efSearch para HNSW y nprobe para IVF pueden ajustarse durante la búsqueda."""
    if ef_search is not None and hasattr(index, "hnsw"):
        index.hnsw.efSearch = ef_search
    if nprobe is not None and hasattr(index, "nprobe"):
        index.nprobe = nprobe
    faiss.normalize_L2(query_embedding.reshape(1, -1))
    distances, indices = index.search(query_embedding.reshape(1, -1), k)
    return distances[0], indices[0]

Paso 5. Transferencia de función: votación por mayoría de términos GO

Se agregan las anotaciones GO de los k vecinos más cercanos para predecir los términos GO de la secuencia de consulta.

from collections import Counter


def transfer_go_annotations(
    neighbor_accessions: list[str],
    neighbor_distances: list[float],
    db_path: Path,
    min_votes: int = 2,
    distance_threshold: float = 0.75,
) -> dict[str, float]:
    """Agrega los términos GO de los vecinos mediante votación mayoritaria ponderada por distancia.

    Returns: {go_id: confidence_score} diccionario.
    """
    weighted_votes: Counter = Counter()
    total_weight = 0.0
    for acc, dist in zip(neighbor_accessions, neighbor_distances):
        if dist < distance_threshold:
            continue
        weight = dist  # Utilizar la propia similitud coseno como peso
        total_weight += weight
        for go_id in lookup_go_terms(acc, db_path):
            weighted_votes[go_id] += weight
    if total_weight == 0:
        return {}
    # Devolver solo los términos GO que aparecen al menos min_votes veces
    return {
        go_id: score / total_weight
        for go_id, score in weighted_votes.items()
        if score >= min_votes * (total_weight / len(neighbor_accessions))
    }

Paso 6. Benchmarking de BLAST (curva PR)

El estándar de oro para comparar la capacidad de detección de homólogos distantes son las bases de datos basadas en estructuras, como SCOP y CATH. En este caso, se toma como referencia el conjunto de pruebas CAFA-3 [4].

from sklearn.metrics import precision_recall_curve, auc


def evaluate_pr_curve(
    query_pairs: list[tuple[str, str]],  # (query_acc, target_acc)
    ground_truth: dict[str, set[str]],   # query_acc → set of true homolog acc
    method: callable,                     # method(query_acc) → [(target_acc, score)]
) -> dict:
    """Cuantificar la capacidad de detección mediante la curva PR."""
    y_true, y_score = [], []
    for query_acc, target_acc in query_pairs:
        preds = method(query_acc)  # [(acc, score), ...]
        scored = {acc: score for acc, score in preds}
        y_true.append(1 if target_acc in ground_truth.get(query_acc, set()) else 0)
        y_score.append(scored.get(target_acc, 0.0))
    precision, recall, _ = precision_recall_curve(y_true, y_score)
    return {
        "auprc": auc(recall, precision),
        "recall_at_p90": max(
            (r for p, r in zip(precision, recall) if p >= 0.9), default=0.0
        ),
    }

Pipeline integrado

def full_pipeline(
    query_sequence: str,
    embedder: ESM3Embedder,
    index: faiss.Index,
    accession_map: list[str],  # orden del índice → accession
    go_db: Path,
    k: int = 20,
) -> dict:
    """Búsqueda y transferencia de función para una nueva secuencia."""
    query_emb = embedder.embed(query_sequence).numpy()
    distances, indices = search(index, query_emb, k=k, ef_search=128)
    neighbor_accs = [accession_map[i] for i in indices]
    go_predictions = transfer_go_annotations(neighbor_accs, distances.tolist(), go_db)
    return {
        "query_length": len(query_sequence),
        "top_neighbors": list(zip(neighbor_accs, distances.tolist())),
        "predicted_go_terms": go_predictions,
    }

Rendimiento, costo y casos de fallo conocidos

Referencias de rendimiento (citación de benchmarks públicos)

Aunque Foldseek sigue siendo muy potente, requiere un cálculo previo de predicción estructural (como AlphaFold), mientras que los embeddings de ESM tienen la ventaja de permitir la búsqueda inmediata utilizando solo la secuencia.

Referencias de medición real por índice (cálculos basados en benchmarks oficiales de FAISS)

Recomendación: Para servicios en tiempo real, usar HNSW; para análisis por lotes, usar IVF; para estándares de benchmark, usar Flat.

Costo esperado de reproducción para el estudiante

• Costo de API: 0 (totalmente local).
• Extracción de embeddings: con una RTX 4060 de 8GB, aproximadamente 2～4 horas para 500k secuencias (40～60 secuencias por segundo).
• Construcción del índice FAISS: 20～30 minutos en CPU para HNSW M=32.
• Búsqueda: 3～6ms por consulta.

3 casos de fallo conocidos (recopilación de la comunidad y literatura)

1. Pérdida de información de dominios debido a la truncación de secuencias largas (> 1024 AA)

2. Al incrustar solo los primeros 1024 AA de una proteína multidominio, los dominios posteriores no se reflejan en el embedding Síntoma: Al incrustar solo los primeros 1024 AA de una proteína multidominio, los dominios posteriores no se reflejan en el embedding. Causa: Límite de longitud máxima de contexto de ESM. Solución: (a) Incrustación por dominio seguida de promedio o concatenación, (b) generar múltiples embeddings mediante una ventana deslizante (sliding window) y aplicar max-pooling. Fuente: HuggingFace ESM Community — GitHub Issues [7].

3. Un valor de nprobe insuficiente en el índice IVF omite homólogos remotos Síntoma: Mantener nprobe en su valor predeterminado (1) provoca una caída abrupta en el recall, ya que es probable que los homólogos remotos se encuentren en otros clústeres. Causa: IVF solo visita los clústeres cercanos al cuantificador, pero los homólogos distantes tienen linajes evolutivos distintos y se asignan a otros clústeres. Solución: Establecer nprobe en al menos el 1% de nlist (50～100). Para servicios en producción, se recomienda una lógica de nprobe adaptativo. Fuente: FAISS Wiki — "Guidelines for choosing an index" [8].

4. La falta de normalización de los embeddings causa confusión entre similitud de coseno y L2 Síntoma: Al usar IndexFlatIP sin normalizar los embeddings, el producto interno (inner product) se desvía de la similitud de coseno, distorsionando la comparación entre proteínas con diferentes longitudes de secuencia. Causa: La norma de los embeddings de ESM varía según la longitud y composición de la secuencia. Solución: faiss.normalize_L2() Es obligatorio llamar a la normalización tanto antes de añadir al índice como antes de la consulta. Fuente: FAISS FAQ — "Normalization for cosine similarity" [9].

Ideas de ampliación

• Integración con la base de datos vectorial Chroma: Utilizar Chroma o Qdrant en lugar de FAISS cuando se requiera filtrado por metadatos (por ejemplo, buscar solo un taxón específico o un campo GO particular).
• Búsqueda híbrida: Ensamble de KNN de embeddings + puntuación HSP de BLAST. Ventajoso en situaciones que exigen alta precisión.
• Combinación de embeddings estructurales: Integrar el track de estructura de ESM3 con el alfabeto 3Di de Foldseek para una búsqueda dual de secuencia-estructura.
• Servicio de interfaz de usuario: Construir un servicio interno de laboratorio con FastAPI + Streamlit. Al ingresar una nueva secuencia, devuelve inmediatamente el top-10 de secuencias similares y la predicción de GO.

Próximo episodio

• Episodio 08 docking-hybrid-diffusion: Listado de proteínas diana candidatas mediante búsqueda de embeddings → Validación de docking con DiffDock-Glide.
• Episodio 11 structure-affinity-boltz: Búsqueda de proteínas similares mediante embeddings y predicción de la afinidad de unión con Boltz-2 (pipeline de reposicionamiento de fármacos).
• Episodio 12 single-cell-perturbation: Selección de candidatos diana para la perturbación génica unicelular mediante similitud de embeddings.
• Episodio 13 protein-design-multimodal: Uso del track de estructura/función de ESM3, en correspondencia con el track de secuencia de este episodio.
• Episodio 14 bio-mcp-agent: Exposición del pipeline de búsqueda de este episodio como herramienta MCP para su uso por agentes autónomos.

Referencias

1. Rives A, Meier J, Sercu T, et al. "Biological structure and function emerge from scaling unsupervised learning to 250 million protein sequences." PNAS 2021. https://www.pnas.org/doi/10.1073/pnas.2016239118
2. Rao R, Meier J, Sercu T, Ovchinnikov S, Rives A. "Transformer protein language models are unsupervised structure learners." ICLR 2021. https://openreview.net/forum?id=fylclEqgvgd
3. van Kempen M, Kim S, Tumescheit C, et al. "Fast and accurate protein structure search with Foldseek." Nature Biotechnology 2024. https://www.nature.com/articles/s41587-023-01773-0
4. CAFA (Critical Assessment of Function Annotation): https://www.biofunctionprediction.org/cafa/
5. Hayes T, Rao R, Akin H, et al. "Simulating 500 million years of evolution with a language model." bioRxiv 2024 (ESM3). https://www.biorxiv.org/content/10.1101/2024.07.01.600583v1
6. UniProt-GOA: https://www.ebi.ac.uk/GOA/downloads
7. HuggingFace ESM Community Discussion: https://huggingface.co/facebook/esm2_t33_650M_UR50D/discussions
8. FAISS Wiki — Guidelines for choosing an index: https://github.com/facebookresearch/faiss/wiki/Guidelines-to-choose-an-index
9. FAISS FAQ — Normalization: https://github.com/facebookresearch/faiss/wiki/FAQ
10. EvolutionaryScale ESM3 GitHub: https://github.com/evolutionaryscale/esm
11. FAISS GitHub: https://github.com/facebookresearch/faiss
12. UniProt Swiss-Prot: https://www.uniprot.org/
13. Biopython: https://biopython.org/
14. SCOPe (Structural Classification of Proteins extended): https://scop.berkeley.edu/
15. CATH (Class · Architecture · Topology · Homology) database: https://www.cathdb.info/