Búsqueda de motivos: encontrar sitios de unión de factores de transcripción en secuencias promotoras

Al finalizar este tema

Podrán crear una herramienta que busque simultáneamente múltiples sitios de unión de factores de transcripción en secuencias promotoras, combinando el Trie y la búsqueda binaria aprendidos en el libro de texto. Comprenderán mediante código la estructura de datos subyacente a los algoritmos prácticos de coincidencia de múltiples patrones como Aho-Corasick.

Este artículo es un ejemplo educativo general. La búsqueda de motivos en la práctica utiliza modelos mucho más sofisticados, como matrices de peso de posición (PWM) y modelos de fondo.

"500 motivos × secuencia de 100 kb": la trampa del enfoque ingenuo

Supongamos que intentan encontrar todos los sitios de unión de factores de transcripción conocidos (cada uno de 6 a 12 pb) en una secuencia promotora de ratón de 100 kb.

Enfoque ingenuo:

def naive_search(sequence: str, motifs: list[str]) -> list[tuple[int, str]]:
    hits = []
    for motif in motifs:
        for i in range(len(sequence) - len(motif) + 1):
            if sequence[i:i+len(motif)] == motif:
                hits.append((i, motif))
    return hits

El problema de este enfoque es que itera sobre cada motivo. La complejidad temporal es O(número de motivos × longitud de la secuencia) = O(500 × 100.000) = 5 × 10⁷. En Python, esto tarda unos segundos.

¿Y si la secuencia fuera el genoma humano de 3 mil millones de pares de bases? 5 × 10⁵ × 3 × 10⁹ = 1,5 × 10¹⁵. No terminaría antes del fin del universo.

La solución correcta es usar un trie. Combinamos los 500 motivos en un único árbol y recorremos la secuencia una sola vez, siguiendo el árbol para detectar las coincidencias. La complejidad temporal se reduce a O(longitud de la secuencia + número de coincidencias).

De la caja negra a los componentes: explorando el trie

Componente 1: Estructura de datos Trie

Un trie es un árbol en el que cada nodo representa un carácter y la ruta desde la raíz hasta una hoja representa una cadena.

Si construimos un trie para el conjunto de motivos {"TATA", "TATT", "GCGC"}:

root
        /  \
       T    G
       |    |
       A    C
       |    |
       T    G
      / \   |
     A*  T* C*

* indica el final de un motivo completo. Si tres motivos comparten la parte inicial (TA), esa porción se almacena una sola vez en el árbol.

Implementación en Python:

class TrieNode:
    def __init__(self) -> None:
        self.children: dict[str, TrieNode] = {}
        self.matches: list[str] = []


class Trie:
    def __init__(self) -> None:
        self.root = TrieNode()

    def insert(self, pattern: str) -> None:
        node = self.root
        for char in pattern:
            if char not in node.children:
                node.children[char] = TrieNode()
            node = node.children[char]
        node.matches.append(pattern)

Componente 2: Búsqueda mediante recorrido en el árbol Trie

Ahora, se recorre el árbol Trie desde cada posición de la secuencia para verificar las coincidencias.

def search(trie: Trie, sequence: str) -> list[tuple[int, str]]:
    hits: list[tuple[int, str]] = []

    for i in range(len(sequence)):
        node = trie.root
        j = i
        while j < len(sequence) and sequence[j] in node.children:
            node = node.children[sequence[j]]
            for match in node.matches:
                hits.append((i, match))
            j += 1

    return hits

Desde cada posición i, se recorre el árbol hasta la profundidad máxima (longitud del motivo más largo). Si la longitud máxima de 500 motivos es 12, se realizan hasta 12 consultas por posición. La complejidad total es O(longitud de la secuencia × longitud máxima del motivo).

La optimización práctica (Aho-Corasick) añade enlaces de fallo para calcular de antemano la posición de reinicio. En este tutorial, solo establecemos los conceptos básicos utilizando un trie estándar.

Comparación con la búsqueda binaria

Para comprender el enfoque del trie, resulta claro compararlo con una alternativa: arreglo ordenado + búsqueda binaria.

Enfoque basado en un arreglo de motivos ordenados:

def sorted_search(sequence: str, sorted_motifs: list[str], max_motif_len: int) -> list[tuple[int, str]]:
    from bisect import bisect_left

    hits = []
    for i in range(len(sequence)):
        for length in range(1, max_motif_len + 1):
            substring = sequence[i:i+length]
            idx = bisect_left(sorted_motifs, substring)
            if idx < len(sorted_motifs) and sorted_motifs[idx] == substring:
                hits.append((i, substring))
    return hits

Complejidad temporal de este enfoque: O(longitud de la secuencia × longitud máxima del motivo × logaritmo del número de motivos).

Aunque solo hay un factor logarítmico adicional en comparación con el trie, en la práctica, el trie es más rápido. La razón es la localidad de caché. El recorrido del trie implica accesos a memoria contiguos para los nodos relacionados, mientras que la búsqueda binaria realiza accesos aleatorios dispersos por el array.

La búsqueda binaria es ventajosa frente al trie cuando el conjunto de motivos cambia con poca frecuencia y el tiempo de carga inicial es crítico. La construcción del trie es lenta, pero su búsqueda es rápida. Por otro lado, la búsqueda binaria solo requiere ordenamiento, por lo que su carga inicial es rápida. Se trata de un compromiso.

