Bioinfo 02a - Departamento de Informática

Transcript Bioinfo 02a - Departamento de Informática

IV Alineamiento múltiple de secuencias

Andrés Moreira Departamento de Informática UTFSM

Alineamiento múltiple de secuencias

Alineamiento múltiple: •se alinea un grupo de secuencias •(no es lo mismo que alinearlas a todas de a pares!).

•El alineamiento múltiple busca identificar letras homólogas en secuencias homólogas. •En general se aplica sobre grupos de secuencias que ya se saben homólogas.

Aún en sus variantes más simples, es NP-duro.

Alineamiento múltiple de secuencias

El alineamiento múltiple ayuda a: Hacer árboles filogenéticos, y estudiar la evolución a nivel molecular.

Alineamiento múltiple de secuencias

El alineamiento múltiple ayuda a: Hacer árboles filogenéticos, y estudiar la evolución a nivel molecular.

Preguntas del tipo: •¿De dónde salió la secuencia X?

•¿Qué filogenia tienen las secuencias X, Y, Z?

•¿Este genoma, contiene algún gen que haya llegado por transferencia horizontal?

•¿Hay genes descendientes una duplicación?

•...etc

Alineamiento múltiple de secuencias

El alineamiento múltiple ayuda a: Caracterizar familias de proteínas  identificación de regiones conservadas, funcionalmente importantes (sitios activos de enzimas, sitios de asociación con otras moléculas...)  inferencia de función (cuando alguien entra a la familia)  inferencia de estructura (global y de regiones específicas) para miembros nuevos de la familia

Alineamiento múltiple de secuencias

El alineamiento múltiple ayuda a: Caracterizar familias de proteínas  También sirve para buscar miembros nuevos, buscando homólogos en una BD, que sigan el mismo “perfil”.

A C G T

Perfiles y logos

Alineamiento y perfil: A A A T A 0.8

0 0 0.2

C T C C T 0 0.6

0 0.4

A A G G C 0.4

0.2

0.4

0 T A C G G 0.2

0.2

0.4

0.2

A A A A A 1 0 0 0 T T T T G 0 0 0.2

0.8

Perfiles y logos

Sequence Logo

A C G T 0.8

0 0 0.2

0 0.6

0 0.4

0.4

0.2

0.4

0 0.2

0.2

0.4

0.2

1 0 0 0 0 0 0.2

0.8

Perfiles y logos

En cada columna tenemos una densidad de probabilidad; podemos calcular su entropía:  Si

p ik

es la frecuencia de la letra k en la columna i, entonces

(

)

= −

k p ik

log

p ik

•Es una forma de evaluar la calidad del alineamiento en esa columna: H(i)=0  todas iguales (la misma letra) H(i) máximo  cero información!

•Se suelen tomar los gaps como una letra más, aunque hay variaciones.

Perfiles y logos

En el logo, se escala cada columna por su entropía (con respecto al máximo, o mejor dicho, a su diferencia con respecto al máximo: columnas más altas reducen más la incerteza).

 Más útil que 

Perfiles y logos

También se hacen de proteínas:

Perfiles y logos

En ese caso se colorean los aminoácidos por grupos: •Hidrofóbicos: A, V, I, L, M, F, P, W •Básicos: K, R, H •Acídicos: D, E •Polares: S, T, N, Q, G, Y, C

¿Cómo evaluar?

¿Cómo elegimos el alineamiento?

Opción 1: el que construye un experto, a ojo. Se usan. Y a veces son excelentes.

Opción 2: el que encontramos con un algoritmo. Pero en ese caso necesitamos una función a maximizar.

 ¿Qué es un buen alineamiento múltiple?

¿Cómo evaluar?

X En rigor, quisiéramos buscar el árbol filogenético que nos diera una historia más probable...

U V C D W E A B

(

)

(

)

(

)

(

)

(

)

(

)

(

)

(

)

(

) (nos deshacemos de lo desconocido, sumándolo)

, 

w x

(

)

(

)

(

)

(

)

(

)

(

)

(

)

(

)

(

) ... pero es impracticable maximizar eso.

 Falta demasiada información (e.g., las ramas pueden ser de distintos largos), y además es demasiado cómputo.

¿Cómo evaluar?

Lo más común: se dan puntajes por columnas del alineamiento, y se busca optimizar la suma de esos puntajes (es la extensión más obvia de lo que se hace en el alineamiento de a pares).

Se usan varias opciones. Entre otras: •Entropía •“Suma de a pares”: suma el puntaje de cada par de letras de la columna (calculado como al alinear de a pares).

¿Cómo evaluar?

Suma de a pares:

...A...

...S...

...T...

AA: 4, AS: 1, AT:0 AS: 1, AT: 0 ST: 1 Puntaje: 4+1+0+1+0+1 = 7

Desventaja de la suma de a pares: •No tiene detrás un modelo probabilista razonable.

Desventaja de la entropía (y de la suma de a pares): •No tiene detrás un modelo de la filogenia!

Se parcha. Se ponen pesos para distintos niveles de semejanza (en alguna medida, recoge filogenia).

Prog. dinámica multidimensional

Ok, escogemos una forma de poner puntaje, por columna. ¿Cómo lo optimizamos?

