¿Por qué se necesita un HMM?

En la alineación aprendida anteriormente, se correspondían directamente los caracteres de dos secuencias. Sin embargo, en las secuencias reales hay estados ocultos que no se pueden observar directamente. ¿Es este lugar un exón o un intrón? ¿Esta región es rica en GC o rica en AT? ¿Esta posición está dentro o fuera de una isla CpG? Para responder a estas preguntas, es necesario extraer la disposición de los estados que generaron los caracteres, en lugar de centrarse únicamente en los propios caracteres.

El Modelo Oculto de Markov (Hidden Markov Model, HMM) aborda precisamente este problema. Dada una secuencia observada, se remonta mediante un modelo probabilístico a la disposición de estados subyacente. En esta sección se definen las tres probabilidades del HMM y se deriva el algoritmo de Viterbi para encontrar dicha disposición de estados en una cuadrícula. La similitud con la cuadrícula de Needleman-Wunsch es sorprendente.

Las tres probabilidades del Modelo Oculto de Markov

Un HMM se define mediante los siguientes cinco componentes:

• Conjunto de estados S — por ejemplo, rico en GC (H), rico en AT (L)
• Conjunto de observaciones V — por ejemplo, el alfabeto del ADN A, C, G, T
• Probabilidad inicial pi — probabilidad de cada estado al inicio de la secuencia
• Probabilidad de transición a — probabilidad de pasar del estado i al estado j
• Probabilidad de emisión b — probabilidad de que aparezca el carácter v dado el estado s

Existen dos supuestos clave:

1. Propiedad de Markov: el siguiente estado depende únicamente del estado actual (el historial pasado es irrelevante).
2. Independencia de las observaciones: el carácter observado depende únicamente del estado actual.

Gracias a estos dos supuestos, la probabilidad conjunta de la disposición de estados se puede descomponer mediante multiplicación. La probabilidad conjunta de la disposición de observaciones y la disposición de estados es igual al producto de la probabilidad inicial por la primera emisión, multiplicado por (transición por emisión) hasta el tiempo T.

Ejemplo sencillo — Detector de regiones ricas en GC vs ricas en AT

Simplifiquemos el problema de detectar islas CpG en un genoma bacteriano. Dos estados y cuatro caracteres en el alfabeto de observaciones.

• Estados: H (rico en GC, alto) · L (rico en AT, bajo)
• Observaciones: A, C, G, T

Tabla de probabilidades configurada arbitrariamente (en la realidad se estima a partir de datos; esto se aborda en la sección M20 sobre Baum-Welch).

Probabilidades iniciales:

Probabilidades de transición a:

El estado H tiene una fuerte inercia para mantenerse como H, y el estado L también tiende a mantenerse como L. Los cambios de estado implican penalizaciones en cada paso.

Probabilidades de emisión b:

El estado H emite más C y G, mientras que el estado L emite más A y T. Esto coincide cualitativamente con las estadísticas de las regiones ricas en GC.

Supongamos que la secuencia observada es GCAT. El número de arreglos de estados posibles para generar esta secuencia es $2^4$, es decir, 16. ¿Cuál es el arreglo con la probabilidad más alta?

Algoritmo de Viterbi — Arreglo de estados de máxima probabilidad

La pregunta que plantea Viterbi es la siguiente: Dada una secuencia de observaciones, ¿cuál es el arreglo de estados cuya probabilidad conjunta es máxima?

El número de arreglos de estados posibles crece exponencialmente con la longitud de la secuencia elevada a la potencia del número de estados. Un enfoque exhaustivo es inviable. En su lugar, reformulamos el problema como uno de encontrar la ruta óptima en una cuadrícula.

Definición de la cuadrícula: Las filas representan el tiempo $t$ y las columnas los estados $s$. La celda $(t, s)$ almacena la probabilidad de la subruta óptima que termina en el estado $s$ en el instante $t$.

Ecuación recursiva (descrita en lenguaje natural):

• $V(t, s) = \max_{u} [ V(t-1, u) \times a(u, s) ] \times b(s, o_t)$
• Condición inicial: $V(1, s) = \pi(s) \times b(s, o_1)$

Esta ecuación recursiva es análoga a la de Needleman-Wunsch. Las únicas diferencias son dos:

1. En lugar de sumas, se utilizan multiplicaciones.
2. Los valores de las celdas de la cuadrícula no son puntuaciones de alineación, sino probabilidades de subarreglos de estados.

