Alineación de secuencias: alineación exacta de dos ADN mediante programación dinámica

Al finalizar este tema

Podrán construir una herramienta de alineación que, al combinar la programación dinámica, la comparación entre recursión e iteración y las matrices 2D aprendidas en el libro de texto, determine con precisión dónde ocurren inserciones, deleciones y sustituciones al alinear exactamente dos secuencias de ADN. Comprenderán mediante código los principios internos del alineamiento que BLAST realiza aparentemente como por arte de magia.

Este texto es un ejemplo educativo general. Dado que el análisis de variantes es una tarea cotidiana en la bioinformática, se ha utilizado como tema central.

"¿Por qué no coinciden las posiciones de las mutaciones?" — La trampa de la comparación ingenua

Supongamos que están comparando una subsecuencia silvestre del gen BRCA1 con una secuencia mutante secuenciada a partir de un paciente.

tipo silvestre:   ATGCTAGCATGCA
mutante: ATGCAGCATGCA

La comparación más sencilla consiste en alinear los elementos uno a la vez, según su posición.

def compare_naive(seq1: str, seq2: str) -> list[tuple[int, str, str]]:
    diffs = []
    for i in range(min(len(seq1), len(seq2))):
        if seq1[i] != seq2[i]:
            diffs.append((i, seq1[i], seq2[i]))
    return diffs

Al aplicar esto a las dos secuencias mencionadas anteriormente, obtenemos lo siguiente.

posición 4: T → A
posición 5: A → G
posición 6: G → C
posición 7: C → A
posición 8: A → T
posición 9: T → G
posición 10: G → C
posición 11: C → A

Es extraño. Después de la posición 4, casi todas las posiciones son diferentes. ¿Es posible que haya tantas mutaciones?

Si se comparan visualmente ambas secuencias, se puede ver la respuesta: en el tipo salvaje, simplemente se ha eliminado (deleción) la T en la posición 4.

tipo silvestre:   ATGC[T]AGCATGCA
mutante: ATGC AGCATGCA
                ↑ A partir de aquí, el desplazamiento aumenta en uno

La eliminación de un solo carácter desplazó todas las posiciones siguientes, dando la impresión de que se produjeron múltiples sustituciones. Una comparación simple no detecta estas eliminaciones e inserciones.

La solución a este problema es el alineamiento de secuencias. Se insertan tantas brechas (gaps, - caracteres) como sea necesario entre las dos secuencias para corregir las discrepancias.

Tipo silvestre:   ATGCTAGCATGCA
Mutante: ATGC-AGCATGCA   ← hueco en la posición 4

Ahora, las posiciones posteriores a la 5 coinciden perfectamente. Queda claro que la verdadera variante es una deleción.

El problema es decidir dónde insertar este espacio. Cuando dos secuencias tienen alrededor de 100 caracteres, el número de combinaciones posibles para insertar espacios crece exponencialmente. ¿Cómo podemos encontrar la combinación que produce el alineamiento "más natural"? Este artículo construye la respuesta desde cero.

Primero, veamos el producto final (ejecutemos la caja negra)

La herramienta que construiremos tomará dos secuencias y devolverá el alineamiento óptimo.

aligned1, aligned2, score = align(seq1, seq2)

print(aligned1)   # 'ATGCTAGCATGCA'
print(aligned2)   # 'ATGC-AGCATGCA'
print(score)      # 8 (ejemplo)

=== Resultado del alineamiento ===
Wild:   ATGCTAGCATGCA
Mut:    ATGC-AGCATGCA
Score:  8

=== Comparación ingenua sin alineamiento ===
Wild:   ATGCTAGCATGCA
Mut:    ATGCAGCATGCA
Posiciones de discrepancia: 8 (casi todas)

=== Mutaciones reales después del alineamiento ===
Posición 4: Deleción de T — 1 mutación real

Aunque las dos secuencias son idénticas, la interpretación cambia por completo antes y después de la alineación. La alineación es el primer paso y una etapa crucial en la interpretación de variantes. A continuación, analizaremos los principios que sustentan la creación de esta alineación.

¿De qué componentes está compuesta esta herramienta (desglose de componentes)?