Ejemplo de uso práctico

Verifiquemos con un escenario sencillo.

# Motivos de unión de factores de transcripción conocidos (ejemplo)
motifs = [
    "TATAAA",       # TATA box
    "CAAT",         # CAAT box
    "GGGCGG",       # GC box
    "CACGTG",       # E-box
    "TGACTCA"       # AP-1
]

# Secuencia promotora virtual
promoter = "GCTATAAACCAATGGGCGGATGCACGTGCCCTGACTCAAG"

# Construcción del trie
trie = Trie()
for motif in motifs:
    trie.insert(motif)

# Búsqueda
hits = search(trie, promoter)
for pos, motif in sorted(hits):
    print(f"Posición {pos}: {motif}")

# Salida:
# Posición 2: TATAAA
# Posición 9: CAAT
# Posición 12: GGGCGG
# Posición 21: CACGTG
# Posición 28: TGACTCA

Se detectaron los cinco motivos en un único recorrido de la secuencia.

Relleno: los dos espacios en blanco que deben completar

Espacio en blanco 1: Compatibilidad con los códigos IUPAC

Los sitios reales de unión de los factores de transcripción presentan variabilidad. Se representan mediante códigos IUPAC (R = A/G, Y = C/T, W = A/T, etc.).

IUPAC = {
    "A": {"A"}, "C": {"C"}, "G": {"G"}, "T": {"T"},
    "R": {"A", "G"}, "Y": {"C", "T"}, "S": {"C", "G"},
    "W": {"A", "T"}, "K": {"G", "T"}, "M": {"A", "C"},
    "B": {"C", "G", "T"}, "D": {"A", "G", "T"},
    "H": {"A", "C", "T"}, "V": {"A", "C", "G"},
    "N": {"A", "C", "G", "T"}
}


def search_with_iupac(trie: Trie, sequence: str) -> list[tuple[int, str]]:
    """
    El motivo puede contener códigos IUPAC.
    Ejemplo: "TATA**W**A" coincide tanto con TATAAA como con TATATA.
    """
    # TODO: Al recorrer el trie, si la clave 'children' de cada nodo es un código IUPAC, coincidir con el conjunto expandido
    pass

Pista: Mantenga la estructura del árbol tal cual y realice la búsqueda coincidiendo sequence[j] con las claves de cada nodo hijo en el conjunto IUPAC ampliado.

Espacio vacío 2: Búsqueda automática de complementos inversos

Dado que el ADN es de doble cadena, los motivos también pueden aparecer como complementos inversos.

def reverse_complement(seq: str) -> str:
    complement = {"A": "T", "T": "A", "G": "C", "C": "G"}
    return "".join(complement.get(b, b) for b in reversed(seq))


def search_both_strands(motifs: list[str], sequence: str) -> list[tuple[int, str, str]]:
    """
    Insertar cada motivo y su complemento inverso en el trie y buscar.
    Devuelve: (posición, motivo, "forward" o "reverse")
    """
    # TODO: Al construir el trie, insertar tanto la dirección directa como la complementaria inversa de cada motivo
    # Registrar el motivo original y la dirección de cada entrada
    pass

Pista: TrieNode.matches se expande a list[tuple[str, str]] (motivo original, dirección).

Reflexión: diferencias con la búsqueda práctica de motivos

Matriz de ponderación de posición (PWM): en la práctica, la mayoría de los sitios de unión de factores de transcripción no coinciden mediante una coincidencia estricta de cadenas, sino mediante una coincidencia probabilística. Se asigna una probabilidad a cada nucleótido en cada posición para calcular la puntuación de las secuencias candidatas. Los motivos de bases de datos como JASPAR utilizan este formato.

Modelo de fondo: en la práctica, se calcula la probabilidad de coincidencia aleatoria utilizando un modelo de fondo, y solo se informan las coincidencias estadísticamente significativas. FIMO (suite MEME) es la herramienta estándar.

Aho-Corasick: si se añaden enlaces de fallo al trie, el tiempo se vuelve verdaderamente lineal: O(longitud de la secuencia + número de coincidencias). Es el estándar para la coincidencia múltiple de patrones en la práctica.

Escala del genoma completo: para buscar en el genoma humano (3 mil millones de pares de bases), las estructuras de datos indexadas, como los arreglos de sufijos y el índice FM, son alternativas. Son el núcleo de los alineadores como BWA y Bowtie.

Proyecto de expansión

1. Carga de datos JASPAR: descarga 100 PWM de factores de transcripción humanos de la base de datos JASPAR real y explora los promotores con tu herramienta.

2. Visualización: dibuja los resultados de la búsqueda como un mapa de promotores. Utiliza matplotlib o pistas de estilo IGV.

3. Expansión de Aho-Corasick: añade enlaces de fallo para completar el mecanismo de coincidencia en tiempo lineal real.

Mapa de componentes de este capítulo

• [F] Trie: fusiona múltiples cadenas en un árbol. Nodo = carácter, camino = cadena.
• [F] Búsqueda binaria: acceso alternativo con una matriz ordenada + bisect. Comprende el equilibrio entre el trie y la búsqueda binaria.
• [W] E/S de archivos: análisis de FASTA, etc. (se proporciona en el script completo).

[F] = implementado por ti / [W] = se proporciona con código completo.