¿Por qué abordar este capítulo primero?

El alineamiento de secuencias es el alfabeto de la bioinformática. Tanto BLAST como BWA, e incluso el AAM (Alineamiento Múltiple de Secuencias) de AlphaFold, son extensiones de los dos problemas de alineamiento de secuencias que tratamos en este capítulo. Sin embargo, al buscar materiales en coreano dentro de Corea del Sur, la mayoría dice simplemente "con usar BLAST basta". No hay prácticamente ningún material que pueda responder por qué BLAST es tan rápido o por qué BLOSUM utiliza esos números específicos.

En este capítulo, no nos limitamos a usar herramientas; derivamos las ecuaciones de recurrencia, rellenamos manualmente la cuadrícula e implementamos desde cero en Python. La tabla que rellenemos hoy se convertirá en la expansión de semillas de BWA que manejaremos mañana y en la columna vertebral del Viterbi de HMM que aprenderemos pasado mañana. No nos apresuremos.

Una tarde de 1970 en Princeton

Saul Needleman y Christian Wunsch publicaron su artículo en J Mol Biol en marzo de 1970. La pregunta que plantearon es hoy demasiado familiar para nosotros.

"Tenemos dos secuencias de proteínas. ¿Qué tan similares son?"

Ahora parece obvio, pero en aquel entonces ni siquiera existía una definición de qué significaba "qué tan similares". ¿Cómo debemos alinearlas para obtener el mejor resultado? Por ejemplo, dados dos conjuntos de secuencias GATTACA y GCATGCU, alguien podría emparejarlas carácter por carácter desde el principio, otro desde el final, y otro abriendo huecos en el centro. ¿Cuál es la respuesta correcta?

La respuesta dada por Needleman y Wunsch es sorprendentemente simple: la que tiene la puntuación más alta entre todos los alineamientos posibles. El problema radica en que el número de "todos los alineamientos posibles" explota exponencialmente con respecto a las longitudes de las secuencias $m$ y $n$. Resolvieron este problema mediante programación dinámica. Dado que Bellman nombró la PD en 1953, esto llegó al campo de la bioinformática apenas 17 años después.

Función de puntuación · Cuadrícula · Ecuación de recurrencia

¿Cómo se asigna la puntuación?

Al alinear dos secuencias ocurren tres tipos de eventos:

1. Coincidencia (match): correspondencia entre caracteres iguales → +s_match (por ejemplo, +1)
2. No coincidencia (mismatch): correspondencia entre caracteres diferentes → -s_mismatch (por ejemplo, -1)
3. Hueco (gap): un espacio vacío en uno de los lados → -d (por ejemplo, -2)

Supongamos que tenemos el siguiente alineamiento:

G A T T A C A
G - A T T A C

3 coincidencias (G, A, T), 3 discrepancias (A-A, T-T, A-C), y 1 hueco (--). Al asignar +1 por coincidencia, -1 por discrepancia y -2 por hueco, la puntuación total es 3 × 1 − 3 × 1 − 1 × 2 = −2.

Rellenemos la cuadrícula

El punto clave es definir de la siguiente manera la puntuación del alineamiento óptimo F(i, j) para las secuencias X = x₁x₂…xₘ e Y = y₁y₂…yₙ:

F(i, j) = max {
    F(i-1, j-1) + s(xᵢ, yⱼ)   ← diagonal: coincidencia / desacuerdo
    F(i-1, j  ) - d            ← desde arriba: carácter en X, gap en Y
    F(i  , j-1) - d            ← desde la izquierda: carácter en Y, gap en X
}

Aquí, s(xᵢ, yⱼ) es la puntuación de coincidencia si los dos caracteres son iguales, o la puntuación de desacuerdo si no lo son.

Las condiciones de frontera son las siguientes:

• F(0, 0) = 0
• F(i, 0) = −i · d
• F(0, j) = −j · d

La primera fila y la primera columna representan simplemente la acumulación continua de huecos.

¿Por qué esta ecuación de recurrencia es correcta?

Si F(i, j) representa el mejor alineamiento entre los primeros i caracteres de X y los primeros j caracteres de Y, entonces su última columna debe corresponder a una de tres posibilidades: diagonal (ambos son caracteres), arriba (solo un carácter de X) o izquierda (solo un carácter de Y). Si ya conocemos las soluciones óptimas de los subproblemas anteriores F(i-1, j-1), F(i-1, j) y F(i, j-1), entonces el máximo entre esos tres valores constituye la solución óptima actual. Esta es la subestructura óptima, la razón por la cual la programación dinámica es aplicable.

Relleno manual

Rellenemos la cuadrícula con X = GATTACA, Y = GCATGCU, coincidencia +1, desacuerdo −1 y hueco −2. Para ahorrar espacio, lo acortaremos a X = GAT, Y = GAC.

• F(1,1): G vs G — coincidencia → F(0,0)+1 = 1. Gana la diagonal.
• F(2,2): A vs A — coincidencia → F(1,1)+1 = 2.
• F(3,3): T vs C — desacuerdo → max(F(2,2)−1, F(2,3)−2, F(3,2)−2) = max(1, −2, −2) = 1.

El valor en la esquina inferior derecha, +1, es la puntuación óptima. A partir de aquí, siguiendo las flechas hacia atrás (traceback), se obtiene el alineamiento real.

G A T
G A C

Partida 2 + Desajuste 1 = 2 − 1 = 1. Coincide con el valor de la cuadrícula.

Implementación completa en Python en 20 líneas