alineador de secuencias
   ┌──────────────────────────────────────────────────┐
   │  [Entrada] Cargar dos secuencias ─── Componente: entrada/salida de cadenas/archivos  │  ← Proporcionado (herramienta)
   │              │                                      │
   │              ▼                                      │
   │  [Paso 1] Llenar la matriz de puntuación                          │
   │        Componente: matriz 2D + programación dinámica                    │  ← Crear directamente ★
   │              │                                      │
   │              ▼                                      │
   │  [Paso 2] Restaurar la ruta óptima mediante retroceso                     │
   │        Componente: selección recursiva/iterativa                           │  ← Crear directamente ★
   │              │                                      │
   │              ▼                                      │
   │  [Salida] Dos secuencias alineadas + puntuación                          │
   └──────────────────────────────────────────────────┘

📌 If these concepts are new (links above)

• Dynamic programming — avoiding exponential complexity by storing subproblems
• Recursion vs iteration — trade-off between elegance and efficiency
• 2D array — values stored in a grid

The new concepts you will implement are dynamic programming, 2D matrix, and recursion/iteration choice (3 items). Sequence handling is provided as a ready-made tool. Just 3 items — within cognitive limits.

Step 1: Prepare sequences (provided)

First, prepare two sequences to align. In practice, you would read them from a FASTA file, but for browser-based practice, they are defined as strings.

# Secuencia parcial BRCA1 silvestre (ejemplo educativo)
wild = "ATGCTAGCATGCA"

# Secuencia mutante (deleción de T en la posición 4)
mut = "ATGCAGCATGCA"

assert len(wild) == 13
assert len(mut) == 12
assert all(c in "ACGT" for c in wild)
assert all(c in "ACGT" for c in mut)

El objetivo es alinear estas dos secuencias superponiéndolas e insertando espacios - de manera adecuada.

Paso 2: Definir las reglas de puntuación (solución proporcionada)

Primero, debemos definir cómo medir la "calidad" de la alineación. El método estándar es un sistema de puntuación.

MATCH = 1        # +1 si se alinean bases iguales
MISMATCH = -1    # -1 si se alinean bases diferentes
GAP = -2         # -2 por cada inserción de hueco

def score_pair(a: str, b: str) -> int:
    """Puntuación del alineamiento entre dos bases (o huecos)."""
    if a == b:
        return MATCH
    return MISMATCH

assert score_pair("A", "A") == 1
assert score_pair("A", "T") == -1

Esta regla de puntuación determina la preferencia del algoritmo de alineación.

• Si la penalización por huecos es baja (por ejemplo, -0.5), se prefiere una alineación con muchos huecos.
• Si la penalización por huecos es alta (por ejemplo, -5), se toleran más desajustes que huecos.

En la práctica de la bioinformática, se utilizan matrices de puntuación sofisticadas como BLOSUM/PAM. Aquí, usamos una regla simple con fines didácticos.

Paso 3 de creación: comenzar con la recursión y experimentar la explosión ★ (recursión vs. iteración)

✍️ Sección para completar directamente. Componente = recursión. Primero, calcula la puntuación de alineación mediante recursión y comprende por qué no es práctico.

Expresar el problema de la puntuación de alineación con recursión es sorprendentemente simple. Comenzando desde el final de ambas secuencias, se prueban tres opciones en cada posición mediante recursión:

• Alinear el último carácter de ambas secuencias → score_pair(s1[-1], s2[-1]) + align(s1[:-1], s2[:-1])
• Alinear el último carácter de s1 con un hueco → GAP + align(s1[:-1], s2)
• Alinear el último carácter de s2 con un hueco → GAP + align(s1, s2[:-1])

Se elige la opción con la puntuación más alta. Implementado con recursión, esto se vería así:

def align_recursive(s1: str, s2: str) -> int:
    """Calcula recursivamente solo la puntuación del alineamiento óptimo (sin retroceso)."""
    if not s1:
        return len(s2) * GAP
    if not s2:
        return len(s1) * GAP
    return max(
        align_recursive(s1[:-1], s2[:-1]) + score_pair(s1[-1], s2[-1]),
        align_recursive(s1[:-1], s2) + GAP,
        align_recursive(s1, s2[:-1]) + GAP,
    )