Primera idea: extender el método de programación dinámica usado para alinear de a pares.

 Así como antes avanzábamos por cada una de las secuencias, moviéndonos en un plano, ahora avanzamos por k secuencias.

0 1 A A - 1 - A A 2 T T T 3 G - G 4 C C C

(0,0,0) (1,2,1) (3,3,3)

  

(1,1,0)



(2,3,2)



(4,4,4)

Prog. dinámica multidimensional

s i,j,k

Prog. dinámica multidimensional

= max

s i-1,j-1,k-1 s i-1,j-1,k s i-1,j,k-1 s i,j-1,k-1 s i-1,j,k s i,j-1,k s i,j,k-1 + +



(v i , w j , u k )



(v i , w j , _ ) +



(v i , _, u k ) +



(_, w j , u k ) +



(v i , _ , _) (_, w j , _) +



(_, _, u k )

Prog. dinámica multidimensional

El orden es O(2 k n k ), para k secuencias de largo n.

Para una proteína típica, n=100.

 Más allá de k=4,5, se vuelve inservible.

Para k=7, n=100, requeriría varios terabytes de RAM.

 No se usa!!

Otras ideas...

Otra idea: alinear todos los pares de secuencias, y combinar todos esos alineamientos en uno común.

Después de todo, un alineamiento múltiple induce alineamientos de a pares: (se proyecta el camino 3d sobre los lados del cubo)

Otras ideas...

x: y: z: AC-GCGG-C AC-GC-GAG GCCGC-GAG

induce:

x: y: ACGCGG-C; x: ACGC-GAC; z: AC-GCGG-C; y: GCCGC-GAG; z: AC-GCGAG GCCGCGAG

...así que cabría la esperanza de hacer el proceso inverso.

Otras ideas...

Pero: •Alineamiento óptimo múltiple no induce óptimos de a pares.

•Alineamientos óptimos de a pares son difíciles de combinar en un buen alineamiento múltiple.

•Alineamientos de a pares pueden ser inconsistentes (es decir, puede no existir uno múltiple que los induzca).

 no va por ahí la cosa.

Alineamiento progresivo

Algo que sí se usa (y mucho): alineamiento progresivo.

Idea: alineo primero dos secuencias. Luego alineo el resultado de eso con una tercera secuencia. Y así sucesivamente.

¿En qué orden?

Una opción sería tomar las dos más cercanas (evaluando la similaridad de a pares). Luego, la que se acerque más al resultado de eso. Y así, sucesivamente. No se hace.

Alineamiento progresivo

CLUSTAL: •Calcula la similaridad de cada par (usando NW).

C E D B A •Construye un “árbol guía”.

•Alinea primero lo más cercano, y va subiendo por el árbol.

Ojo: se requiere alinear secuencias con alineamientos (y alineamientos con alineamientos!).

A U V B C X D W E

Alineamiento progresivo

v 1 v 2 v 3 v 4

El arbol guía se construye con neighbour joining (un algoritmo que veremos al hablar de filogenia).

v 1

v 2 v 3

.17 .87 .28 -

v 4

.59 .33 .62 -

v v v v 1 3 4 2

Calcular: v 1,3 = alineamiento(v 1 , v 3 ) v 1,3,4 v 1,2,3,4 = alineamiento((v = alineamiento((v 1,3 ),v 1,3,4 4 ) ),v 2 )

Alineamiento progresivo

Pueden aparecer nuevos gaps en el camino.

S1 S3 S1 S3 S4 S2 S2 S4 S1 S3 S2 S4

Gap nuevo para optimizar el alineamiento de (S2,S4) con (S1,S3)

Alineamiento progresivo

¿Cómo asignar puntajes al alinear una letra contra un perfil, o un perfil contra otro?

•Lo sencillo: escribir la nueva columna completa (con lo que ya llevo alineado en cada lado para esa posición), y evaluar el puntaje de la columna.

•Alternativa: comparar sólo con el perfil: No será A contra AC, sino A contra (0.5,0.5,0,0). Esto es útil cuando hay muchas secuencias, y cuando se usa un perfil para buscar en una base de datos.

Alineamiento progresivo

Pero ojo: si se compara una letra y contra un perfil p (en que p(x) es la frecuencia de la letra x ), y estamos usando la matriz de sustitución Q, entonces el puntaje NO es 

x p