La razón por la que la programación dinámica (DP) es válida es exactamente la misma: el arreglo óptimo que termina en el estado $s$ en el instante $t$ consiste en el arreglo óptimo que termina en algún estado $u$ en el instante $t-1$, al cual se le añade una transición de $u$ a $s$ y una emisión del carácter $o_t$ desde $s$. Esto refleja la propiedad de subestructura óptima.

Llenemos la cuadrícula manualmente

Para las observaciones GCAT, llenamos la cuadrícula usando las tablas de probabilidad anteriores.

t=1, observación G:

• $V(1, H) = 0.5 \times 0.4 = 0.20$
• $V(1, L) = 0.5 \times 0.1 = 0.05$

t=2, observación C:

• $V(2, H) = \max(0.20 \times 0.7, 0.05 \times 0.4) \times 0.4 = 0.14 \times 0.4 = 0.056$
• $V(2, L) = \max(0.20 \times 0.3, 0.05 \times 0.6) \times 0.1 = 0.06 \times 0.1 = 0.006$

Registramos que ambas celdas provienen del estado H (puntero hacia atrás).

t=3, observación A:

• $V(3, H) = \max(0.056 \times 0.7, 0.006 \times 0.4) \times 0.1 = 0.0392 \times 0.1 = 0.00392$
• $V(3, L) = \max(0.056 \times 0.3, 0.006 \times 0.6) \times 0.4 = 0.0168 \times 0.4 = 0.00672$

Aquí, por primera vez, la celda L supera a la celda H. Esto se debe a que A se emite mucho mejor desde el estado L.

t=4, observación T:

• V(4, H) = max(0.00392 × 0.7, 0.00672 × 0.4) × 0.1 = 0.002744 × 0.1 = 0.0002744
• V(4, L) = max(0.00392 × 0.3, 0.00672 × 0.6) × 0.4 = 0.004032 × 0.4 = 0.0016128

Resumen en tabla:

La probabilidad máxima es V(4, L) = 0.0016128. Al retroceder mediante los backpointers, la secuencia de estados óptima es H H L L. Las dos primeras letras de la observación GCAT se interpretan como generadas por el estado rico en GC, y las dos últimas por el estado rico en AT.

Este es el resultado de Viterbi: hemos extraído la secuencia de estados más probable oculta detrás de las observaciones mediante la ruta óptima en la cuadrícula.

Escalado logarítmico — un truco esencial en la práctica

En el ejemplo anterior, la longitud de la secuencia era 4, por lo que las multiplicaciones fueron seguras. En secuencias reales, que pueden tener miles o decenas de miles de posiciones, seguir multiplicando probabilidades hace que los valores converjan a 0, provocando un desbordamiento por debajo del rango (underflow).

La solución es el escalado logarítmico: convierte las multipliciones en sumas.

• log V(t, s) = para cada estado anterior u, max de [ log V(t-1, u) más log a(u, s) ] y luego suma log b(s, o_t)

log(0) se trata como menos infinito (en Python float('-inf')). Este truco es común en HMM y en la mayoría de las programaciones dinámicas probabilísticas.

Implementación en Python (30 líneas)

import math

def viterbi(obs, states, start_p, trans_p, emit_p):
    T = len(obs)
    N = len(states)
    log = math.log

    V = [[float('-inf')] * N for _ in range(T)]
    back = [[0] * N for _ in range(T)]

    for s in range(N):
        V[0][s] = log(start_p[s]) + log(emit_p[s][obs[0]])

    for t in range(1, T):
        for s in range(N):
            best_prev, best_score = 0, float('-inf')
            for u in range(N):
                score = V[t-1][u] + log(trans_p[u][s])
                if score > best_score:
                    best_score, best_prev = score, u
            V[t][s] = best_score + log(emit_p[s][obs[t]])
            back[t][s] = best_prev

    path = [0] * T
    path[T-1] = max(range(N), key=lambda s: V[T-1][s])
    for t in range(T-2, -1, -1):
        path[t] = back[t+1][path[t+1]]

    return [states[s] for s in path], V[T-1][path[T-1]]

states = ['H', 'L']
obs_map = {'A': 0, 'C': 1, 'G': 2, 'T': 3}
obs = [obs_map[c] for c in 'GCAT']
start_p = [0.5, 0.5]
trans_p = [[0.7, 0.3], [0.4, 0.6]]
emit_p  = [[0.1, 0.4, 0.4, 0.1],
           [0.4, 0.1, 0.1, 0.4]]