# Verificación: solo con ejemplos cortos
assert align_recursive("AA", "AA") == 2       # dos coincidencias
assert align_recursive("AT", "A") == 1 + GAP  # coincidencia + hueco = 1 - 2 = -1

Esto proporciona la respuesta correcta. Sin embargo, no se puede utilizar con secuencias largas. Si ejecuta el siguiente código, comprenderá la razón.

import time

# Con una longitud de unos 15, se entra en el infierno de la recursión
short1 = "ATGCTAGCATGCAAG"
short2 = "ATGCAGCATGCAAG"

t0 = time.time()
score = align_recursive(short1, short2)
elapsed = time.time() - t0
print(f"Recursión con longitud 15: {elapsed:.2f}s")
# Tardará unos segundos. Para longitud 20 tardará minutos, para 25 horas

La razón es que la recursión repite el cálculo de los mismos subproblemas cientos de millones de veces. El valor align_recursive("ATGCT", "ATGC") se calcula repetidamente en múltiples ramas del árbol recursivo. Cada vez, desde el principio.

Este es el poder de la complejidad temporal O(3^n). En este punto, entra en juego la programación dinámica.

🤔 Indicación para la autoexplicación Dibuje align_recursive("AT", "A") como un árbol recursivo. ¿Cuántos subproblemas aparecen repetidamente? (Basta con expandir hasta longitudes 3 y 4 para que las superposiciones sean evidentes).

Paso 4 de construcción: guardar la respuesta en una cuadrícula ★ (matriz 2D + programación dinámica)

✍️ Sección para completar manualmente. Componente = matriz 2D + programación dinámica. Almacenar las respuestas de los subproblemas en una cuadrícula evita cálculos repetidos.

La idea es la siguiente: guardamos la respuesta de align(s1[:i], s2[:j]) en la cuadrícula dp[i][j]. Así, no será necesario calcular el mismo subproblema dos veces.

🔎 ¿Qué es la programación dinámica? (draw.io — dynamic-programming) Cuando la respuesta a un problema grande se construye a partir de las respuestas de problemas más pequeños, guardar las respuestas de los problemas pequeños evita tener que calcularlas dos veces. La explosión exponencial de la recursión se convierte en tiempo polinómico con una sola técnica de memorización. El enfoque de llenado de cuadrícula se conoce específicamente como "programación dinámica ascendente".

El tamaño de la cuadrícula es (len(s1)+1) × (len(s2)+1). Una fila/columna adicional sirve como condición base para el alineamiento con cadenas vacías.

def build_score_matrix(s1: str, s2: str) -> list[list[int]]:
    \"\"\"dp[i][j] = puntuación del alineamiento óptimo entre s1[:i] y s2[:j].\"\"\"
    n, m = len(s1), len(s2)
    dp = [[0] * (m + 1) for _ in range(n + 1)]

    # Caso base: alinear con cadenas vacías solo implica huecos
    for i in range(1, n + 1):
        dp[i][0] = dp[i-1][0] + GAP
    for j in range(1, m + 1):
        dp[0][j] = dp[0][j-1] + GAP

    # Rellenar la cuadrícula (Bottom-Up)
    for i in range(1, n + 1):
        for j in range(1, m + 1):
            match = dp[i-1][j-1] + score_pair(s1[i-1], s2[j-1])
            delete = dp[i-1][j] + GAP    # carácter de s1 y hueco
            insert = dp[i][j-1] + GAP    # hueco y carácter de s2
            dp[i][j] = max(match, delete, insert)
    return dp

dp = build_score_matrix(wild, mut)

# Verificación: la celda final debe coincidir con la respuesta recursiva
assert dp[len(wild)][len(mut)] == align_recursive(wild, mut)
# Verificación del tamaño
assert len(dp) == len(wild) + 1
assert len(dp[0]) == len(mut) + 1
# La columna izquierda es la acumulación pura de gaps
assert dp[3][0] == 3 * GAP

Se obtiene el mismo resultado, pero la velocidad es muy diferente.

import time

# Antes, ¿cuánto tiempo tardó la recursión para una longitud de 15?
t0 = time.time()
score = build_score_matrix(short1, short2)[len(short1)][len(short2)]
elapsed = time.time() - t0
print(f"Longitud 15 programación dinámica: {elapsed*1000:.2f}ms")
# Nivel de milisegundos