( pues ese Q se obtuvo como Q(a,b)=log q (a,b), donde q (a,b) representaba la probabilidad de que el alineamiento fuera no aleatorio. Para que siga teniendo sentido sumar, hay que usar log y para dos perfiles, 

x x

(

log (



y y

)

(

) (

)

2 (

)

(

)

Alineamiento progresivo

Clustal es lejos lo más usado hoy en día por los biólogos para alinear secuencias.

http://www.clustal.org/ Viene en dos sabores: ClustalW (línea de comando), ClustalX (+GUI). La GUI no es algo menor: alinear es un arte.

La W viene de “weighted”: añade pesos en distintos niveles del alineamiento (según la distancia a la raíz), quitándole importancia excesiva a grupos de secuencias muy similares.

Alineamiento progresivo

Otras cosas que agrega: •Usa matrices de sustitución distintas, según nivel de semejanza entre secuencias.

•La penalización de gaps también la va variando: •Se penalizan menos allí donde ya hay gaps •Se penalizan más cerca de donde ya hay gaps.

•Además los modifica con consideraciones bioquímicas.

•Hay posibilidades para que el usuario predetermine el árbol guía, o para realinear tramos de la secuencia.

Alineamiento progresivo

Problemas del alineamiento progresivo: •Se casa mucho con los alineamientos con que parte; esos gaps y alineación se acarrean después. Si las secuencias no son muy parecidas, puede ser malo.

Respuesta: métodos iterativos, y otros.

•Es global! La semejanza podría haber estado en trozos de las secuencias.

Alineamiento progresivo

Lo que debiera ser Lo que da Clustal Respuestas a eso: 1. Cortar antes los tramos similares y alinearlos.

2. Usar algoritmos que se preocupan de eso.

Alineamiento iterativo

Iterar: para deshacerse de la dependencia de lo inicial.

En Clustal es posible iterar de varias formas: •Usar el alineamiento para construir un árbol filogenético, y usar eso como árbol guía en un nuevo alineamiento. [Eso requiere pasar por otro software].

•Realinear las secuencias ya alineadas (con sus gaps). +Hacerlo a partir de varios sets de parámetros distintos. Ir viendo qué zonas son importantes. Realinear a mano. Etc. (El manual recomienda generar al menos 10 alineamientos).

Alineamiento iterativo

Otros métodos iterativos: •DIALIGN: pone énfasis en alineamientos locales entre las secuencias.

•MUSCLE: va revisando el árbol guía

Alineamiento iterativo: MUSCLE

MUSCLE, idea base: va revisando el árbol guía.

1) Parte como Clustal 2) Sobre la marcha va revisando: acaso el árbol guía ha cambiado, y acaso alguna parte se puede realinear mejor. También considera (testea) otros cambios al árbol.

¿Ha cambiado el árbol?

• Calcula la similaridad de los alineamientos de a pares inducidos por el actual alineamiento múltiple.

• Si eso da un árbol distinto al que está usando, realínea según eso, y ve si el puntaje general mejora.

Alineamiento : DIALIGN II

Énfasis en los alineamientos locales.

•Busca buenos (=buen P-value) alineamientos locales (de a pares).

•Los va re-evaluando según qué tan compatibles son unos con otros...

•Hasta que obtiene un buen set de alineamientos locales ordenados y coherentes. Y con eso, alínea .

Alineamiento : DIALIGN II

Es especialmente útil para: •similitudes que son sólo locales •secuencias muy lejanas Y es una gracia el no tener que tratar los “gaps” de manera ad-hoc.

Alineamiento iterativo: alg. evolucionarios

Métodos evolucionarios: Lo que se muta es la ubicación de los gaps (un alineamiento no es otra cosa que una distribución de gaps!).

• Algoritmos genéticos: SAGA • Simulated Annealing: MSASA Variante “hill climbing”: muta en direcciones provechosas.

Alineamiento iterativo: alg. evolucionarios

Cruce de alineamientos Una mutación ...etc

Alineamiento: nota al margen

•Hay muuuchos algoritmos y software de MSA.

•Casi toda técnica de optimización conocida por el hombre, se ha aplicado. Para bien o para mal.

•Clasificar los algoritmos es difícil; combinan características de varios lados (estocásticos/deterministas, progresivos/simultaneos, exactos/heurísticos, iterativos o no...), y mutan de año en año.

•Aquí daré sólo las ideas principales de algunos; meternos al detalle sería interminable.

Alineamiento: nota al margen

•Esa diversidad, y los avances permanentes, son importantes. El alineamiento (de a pares y múltiple) es ubicuo en bioinformática.

Dos reviews instructivos, legibles y no muy largos: •Multiple sequence alignment, R.C. Edgar & S. Batzoglou (2006), Curr Op Struct Biol 16(3) 368-373; doi:10.1016/j.sbi.2006.04.004

 Disponible en http://ai.stanford.edu/~serafim/Publications/2006_MSA_COSB.pdf

•Recent evolutions of multiple sequence alignment algorithms, C. Notredame (2007), PLoS Comput Biol 3(8) e123, doi:10.1371/journal.pcbi.0030123

 Disponible ahí.

Alineamiento: T-Coffee

Idea: Generar una biblioteca de alineamientos (globales y locales) de a pares entre las secuencias. Tratar de hacer el alineamiento múltiple que sea consistente con esa biblioteca.

•Primero genera una biblioteca, haciendo alineamientos de a pares, global y local, entre todas las secuencias.

Alineamiento: T-Coffee

•Asigna puntaje a los matches de esos alineamientos. En un proceso llamado extensión de la biblioteca, usa terceras secuencias para modificar los puntajes en un par dado.

•Al final tiene, sobre los alineamientos de a pares, puntajes que reflejan la compatibilidad con otros alineamientos.

•Y usa eso para hacer la prog. din., en alineamiento progresivo.

Alineamiento: T-Coffee

+ Presta más atención a homologías locales que Clustal.