path, log_prob = viterbi(obs, states, start_p, trans_p, emit_p)
print(path, math.exp(log_prob))

Con 30 líneas está hecho. Incluso si el número de estados es 10 y la longitud de la secuencia es 10,000, el tamaño de esta cuadrícula sigue siendo de 100,000 celdas. Sigue calculándose instantáneamente.

Complejidad

• Tiempo: T por N al cuadrado. T es la longitud de la secuencia y N es el número de estados. Esto se debe a que, para cada estado en cada instante, se recorren todos los estados del instante anterior.
• Espacio: T por N. Para toda la cuadrícula y los backpointers.

Es exactamente la misma clase que Needleman-Wunsch. Simplemente, el número de estados de HMM reemplaza al tamaño del alfabeto.

Mapeo de CS — Tablero probabilístico de DP en cuadrícula

El problema de la ruta óptima en la cuadrícula DP tratado en DryBench se reproduce aquí tal cual.

• Coordenadas de la cuadrícula = (instante t, estado s)
• Ruta = un arreglo de estados
• Función de puntuación = suma acumulada de probabilidades logarítmicas
• Subestructura óptima = el arreglo óptimo hasta t-1 seguido de la última transición y emisión
• Retroceso (Traceback) = restauración de la ruta óptima mediante backpointers

Si Needleman-Wunsch era una cuadrícula secuencia-secuencia, Viterbi es una cuadrícula instante-estado. Solo han cambiado la multiplicación por suma y las puntuaciones de alineación por probabilidades logarítmicas; el esqueleto es completamente idéntico.

Esta similitud no es casualidad. DP cambia de ropa —multiplicación probabilística, optimización mediante sumas, caminos más cortos— pero su esqueleto es uno solo. Dominarlo en la práctica hace que todos los futuros algoritmos de modelos probabilísticos (Forward-Backward, Baum-Welch, Beam Search, etc.) se vean como ramas que brotan de la misma raíz.

Evaluación en Rosalind

Abra el problema HMMV de Rosalind. Se proporcionan los parámetros del HMM y las observaciones, y se solicita el arreglo de estados de Viterbi. Pegar el código anterior tal cual pasará la evaluación. La experiencia de que sus propias 30 líneas superen la evaluación real profundiza más la comprensión que usar herramientas cien veces.

Ramas hacia el siguiente capítulo

• Siguiente capítulo (M19): Algoritmo Forward-Backward — Si Viterbi extrae "una única ruta con mayor probabilidad", Forward-Backward extrae "la probabilidad de cada estado en cada instante". La relación entre estos dos algoritmos es el núcleo de esta serie.
• Dos capítulos después (M20): Baum-Welch EM — Ahora hemos establecido la tabla probabilística arbitrariamente, pero en la práctica debemos aprender esta tabla únicamente a partir de los datos observados. El bucle de autoaprendizaje del algoritmo EM.
• Tres capítulos después (M21): Profile HMM y detección de genes — Una sección práctica que extiende el modelo simple de 2 estados para detectar estructuras génicas (exones, intrones e intergénicos).

Para profundizar más

El texto principal es una narrativa reconstruida por BPD. Si ya comprende los principios, puede profundizar con las siguientes clases ortodoxas.

• Rabiner (1989), A tutorial on hidden Markov models and selected applications in speech recognition, Proceedings of the IEEE — Un tutorial que captura la esencia de los modelos ocultos de Markov. Sigue siendo una referencia máxima hoy en día.
• MIT OCW 7.91J — Conferencias de HMM del profesor Christopher Burge (con subtítulos completos). Se estudian las inducciones de retícula mediante pizarra.
• Durbin, Eddy, Krogh, Mitchison Biological Sequence Analysis — El libro de texto estándar para aplicar HMM a secuencias biológicas.
• Rastreado de Textos de Bioinformática de Rosalind — Resolver los tres problemas HMMV, HMMFB y BW en orden hará que esta serie se asiente en tu práctica.

Antes de encontrarnos en el próximo capítulo, resuelva manualmente estos ejemplos una vez. Escriba las ocho multiplicaciones en papel y verifique que los resultados coincidan con la tabla anterior. Esta experiencia es lo que hace que el resto de los HMM se asiente en su práctica.