Lo que antes tardaba segundos con la recursión, ahora se resuelve en milisegundos. Esta es la potencia de O(n×m). Cada celda de la cuadrícula se calcula solo una vez.

🔎 ¿Por qué la cuadrícula es la solución (Drawee — matrix-2d-array) Los subproblemas que dependen de dos índices (i, j) se representan naturalmente en una cuadrícula 2D. Cada celda representa un subproblema y las flechas entre celdas indican las dependencias. Cuando esta estructura es visible, resulta más fácil abordar los problemas de programación dinámica.

🤔 Indicación para la autoexplicación Para calcular dp[i][j], solo necesitamos dp[i-1][j-1], dp[i-1][j] y dp[i][j-1]. ¿Es necesario almacenar toda la cuadrícula? En realidad, con solo la fila anterior, podemos rellenar la siguiente, por lo que podemos obtener solo la puntuación con un espacio de O(min(n,m)). ¿Cuándo no sería posible esta optimización? (Pista: retroceso)

Paso 5: Reconstruir la solución rastreando la cuadrícula ★

✍️ Sección para completar manualmente. El elemento clave es el retroceso iterativo. Comenzamos en la última celda de la cuadrícula y rastreamos el camino inverso para ver cómo se obtuvo ese valor.

Conocer solo la puntuación y conocer la alineación exacta son cosas diferentes. Para reconstruir la alineación real (dónde se insertan los huecos), debemos realizar un retroceso desde el final hasta el inicio de la cuadrícula.

En cada celda dp[i][j], retrocedemos para ver cuál de las tres opciones elegimos:

• Si elegimos dp[i-1][j-1] + score_pair(s1[i-1], s2[j-1]) → movimiento diagonal (coincidencia/no coincidencia)
• Si elegimos dp[i-1][j] + GAP → movimiento hacia arriba (carácter de s1 y hueco)
• Si elegimos dp[i][j-1] + GAP → movimiento hacia la izquierda (hueco y carácter de s2)

def traceback(dp: list[list[int]], s1: str, s2: str) -> tuple[str, str]:
    aligned1: list[str] = []
    aligned2: list[str] = []
    i, j = len(s1), len(s2)

    while i > 0 or j > 0:
        current = dp[i][j]
        if i > 0 and j > 0 and current == dp[i-1][j-1] + score_pair(s1[i-1], s2[j-1]):
            aligned1.append(s1[i-1])
            aligned2.append(s2[j-1])
            i -= 1
            j -= 1
        elif i > 0 and current == dp[i-1][j] + GAP:
            aligned1.append(s1[i-1])
            aligned2.append("-")
            i -= 1
        else:
            aligned1.append("-")
            aligned2.append(s2[j-1])
            j -= 1

    return "".join(reversed(aligned1)), "".join(reversed(aligned2))

a1, a2 = traceback(dp, wild, mut)

# Verificación: longitud idéntica · restauración de la secuencia original al eliminar gaps
assert len(a1) == len(a2)
assert a1.replace("-", "") == wild
assert a2.replace("-", "") == mut
# El punto de deleción que insertamos debe ser encontrado realmente
assert "-" in a2                 # Debe haber un gap en el lado de la mutación

Esta función se implementó mediante iteración, pero la misma lógica también se puede implementar con recursión. Cuando el tamaño de la cuadrícula es muy grande, la iteración es más adecuada para evitar desbordamientos de pila, mientras que la recursión puede mejorar la legibilidad del código en algunos casos. La elección depende del tamaño del problema.

🔎 ¿Cuándo usar recursión e iteración? (drivere — recursion-vs-iteration) La recursión es elegante para expresar problemas de forma natural. La iteración es más segura porque no utiliza la pila y, a menudo, es más rápida. Problemas pequeños y necesidad de elegancia → recursión / Problemas grandes y necesidad de estabilidad → iteración. Esta herramienta tiene tres pasos (llenado de la cuadrícula) implementados con iteración, y el paso cuatro (retroceso) también se implementó con iteración. Ambos podrían implementarse con recursión, pero al aumentar la longitud de la secuencia, aumenta el riesgo de desbordamiento de pila, por lo que se optó por la iteración.