+ Permite incluir información de otras fuentes! A la biblioteca primaria se le pueden agregar alineamientos obtenidos con Clustal, Dialign, información biológica, lo que sea.

Es de lo mejorcito que ha aparecido.

Alineamiento: POA

POA: Partial Order Alignment Idea: no resumir los alineamientos en perfiles, sino que dibujar un grafo. Y alinear vía programación dinámica, pero contra ese grafo.

Alineamiento: POA

La idea no es mala... Y permite conservar la información de un alineamiento, en una representación compacta.

Alineamiento: POA

Alínea permitiendo agregar caminos alternativos por sobre la matriz 2d.

Y ayuda a visualizar estructuras compartidas parcialmente:

Alineamiento: POA

Problemas: •Se vuelve poco práctico si hay demasiada bifurcación: con sets muy grandes de secuencias, o secuencias evolutivamente distantes.

•Con frecuencia lo que queremos es generalizar: discernir los motivos que caracterizan una familia de proteínas, o una secuencia de regulación en el DNA, etc.

Alineamiento: benchmarks

Como ya dijimos, los benchmarks aquí son humanos: sets curados a mano, a partir de amplia información filogenética, estructural, etc.

El más conocido es BaliBase: •Son “correctos” •Tiene sets “difíciles”, agrupados según naturaleza de la dificultad.

•Lo más difícil es donde hay mucha distancia evolutiva: la twilight zone!

http://www.people.virginia.edu/~wrp/papers/ismb2000.pdf

Alineamiento: benchmarks

BaliBase, Thompson et al, NAR, 1999, PROBLEMA Descripción 1. Grupo de secuencias cercanas 2. Una secuencia lejana 3. Dos grupos distantes 4. Inserción interna 5. Inserción en un extremo

Alineamiento: benchmarks

PROBLEMA Estrategia ClustalW, T-coffee, MSA T-Coffee T-Coffee Dialign T-Coffee Dialign T-Coffee

Alineamiento: benchmarks

...aunque, de otra fuente: Pero parece ser que: Dialign para cosas lejanas, T-Coffee para lo demás.

Como regla a primera vista, claro. Toda la fauna que anda afuera tiene razón de ser, según el caso. Y en último término, T-Coffee puede combinar.

Alineamiento: benchmarks

Lista de algunos:

MAFFT Progressive www.biophys.kyoto-u.jp/katoh POA Progressive/Simulataneous www.bioinformatics.ucla.edu/poa MUSCLE Progressive/Iterative www.drive5.com/muscle/

Alineamiento: benchmarks

En todo caso, la mayoría de los biólogos sigue usando Clustal.

Alineamiento: benchmarks

En todo caso, la mayoría de los biólogos sigue usando Clustal.

Volvamos a eso de generalizar. Extracción de motivos  patrones.

Motivos

Motivo ( motif ): un patrón que puede estar presente en una secuencia.

Puede ser más o menos estricto. Puede reflejar un sitio activo de una proteína, o un punto de interacción entre moléculas.

Motivos

Con los motivos se necesita: •Un lenguaje para describirlos: para saber que tipos de motivos se pueden encontrar.

•Una forma de evaluarlos: una función objetivo, que defina los motivos que son “interesantes”.

•Algoritmos para determinarlos, y para localizarlos.

  Determinarlos: extraer un motivo a partir de un conjunto de secuencias.

Localizarlos: dado un motivo, saber reconocerlo.

Motivos

Lenguaje para describirlos: Se usan varios.

•Expresiones regulares (no necesito explicárselas) •Perfiles (ya vistos, es la matriz de frecuencias x ubicación) •Modelos markovianos ocultos (HMM)

Cadenas de Markov

En una cadena de Markov, el sistema puede estar en una serie de posibles estados (en un momento dado), y entre un instante y el siguiente cambia de estado de acuerdo a una cierta matriz de probabilidad.

A B Componentes: •Set de estados.

•Una distribución inicial de probabilidad sobre los estados •La matriz de transición.

Cadenas de Markov

Estados: día soleado, lluvioso, nublado.

Matriz de transición: probabilidades para hoy, dado el estado de ayer Una distribución inicial.

Cadenas de Markov

Hay un montón de teoría sobre cadenas de Markov; se usan mucho. Aquí nos interesa un tipo un poco más complicado de proceso Markoviano: 1. El grafo será acíclico: no hay “vuelta atrás” en los estados.

2. En cada paso, no observamos el estado, sino un estado “de salida”, que depende del estado en que estemos, con cierta probabilidad.

Como no observamos directamente el estado en que estamos, se habla de modelos ocultos de Markov.

HMM

"Hidden Markov Model" (HMM).

Ejemplo (no acíclico): un casino tiene un dado bueno, y un dado cargado (con P(6)=½, el resto uniforme). Cambia del bueno al cargado con probabilidad 0.05, y del cargado al bueno con probabilidad 0.1.

0.95

0.05

0.90

Bueno 0.1

Cargado Eso sería como cadena de Markov. Pero supongamos que no sabemos que dado se está usando. Sólo vemos el resultado de las tiradas.

0.95

HMM

0.05

0.90

Bueno 0.1

Cargado HMM: la variable que observamos depende de estados internos que no vemos.

0.95