def needleman_wunsch(x, y, match=1, mismatch=-1, gap=-2):
    m, n = len(x), len(y)
    F = [[0] * (n + 1) for _ in range(m + 1)]

    # Condiciones de frontera
    for i in range(m + 1):
        F[i][0] = i * gap
    for j in range(n + 1):
        F[0][j] = j * gap

    # Llenado de la cuadrícula
    for i in range(1, m + 1):
        for j in range(1, n + 1):
            s = match if x[i-1] == y[j-1] else mismatch
            F[i][j] = max(
                F[i-1][j-1] + s,
                F[i-1][j]   + gap,
                F[i][j-1]   + gap,
            )
    return F[m][n], F

score, table = needleman_wunsch("GATTACA", "GCATGCU")
print(score)  # Puntuación óptima

Esto es todo. 20 líneas. El EMBOSS needle utilizado en la práctica también se apoya finalmente sobre este esqueleto.

Adjuntar el rastro de retroceso

Si desea obtener realmente la alineación, añada una función de rastro de retroceso.

def traceback(F, x, y, match=1, mismatch=-1, gap=-2):
    i, j = len(x), len(y)
    ax, ay = "", ""
    while i > 0 or j > 0:
        s = match if (i > 0 and j > 0 and x[i-1] == y[j-1]) else mismatch
        if i > 0 and j > 0 and F[i][j] == F[i-1][j-1] + s:
            ax = x[i-1] + ax; ay = y[j-1] + ay; i -= 1; j -= 1
        elif i > 0 and F[i][j] == F[i-1][j] + gap:
            ax = x[i-1] + ax; ay = "-" + ay; i -= 1
        else:
            ax = "-" + ax; ay = y[j-1] + ay; j -= 1
    return ax, ay

Se reconstruyen las coincidencias, gaps y desacuerdos retrocediendo desde la esquina inferior derecha para determinar la dirección de origen.

Evaluación con Rosalind

Rosalind GLOB entra en el conjunto de datos, descarga y pega el código anterior tal cual para ejecutarlo, luego envía un único número entero como resultado. Puedes verificar la respuesta correcta en 5 minutos. La experiencia de ver que tu propio algoritmo pasa realmente es más poderosa que ver tutoriales 100 veces.

Complejidad y la barrera de la práctica

Needleman-Wunsch tiene tiempo O(mn) · espacio O(mn). Si las longitudes de las dos secuencias son solo 10,000 cada una, hay que llenar 10⁸ celdas. Comparar directamente el genoma humano (3 × 10⁹) con otra secuencia genómica de esta manera es irreal.

Por lo tanto, la práctica se divide en dos caminos.

1. Alineamiento local: Extrae solo las partes similares de las dos secuencias, no todo el conjunto. La siguiente entrega cubrirá Smith-Waterman (M07), que es ese camino.
2. Heurística: En lugar de llenar la cuadrícula, encuentra rápidamente una solución aproximada mediante el método de semilla (seed) + extensión. La entrega tres después cubrirá BLAST (M10), que es ese camino.

Herramientas de mapeo como BWA también llenan la cuadrícula de Needleman-Wunsch solo localmente y solo en las partes necesarias. Es decir, todas las herramientas de alineamiento que encontrarás en el futuro se tratan de cómo recortar e insertar eficientemente esta cuadrícula.

Mapeo de CS — El entero de DP

Recuerda el problema de la ruta óptima en la cuadrícula DP tratado en DryBench. Needleman-Wunsch es la versión de cadenas de caracteres de ese problema.

• Coordenadas de la cuadrícula = Combinación de subcadenas de las dos secuencias
• Ruta = Un alineamiento
• Función de puntuación = Costo de herencia de cada celda
• Estructura óptima parcial = Reutilizar resultados anteriores mientras se llena la cuadrícula
• Traceback = Restaurar la ruta óptima

Si aprendiste la distancia de edición (distancia de Levenshtein) en CS, es exactamente eso. La distancia de edición asigna +1 a cada inserción, eliminación y sustitución, buscando el mínimo; mientras que Needleman-Wunsch busca maximizar recompensando los acordes desde una perspectiva evolutiva y penalizando los gaps e incompatibilidades. Solo invirtiendo los signos, es exactamente el mismo algoritmo.

Rama hacia la siguiente entrega

• Siguiente entrega (M07): Smith-Waterman — Expansión a alineación local. El momento en que se cortan los valores negativos a 0 en la cuadrícula cambia completamente el problema.
• Dos entregas después (M08): Penalización de gaps afines — Debe haber una penalización diferente para abrir un gap y continuar uno, para ajustarse mejor a los eventos evolutivos reales.
• Tres capítulos después (M09): BLOSUM/PAM — match=+1, mismatch=-1 es demasiado brusco. Los puntajes se determinan mediante probabilidades evolutivas de los aminoácidos.

Si desea profundizar más

El texto principal ha sido reescrito directamente por BPD. Si ya comprende los principios, puede avanzar con los cursos especializados a continuación.

• MIT OCW 7.91J — Conferencias de Christopher Burge sobre Alineamiento Local (BLAST) y Estadística (subtítulos completos, traducción automática superior al 95%). Aborda la inducción de la cuadrícula con mayor rigor matemático.
• Harvard STAT115 — Conferencias de Xiaole Shirley Liu sobre Alineamiento Temprano de Secuencias (uno contra uno) (subtítulos completos, excelente traducción automática). Muestra cómo llenar la cuadrícula paso a paso mediante pizarra.
• Stanford CS262 — Lista de reproducción del profesor Gill Bejerano (subtítulos completos). Destaca por su interpretación desde la perspectiva de la informática sobre el problema de la ruta óptima.
• Libro de referencia: Compau & Pevzner Bioinformatics Algorithms — Colección gratuita de problemas integrada con Rosalind. Permite continuar directamente con las implementaciones en Python de este capítulo.

Al resolver GLOB, LOCA y GAFF en Rosalind utilizando Python, estará listo para abordar el siguiente capítulo. ¡Debe poner manos a la obra!