🤔 Indicación autoexplicativa En traceback, seleccionamos uno de los tres candidatos. Si las puntuaciones de dos candidatos son iguales (empate), ¿cuál deberíamos elegir? ¿Podría cambiar el resultado real de la ordenación? (Pista: puede haber múltiples ordenaciones óptimas).

Unificar las piezas: clase de ordenación completa

Unimos ambas funciones en una sola herramienta.

class Aligner:
    def __init__(self, match=1, mismatch=-1, gap=-2):
        self.match = match
        self.mismatch = mismatch
        self.gap = gap

    def _score(self, a: str, b: str) -> int:
        return self.match if a == b else self.mismatch

    def align(self, s1: str, s2: str) -> tuple[str, str, int]:
        n, m = len(s1), len(s2)
        dp = [[0] * (m + 1) for _ in range(n + 1)]
        for i in range(1, n + 1):
            dp[i][0] = dp[i-1][0] + self.gap
        for j in range(1, m + 1):
            dp[0][j] = dp[0][j-1] + self.gap
        for i in range(1, n + 1):
            for j in range(1, m + 1):
                dp[i][j] = max(
                    dp[i-1][j-1] + self._score(s1[i-1], s2[j-1]),
                    dp[i-1][j] + self.gap,
                    dp[i][j-1] + self.gap,
                )
        # Retroceso
        a1: list[str] = []
        a2: list[str] = []
        i, j = n, m
        while i > 0 or j > 0:
            if i > 0 and j > 0 and dp[i][j] == dp[i-1][j-1] + self._score(s1[i-1], s2[j-1]):
                a1.append(s1[i-1]); a2.append(s2[j-1]); i -= 1; j -= 1
            elif i > 0 and dp[i][j] == dp[i-1][j] + self.gap:
                a1.append(s1[i-1]); a2.append("-"); i -= 1
            else:
                a1.append("-"); a2.append(s2[j-1]); j -= 1
        return "".join(reversed(a1)), "".join(reversed(a2)), dp[n][m]


aligner = Aligner()
a1, a2, score = aligner.align(wild, mut)
print(a1)
print(a2)
print(f"score = {score}")

# Verificación: los resultados de la clase coinciden con los de la función
dp = build_score_matrix(wild, mut)
assert score == dp[len(wild)][len(mut)]

Esta herramienta es una versión simplificada del algoritmo Needleman-Wunsch que estamos analizando hoy. Las herramientas prácticas (como BLAST, needle de EMBOSS, etc.) simplemente añaden matrices de puntuación más complejas, puntuaciones separadas para la apertura y extensión de huecos, opciones de alineamiento local y otros detalles. La estructura básica es exactamente la misma que acabas de crear.

Análisis exhaustivo del rendimiento: ¿por qué es tan rápido?

import time