Si aparece ¿cuál es la secuencia de estados internos más probable?

312453666641 1:1/6 2:1/6 3:1/6 4:1/6 5:1/6 6:1/6 Bueno 0.05

0.10

1:1/10 2:1/10 3:1/10 4:1/10 5:1/10 6:1/2 Cargado 0.90

 BBBBBBCCCCCC

HMM

Esquema general de un HMM (hidden markov model): Componentes: •Set de estados ocultos

(aquí, x

“alfabeto”

1 , x 2 , x (aquí, y 1 3 , y )

•Set de estados de salida, o

2 , y 3 )

•Distribución inicial •Matriz de transición entre estados ocultos)

(P ij (Q )

•Matriz de output (prob. de estados de salida, dados los

ij )

HMM

Problemas usuales con HMM: •Dado un HMM y una secuencia observada, ¿cuál es la probabilidad de que la secuencia sea producida por ese HMM?

•Dado un HMM y una secuencia observada, ¿cuál es la secuencia más probable de estados internos que la genera?

•Dado un conjunto de secuencias observadas, ¿cuál es el mejor HMM?

Algoritmos: forward, Viterbi, Baum-Welch

Aplicando a secuencias genéticas: Los estados internos son las columnas “relevantes” de un alineamiento.

HMM

Si no hay gaps, entonces se reduce a un perfil: M 1 M 2 M 3 M 4

“Motif HMM”

M 5

HMM

Los “problemas usuales” se convierten en: •Dado un perfil HMM, esta secuencia, ¿es probable que pertenezca al club?

•Dado un perfil HMM, esta secuencia, ¿cómo la alíneo con las previas?

•Dado este conjunto de secuencias, ¿cómo lo resumo con un perfil tipo HMM?

Todas son importantes!

ACA - - - ATG TCA ACT ATC ACA C - - AGC AGA - - - ATC ACC G - - ATC inserción Prob. de output

HMM

Prob. de transición Sec. de consenso: ACAC - - ATC P (ACACATC) = 0.8 x 1 x 0.8 x 1 x 0.8 x 0.6 x 0.4 x 0.6 x 1 x 1 x 0.8 x 1 x 0.8 = 4.7 x 10 -2

HMM

Hay que distinguir entre la secuencia observada, y el camino que se usa para generarla (que no se ve!).

Para una secuencia observada, pueden haber varios caminos alternativos.

Si conocemos la secuencia observada a=a 1 ,...a

L , y también el camino, s=s 1 ,...,s L (supondremos que el carácter vacío es parte del alfabeto, para tener la misma longitud), entonces la probabilidad conjunta es

(

) 

(

1 ) 

  2 ...

L P

(

a i

O con la notación de antes,

s i

)

(

s i s i

 1 )

(

) 

Q s

1 ,

1 

  2 ...

L Q s i

a i P s i

 1 ,

s i

HMM

Pero si no conocemos los estados internos (s), entonces sólo podemos calcular

(

)  

s P

(

) donde sumamos sobre todos los s que pueden generar a.

Nota: se suele agregar un estado de inicio y uno de término, ambos “mudos” (o bien, con outputs especiales, “^” y “$”, por ejemplo).

^ ACG $ 0 1 2 3

Este es un caso fácil: cada a tiene un único s.

HMM

ACA - - - ATG TCA ACT ATC ACA C - - AGC AGA - - - ATC ACC G - - ATC P(ACA -- ATG) = 3.3x10

-2 P(TCA ACT ATC) = 0.0075x10

-2 P(ACA C- AGC) = 1.2x10

-2 P(AGA -- ATC) = 3.3x10

-2 P(ACC G- ATC) = 0.59x10

-2 4.9

3.0

5.3

4.9

4.6

Consenso: P(ACA C- ATC) = 4.7x10

-2 P(ACA -- ATC) = 13.1x10

-2 6.7

6.3

Excepcional: P(TGC T- AGG) = 0.0023x10

-2 -0.97

Usamos log odds para tener cifras más fáciles. Odds: prob. / prob. del caso nulo.

Nulo en este caso: secuencia aleatoria de largo L  P=(0.25) L Log Odds(a) = log(P(a)/(0.25) L ) Valor positivo  Más probable que modelo nulo!

HMM: Viterbi

Dada a y un HMM M, ¿cuál es la s más probable?

Respuesta: algoritmo de Viterbi.

 Programación dinámica, similar a Smith-Waterman.

 También rellenamos una matriz.

 En un eje ponemos la secuencia observada, en el otro los estados ocultos.

 En la posición (x,y) la matriz contiene la máxima probabilidad de un camino que termine en el estado oculto x-ésimo con la observación a y .(pasando por a 1 ,...a

y-1 ).

HMM: Viterbi

^ ACG

( 0 , 0 )  1 ,

( 0 ,

)  0 

 0

(

 1 ) 

(

 1 0 1 2

)  max

[

(

)

(



)] ...y vamos guardando el argmax, para poder reconstruir la ruta al llegar al final, desde (3,$).

$ 3

x y

0 (inicio) 1 2 3 (fin)

0 0 1 0 Secuencia observada A C A C G V(2,c) $

HMM: Forward

Dada a, y un HMM M, ¿cuánto vale P( a/M )?

Respuesta: algoritmo Forward.

 Idéntico a Viterbi, pero reemplazamos “max” por “sum”.

 En la posición (x,y) la matriz contiene la suma de las probabilidades de los caminos que terminan en el estado oculto x-ésimo con la observación a y .(pasando por a 1 ,...a

y-1 ).

( 0 , 0 )  1 ,

( 0 ,

)

(

  1 ) 0  

k P

(

  1 0

)  

w F

(

)

(



)

HMM: Decodificación posterior

Posterior decoding contesta a la pregunta: dada a y un HMM M que generó a, ¿cuál es la probabilidad de que el estado oculto x haya generado la letra a y ?

Esto sirve, por ejemplo, para saber cómo se corresponden los aminoácidos de una proteína con posiciones estructurales conservadas.

Respuesta: 1) Se hace un algoritmo Backward, análogo a Forward. La variable B(x,y) será la suma de probalidades de caminos que parten en el estado x con la observación a y .

2) La respuesta a la pregunta será

(

)

(

)

HMM: Construcción

Queda la pregunta del millón: ¿cómo construímos el HMM?

Hay dos problemas: •Determinar la arquitectura (estados ocultos y sus conexiones) •Determinar los parámetros (probabilidades de transición y de output) En el caso más general, la arquitectura es un arte.

Para alineamientos, se puede usar esta forma estándar. ¿Cuántos nodos “M”?   La cantidad de columnas mayoritariamente “no gap” o el promedio de longitudes

HMM: Construcción

A veces las secuencias ya están alineadas, y queremos un HMM que represente el alineamiento.

En ese caso, si las secuencias no son muchas, puede hacerse “a mano”: contamos las frecuencias de letras, vemos los desvíos... . Donde haya gaps, ponemos insert o delete según la cantidad de gaps (>50%  insert).

La corrección que se suele hacer, para generalizar: agregar probabilidades pequeñas de emisión para letras ausentes.

(Pondero las frecuencias observadas con las frecuencias “a priori”, dando un cierto peso a cada parte).

ACA - - - ATG TCA ACT ATC ACA C - - AGC AGA - - - ATC ACC G - - ATC

HMM: Baum-Welch

Si las secuencias no están alineadas, entonces es más complicado. Recibo un conjunto de secuencias, {a 1 ,..a

N }.

Si supiera cuál es el camino que genera cada una, podría estimar las probabilidades de transición y de emisión contando la frecuencia con que ocurren. Notemos A ij la cantidad de veces que usan la transición i  j, y B ij la cantidad de veces que la letra j se emite en el estado i. Entonces podría estimar ˆ

 

k A i A i

 

k B i B i

HMM: Baum-Welch

Por otro lado, si tuviera los parámetros de la HMM, entonces podría saber la probabilidad de los caminos que generan las secuencias, y por lo tanto la frecuencia con que usarán cada emisión y cada transición. En particular, la probabilidad de que se use una transición de x a w sería

 1

sería

k N

  1

( 1

a k

) largo (

  1

a k F

) (

)

(



)

(

 1

)

(

Y la probabilidad de emitir una letra e desde el estado x ˆ

 1

N k N

  1

( 1

a k

) largo (

  1 ,

a k F

) (

)

(

)

a k y



 1 )

HMM: Baum-Welch

Entonces, si tuviera los parámetros, sabría las probabilidades de los caminos, y si tuviera los caminos, podría estimar los parámetros.

E-step: estimar datos faltantes  algoritmo EM “Expectation Maximization” •Se aplica cuando desconocemos parámetros de un modelo que genera ciertos datos, y tampoco conocemos una parte de eso que se generó.

Parámetros Información faltante M-step: estimar parámetros (maximizando verosimilitud)

HMM: Baum-Welch

En este caso, los parámetros desconocidos son los del HMM, y de lo que se genera, conocemos la secuencia pero no el camino.

•Partimos con alguna “adivinanza” (posiblemente uniforme) sobre los parámetros.

E-step: estimar datos faltantes: Parámetros: P y Q ˆ

, Información faltante: Caminos para {a k } ˆ

•Iteramos hasta converger.

M-step: estimar parámetros: ˆ

, ˆ

HMM: Construcción

Entonces, para construir un HMM con secuencias no alineadas: •Definir arquitectura (longitud); por lo general se usa la longitud promedio.

•Aplicar Baum-Welch. Ojalá con varias condiciones de inicio distintas, para evitar óptimos locales.

•Una vez obtenido el HMM, alinear las secuencias encontrando para c/u su alineación con el HMM (vía Viterbi).

HMM presentes localmente

En principio, una secuencia se alínea globalmente contra un HMM. Pero podemos permitir alineamientos locales!

Q  emite letras no alineadas.

Viterbi sigue funcionando.

HMM

Sirven para: •Representar un perfil, de un alineamiento nuevo, o una familia de secuencias o subsecuencias (lugares en DNA, familias de proteínas, dominios funcionales en proteínas...)  Se puede pensar en los HMM como secuencias “generalizadas”, con probabilidades para las variaciones respecto al patrón.

•Crear alineamientos (usando Baum-Welch). Pero sólo cuando hay hartas secuencias; tienen demasiados parámetros para ser entrenados con poca información.