def bench(n: int, m: int):
    s1 = "AT" * (n // 2)
    s2 = "AG" * (m // 2)
    t0 = time.time()
    build_score_matrix(s1, s2)
    return time.time() - t0

for size in [10, 50, 100, 200]:
    t = bench(size, size)
    print(f"Longitud {size:4d}: {t*1000:.2f}ms  ({size*size} celdas)")
assert bench(100, 100) < bench(200, 200)

Los valores numéricos varían según la máquina, pero la dirección es siempre la misma. Es directamente proporcional al número de celdas (n×m). La diferencia clave es que, mientras que la versión recursiva crece exponencialmente a 3^n, la versión de programación dinámica crece a n×m.

🔎 Resumen con notación Big-O (Draper — big-o-notation)

• Versión recursiva ingenua: O(3^n) — exponencial. Ya es inviable para secuencias de longitud 30.
• Versión de programación dinámica: O(n·m) — polinomial. Incluso para secuencias de 10.000 elementos, tarda solo unos segundos.
• Almacenar una sola celda en la cuadrícula transforma el crecimiento exponencial en polinomial. Esta es la magia de la programación dinámica.

Existen otros enfoques (reflexión multietapa)

• Alineamiento global vs. local: Lo que hemos creado es Needleman-Wunsch (global), que alinea toda la secuencia de principio a fin. En la práctica, a menudo se utiliza Smith-Waterman (local), que busca fragmentos de una secuencia dentro de otra. La estructura básica del algoritmo es casi idéntica; solo difieren el sistema de puntuación y las condiciones iniciales. Cuándo usar cada uno: para secuencias filogenéticas evolutivas = global / para encontrar dominios dentro de un gen = local.
• Heurísticas de la familia BLAST: Nuestro algoritmo de programación dinámica escala proporcionalmente a n·m (las longitudes de las dos secuencias), lo que sigue siendo elevado para genomas de gran tamaño (miles de millones de pares de bases). BLAST obtiene soluciones aproximadas rápidamente sacrificando la solución exacta (semilla → extensión). Compromiso: 100% óptimo vs. velocidad práctica.
• Refinamiento de las matrices de puntuación: Nosotros utilizamos una regla binaria de coincidencia/no coincidencia, pero en la práctica se utilizan matrices como BLOSUM62/PAM250 que asignan diferentes puntuaciones a cada par de residuos. Esto refleja la similitud química de las secuencias de proteínas.
• Penalización de huecos afín: Nosotros utilizamos un coste fijo de -2 por cada hueco, pero en realidad se aplica una regla doble en la que el coste inicial para abrir un hueco es alto y el coste para extenderlo es bajo (por ejemplo, apertura -10, extensión -1). Esto refleja la observación evolutiva de que un único hueco largo es más natural que varios huecos cortos.

Punto clave: "Almacenar los posibles subproblemas en una cuadrícula para evitar cálculos repetidos". Este principio se reutiliza innumerables veces más allá del alineamiento de secuencias, como en la búsqueda de rutas, la distancia de edición y la recuperación de información. La estructura básica que acaban de crear abre muchas posibilidades.

Próximos pasos (enlaces de salida en la parte inferior)

• Al rellenar la cuadrícula, ¿por qué se visitan las celdas en este orden? → Dependencias y orden de visita
• Optimización de la memoria cuando las secuencias son muy largas → Técnica de compresión de programación dinámica 1D
• Combinación con los índices creados anteriormente → Aplicaciones Indexación de bases de datos de secuencias

Práctica guiada (problema independiente)

1. Escenario de inserción: Alinee la secuencia silvestre ATGCATGCA con la mutante ATGCXATGCA (X = inserción de una base adicional) y verifique mediante aserciones que la posición de inserción se haya determinado correctamente.
2. Ajuste de puntuación: ¿Cómo cambian los resultados si utiliza Aligner(gap=-5)? ¿Y si gap=-0.5? Resuma cómo la penalización por huecos determina el "carácter" del alineamiento.
3. Alineamientos óptimos múltiples: Cree un par de secuencias que tenga dos alineamientos óptimos y escriba una versión extendida de traceback que devuelva todos los alineamientos óptimos.
4. Desafío: Memoria O(min(n,m)): Si solo se necesita la puntuación, basta con mantener solo la fila anterior en lugar de toda la cuadrícula. Implemente esta versión optimizada. (Esto implica renunciar al retroceso, pero es útil para secuencias muy largas.)

Resumen

Hemos resuelto el problema de "encontrar el alineamiento óptimo entre dos secuencias" mediante tres componentes clave.

• La programación dinámica transformó una explosión exponencial (O(3^n)) en tiempo polinómico (O(n·m)).
• La matriz 2D sirvió como la cuadrícula que almacena las respuestas a los subproblemas.
• La elección entre recursión e iteración garantizó la estabilidad tanto en el llenado de la cuadrícula como en el retroceso.

Cuando BLAST arroja resultados de alineamiento como si fuera magia, ahora podrá ver lo que hay detrás. No es magia, sino simplemente aplicar las secuencias a la cuadrícula de programación dinámica aprendida en los libros de texto.

Este texto es un ejemplo educativo general. Las herramientas prácticas de alineamiento de secuencias (BLAST/EMBOSS/Bowtie, etc.) incorporan matrices de puntuación sofisticadas, penalizaciones dobles por huecos, heurísticas de "seed-and-extend" y paralelización, entre otros elementos. La versión detallada de estos componentes puede ser añadida por ustedes sobre este esqueleto o delegada a herramientas validadas.