HMM

Sirven para: •Alinear una secuencia nueva c/r al HMM (usar Viterbi).

•Clasificar una secuencia nueva: usando Forward evaluamos la probabilidad de que pertenezca a distintas familias (cada una caracterizada por su HMM).

•Buscar nuevos miembros (idem, Forward, pero varío la query: es el HMM, no la secuencia).

HMM: software

•SAM (Sequence Alignment and Modeling System): el nombre no menciona HMM, pero se trata de eso.

•Meta-MEME: “is a software toolkit for building and

using motif-based hidden Markov models of

biological sequences”. La parte “MEME” es “Multiple EM for Motif Elicitation”.

•HMMER: búsqueda en base de datos mediante HMM.

HMM: datos

Pfam: base de datos de dominios y familias de proteínas.

•Un MSA curado a mano se usa como “semilla”  Pfam-A .

•Con eso se crea un HMM.

•Se usa el HMM para agregar otros miembros de manera automática (el umbral de puntaje se escoge a mano)  Pfam-B .

•Contiene más de 9000 familias; el 75% de las proteínas conocidas tienen al menos un match en Pfam.

+ datos

Bases de datos relacionadas: hay varias que también agrupan dominios y familias.

•ProSite: Curada a mano, altamente confiable. Usa expresiones regulares, y ahora también perfiles. •ProDom: Usa secuencias de consenso para la representación. Generada haciendo clustering sobre secuencias de SwissProt y TrEmbl.

•InterPro: integra a varias (Pfam, ProSite, ProDom...)

Algo más sobre motivos

•Hay softwares que se dedican a buscar motivos conocidos en una secuencia (consultan ProSite, Pfam, etc.).

•Otro problema es encontrar un motivo nuevo en un conjunto de secuencias dado. Una opción es hacer un alineamiento múltiple, y ver las regiones mejor alineadas.

•Eso falla si las secuencias sólo comparten el motivo, y lo tienen en lugares distintos. Entonces podríamos pensar en algoritmos locales, como Dialign.

•Hay algoritmos más específicos, vía EM, y Gibbs sampling.

Idea del Gibbs sampling

•N secuencias, largos posiblemente distintos.

•Una longitud L para el motivo desconocido.

•Una adivinanza inicial, con una ubicación para el motivo en cada secuencia. Eso da un MSA (aunque pueda ser pésimo).

Itero los siguientes pasos: •elimino una de las ubicaciones •escojo una nueva ubicación en esa secuencia, de manera probabilista

Idea del Gibbs sampling

Para escoger la nueva ubicación, le asigno a cada posible nueva ubicación una probabilidad que refleja qué tan bien esa ubicación le vendría al MSA actual de las otras N-1 ubicaciones:  es el producto, sobre i=1..L, de la frecuencia de la letra candidata en esa posición del MSA, dividido por la frecuencia “background”.

Puede agregarse un paso que haga variar L, si conviene.

Idea del Gibbs sampling

5’- TCTCTCTCCACGGCTAATTAGGTGATCATGAAAAAA TGAAAAATTC ATGAGAAAAGAGTCA GACATCGAAA CATACAT

…HIS7 …ARO4

5’- CACATCCAACGAATCACCTCACCGTTATCGTGACTCACTTTCTTTCGCATC GCCGAAGTGC CATAAAAAATATTTTTT 5’- TGCGAACAAAA GAGTCATTAC AACGAGGAAATAGAAGAAAATGAAAAATTTTCGACAAAATGTATAGTCATTTCTATC

…ILV6 …THR4

5’- ACAAAGGTACCTTCCTGGCCAATCTCACAGATTTAATATA

GTAAATTGTC

ATGCATATGACTCATCCCGAACATGAAA 5’- ATTGATTGACTCATTTTCCTCTGACTACTACCAGTTCAAAATGTTAGAGAAAAATAGAAAAGCAGAAAAAATAAATAA 5’- GGCG

CCACAGTCCG

CGTTTGGTTATCCGGCTGACTCATTCTGACTCTTTTTTGGAAAGTGTGGCAT

GTGCTTCACA

…ARO1 …HOM2 …PRO3

**********

(En este ejemplo se le da aún más libertad a las ubicaciones: pueden estar en cualquier hebra, y en cualquier secuencia –de modo que puede haber dos en una secuencia, y ninguna en otra.)

Idea del Gibbs sampling

5’-

Agrego?

TCTCTCTCCA

CGGCTAATTAGGTGATCATGAAAAAATGAAAAATTCATGAGAAAAGAGTCA

GACATCGAAA

CATACAT 5’- ATGGCAGAATCACTTTAAAACGTGGCCCCACCCGCTGCACCCTGTGCATTTTGTACGTTACTGCGAAATGACTCAACG 5’- CACATCCAACGAATCACCTCACCGTTATCGTGACTCACTTTCTTTCGCATC

GCCGAAGTGC

CATAAAAAATATTTTTT 5’- TGCGAACAAAA

GAGTCATTAC

AACGAGGAAATAGAAGAAAATGAAAAATTTTCGACAAAATGTATAGTCATTTCTATC

…HIS7 …ARO4 …ILV6 …THR4

5’- ACAAAGGTACCTTCCTGGCCAATCTCACAGATTTAATATA

GTAAATTGTC

ATGCATATGACTCATCCCGAACATGAAA 5’- ATTGATTGACTCATTTTCCTCTGACTACTACCAGTTCAAAATGTTAGAGAAAAATAGAAAAGCAGAAAAAATAAATAA 5’- GGCG

CCACAGTCCG

CGTTTGGTTATCCGGCTGACTCATTCTGACTCTTTTTTGGAAAGTGTGGCAT

GTGCTTCACA

…ARO1 …HOM2 …PRO3

TGAAAAATTC GACATCGAAA GCACTTCGGC GAGTCATTAC GTAAATTGTC CCACAGTCCG TGTGAAGCAC Mejora?

TCTCTCTCCA GACATCGAAA GCACTTCGGC GAGTCATTAC GTAAATTGTC CCACAGTCCG TGTGAAGCAC

Idea del Gibbs sampling

5’- TCTCTCTCCACGGCTAATTAGGTGATCATGAAAAAATGAAAAATTCATGAG

AAAAGAGTCA

GACATCGAAACATACAT

…HIS7

5’- ATGGCAGAATCACTTTAAAACGTGGCCCCACCCGCTGCACCCTGTGCATTTTGTACGTTACTGCG

AAATGACTCA

ACG

…ARO4

5’- CACATCCAACGAATCACCTCACCGTTATCG

TGACTCACTT

TCTTTCGCATCGCCGAAGTGCCATAAAAAATATTTTTT 5’- TGCGAAC

AAAAGAGTCA

TTACAACGAGGAAATAGAAGAAAATGAAAAATTTTCGACAAAATGTATAGTCATTTCTATC

…ILV6 …THR4

5’- ACAAAGGTACCTTCCTGGCCAATCTCACAGATTTAATATAGTAAATTGTCATGCATA

TGACTCATCC

CGAACATGAAA 5’- ATTGAT

TGACTCATTT

TCCTCTGACTACTACCAGTTCAAAATGTTAGAGAAAAATAGAAAAGCAGAAAAAATAAATAA 5’- GGCGCCACAGTCCGCGTTTGGTTATCCGGC

TGACTCATTCTGACTCTTTT

TTGGAAAGTGTGGCATGTGCTTCACACA

…ARO1 …HOM2 …PRO3

AAAAGAGTCA AAATGACTCA AAGTGAGTCA AAAAGAGTCA GGATGAGTCA AAATGAGTCA GAATGAGTCA AAAAGAGTCA

(Es el motivo de un sitio de regulación en genes de S. cerevisiae)

Gibbs sampling, EM

The Gibbs Motif Sampler Homepage: http://bayesweb.wadsworth.org/gibbs/gibbs.html

 software bajable, y también servicio online Otro servidor: http://www.cbs.dtu.dk/biotools/EasyGibbs/

Idea con EM

Los métodos EM no escogen ubicación ni hacen un MSA: •su parámetro a iterar es un motivo •le asignan una probabilidad a cada ubicación, según ese motivo •recalculan el motivo como un perfil de todas las ubicaciones, ponderadas por esas probabilidades Hay diversas variaciones sobre estas ideas (EM y Gibbs).

Dominios, motivos, familias

•Motivos: patrones locales. Pueden ser sitios activos de proteínas, dominios funcionales, regiones reconocibles de DNA, etc.

•Dominios: los segmentos funcionales en las proteínas (que son relativamente modulares). Un dominio suele estar caracterizado por un motivo.

Dominios, motivos, familias

•Familia: una familia de proteínas comparte la mayoría de sus dominios; puede haber variaciones. Y dentro de “superfamilias” de proteínas, más variación aún.

 Los motivos salen de MSA, pero también son útiles para construir MSA, si detectamos primero motivos conocidos en las secuencias, y los usamos para anclar el alineamiento.

Proteínas con múltiples dominios

Alinear bien proteínas multi-dominio es un problema. Especialmente si lo que queremos es no sólo identificar dominios, sino relación filogenética.

Proteínas con múltiples dominios

En este caso, las proteínas alineadas en A tienen ancestro común, pero no las alineadas en B. Sin embargo las de B resultan similares por culpa de los dominios IG compartidos (IG es un dominio “promiscuo”, que suele insertarse en proteínas). ¿Cómo distinguir?

Proteínas con múltiples dominios

Una aproximación:

Sequence Similarity Network Reveals Common Ancestry of Multidomain Proteins

N. Song et al, PLoS Comput Biol. 2008 4(5): e1000063.

10.1371/journal.pcbi.1000063

Estudian el grafo de las secuencias similares, dentro de una base de datos, a dos secuencias dadas. Según la forma de ese grafo, se distingue la ascendencia común vs. la homología sólo de un trozo.

Proteínas con múltiples dominios

Los vecinos comunes son más que los vecinos únicos, en el caso de verdadera homología de las secuencias.

En el caso en que sólo hay homología de algunos dominios, resulta al revés.

Definen un coeficiente de “correlación de vecindades”, que discrimina bien entre homólogos verdaderos y los que no.