M´aster en T´ecnicas Estad´ısticas - Proyecto Fin de M´aster

Cadenas de Markov en la Investigaci´ on del Genoma Autora: Sara Prada Alonso ´ Directora de proyecto: Mar´ıa de los Angeles Casares de Cal

8 de Julio de 2013

Agradecimientos En primer lugar, me gustar´ıa expresar mi m´as sincero agradecimiento a mi di´ rectora de proyecto, Mar´ıa de los Angeles Casares de Cal, por haber confiado en m´ı como alumna para realizar el presente trabajo, y haberme dado todas las facilidades para poder desarrollarlo a lo largo de estos meses. Agradezco encarecidamente su ayuda, inter´es, dedicaci´on y labores de supervisi´on. Del mismo modo, quiero darle las gracias a Antonio G´omez Tato por ayudarme a comprender mejor la materia de estudio y sus aplicaciones, especialmente en el campo computacional. Le agradezco igualmente y de forma encarecida su labor de supervisi´on y soporte en la b´ usqueda de cierta informaci´on. Todos las materias cursadas durante el M´aster me han aportado los conocimientos necesarios para comprender mejor y m´as r´apidamente nuevos temas a abordar, como el que presento en esta memoria. Gracias a mis profesores y a la formaci´on personal y profesional que me han aportado. Por u ´ltimo, agradecer a mi familia y amigos su apoyo incondicional en este tiempo de trabajo.

i

ii

AGRADECIMIENTOS

Resumen Los modelos de Markov ocultos est´an siendo utilizados actualmente para modelizar familias de prote´ınas en la b´ usqueda autom´atica de genes, y en el alineamiento m´ ultiple de secuencias, entre otras aplicaciones. Existen varios algoritmos para estimar los par´ametros de estos modelos, a partir de los datos. Se trata de hacer una revisi´on de los mismos y aplicarlo a datos reales.

iii

iv

RESUMEN

´Indice general 1. Introducci´ on general

1

2. Nociones b´ asicas de biolog´ıa

5

3. Introducci´ on a las Cadenas de Markov 3.1. Introducci´on a las Cadenas de Markov Finitas . . . . 3.2. Probabilidades de transici´on y Matriz de transici´on . 3.3. Distribuci´on inicial . . . . . . . . . . . . . . . . . . . 3.4. Distribuci´on de probabilidad en la etapa n-´esima . . . 3.5. Orden de una Cadena de Markov . . . . . . . . . . . 3.6. Clasificaci´on de estados . . . . . . . . . . . . . . . . . 3.7. Distribuciones estacionarias . . . . . . . . . . . . . . 3.8. Ejemplo: una cadena de Markov como modelo para el 3.8.1. Ejemplo . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . ADN . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

13 13 14 16 16 17 17 19 20 20

4. An´ alisis de las secuencias del ADN 4.1. Contrastes de independencia . . . . . . . . . . . 4.2. Contraste de bondad de ajuste . . . . . . . . . . 4.3. Modelizaci´on de “se˜ nales” en el ADN . . . . . . 4.3.1. Matrices de Peso. Independencia . . . . 4.3.2. Dependencias de Markov . . . . . . . . . 4.4. An´alisis de Patrones . . . . . . . . . . . . . . . 4.5. B´ usqueda de repeticiones en secuencias de ADN 4.5.1. Repeticiones en tandem . . . . . . . . . 4.5.2. Repeticiones dispersas . . . . . . . . . . 4.5.3. Estudio de repeticiones . . . . . . . . . . 4.5.4. El problema de la detecci´on de patrones 4.6. Aplicaci´on: las islas CpG . . . . . . . . . . . . . 5. Los 5.1. 5.2. 5.3.

modelos de Markov ocultos Introducci´on . . . . . . . . . . . El algoritmo de Viterbi . . . . . Aplicaciones: b´ usqueda de genes 5.3.1. GENSCAN . . . . . . . 5.4. Caso pr´actico: las islas CpG . .

. . . . .

v

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

23 24 28 29 30 31 33 34 35 35 35 36 37

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

43 43 49 54 56 59

vi

´INDICE GENERAL 5.5. Otras aplicaciones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 5.5.1. Modelos de Markov ocultos en Bioinform´atica: estudio de las prote´ınas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 5.5.2. Modelos de Markov ocultos en filogenia molecular . . . . . . . 70

Cap´ıtulo 1 Introducci´ on general La aparici´on de los secuenciadores autom´aticos y sus desarrollos posteriores ha puesto a disposici´on de los investigadores ingentes cantidades de secuencias de ADN que necesitan ser analizadas. Una de las primeras tareas consiste en identificar las diferentes “regiones” de la secuencia y su anotaci´on posterior (estimar su funci´on en comparaci´on con lo encontrado en otras especies). Es conocido que la mayor parte del genoma humano es “no codificante” (no da lugar a un gen) aunque est´a cada vez m´as claro que esas regiones, anta˜ no denominadas “ADN basura”, contienen informaci´on relevante y son determinantes a la hora de comprender el comportamiento de la c´elula. La Bioinform´atica se encarga, entre otras cosas, del dise˜ no de herramientas computacionales (utilizando t´ecnicas matem´aticas) para el discernimiento de la estructura del ADN y la identificaci´on de las diferentes partes del mismo. Los m´etodos empleados son muy diversos, pero entre ellos destaca el uso de modelos probabil´ısticos m´as o menos complejos para intentar encontrar y clasificar las diferentes regiones, como pueden ser genes, pseudogenes, micro ARNs, regiones promotoras, intrones, exones, regiones CpG, etc. Entre los modelos probabil´ısticos m´as usados est´an los modelos de Markov ocultos, originalmente introducidos en el an´alisis de textos, pero que est´an siendo utilizados con ´exito en Bioinform´atica y Filogenia Molecular, por citar dos disciplinas relacionadas o que utilizan informaci´on de secuencias de ADN. La mol´ecula de ADN se codifica por una secuencia de letras (A, C, G, T) que resume su estructura qu´ımica. Un primer modelo probabil´ıstico consistir´ıa en utilizar la distribuci´on multinomial, es decir, que la probabilidad de aparici´on de una secuencia de longitud n seguir´ıa una distribuci´on multinomial con par´ametros n, pA , pC , pG y pT (pA + pC + pG + pT = 1), lo cual significar´ıa que la aparici´on de un nucle´otido en una posici´on determinada ser´ıa independiente de las dem´as. Pero ´este no ser´a un buen modelo, ya que no tiene en cuenta la informaci´on que aportan los nucle´otidos adyacentes. Para incorporar esa informaci´on, puede utilizarse un modelo de Markov que determine la probabilidad de aparici´on de un nucle´otido en una posici´on i + 1 en funci´on de lo que ocurre en la posici´on i de la cadena (o incluso en 1

2

´ GENERAL CAP´ITULO 1. INTRODUCCION

las posiciones i − k, i − (k + 1), . . . , i. Este modelo funcionar´ıa si toda la secuencia se comportase de manera semejante, pero no es as´ı; existen “regiones” diferentes a lo largo del genoma, por ejemplo, regiones ricas en dinucle´otidos CG, llamadas islas CpG. En esos casos, la probabilidad de aparici´on de un nucle´otido en una posici´on (sitio) i de la cadena ya no s´olo depende del nucle´otido o nucle´otidos anteriores sino que tambi´en influye el “estado” o la regi´on que lo precede. Esos estados est´an ocultos y de ah´ı la necesidad de modelizar o utilizar modelos de probabilidad como los “modelos de Markov ocultos”. Desde sus primeras aplicaciones hasta la actualidad, los modelos de Markov ocultos (HMM: Hidden Markov Models) han ido aumentando su complejidad. Veremos someramente uno de los que se utilizan para la determinaci´on o b´ usqueda autom´atica de genes mediante uno de los programas m´as utilizados, GENSCAN, cuyo estudio queda fuera del objetivo de esta memoria. Como ejemplo de aplicaci´on de los modelos de Markov ocultos estudiaremos el problema de la determinaci´on de las regiones ricas en CG (islas CpG). Debido al proceso biol´ogico de la metilaci´on (adici´on de un grupo metilo (-CH3) a una mol´ecula), existe una probabilidad bastante alta de que esta metilaci´on de C mute en una T , con la consecuencia de que, en general, los dinucle´otidos CpG son m´as raros en el genoma de lo que podr´ıa esperarse de las probabilidades independientes de C y G. Debido a importantes razones biol´ogicas, el proceso de metilaci´on se suprime en determinadas extensiones cortas del genoma, como alrededor de los “promotores” o regiones de inicio de muchos genes. En estas regiones, observamos m´as dinucle´otidos CpG que en ning´ un otro lugar, y de hecho m´as nucle´otidos C y G en general. Tales regiones se denominan islas CpG y, en contraste, el resto del genoma es el “oc´eano”. Normalmente est´an compuestas por cientos de bases de longitud. Una primera aproximaci´on al problema anterior consiste en ver si una regi´on determinada es o no una isla CpG. Para ello, basta utilizar un test de raz´on de verosimilitudes. Una etapa posterior consistir´a en utilizar modelos de Markov ocultos muy sencillos, como puede ser el del “casino deshonesto”, que no tiene en cuenta la informaci´on de los nucle´otidos precedentes. Existen modelos de Markov ocultos sofisticados tambi´en para atacar el problema, uno de ellos ha aparecido recientemente en el art´ıculo de Irizarry et al. ([24]). Nosotros nos detendremos en la aplicaci´on del modelo de Markov oculto m´as utilizado actualmente para estudiar el tema. Una vez escogido el modelo y estimados sus par´ametros mediante una secuencia o secuencias “de entrenamiento”, se necesita un algoritmo para aplicarlo y determinar las diferentes regiones de una nueva secuencia. Uno de los algoritmos m´as utilizados es el algoritmo de Viterbi, que encuentra el “camino o´ptimo” mediante t´ecnicas de programaci´on din´amica. Como ejemplo, hemos estudiado el uso de los modelos de Markov ocultos en el cromosoma 10 del genoma del “Danio rerio” o “pez cebra” (uno de los “genomas-modelo” estudiados en la actualidad). Hemos descargado la secuen-

3 cia del cromosoma 10 as´ı como las regiones ya reconocidas como islas, tal como est´an en las bases de datos dedicadas a estas a´reas; y utilizando el algoritmo de Viterbi, programado en “R”, hemos estimado la aparici´on de regiones “CpG”. Terminamos esta memoria con un par de notas sobre algunos de los casos m´as recientes del uso en Bioinform´atica de los modelos de Markov ocultos.

4

´ GENERAL CAP´ITULO 1. INTRODUCCION

Cap´ıtulo 2 Nociones b´ asicas de biolog´ıa La mol´ecula de ADN de cada organismo est´a formada por dos hebras de nucle´otidos. Existen cuatro tipos distintos de nucle´otidos, cuyas bases son: Adenina, Timina, Guanina y Citosina. La Adenina s´olo puede unirse a la Timina y la Guanina a la Citosina, por eso decimos que la Adenina es complementaria de la Timina; y la Citosina es complementaria de la Guanina. Cada nucle´otido de una cadena est´a unido al nucle´otido que se encuentra enfrente en la otra cadena, y las dos hebras est´an plegadas formando una “doble h´elice”. Los cuatro nucle´otidos descritos forman largas secuencias de un “alfabeto” de cuatro letras A, T , G y C (para Adenina, Timina, Guanina y Citosina; respectivamente), llamado c´odigo gen´etico. Estas secuencias experimentan cambios dentro de una poblaci´on a lo largo de varias generaciones, como las mutaciones aleatorias que surgen y se mantienen en la poblaci´on. Por lo tanto, dos secuencias bastante diferentes pueden derivar de un antepasado com´ un. Supongamos que tenemos dos peque˜ nas secuencias de ADN como las mostradas a continuaci´on, quiz´as de dos especies distintas. Las letras en negrita indican parejas de nucle´otidos que son el mismo en ambas secuencias. GGAGACT GT A GAC AGC T AAT GC T A TA GAAC GCC C T A GC C AC GAGC C C T T A TC Deseamos medir si las dos secuencias muestran una similitud significativa, para evaluar, por ejemplo, si provienen de un ascendiente com´ un. Supongamos una variable aleatoria X que mida el n´ umero de coincidencias en un conjunto de n = 26 observaciones, o letras de nuestro alfabeto en este caso. Si las secuencias son generadas aleatoriamente, teniendo las cuatro letras A, T , G y C igual probabilidad p = 41 de ocurrir en cualquier posici´on, y consideramos que una coincidencia entre observaciones de ambas secuencias es un “´exito”, tendr´ıamos que X es una variable aleatoria que sigue una distribuci´on Binomial de par´ametros n = 26 y p = 41 .

5

´ CAP´ITULO 2. NOCIONES BASICAS DE BIOLOG´IA

6

X v B(n, p) As´ı, las dos secuencias deber´ıan coincidir en un cuarto de las posiciones. Las dos secuencias anteriores coinciden en 11 de 26 posiciones, ¿c´omo de improbable es este resultado si las secuencias son generadas aleatoriamente? La teor´ıa de la probabilidad junto con lo expuesto en el p´arrafo anterior muestran que, asumiendo probabilidades iguales para A, T , G y C en cualquier posici´on, e independencia de todos los nucle´otidos involucrados; la probabilidad de que haya 11 o m´as coincidencias en una comparaci´on de secuencias de longitud 26 es aproximadamente de 0.04, tal y como se indica a continuaci´on P (X ≥ 11) = 1 − P (X < 11) = = 1 − (P (X = 0) + P (X = 1) + . . . + P (X = 10)) = 0.0400845, siendo 26−i    i  1 26 1 1− , P (X = i) = 4 4 i

i = 0, 1, 2, . . . , 26.

Por lo tanto, nuestra observaci´on de 11 coincidencias muestra evidencias significativas de que esto no s´olo es cosa del azar. Estudiaremos el modelo que se encuentra bajo las secuencias de nucle´otidos, y que da lugar con una cierta probabilidad a unas u otras cadenas. A ´este lo llamaremos Modelo de Markov Oculto. Mostramos a continuaci´on con m´as detalles algunas nociones b´asicas de biolog´ıa que se necesitar´an para comprender mejor las explicaciones posteriores ([40]). El ´acido desoxirribonucleico (ADN) es la informaci´on macromolecular b´asica de la vida. Consiste en un pol´ımero de nucle´otidos, en el cual cada nucle´otido est´a compuesto por un az´ ucar (la desoxirribosa) y un grupo fosfato, conectado a una base nitrogenada de uno de los cuatro posibles tipos: adenina, guanina, citosina o timina (abreviadas como A, G, C y T , respectivamente). Los nucle´otidos adyacentes en una sola hebra de ADN est´an conectados por un enlace qu´ımico entre el az´ ucar de uno de ellos y el grupo fosfato del siguiente. La estructura cl´asica de “doble h´elice” del ADN se forma cuando dos hebras de ADN forman enlaces hidrogenados entre sus bases nitrogenadas, resultando en la ya familiar estructura de escalera. Bajo condiciones normales, estos enlaces hidrogenados se forman solo entre particulares pares de nucle´otidos (en referencia a una de las parejas de bases): adenina se empareja solamente con timina, y guanina con citosina. Dos hilos de ADN son complementarios si la secuencia de bases en cada uno es tal que se emparejan correctamente a lo largo de toda la longitud de ambas cadenas. La secuencia en la cual las diferentes bases ocurran en una particular hebra de ADN representa la informaci´on gen´etica cifrada en esa hebra. En virtud de la especificidad del emparejamiento de nucle´otidos, cada una de las dos hebras de cualquier mol´ecula de ADN contiene toda la informaci´on presente en la otra. Hay tambi´en una polaridad qu´ımica en las cadenas

7 de polinucle´otidos de modo que la informaci´on contenida en las bases A, G, C y T es sintetizada y decodificada en cualquier direcci´on. La mol´ecula de ADN puede ser circular o lineal, y puede estar compuesta de 10,000 mil millones de nucle´otidos en una cadena larga.

Figura 2.1: Situaci´on del ADN dentro de la c´elula. En la c´elula, el ADN se organiza en cromosomas, que se encuentran en el n´ ucleo de ´esta. Se trata de una sola pieza de espiral de ADN que contiene muchos genes, elementos reguladores y otras secuencias de nucle´otidos. Los cromosomas var´ıan ampliamente entre los diferentes organismos. Las c´elulas humanas contienen 23 pares de cromosomas, un miembro de cada par heredado de la madre y el otro del padre. Los dos cromosomas en un par son pr´acticamente id´enticos, con la excepci´on del cromosoma sexual, para el cual hay dos tipos, X e Y (XX para la mujer y XY para el hombre, siendo la pareja la que determina el sexo). Cada c´elula del cuerpo contiene copias id´enticas del conjunto entero de 23 pares de cromosomas. Como curiosidad, se˜ nalar que si se desenrollaran y pusieran en fila los cromosomas en cada una de las c´elulas (un humano adulto tiene entre 10 y 50 billones de c´elulas) la longitud total de ADN ser´ıa de unos 2 metros. Si se sumara la longitud del ADN de todas las c´elulas de una sola persona, se podr´ıa rodear la circunferencia terrestre 500,000 veces. El conjunto total del ADN de un organismo es el genoma; que contiene m´as de tres billones de pares de bases. El genoma humano o secuencia completa de ADN de un organismo constituye la informaci´on gen´etica heredable del n´ ucleo celular. Un gen es un trozo de ADN que lleva la informaci´on para que se fabrique una prote´ına, que puede hacerse de m´ ultiples formas. Los genes est´an en el n´ ucleo de las c´elulas y las prote´ınas que codifican, que son las que controlan los caracteres, se fabrican

8

´ CAP´ITULO 2. NOCIONES BASICAS DE BIOLOG´IA

en el citoplasma de la c´elula (parte entre en n´ ucleo y la membrana que delimita la c´elula). Existen, por tanto, mol´eculas intermediarias que pueden “copiar” un trozo de la cadena de ADN, atravesar la membrana del n´ ucleo y ya en el citoplasma traducir la informaci´on almacenada en el ADN. Son las mol´eculas de ARN mensajero. Un cromosoma humano est´a compuesto principalmente del llamado ADN basura o ADN interg´enico, cuya funci´on no est´a del todo comprendida. El ADN interg´enico es ADN que no codifica prote´ınas. Ejemplos de ´este son los intrones, que son segmentos internos dentro de los genes y que se eliminan a nivel de ARN; o los pseudogenes, inactivados por una inserci´on o supresi´on. Intercalados con estas ´areas ´ de ADN est´an los genes. Estos se organizan en exones, que son las secuencias que ser´an eventualmente utilizadas por la c´elula, alternados con los intrones, que ser´an desechados. Actualmente se piensa que el genoma humano contiene aproximadamente entre 30,000 y 40,000 genes. La informaci´on en estos genes se transcribe en ARN (´acido ribonucleico), y en muchos casos finalmente en prote´ınas.

Figura 2.2: Situaci´on del gen dentro del cromosoma. La estructura de los genes y el mecanismo de expresi´on de ´estos en organismos eucariotas es mucho m´as complicado que en los procariotas. En organismos eucariotas, la regi´on de ADN codificante en prote´ına no es continua normalmente. Esta regi´on est´a compuesta por extensiones de exones e intrones. Durante la transcripci´on, ambos son transcritos en ARN, en su orden lineal. Despu´es de eso, tiene lugar un proceso llamado “de empalme” o “de ayuste” (splicing), en el cual las secuencias de intrones se extraen y se descartan de la secuencia de ARN. Los segmentos restantes de ARN, los que corresponden a los exones, se ligan para formar la hebra madura de ARN. Un gen multi-ex´on t´ıpico tiene la siguiente estructura: empieza con la regi´on promotora, a la que le sigue una regi´on transcrita pero no codificante llamada

9 regi´on no traducida 5’ (5’UTR: 5’ Untranslated Region). Le sigue el ex´on inicial que contiene el cod´on de inicio (secuencia de grupos de tres nucle´otidos que marcan el inicio del ex´on). Siguiendo al ex´on inicial, existen series alternadas de intrones y exones internos, seguidos por el ex´on final, que contiene el cod´on de terminaci´on o final (secuencia de grupos tres nucle´otidos que “se˜ nalan” el final del ex´on). Despu´es, nos encontramos otra regi´on no codificante llamada regi´on no traducida 3’ (3’UTR). Terminando el gen eucariota, hay una se˜ nal de poliadenilaci´on (Poli(A): el nucle´otido Adenina se repite varias veces). Los l´ımites ex´on-intr´on (es decir, los lugares de empalme) est´an se˜ nalados por cortas secuencias espec´ıficas (2 pares de bases de longitud). Una 5’UTR (3’UTR) terminada en un intr´on (ex´on) se llama sitio donor (donor site), y una 3’UTR (5’UTR) que termine con un intr´on (ex´on) se llama sitio aceptor (acceptor site).

Figura 2.3: Estructura de un gen. El primer paso en el proceso de conversi´on del gen a prote´ına es la transcripci´on, la creaci´on de la mol´ecula de ARN utilizando la secuencia de ADN del gen como plantilla. La transcripci´on se inicia en las secuencias no codificantes o promotoras, localizadas inmediatamente antes del gen. Como el ADN, el ARN se compone de series de nucle´otidos, pero con varias diferencias importantes: el ARN est´a compuesto por una sola hebra, que contiene el “az´ ucar ribosa”, y sustitutos de la base nitrogenada Uracilo por Timina. El Uracilo es una de las cuatro bases nitrogenadas que forman parte del ARN y en el c´odigo gen´etico se representa con la letra U. Despu´es de la transcripci´on, que incluye la eliminaci´on de intrones, el ARN ir´a a diferentes

10

´ CAP´ITULO 2. NOCIONES BASICAS DE BIOLOG´IA

destinos dentro de la c´elula. De inter´es particular es el ARNm (ARN mensajero), que ser´a convertido en prote´ına. Una prote´ına est´a compuesta de una secuencia de amino´acidos, que se unir´an en el orden indicado en el “mensaje” que se encuentra en el ARN mensajero. Hay 20 amino´acidos que forman las prote´ınas. Cada uno de esos amino´acidos est´a representado por una o m´as secuencias de tres ARN nucle´otidos conocidos como un cod´on; por ejemplo, la secuencia de ARN AAG codifica el amino´acido lisina. La lisina es uno de los 10 amino´acidos esenciales para los seres humanos. La combinaci´on de cuatro posibles nucle´otidos en grupos de tres da lugar a V R4,3 = 43 o 64 codones, es decir, que la mayor´ıa de los amino´acidos est´an codificados por m´as de un cod´on. El ribosoma realiza el cambio del ARNm en prote´ına. Los ribosomas son complejos macromoleculares de prote´ınas y a´cido ribonucleico (ARN) que se encuentran en el citoplasma y se encargan de sintetizar prote´ınas a partir de la informaci´on gen´etica que les llega del ADN transcrita en forma de ARN mensajero (ARNm). El ribosoma empareja cada cod´on en la secuencia de ARN con el apropiado amino´acido, y despu´es a˜ nade el amino´acido sobre la prote´ına en creaci´on. El proceso de cambio est´a mediado por dos tipos especiales de cod´on: el cod´on de inicio se˜ nala la ubicaci´on en el ARN molecular donde la traducci´on (o s´ıntesis de prote´ınas) debe comenzar, mientras que los codones finales se˜ nalan el lugar donde la traducci´on debe terminar. Una vez que se compone la secuencia de amino´acidos, se “monta” una particular prote´ına, ´esta se disocia del ribosoma y se pliega en una espec´ıfica forma tridimensional. La funci´on de la prote´ına depende fundamentalmente de su estructura tridimensional y su secuencia de amino´acidos. Las prote´ınas pasan a llevar a cabo una variedad de funciones en la c´elula, cubriendo todos los aspectos de las funciones celulares del metabolismo. Actualmente, estas funciones celulares han sido asignadas a solamente una peque˜ na proporci´on de los genes, incluso en los modelos de organismos mejor conocidos. Para poder asignar funciones a los genes restantes, es u ´til examinar “expresiones patrones” (o repeticiones de secuencias de nucle´otidos) de estos genes en varios tejidos. La tecnolog´ıa “microarray” desarrollada en los a˜ nos pasados, permite la medida de los niveles de ARNm para decenas de miles de genes simult´aneamente. Esto proporciona una eficiente y buena manera de determinar la “expresi´on patr´on” (o secuencia repetida) de los genes en muchos tipos diferentes de tejidos, pero al mismo tiempo proporciona nuevos retos en informaci´on catalogada y an´alisis estad´ıstico.

11

Figura 2.4: Proceso de creaci´on de la prote´ına.

12

´ CAP´ITULO 2. NOCIONES BASICAS DE BIOLOG´IA

Cap´ıtulo 3 Introducci´ on a las Cadenas de Markov 3.1.

Introducci´ on a las Cadenas de Markov Finitas

En este cap´ıtulo abordaremos una breve introducci´on sobre procesos estoc´asticos, en particular sobre las cadenas de Markov en tiempo discreto y finitas, repasando sus caracter´ısticas y propiedades. El objetivo es introducir el material necesario para la comprensi´on de los modelos que desarrollaremos posteriormente. Un proceso estoc´astico es un conjunto de variables aleatorias (estoc´asticas) {Xt , t ∈ T }. El par´ametro o ´ındice t se interpreta generalmente como “el tiempo”, y los posibles valores de la v.a. Xt se interpretan como los posibles “estados” del proceso en el tiempo t. El conjunto de par´ametros T puede ser continuo o discreto. Si T es discreto diremos que el proceso es de par´ametro o de tiempo discreto. Si T es un intervalo de la recta real, entonces diremos que el proceso es de par´ametro o tiempo continuo. Cada una de las variables aleatorias del proceso tiene su propia funci´on de distribuci´on de probabilidad y, entre ellas, pueden estar correlacionadas o no. Cada variable o conjunto de variables sometidas a influencias o impactos aleatorios constituye un proceso estoc´astico. T´ıpicamente, describen alg´ un fen´omeno que evoluciona en el tiempo o en el espacio. Dentro de los procesos estoc´asticos, encontramos los procesos de Markov, que son aquellos procesos discretos en los que la evoluci´on en el tiempo (generaciones) o en el espacio (secuencias biol´ogicas) s´olo depende del estado actual y no de los anteriores. Un proceso de Markov, llamado as´ı por el matem´atico ruso Andr´ei Markov, es un fen´omeno aleatorio dependiente del tiempo para el cual se cumple una propiedad espec´ıfica: la propiedad de Markov. En una descripci´on com´ un, la condici´on de Markov, o “sin memoria”, exige que la propiedad de que la cadena de Markov que se encuentre en un estado j en el instante (n + 1) dependa u ´nicamente del estado 13

´ A LAS CADENAS DE MARKOV CAP´ITULO 3. INTRODUCCION

14

en que se encontraba en el instante n, y que esto se cumpla para cualquier etapa en la que se encuentre la cadena. Frecuentemente el t´ermino cadena de Markov se utiliza para dar a entender que un proceso de Markov tiene un espacio de estados discreto (finito o numerable), es decir, es un proceso de Markov en tiempo discreto. Las Cadenas de Markov son muy u ´tiles en Bioinform´atica, ya que se utilizan en la modelizaci´on de los procesos evolutivos de las cadenas de ADN. De manera formal, dir´ıamos que Definici´ on 3.1 Una Cadena de Markov es una sucesi´on de variables aleatorias {Xn , n ∈ N} que cumple: 1. Cada variable Xn toma valores en un conjunto finito o numerable E, que se denomina espacio de estados. 2. La sucesi´on de variables verifica la condici´on o propiedad de Markov o markoviana: P (Xn+1 = in+1 |X0 = i0 , . . . , Xn = in ) = P (Xn+1 = in+1 |Xn = in ),

(3.1)

donde i0 , . . . , in denotan los estados en los que se encuentra la cadena en cada etapa.

3.2.

Probabilidades de transici´ on y Matriz de transici´ on

Las probabilidades condicionadas P (Xn+1 = j|Xn = i),

(3.2)

son conocidas como las probabilidades de transici´on. La condici´on de homogeneidad en el tiempo, seg´ un la cual la probabilidad de pasar de i a j es independiente de la etapa en la que se encuentre la cadena, hace que la probabilidad de transici´on de i a j sea ∀n, m ∈ N ; i, j ∈ E (3.3) que indica la probabilidad de pasar al estado j desde el estado i, en cualquier etapa. Una Cadena de Markov en tiempo discreto (CMTD) que cumpla esta condici´on de homogeneidad en el tiempo (cadena homog´enea) se dice que tiene probabilidades de transici´on estacionarias. pij = P (Xn+1 = j|Xn = i) = P (Xm+1 = j|Xm = i),

Supongamos que el espacio de estados E que manejamos, es decir el conjunto de posibles valores que puede tomar el proceso en cada una de sus etapas, tiene un n´ umero finito k de estados. En este caso, decimos que la cadena de Markov es finita.

´ Y MATRIZ DE TRANSICION ´ 3.2. PROBABILIDADES DE TRANSICION As´ı, las probabilidades de transici´on entre estados, matriz de transici´on P de tama˜ no k × k. Dada la matriz es de la forma  p11 p12 · · ·  p21 p22 · · ·  P = (pij ) =  .. .. . .  . . . pk1 pk2 · · ·

15

pij , se combinan formando una homogeneidad de la cadena, la p1k p2k .. .

   . 

(3.4)

pkk

Los elementos de esta matriz P verifican que pij > 0, ∀i, j ∈ E. Adem´as, dado que el estado es i en el instante n yX que el proceso debe estar en alg´ un estado en el instante n + 1, se verifica que pij = 1, ∀i ∈ E; es decir, las sumas por j∈E

filas son uno. Esta propiedad es caracter´ıstica de las matrices estoc´asticas, a las que pertenecen las matrices de transici´on. Definici´ on 3.2 Sea P una matriz cuadrada k × k. Llamaremos a P matriz estoc´astica si 1. Todas sus entradas son positivas, pij ≥ 0, para todo i, j. 2. La suma de sus coeficientes en cada fila vale 1, es decir k X

pij = 1,

∀i = 1, . . . , k

(3.5)

j=1

y doblemente estoc´astica si adem´as la suma de sus coeficientes en cada columna tambi´en vale 1, o sea k X

pij = 1,

∀j = 1, . . . , k

(3.6)

i=1

Se tiene que las potencias de una matriz estoc´astica tambi´en lo son ya que, si U es el vector columna k × 1 cuyos elementos son todos unos, y P 0 es otra matriz estoc´astica, entonces 1. P P 0 tiene todos los elementos ≥ 0, pues ambas los tienen por definici´on de matriz estoc´astica. 2. (P P 0 )U = P (P 0 U ) = P U = U , la suma de los elementos de cada fila es igual a la unidad ya que P U = U . La primera igualdad se obtiene por la propiedad asociativa de la multiplicaci´on para las matrices, y la segunda y tercera por definici´on de matriz estoc´astica para P 0 y P , respectivamente. Por lo tanto, podemos decir que cada fila de la matriz de transici´on es una distribuci´on de probabilidad.

´ A LAS CADENAS DE MARKOV CAP´ITULO 3. INTRODUCCION

16

3.3.

Distribuci´ on inicial

Una cadena de Markov est´a completamente determinada por la distribuci´on de probabilidad del estado inicial X0 y por las probabilidades de transici´on pij . La distribuci´on inicial de la cadena se expresa en forma de vector, en el que cada componente indica la probabilidad de que la cadena se encuentre en el estado i en el instante inicial. De esta forma, se conoce el punto de partida del proceso. Se expresa como (0) (0) (0) P (0) = (p1 , . . . , pi , . . . , pk ), (3.7) (0)

donde pi = P (X0 = i), i ∈ E. (0)

Adem´as, la distribuci´on inicial cumple: pi ≥ 0 y

X

(0)

pi = 1.

i∈E

3.4.

Distribuci´ on de probabilidad en la etapa n´ esima

Una vez definida la cadena de Markov en tiempo discreto, se puede obtener la distribuci´on marginal de Xn , es decir, la distribuci´on de la cadena en la etapa n´esima. Si se denota la probabilidad de que en la etapa n la cadena de Markov se en(n) cuentre en el estado i por pi = P (Xn = i), para cada etapa se tiene un vector (n) (n) (n) P (n) = (p1 , . . . , pi , . . . , pk ) que representa la probabilidad de que la cadena se encuentre en cada uno de los posibles estados en la etapa n. Para poder obtener esta distribuci´on se deber´a calcular previamente la matriz de probabilidades de transici´on en n etapas. Partiendo de que se conoce la probabilidad de transici´on en una etapa dada, el siguiente paso es obtener la probabilidad de transici´on en n etapas, (n)

pij = P (Xm+n = j|Xm = i),

∀n, m ∈ N ; i, j ∈ E.

´ Estas se obtienen calculando la potencia n-´esima de la matriz de transici´on P . Este resultado se obtiene a partir de la ecuaci´on de Chapman-Kolmogorov para cadenas de Markov en tiempo discreto X (m) (n) (m+n) pij = pig pgj , ∀n, m ∈ N ; i, j ∈ E. (3.8) g∈E

Intuitivamente, para pasar del estado i al j en (m + n) etapas, se debe pasar por un estado g en m etapas y despu´es ir desde g hasta j en las n etapas restantes. La condici´on de Markov implica que las dos partes de la transici´on de i a j son

3.5. ORDEN DE UNA CADENA DE MARKOV

17

(m+n)

independientes, por lo que se puede escribir pij como el producto de las proba(m+n) (m) (n) bilidades de transici´on y por lo tanto P = P P , donde P (m+n) es la matriz (m+n) de transici´on de la etapa (m + n) y tiene como elementos (pij ). Tomando n = 1 en la ecuaci´on (3.8) se obtiene X (m) (m+1) pij = pig pgj , ∀n, m ∈ N ; i, j ∈ E.

(3.9)

g∈E

Se tiene entonces que P (1) = P , P (2) = P (1) P (1) = P 2 , y as´ı sucesivamente hasta que obtenemos P (n) = P (n−1) P = P (n−2) P 2 = . . . = P (0) P (n) , y con esta f´ormula tenemos la distribuci´on de la cadena en la etapa n. Por lo tanto, la probabilidad de (n) transici´on en n etapas, pij = P (Xm+n = j|Xm = i) se obtiene del elemento (i, j) de la n-´esima potencia de la matriz de transici´on P ; ∀n, m ∈ N ; i, j ∈ E.

3.5.

Orden de una Cadena de Markov

El orden de una cadena de Markov establece el n´ umero de estados anteriores de los cuales depende la probabilidad de un estado, en un instante dado del proceso. As´ı, dado E = {E1 , . . . , Ek }, en una cadena de primer orden tendremos P (Xn+1 = in+1 |X0 = i0 , X1 = i1 , . . . , Xn = in ) = P (Xn+1 = in+1 |Xn = in ), (3.10) y en una cadena de orden dos P (Xn+1 = in+1 |X0 = i0 , . . . , Xn−1 = in−1 , Xn = in ) = = P (Xn+1 = in+1 |Xn−1 = in−1 , Xn = in ),

(3.11)

para i ∈ E y n ∈ N. An´alogamente se definir´ıa una cadena de Markov de orden mayor que dos.

3.6.

Clasificaci´ on de estados

Dentro de los estados de una cadena de Markov en tiempo discreto, se puede hacer una clasificaci´on de ´estos seg´ un las transiciones permitidas entre ellos y las probabilidades de pasar de un estado a otro. Para realizar dicha clasificaci´on es necesario primero definir los conceptos de probabilidad y tiempos de primera pasada. Definici´ on 3.3 Probabilidad de primera pasada. Es la probabilidad de que, empezando en i, la cadena pase por primera vez por el estado j en la etapa n. Se denota por (n) fij = P (Xn = j, Xr 6= j, ∀r < n|X0 = i). (3.12)

18

´ A LAS CADENAS DE MARKOV CAP´ITULO 3. INTRODUCCION

Definici´ on 3.4 Probabilidad de pasada. Se trata de la probabilidad de que la cadena llegue alguna vez al estado j partiendo del estado i. fij =

∞ X

(n)

fij = P (∃ n; Xn = j|X0 = i).

(3.13)

n=1

Definici´ on 3.5 Tiempo de primera pasada. Es el n´ umero de etapa en la que la cadena llega por primera vez al estado j cuando parte del estado i. Se determina por la siguiente variable aleatoria: Nij = {no de la primera etapa en la cual la cadena est´a en j partiendo de i}. Se describe entonces (n)

fij = P (Nij = n), fij =

∞ X

(n)

fij = P (Nij < ∞).

(3.14)

n=1

Se establece la siguiente relaci´on entre las probabilidades de transici´on en n (n) (n) etapas, pij , y las probabilidades de primera pasada ,fij , (n)

(1) (n−1)

pij = fij pjj

(2) (n−2)

+ fij pjj

(n−1)

+ . . . + fij

(n)

pjj + fij .

(3.15)

Se puede hacer entonces la siguiente clasificaci´on del espacio de estados E de una cadena de Markov en tiempo discreto 1. Estado recurrente: Un estado j ∈ E se dice que es recurrente si fjj = 1, es decir, es seguro que la cadena va a volver al estado j una vez que ya ha llegado a ´el en alguna etapa. De otra forma, es recurrente, si y solo si el n´ umero de veces que se espera que la cadena pase por ´el cuando ya parte de ´el es infinito. 2. Estado transitorio: Se dice que un estado j ∈ E es transitorio si no es recurrente, es decir, fjj < 1. En este caso se espera un n´ umero finito de visitas a dicho estado. 3. Estado que comunica con otro estado: Se dice que un estado i ∈ E comunica con otro estado j ∈ E, si fij > 0, es decir, existe la posibilidad de que la cadena llegue al estado j partiendo del estado i. Se cumple tambi´en si alguna potencia de la matriz de transici´on otorga probabilidad no nula a la transici´on entre dos estados. Se dice que una cadena de Markov es irreducible si todos sus estados comunican entre s´ı. 4. Estado que intercomunican. Dos estados, i, j ∈ E intercomunican si i comunica con j y j comunica con i. 5. Estado ef´ımero: Se dice que un estado j ∈ E es ef´ımero si pij = 0, ∀i ∈ E. No se puede llegar a ´el desde ning´ un otro, s´olo se puede salir desde ´el hacia cualquier otro.

3.7. DISTRIBUCIONES ESTACIONARIAS

19

6. Estado absorbente: Un estado j ∈ E se dice absorbente si es imposible abandonarlo, es decir, pjj = 1. Una vez que se alcanza, la cadena s´olo puede mantenerse en ´el. Adem´as, una cadena de Markov en tiempo discreto se dice que es absorbente si tiene al menos un estado absorbente y desde cada uno de los estados es posible alcanzar alguno de los absorbentes en un n´ umero finito de etapas. Esta clase de cadenas de Markov se utilizan en Bioinform´atica para, por ejemplo, encontrar la distancia media entre sucesivas ocurrencias de una “palabra” (como AAG) en una peque˜ na secuencia de ADN. 7. Estado peri´odico: Un estado i ∈ E se dice que es peri´odico con per´ıodo l > 1 si l es el menor n´ umero tal que todas las secuencias de transiciones que parten del estado i y regresan al estado i tienen una longitud m´ ultiplo de l. Si un estado no es peri´odico se llama aperi´odico. Una cadena de Markov se dice aperi´odica si todos sus estados son aperi´odicos. Asumiremos que todas las cadenas de Markov que veremos en adelante son finitas, aperi´odicas e irreducibles.

3.7.

Distribuciones estacionarias

Supongamos que una cadena de Markov tiene una matriz de transici´on P y que en el instante t la probabilidad de que el proceso est´e en el estado j es ϕj , j = 1, 2, . . . , k. Esto implica que la probabilidad de que en el instante t + 1 el proceso est´e en el k X ϕr prj . Supongamos que, para todo j, estas dos probabilidades son estado j es r=1

iguales, por lo tanto ϕj =

k X

ϕr prj ,

j = 1, 2, ..., k.

(3.16)

r=1

En este caso decimos que la distribuci´on de probabilidad (ϕ1 , ϕ2 , ..., ϕk ) es estacionaria; esto es, no ha cambiado entre los instantes t y t + 1, y por lo tanto nunca cambiar´a. Definici´ on 3.6 Una distribuci´on Φ = {ϕj }j∈E sobre E se dice estacionaria respecto de una cadena de Markov con matriz de transici´on PX si verifica que ΦP = Φ. De la ϕj = 1 y se verifica misma forma, Φ es una distribuci´on estacionaria si j∈E

ϕj =

X i∈E

ϕi pij ,

∀j ∈ E

(3.17)

´ A LAS CADENAS DE MARKOV CAP´ITULO 3. INTRODUCCION

20

Los valores de la distribuci´on estacionaria se pueden interpretar como la proporci´on final de tiempo que la cadena ha pasado en cada estado a lo largo de su evoluci´on, con probabilidad 1. Se puede decir que es tambi´en la proporci´on, a largo plazo, de etapas en las que la cadena se encuentra en el estado i a lo largo de su evoluci´on, si ha partido de i o de otro estado recurrente que intercomunica con i. Para una cadena de Markov finita, aperi´odica e irreducible, existe una u ´nica distribuci´on estacionaria. Si una cadena de Markov es finita, aperi´odica e irreducible, entonces si n crece, P (n) se aproxima una matriz k ×k con todas las filas iguales al vector (ϕ1 , ϕ2 , ..., ϕk ), que es la distribuci´on estacionaria de la cadena de Markov, siendo P (n) la matriz de las probabilidades de transici´on al cabo de n pasos. La forma de esta matriz muestra que no importa cual sea el estado inicial, o cual sea la distribuci´on de probabilidad del estado inicial, y que la probabilidad de que n pasos despu´es el proceso est´e en el estado j se aproxima cada vez m´as, cuando n → ∞, al valor ϕj .

3.8.

Ejemplo: una cadena de Markov como modelo para el ADN

Es improbable que una secuencia aleatoria compuesta por las “letras” A, C, G y T sea un buen modelo para el patr´on de nucle´otidos en una secuencia g´enica, como ya vimos en el cap´ıtulo anterior. Una cadena de Markov con {A, C, G, T } como estados podr´ıa ser una mejor aproximaci´on a la realidad: las probabilidades para el nucle´otido en la posici´on n + 1 dependen del nucle´otido en la posici´on n (sin embargo, en la realidad las dependencias son m´as complejas). En este caso, se sustituye el concepto del tiempo t por la “posici´on n” en la secuencia de ADN. Si el espacio de estados es E = {A, C, G, T }, entonces la matriz de transici´on ser´a de la forma   pAA pAC pAG pAT  pCA pCC pCG pCT   P = (3.18)  pGA pGC pGG pGT  , pT A pT C pT G pT T con pij = P (Xn+1 = j|Xn = i), para n ≥ 1, donde i, j ∈ E. Adem´as, la distribuci´on estacionaria ser´a (para k = 4 estados) ϕ0 = (ϕA , ϕC , ϕG , ϕT ).

3.8.1.

Ejemplo

Supongamos que tenemos la secuencia de ADN CCGAT, de longitud 5. Vamos a hallar la probabilidad de la secuencia dado el modelo de Markov de primer orden caracterizado por (0)

(0)

(0)

(0)

1. la distribuci´on inicial P (0) = (pA , pC , pG , pT ) = (0.2, 0.1, 0.1, 0.6).

3.8. EJEMPLO: UNA CADENA DE MARKOV COMO MODELO PARA EL ADN21 2. la matriz de transici´on 

0.10  0.35 P =  0.30 0.60

0.80 0.10 0.20 0.10

0.05 0.10 0.20 0.25

 0.05 0.45  . 0.30  0.05

(3.19)

Se trata de una cadena de Markov finita, aperi´odica e irreducible. El espacio de estados E es en nuestro caso E = {A, C, G, T }, finito con k = 4. La cadena es irreducible ya que todos los estados se comunican entre s´ı (las probabilidades de transici´on en cualquier etapa ser´an no nulas). Adem´as, no existe ning´ un estado (n) (n) peri´odico, es decir, no existe l tal que pii = 0 para n = l, 2l, 3l, . . ., o pii 6= 0 para n = l, 2l, 3l, . . .; cualquiera que sea la etapa n. Sean las variables aleatorias Xn para n = 0, 1, 2, 3, 4, que indican el tipo de nucle´otido j en la posici´on n de la secuencia CCGAT, ∀j ∈ E. X0 es una variable aleatoria que determina el nucle´otido en la posici´on inicial de la secuencia. De esta forma, X1 indica el nucle´otido en la primera posici´on de la secuencia, X2 el nucle´otido en la segunda, y as´ı sucesivamente. Calculemos entonces la probabilidad de la secuencia CCGAT, es decir, tenemos que hallar P (CCGAT ) = P (X0 = C) × P (X1 = C|X0 = C) × P (X2 = G|X1 = C) × × P (X3 = A|X2 = G) × P (X4 = T |X3 = A) = (0)

= pC × pCC × pCG × pGA × pAT (0)

Para empezar, tenemos que pC = P (X0 = C) = 0.1, teniendo en cuenta la distribuci´on inicial de la cadena. (1)

Por otro lado, pCC = pCC = P (X1 = C|X0 = C) = 0.1, observando la matriz de transici´on P indicada anteriormente. De la misma forma, procedemos con todas las probabilidades de transici´on indicadas en la matriz P , y as´ı obtenemos que la probabilidad de la secuencia de ADN CCGAT es, en este caso, (0)

P (CCGAT ) = pC × pCC × pCG × pGA × pAT = = 0.1 × 0.1 × 0.1 × 0.3 × 0.05 = 0.000015. Si consider´asemos que esta secuencia se distribuye de forma iid (independiente e id´enticamente distribuida), siendo la probabilidad de cada nucle´otido la misma que la de la distribuci´on inicial de la cadena de Markov indicada, la probabilidad de la secuencia ser´ıa P (CCGAT ) = pC × pC × pG × pA × pT = = 0.1 × 0.1 × 0.1 × 0.2 × 0.6 = 0.00012, que es mayor que la probabilidad calculada teniendo en cuenta que la secuencia sigue un modelo de Markov de primer orden como el anterior. Se observa una muy leve diferencia debido al tama˜ no peque˜ no de la secuencia.

22

´ A LAS CADENAS DE MARKOV CAP´ITULO 3. INTRODUCCION

Cap´ıtulo 4 An´ alisis de las secuencias del ADN Despu´es de haber introducido la teor´ıa necesaria anterior, pasamos ya al an´alisis de secuencias de ADN y al estudio sus caracter´ısticas. Desde que el bacteri´ofago P hi − X174 fue secuenciado en 1977-1978, las secuencias de ADN de cientos de organismos han sido decodificadas y guardadas en bases de datos. Esos datos son analizados para determinar, entre otras cosas, los genes que codifican ciertas prote´ınas. Una comparaci´on de genes en una especie o entre especies puede mostrar similitudes entre funciones de prote´ınas, o relaciones entre especies. Con la creciente cantidad de datos, desde hace mucho se ha vuelto poco pr´actico analizar secuencias de ADN manualmente. Hoy en d´ıa, se utilizan programas inform´aticos para estudiar el genoma de miles de organismos, conteniendo miles de millones de nucle´otidos. Estos programas pueden compensar mutaciones (con bases intercambiadas, borradas o insertadas) en la secuencia de ADN, para identificar secuencias que est´an relacionadas, pero que no son id´enticas. Una variante de este alineamiento de secuencias se usa en el “proceso de secuenciaci´on”. El m´etodo de secuenciaci´on conocido como “shotgun” o “por perdigonada” fue utilizado por el Instituto de Investigaci´on Gen´omica para secuenciar el primer genoma de una bacteria, la Haemophilus influenzae. No da una lista secuencial de nucle´otidos, pero en cambio nos ofrece las secuencias de miles de peque˜ nos fragmentos de ADN (cada uno de aproximadamente 600 a 800 nucle´otidos de largo). Las terminaciones de estos fragmentos se superponen y, cuando son alineados de la manera correcta, constituyen el genoma completo del organismo en cuesti´on. El secuenciamiento “shotgun” proporciona datos de secuencia r´apidamente, pero la tarea de ensamblar los fragmentos puede ser bastante complicada para genomas muy grandes. En el caso del Proyecto del Genoma Humano, se requirieron varios meses de tiempo de procesador para ensamblar los fragmentos. El “shotgun sequencing” es el m´etodo elegido para todos los genomas secuenciados hoy en d´ıa, y los algoritmos de ensamblado gen´omico son un a´rea cr´ıtica de la investigaci´on en Bioinform´atica. No profundizaremos en el an´alisis de este m´etodo debido a su complejidad t´ecnica ([40]).

23

24

´ CAP´ITULO 4. ANALISIS DE LAS SECUENCIAS DEL ADN

Otro aspecto importante de la Bioinform´atica en el an´alisis de secuencias es la b´ usqueda autom´atica de genes. Veremos de forma introductoria c´omo modelizar secuencias de genes en el cap´ıtulo 5, como aplicaci´on de los Modelos de Markov Ocultos. En su nivel m´as elemental, como ya comentamos, la estructura del ADN puede pensarse como largas secuencias de nucle´otidos. Estas secuencias est´an organizadas en secuencias codificantes, o genes, separadas por largas regiones interg´enicas de secuencias no codificantes. La mayor´ıa de los genes eucariotas tienen un nivel adicional de organizaci´on: dentro de cada gen, las secuencias codificantes (exones) son normalmente interrumpidas por tramos de secuencias no codificantes (intrones). Durante la transcripci´on del ADN en ARN mensajero (ARNm), los intrones se eliminan y no aparecen en el u ´ltimo ARNm, el cual se traduce en secuencia de prote´ına. Las regiones interg´enicas y los intrones tienen diferentes propiedades estad´ısticas que los exones. Para capturar estas propiedades en un modelo, podemos construir procedimientos estad´ısticos y comprobar si una parte no caracterizada de ADN es parte de la regi´on codificante del gen. El modelo est´a basado en un conjunto de datos “training” o datos “de entrenamiento” tomados de secuencias ya conocidas. En el modelo m´as sencillo, se asume que los nucle´otidos en posiciones distintas son independientes y est´an id´enticamente distribuidos (iid). Si ´este es el caso, las diferencias entre ADN codificante y no codificante podr´ıan detectarse por las diferencias entre las frecuencias de los cuatro nucle´otidos en los dos casos distintos (intr´on vs ex´on). Estas distribuciones pueden estimarse mediante los datos de entrenamiento. Si no se utilizan datos de entrenamiento como referencia, el an´alisis es m´as complicado. En las aplicaciones pr´acticas que presentaremos en esta memoria s´ı los utilizaremos.

4.1.

Contrastes de independencia

La precisi´on de los procedimientos a llevar a cabo depende de la precisi´on de las suposiciones hechas. Podemos intuir que suponer que hay independencia entre los nucle´otidos suele ser una simplificaci´on excesiva, puesto que puede haber correlaciones, por ejemplo, entre los nucle´otidos debido a su pertenencia a uno u otro cod´on. Por ello, es importante, entre otras cosas, desarrollar un contraste de independencia sobre la secuencia de nucle´otidos. Este contraste est´a basado en el an´alisis de una cadena de Markov. Al utilizar cadenas de Markov en el an´alisis de secuencias, el concepto del “tiempo t” se reemplaza por la “posici´on a en la secuencia”. Nuestras variables aleatorias o sucesos ser´an “los nucle´otidos en posiciones dadas de la secuencia de ADN” y, por tanto, el espacio de estados estar´a compuesto por los cuatro nucle´otidos: E = {A, C, G, T }, de tama˜ no k = 4. Desarrollamos entonces un test 2 de independencia chi-cuadrado χ sobre la secuencia de nucle´otidos que queremos analizar, que contrasta la hip´otesis nula de que los nucle´otidos en posiciones distintas sean independientes, frente a la hip´otesis alternativa de que un nucle´otido en una posici´on dada dependa del nucle´otido en la posici´on anterior. De esta forma, la

4.1. CONTRASTES DE INDEPENDENCIA

25

distribuci´on del contraste es una cadena de Markov en la cual todas las filas de la matriz de transici´on son id´enticas. Por lo tanto, la hip´otesis alternativa la pensamos como una cadena de Markov donde la probabilidad de un nucle´otido en una posici´on dada, dependa del nucle´otido en la posici´on anterior, es decir, una cadena de Markov de orden uno. As´ı, como ya vimos en el ejemplo del cap´ıtulo anterior, bajo la hip´otesis nula de independencia entre nucle´otidos, y dada la secuencia de ADN X = GATTACA, por ejemplo, la probabilidad de tal secuencia bajo el modelo de independencia ser´a P (X) = pG pA pT pT pA pC pA = p3A p1C p1G p2T , donde pi indica la probabilidad del nucle´otido i en la secuencia, ∀i ∈ {A, C, G, T }. Este contraste es de inter´es para evaluar si el modelo de Markov bajo la hip´otesis alternativa describe la realidad significativamente mejor que el modelo de independencia y, por tanto, podr´ıa aumentar la precisi´on de nuestros procedimientos de predicci´on. Adem´as, si lo hace, podr´ıa considerarse incluso el modelo con un nivel de dependencia m´as complejo (una cadena de Markov con mayor orden). Supongamos que la longitud de la secuencia de ADN que queremos analizar (nuestra muestra) es n. El contraste estad´ıstico de independencia es un contraste de asociaci´on en una tabla 4 × 4 de doble entrada como la siguiente, denominada tabla de contingencia.

A C G T Total

A n11 n21 n31 n41 n,1

C n12 n22 n32 n42 n,2

G n13 n23 n33 n43 n,3

T n14 n24 n34 n44 n,4

Total n1, n2, n3, n4, n

Tabla 4.1: Tabla de contingencia de los cuatro nucle´otidos. Tenemos entonces un contraste de independencia con la hip´otesis nula H0 : pij = pi × pj ,

∀i, j ∈ E.

Los datos deben aparecer como se muestra en la tabla 4.1, donde las filas representan el nucle´otido en la posici´on a y las columnas el nucle´otido en la posici´on a + 1. El elemento (i, j) de la tabla representa la frecuencia observada nij , que indica el n´ umero de transiciones Xdel estado o nucle´otido i al estado o nucle´otido j en una etapa, para i, j ∈ E; y nij = n. As´ı, cada uno de los elementos ni, y n,j reprei,j∈E

sentan la suma de la fila i (n´ umero de transiciones desde el estado i) y la columna j (n´ umero de transiciones que llegan al estado j), respectivamente, ∀i, j ∈ E. La idea es realizar el contraste anterior comparando las frecuencias esperadas bajo la hip´otesis nula, Tij = npi pj , con las observadas, Oij = nij . Si las cantidades pi y pj

26

´ CAP´ITULO 4. ANALISIS DE LAS SECUENCIAS DEL ADN

no son conocidas, han de ser estimadas a partir de las frecuencias observadas de la siguiente forma n,j ni, , pbj = , (4.1) pbi = n n y por lo tanto ni, n,j Tij = nb pi pbj = , (4.2) n lo que hace perder (k − 1) + (k − 1), con k = 4, grados de libertad adicionales al estad´ıstico de contraste 2 X(k−1)+(k−1)

k X k X (nij − Tij )2 = ' χ2(k−1)+(k−1) . T ij i=1 j=1

(4.3)

que sigue, bajo la hip´otesis nula, una distribuci´on aproximada chi-cuadrado con tales grados de libertad. De esta forma, rechazaremos H0 si 2 > χ2(k−1)+(k−1),1−α , X(k−1)+(k−1)

y el p-valor correspondiente ser´a menor o igual que α, para α = 0.05 . Por lo tanto, la hip´otesis nula de independencia resulta ser la hip´otesis nula de no asociaci´on en la tabla. Pruebas de este tipo muestran que nucle´otidos contiguos en posiciones de secuencias de ADN son a veces dependientes, y el modelo de una cadena de Markov de primer orden se ajusta a los datos reales significativamente mejor que el modelo de independencia, tanto en regiones de intrones como de exones. Modelos homog´eneos de Markov de orden mayor y modelos no homog´eneos a veces pueden ajustarse mejor a los datos. Podemos ahora extender el an´alisis de independencia realizado a los contrastes de cadenas de Markov de orden superior (mayor que 1), por un procedimiento de raz´on de verosimilitud generalizada. Presentamos este an´alisis en t´erminos de secuencias de ADN tambi´en: queremos saber cu´al es el orden de la cadena de Markov que describe significativamente mejor las secuencias de ADN. La estructura probabil´ıstica de los nucle´otidos en una secuencia de ADN se describe por una cadena de Markov de orden d ≥ 1 si la probabilidad de que cualquier nucle´otido ocurra en una posici´on dada depende de los nucle´otidos de las d anteriores posiciones. Las probabilidades de transici´on para el caso d = 3, por ejemplo, son de la forma P (Xn+3 = in+3 |Xn = in , Xn+1 = in+1 , Xn+2 = in+2 ), donde Xa denota el nucle´otido en la posici´on a del la secuencia, y los ia son tipos espec´ıficos de nucle´otidos, ∀i ∈ E y a = 1, 2, . . . , n. Esta notaci´on muestra que para un d general, habr´a 3×4d probabilidades de transici´on entre distintos nucle´otidos.

4.1. CONTRASTES DE INDEPENDENCIA

27

Supongamos que la hip´otesis nula es que la cadena de Markov es de orden d − 1, y la hip´otesis alternativa que es de orden d, para alg´ un valor de d (lo discutido en la secci´on anterior es an´alogo con d = 1). En t´erminos del contraste anterior,  H0 : el nucle´otido en la posici´on a es independiente del de la posici´on a − d. Ha : el nucle´otido en la posici´on a depende del de la posici´on a − b. para todo b = 1, . . . , d. Aplicamos el test de raz´on de verosimilitud generalizada ([8]), seg´ un el cual Teorema 4.1 Sea X1 , . . . , Xn una muestra aleatoria de una poblaci´on con funci´on de densidad de probabilidad multiparam´etrica dada por f (x, θ), con θ = (θ1 , . . . , θm ). Si θ ∈ Θ0 , en que Θ0 es el correspondiente espacio param´etrico bajo H0 . Entonces el estad´ıstico T (X) = −2log(λ) tiene una distribuci´on asint´otica χ2m−r , donde r es el n´ umero de par´ametros de θ que ha sido completamente especificado en la hip´otesis nula y donde λ tiene la siguiente expresi´on supL0 fn (x, θb0 ) supθ∈Θ0 fn (x, θ) = = . (4.4) λ= b supθ∈Θ fn (x, θ) supL fn (x, θ) En nuestro caso, la funci´on de densidad de probabilidad corresponde a la matriz de transici´on de la cadena de Markov, Pb, donde los elementos de θ representan n cada una de las probabilidades de transici´on estimadas entre nucle´otidos, pbij = niji, , ∀i, j ∈ E. El n´ umero ni, representa la sumaX de todos los procesos que salen del estado i considerado en cada caso, es decir ni, = nij . As´ı, la suma por filas de la matriz j∈E

Pb es l´ogicamente 1. Por lo tanto, la matriz de transici´on Pb estimada por m´axima verosimilitud se puede obtener a partir de la tabla de contingencia de las variables estado anterior y estado actual (tabla 4.1), calculando las distribuciones marginales por filas. El n´ umero de par´ametros r de θ bajo la hip´otesis nula es 3 × 4(d−1) , por definici´on de orden de una cadena de Markov. Por tanto, el n´ umero de par´ametros (m) bajo la hip´otesis alternativa es 3 × 4d . Teniendo en cuenta esto, el estad´ıstico T de la definici´on anterior tendr´a una distribuci´on asint´otica chi-cuadrado con m − r = 3 × 4d − 3 × 4(d−1) = 9 × 4(d−1) grados de libertad. En este caso, la expresi´on de λ ser´a la siguiente Y n pbij 0 ij i,j∈E

λ= Y

n

pbijij

,

(4.5)

i,j∈E

donde pbij 0 denotan las probabilidades de transici´on en una cadena de Markov de orden d − 1, es decir, bajo la hip´otesis nula. En general, se conoce que incluso secuencias de ADN de regiones interg´enicas suelen seguir modelos de cadenas de Markov de orden superior a uno.

28

4.2.

´ CAP´ITULO 4. ANALISIS DE LAS SECUENCIAS DEL ADN

Contraste de bondad de ajuste

Si el contraste de independencia no es aceptado, se aplica otro test de inter´es para saber si la distribuci´on estacionaria de la cadena de Markov es uniforme. De esta forma, si fuese as´ı, los c´alculos llevados a cabo con los datos de entrenamiento para hallar los par´ametros del modelo se facilitar´ıan bastante. Aplicaremos entonces un contraste de bondad de ajuste a tal distribuci´on. La hip´otesis nula ser´a que las distribuciones de probabilidad estacionaria y uniforme son id´enticas, y la hip´otesis alternativa que ambas distribuciones son diferentes. Podemos discutirlo en el caso de una cadena de Markov de primer orden, ya que para cadenas de o´rdenes superiores es an´alogo. Si denotamos por Φ la distribuci´on estacionaria y por U a la funci´on de distribuci´on uniforme, el contraste de hip´otesis planteado se puede escribir de la forma  H0 : Φ ∼ = U. Ha : Φ  U. En el caso de primer orden, la ecuaci´on (3.17) muestra que una condici´on necesaria y suficiente para que la distribuci´on estacionaria sea uniforme es que las probabilidades de transici´on en cada columna de la matriz de transici´on sumen 1. Es el siguiente teorema. Teorema 4.2 Si la matriz de transici´on P de una cadena de Markov con k estados es doblemente estoc´astica, entonces la distribuci´on uniforme U (i) = k1 para todo i ∈ E, es una distribuci´on estacionaria. Demostraci´ on 4.2.1 Observamos que X i∈E

U (i)pij =

1X 1 pij = , k i∈E k

∀j ∈ E

(4.6)

de modo que la distribuci´on uniforme satisface la condici´on U P = U que define una X distribuci´on estacionaria (junto con que, por ser distribuci´on de probabilidad, U (i) = 1). Vemos adem´as, que si la distribuci´on estacionaria es uniforme, nei∈E

cesariamente la matriz P es doblemente estoc´astica, es decir, la suma por filas y columnas de P ha de ser 1. En nuestro caso, suponemos entonces que la distribuci´on uniforme es U = 14 . Esto implica que, por ejemplo, los elementos en la cuarta fila de la matriz de transici´on est´an determinados por los elementos de las primeras tres filas. La probabilidad de cualquier secuencia observada de ADN puede calcularse bajo la hip´otesis nula de que la distribuci´on estacionaria es uniforme. Ya que bajo esta hip´otesis los elementos de cada fila y cada columna de la matriz de transici´on deben sumar 1, hay 9 par´ametros libres (4 × 3 hab´ıa cuando solo las filas sumaban 1, menos 3 que se restan cuando

´ DE “SENALES” ˜ 4.3. MODELIZACION EN EL ADN

29

tambi´en las columnas lo hacen). La probabilidad de la secuencia observada puede maximizarse con respecto a estos par´ametros. Bajo la hip´otesis alternativa, la u ´nica restricci´on es que los elementos de cualquier fila de la matriz de transici´on deben sumar 1, y por lo tanto habr´a 12 par´ametros libres. As´ı, de la misma forma, la probabilidad de la secuencia observada puede ser maximizada respecto a estos par´ametros. De estas dos probabilidades puede calcularse un estad´ıstico de verosimilitud −2logλ como el se˜ nalado en la secci´on anterior, donde λ es λ=

1 supL0 fn (x, θb0 ) supθ∈Θ0 fn (x, θ) 4 = = . = b supθ∈Θ fn (x, θ) supL sup EΦ fn (x, θ)

(4.7)

Seg´ un este contaste, bajo la hip´otesis nula, este estad´ıstico tiene una distribuci´on asint´otica chi-cuadrado con m − r = 12 − 9 = 3 grados de libertad. De esta forma, podemos aplicar un contraste de bondad de ajuste entre ambas distribuciones y comprobar si la distribuci´on estacionaria se identifica significativamente con la distribuci´on uniforme, lo que facilita los c´alculos a la hora de modelizar secuencias de ADN.

4.3.

Modelizaci´ on de “se˜ nales” en el ADN

En el contexto de la gen´omica, la anotaci´on o puntuaci´on es el proceso de marcado de los genes y otras caracter´ısticas biol´ogicas de la secuencia de ADN. El primer sistema software de anotaci´on de genomas fue dise˜ nado en 1995 por Owen White, director de Bioinform´atica de la escuela de medicina de la universidad de Maryland, quien fue miembro del equipo que secuenci´o y analiz´o el primer genoma de un organismo independiente en ser decodificado, la bacteria Haemophilus influenzae (responsable de un amplio rango de enfermedades como meningitis, epiglotitis, neumon´ıa, sepsis y otras de menor gravedad). White construy´o un software para localizar los genes, el ARN de transferencia, y otras caracter´ısticas; as´ı como para realizar las primeras atribuciones de funci´on a esos genes. La mayor´ıa de los actuales sistemas de anotaci´on gen´omica trabajan de forma similar, pero los programas disponibles para el an´alisis del genoma se encuentran en continuo cambio y mejora. En este contexto, se conoce que los genes contienen “se˜ nales” en el ADN, que son secuencias de ADN con un prop´osito espec´ıfico: indicar, por ejemplo, el comienzo y final de la regi´on transcrita, los l´ımites de los exones e intrones, as´ı como otras caracter´ısticas. La maquinaria de la c´elula utiliza estas se˜ nales para reconocer el gen, para editarlo correctamente, y para traducirlo apropiadamente en prote´ına. Si la naturaleza fuese bondadosa, cada se˜ nal consistir´ıa en una u ´nica secuencia de ADN que no aparecer´ıa en ning´ un otro lugar del ADN excepto donde sirviese para su prop´osito espec´ıfico. En la realidad, hay muchas secuencias de ADN que realizan la misma funci´on de se˜ nal; llamamos a ´estas “miembros” de una se˜ nal. Adem´as, los miembros

´ CAP´ITULO 4. ANALISIS DE LAS SECUENCIAS DEL ADN

30

de las se˜ nales tambi´en aparecen aleatoriamente en el ADN no funcional, haciendo dif´ıcil clasificar las se˜ nales funcionales de las no funcionales. En la pr´actica, no todos los miembros de una se˜ nal son conocidos. Nuestro objetivo es utilizar los miembros conocidos para evaluar la probabilidad de que una nueva secuencia no caracterizada de ADN sea tambi´en un miembro de la se˜ nal. Supondremos que los diferentes miembros surgen de antepasados comunes mediante procesos estoc´asticos, por lo tanto, es razonable construir un modelo estad´ıstico de los datos. Algunas se˜ nales requieren solamente modelos simples, mientras que otras necesitan modelos m´as complejos. Cuando el modelo requerido es complejo y los datos son limitados, debemos tener cuidado eligiendo suposiciones simplificadoras en el modelo para utilizar los datos de la manera m´as eficiente posible. Asumimos que todos los miembros de la se˜ nal de inter´es tienen la misma longitud, que denotaremos por n. Esta suposici´on no es demasiado restrictiva, ya que los miembros de muchas se˜ nales tienen la misma longitud, y para los que no la tengan, podemos capturar porciones de ellos, que es normalmente suficiente. Para modelar las propiedades de cualquier se˜ nal debemos tener un conjunto de datos de entrenamiento, esto es, una gran cantidad de datos en la cual los miembros de las se˜ nales sean conocidos. Consideraremos ahora algunos de los modelos de se˜ nales b´asicos que se utilizan en Bioinform´atica.

4.3.1.

Matrices de Peso. Independencia

El contraste de hip´otesis sobre si una secuencia no caracterizada de ADN es miembro de una se˜ nal dada, es m´as f´acil de llevar a cabo cuando el nucle´otido en cualquier posici´on de la se˜ nal es independiente de los nucle´otidos en cualquier otra posici´on de esa se˜ nal. Es, por lo tanto, necesario realizar un contraste de independencia sobre si el nucle´otido en la posici´on a en una se˜ nal es independiente del nucle´otido en la posici´on b, sin necesidad de que a y b sean contiguos. Esto puede hacerse de varias formas. Es natural generalizar el contraste de independencia descrito en la tabla 4.1, que comprueba la independencia en posiciones contiguas. En esta generalizaci´on, “posici´on a” es reemplazada por “posici´on a en la se˜ nal” y “posici´on a + 1” se reemplaza por “posici´on b en la se˜ nal”. As´ı, nij es interpretado como el n´ umero de veces que, en los datos, el nucle´otido i ocurre en la posici´on a de la se˜ nal y el nucle´otido j ocurre en la posici´on b de la se˜ nal. Este contraste se utiliza luego para todos los pares a y b; con i, j ∈ E y a, b = 1, . . . , n. Por tanto, nos planteamos un contraste de independencia chi-cuadrado entre las posiciones de los nucle´otidos en la se˜ nal, es decir la hip´otesis nula es H0 : el nucle´otido en la posici´on a es independiente del de la posici´on b, ∀ a, b = 1, . . . , n. El contraste se desarrolla an´alogamente al test de independen-

´ DE “SENALES” ˜ 4.3. MODELIZACION EN EL ADN

31

cia para cadenas de Markov de orden uno, expuesto en el presente cap´ıtulo (donde las posiciones en la secuencia s´ı eran contiguas). Supongamos que, como resultado de este contraste, puede asumirse independencia. Se construye entonces una matriz 4 × n, donde cada fila corresponde a uno de los 4 nucle´otidos y cada columna a las posibles posiciones de la se˜ nal de longitud n. Su posici´on o entrada (i, a) es la proporci´on de casos en los datos que el nucle´otido i ocurre en la posici´on a de la se˜ nal. Nos referimos a ´esta como una matriz de peso y la denotaremos por M . Un ejemplo para el caso n = 5 se presenta en la siguiente matriz   0.33 0.34 0.19 0.20 0.21  0.31 0.18 0.34 0.30 0.25   M =  0.22 0.27 0.23 0.24 0.21  . 0.14 0.21 0.24 0.26 0.33

(4.8)

Los elementos de la columna a de la matriz dan (de forma aproximada) las probabilidades para los cuatro nucle´otidos (A, G, C y T , respectivamente) en la posici´on a de la se˜ nal, a = 1, . . . , n. La matriz M define la probabilidad P (X|M ) para cualquier secuencia (se˜ nal) X de longitud n. Las matrices de peso se utilizan como una componente de un algoritmo para b´ usqueda de genes que veremos brevemente en el siguiente cap´ıtulo. En general, al modelizar un gen, se necesitan matrices de peso para modelizar ciertas “regiones” de ´este ([40]).

4.3.2.

Dependencias de Markov

Si los nucle´otidos en las posiciones de la se˜ nal no son independientes, una posibilidad es que sus posiciones a lo largo de la se˜ nal sigan un modelo de Markov de primer orden. Bajo estas hip´otesis, los nucle´otidos en cada posici´on tienen una distribuci´on de probabilidad dependiente del nucle´otido en la posici´on anterior. As´ı, estas probabilidades se disponen una matriz de transici´on 4 × 4 cuyo elemento (i, j) es la probabilidad de que el nucle´otido j est´e en la posici´on a, dado el nucle´otido i en la posici´on a − 1. Los elementos de esta matriz son estimados a partir de los datos, como ya hemos visto en el contraste de independencia presentado para secuencias de ADN, no necesariamente se˜ nales. En general, podr´ıan considerarse dependencias en cadenas Markov de mayor orden. Sin embargo, estas matrices de transici´on pueden ser muy grandes a medida que el orden de la cadena de Markov crece y, por lo tanto, la cantidad de datos que se necesitan para satisfacer la estimaci´on puede ser excesiva. Por ello, cuando el conjunto de los datos es limitado, debemos economizar y capturar solamente las dependencias “m´as informativas” en nuestro modelo, que ser´an aquellas posiciones

´ CAP´ITULO 4. ANALISIS DE LAS SECUENCIAS DEL ADN

32

donde encontremos valores mayores de una tabla de contingencia como la tabla 4.1. Un m´etodo para hacer esto lo vemos a continuaci´on. Descomposici´ on de Dependencia M´ axima Puede ser imposible, debido a los datos limitados, obtener satisfactoriamente estimaciones de las distribuciones de probabilidad de una se˜ nal dada, en los casos de cadenas de Markov de orden mayor que uno. Esto motiva la b´ usqueda de un m´etodo de capture aquellas dependencias que sean “m´as informativas”. Describiremos ahora la aproximaci´on, conocida como descomposici´on de dependencia m´axima ([9]), para solucionar este problema. Supongamos que deseamos modelizar una se˜ nal de longitud n. El primer paso es encontrar una posici´on que tenga la “mayor influencia” en las otras. Para ello, construimos una tabla n × n cuyas entradas (a, b) sean los valores del estad´ıstico chicuadrado obtenidos de una tabla de contingencia como tabla 4.1, pero que compare los nucle´otidos de una posici´on fija a con los de otra posici´on b (en lugar de a + 1). Si la hip´otesis de independencia es cierta, esperamos encontrar alrededor de un valor significativo sobre 20, con un error de Tipo I o α del 5 %. As´ı, si algunos de los valores del estad´ıstico chi-cuadrado observados son “un poco” significativos, podr´ıamos concluir que las frecuencias observadas en la muestra no difieren significativamente de las que te´oricamente deber´ıan haberse obtenido bajo independencia, y por tanto no rechazar´ıamos la hip´otesis nula. Si “bastantes” valores son significativos, o si existen valores con “grandes” niveles de significaci´on, entonces podr´ıamos concluir que los nucle´otidos en varias posiciones no son independientes, y rechazar´ıamos la hip´otesis nula al haber encontrado evidencia en el sentido de que las dos variables consideradas est´an relacionadas: nucle´otido en posici´on a con nucle´otido en posici´on b. La tabla siguiente proporciona un ejemplo, con n = 5, en el cual muchos de los valores chi-cuadrado son significativos al nivel del 5 % (indicados por asteriscos). La fila con la suma mayor da una indicaci´on de cu´al es la posici´on que tiene mayor influencia sobre las otras n − 1 posiciones restantes.

1 2 3 4 5

1 2 3 4 5 total 0 34.2* 7.1 37.2* 2.8 81.3 34.2* 0 0.4 72.4* 4.5 111.5 7.1 0.4 0 15.3 98.3* 121.1 37.2* 72.4* 15.3 0 14.2 139.1 2.8 4.5 98.3* 14.2 0 119.8 Tabla 4.2: Ejemplo.

En este ejemplo, la posici´on 4 tiene el mayor total por filas, por lo tanto la toma-

´ 4.4. ANALISIS DE PATRONES

33

remos como la posici´on de mayor influencia y trataremos de evaluar esta influencia. Es este caso decimos que dividimos en la posici´on 4. Construimos primero un modelo que determine las distribuciones para las posiciones 1, 2, 3 y 5 condicionales de la posici´on 4. Para ello, dividimos la secuencia dada en 4 conjuntos de secuencias, cada uno determinado por el nucle´otido en la posici´on 4. Por lo tanto, para cada nucle´otido x tendremos un conjunto Tx constituido por aquellas secuencias donde haya una x en la posici´on 4. Despu´es calculamos px = ndx , para cada x = A, G, C, T , donde nx es el n´ uP mero de miembros de los datos que tienen el nucle´otido x en la posici´on 4, y d = nx . Luego, para cada x = A, G, C, T , calculamos una matriz de peso 4 × 4, Mx , de las secuencias Tx , para las posiciones 1, 2, 3 y 5, de la misma forma que la calculada en el caso de independencia entre nucle´otidos de la se˜ nal. El modelo de dependencia tiene entonces la distribuci´on {pA , pG , pC , pT } junto con las cuatro matrices de peso {MA , MG , MC , MT }. As´ı, para cualquier secuencia X de longitud 5 calcularemos la P (X) como sigue: si el nucle´otido x ocurre en la posici´on 4 de X, la matriz de peso Mx se utiliza para asignar una probabilidad pa a las posiciones a = 1, 2, 3 y 5 de X. Entonces, P (X) = px × p1 × p2 × p3 × p5 . En general, algunos o todos los conjuntos Tx ser´an suficientemente grandes (secuencias de longitud grande) para que podamos repetir el proceso entero recursivamente, dividiendo Tx en una de las posiciones 1, 2, 3 y 5. Es este caso, para cada Tx suficientemente grande construiremos una tabla similar a la 4.2, esta vez 4 × 4, para encontrar la posici´on de esas cuatro que tenga la influencia mayor en las otras tres. Entonces, descompondremos Tx en Txy , para y = A, G, C y T ; continuaremos a Txyz y as´ı sucesivamente, siempre y cuando haya datos suficientes. Una regla que podemos utilizar es parar de buscar posiciones influyentes cuando el conjunto Txyz... tenga menos de 100 secuencias. Cuando ese l´ımite se alcanza, las posiciones restantes se modelan con una matriz de peso, con la que se procede de la forma explicada en la secci´on anterior.

4.4.

An´ alisis de Patrones

Un patr´on es una secuencia de ADN que se repite en determinadas posiciones de otra secuencia m´as larga, correspondiente a distintas partes de la estructura del ADN. Supongamos que estamos interesados en una “palabra”, por ejemplo GAGA. Nos preguntamos las dos siguientes cuestiones de una secuencia de ADN iid de longitud n: ¿cu´al es el n´ umero medio de veces que esta palabra ocurre en un segmento de longitud n? y ¿cu´al es la longitud media entre una ocurrencia de esta palabra y la siguiente?.

34

´ CAP´ITULO 4. ANALISIS DE LAS SECUENCIAS DEL ADN

Este an´alisis asume que los tipos de nucle´otidos en diferentes posiciones son independientes. La suposici´on de independencia se utiliza para poder introducir algunas caracter´ısticas de las propiedades de los patrones en un marco simple. Como existen dependencias entre nucle´otidos en posiciones contiguas, el an´alisis general se extiende al caso de la dependencia de una cadena de Markov. Hay varias razones por las cuales debemos preguntarnos estas cuestiones. Una de ellas es que puede haber alguna raz´on “a priori” para sospechar que la palabra GAGA ocurre significativamente m´as a menudo de lo que deber´ıa si los nucle´otidos son generados en una manera iid. Para comprobar esto, es necesario tener en cuenta aspectos probabil´ısticos de la frecuencia de esta palabra asumiendo iid. Quiz´as sorprendentemente, las frecuencias de otras palabras como GGGG, GAAG, y GAGC son a menudo diferentes de los de la palabra GAGA, incluso bajo iid. Una segunda raz´on ha sido discutida por Bussemaker et al. ([26]). Aqu´ı, el objetivo es descubrir se˜ nales promotoras buscando patrones comunes de ADN sobre regiones de genes. Esto se puede hacer creando un diccionario de palabras de longitudes diferentes, cada una con una asignaci´on de probabilidad. Cualquier palabra que ocurra m´as frecuentemente que la esperada en determinadas regiones de genes, es una candidata para tal se˜ nal. El an´alisis llega a ser complicado, ya que no hay una palabra de inter´es “a priori”, y, en efecto, tampoco tenemos una longitud de palabra que pueda definirse por adelantado y buscarse espec´ıficamente. Para m´as informaci´on sobre el an´alisis de patrones podemos consultar [36], donde se estudia una f´ormula de detecci´on de secuencias patrones calculando sus frecuencias esperadas.

4.5.

B´ usqueda de repeticiones en secuencias de ADN

Describimos ahora las posibles repeticiones de secuencias de ADN presentes en un genoma eucariota y su clasificaci´on de acuerdo con su longitud y forma. Tambi´en abordaremos el problema computacional para detectarlas y los enfoques utilizados. Para m´as informaci´on y detalles v´ease [37]. En un genoma eucariota, los nucle´otidos no se encuentran en igual cantidad, A y T se encuentran en un 60 % (30 % cada uno), y C, G en un 40 % (20 % cada uno). Pero un genoma no s´olo var´ıa en la cantidad de nucle´otidos que contiene, tambi´en presenta una distribuci´on no uniforme de los nucle´otidos a lo largo de la secuencia gen´omica. Incluso existen algunos segmentos con una alta concentraci´on de un par de nucle´otidos (A, T) o (C, G) denominados isocoros. Estos isocoros son objeto de estudio ya que presentan una alta concentraci´on de secuencias codificantes. Oliver et al. ([27]) propusieron un m´etodo para detectar segmentos (isocoros o dominios que

´ 4.5. BUSQUEDA DE REPETICIONES EN SECUENCIAS DE ADN

35

presenten una composici´on similar) denominado segmentaci´on entr´opica. Adem´as de los isocoros, es interesante estudiar repeticiones de secuencias en un genoma ya que, principalmente los organismos eucariotas, contienen un alto n´ umero de secuencias repetidas (de longitud, composici´on y frecuencia variables). Estas repeticiones se pueden encontrar en los genomas de dos formas: “en tandem”, dos o m´as copias consecutivas de un patr´on, o dispersas aleatoriamente a lo largo del genoma.

4.5.1.

Repeticiones en tandem

Estas repeticiones no est´an definidas claramente todav´ıa, no obstante se pueden clasificar de la siguiente forma Clase Long - Nucle´otidos Frec - miles Sat´elites superior a 100 1000 Minisat´elites 9 - 100 100 Microsat´elites 1-8 10 - 1000 Tabla 4.3: Repeticiones en tandem. Los microsat´elites son muy variables, muy abundantes y est´an distribuidos por todo el genoma. Se utilizan principalmente para pruebas de paternidad y construcci´on de mapas gen´eticos.

4.5.2.

Repeticiones dispersas

Las repeticiones dispersas o intercaladas se clasifican de acuerdo con su longitud en cortas, denominadas SINE (short interspersed repeat) y en largas, denominadas LINE (long interspersed ). Estas repeticiones pueden estar formadas por la combinaci´on de nucle´otidos en diferente orden y cantidad. Las repeticiones SINE son fragmentos de ADN cortos repetidos millones de veces y dispersos por todo el genoma. En el genoma humano, se estima que tienen una longitud de 100 a 300 pares de bases y se repiten 1,5 millones de veces (13 % del genoma humano). Las repeticiones LINE son fragmentos de ADN de gran tama˜ no repetidos miles de veces y dispersos por todo el genoma. En el genoma humano se estima que tienen una longitud de 6,000 a 8,000 pares de bases y se repiten 85 mil veces (21 % del genoma humano).

4.5.3.

Estudio de repeticiones

La automatizaci´on en la obtenci´on de secuencias biol´ogicas ha generado un gran n´ umero de estructuras primarias de genomas y prote´ınas, las cuales demandan de herramientas computacionales para su correcta y eficiente anotaci´on. La importancia del estudio de patrones repetidos aparece desde el proceso de secuenciaci´on.

36

´ CAP´ITULO 4. ANALISIS DE LAS SECUENCIAS DEL ADN

Durante este proceso, se detect´o la presencia de bloques de ADN repetidos, lo cual dificult´o el ensamblado del genoma. El proceso de secuenciaci´on autom´atico divide el genoma en miles de fragmentos que luego se amplifican y por u ´ltimo se secuencian ensambl´andolos por superposici´on de sus extremos. Pero antes de unir dos fragmentos, el algoritmo debe determinar si los fragmentos en cuesti´on son producto de la “clonaci´on” (son el mismo fragmento) o se trata de un fragmento repetido en el genoma. Es el procedimiento “shotgun” que ya comentamos al inicio de este cap´ıtulo. En el genoma del hombre, al igual que en el de muchos organimos, se presentan secuencias de nucle´otidos repetidos (patrones o secuencias “motif”). Estos patrones generalmente tienen alguna relaci´on funcional o estructural, ya sea en la secuencia de ADN o en la prote´ına, y se cree que por esta raz´on la evoluci´on los ha preservado. Los patrones repetidos pueden ser u ´tiles en tareas como: validaci´on de m´etodos de clasificaci´on de prote´ınas, asociaci´on directa con alguna funci´on y/o estructura de una prote´ına, predicci´on o identificaci´on de dominios funcionales, o predicci´on de la estructura tridimensional. Espec´ıficamente, para el genoma humano, se estima que el 10 % est´a formando por secuencias repetidas en tandem. Estas repeticiones son probablemente resultado de una replicaci´on inexacta del ADN que con el tiempo se fijaron en el genoma. El estudio de repeticiones en tandem tiene aplicaci´on directa en medicina debido a que algunos patrones repetidos han sido asociados a enfermedades. La investigaci´on en el descubrimiento y caracterizaci´on de motifs para genomas y prote´ınas, ha producido una serie de herramientas inform´aticas, as´ı como la creaci´on de bases de datos bioinform´aticas para administrar la informaci´on generada. Entre las principales bases de datos se encuentran, por ejemplo: BLOCKS, STRBase (sobre repeticiones cortas en tandem), PROSITE (base de datos de familias de prote´ınas), Pfam, etc.

4.5.4.

El problema de la detecci´ on de patrones

La estructura primaria de un genoma se representa computacionalmente como un conjunto de archivos de texto (un archivo para cada cromosoma). Estos archivos contienen la secuencia de nucle´otidos. Cada nucle´otido se representa con una letra may´ uscula A, C, G y T , y si a´ un no se ha determinado el nucl´eotido en alguna posici´on del genoma se utiliza la letra N . El n´ umero de repeticiones a buscar se incrementa exponencialmente con la longitud del patr´on de inter´es. En la tabla inferior se ilustra el n´ umero total de combinaciones posibles de un patr´on de longitud 1 hasta 20. Para un patr´on de longitud n, existen 4n combinaciones posibles, raz´on por la cual este problema presenta un reto computacional de inter´es. La identificaci´on de patrones repetidos en un genoma puede estar dirigida a la b´ usquea de repeticiones dispersas o en tandem, con coincidencia exacta del patr´on o con presencia de “gaps”

´ LAS ISLAS CPG 4.6. APLICACION:

37

Longitud Cantidad Longitud Cantidad 1 4 11 4194304 2 16 12 16777216 3 64 13 67108864 4 256 14 268435456 5 1024 15 1073741824 6 4096 16 4294967296 7 16384 17 17179869184 8 65536 18 68719476736 9 262144 19 27487790644 10 1048576 20 1099511627776 Tabla 4.4: Cantidad de secuencias repetidas por cada longitud de un patr´on. (inserciones o deleciones) y/o mutaciones. El problema se puede plantear as´ı: descubrir(X, w, f, g, lt ) donde la funci´on“descubrir”debe encontrar una secuencia X de longitud |X| = n, el patr´on w de longitud |w| que tenga una frecuencia de aparici´on mayor o igual a f , para descartar secuencias con un bajo n´ umero de repeticiones; y presente un n´ umero de diferencias, incluyendo gaps, menor o igual a g. Si g es igual a cero, el patr´on se buscar´a con una coincidencia exacta. El par´ametro lt indica la longitud de la repetici´on en tandem. Si lt es igual a uno, la b´ usqueda se realizar´a para patrones repetidos de forma dispersa. Un m´etodo para analizar repeticiones deber ser eficiente en tiempo y espacio (lineal), que sea aplicable a la mayor cantidad posible de problemas, y que provea un an´alisis estad´ıstico y visualizaci´on integral de los resultados. Por u ´ltimo, el gran n´ umero de secuencias disponibles que incluye el genoma de varios organismos y la evidente utilidad biol´ogica de descubrir patrones en ellas, hace importante el dise˜ no de herramientas computacionales que ayuden en la caracterizaci´on de patrones de mayor longitud, de cuya existencia se tiene evidencia previa pero su descripci´on y su relevancia biol´ogica no han sido exploradas. Podemos obtener m´as informaci´on sobre el estudio de secuencias patrones en [40].

4.6.

Aplicaci´ on: las islas CpG

Veremos ahora una breve aplicaci´on o ejemplo de la teor´ıa introducida en las secciones anteriores. Se trata de observar de modo pr´actico algunas de las definiciones y estimaciones que hemos visto, en particular para el an´alisis de patrones repetidos en secuencias de ADN.

38

´ CAP´ITULO 4. ANALISIS DE LAS SECUENCIAS DEL ADN

En el genoma humano, donde ocurre el dinucle´otido CG (frecuentemente escrito como CpG para distinguirlo del par de bases C-G en dos cadenas), el nucle´otido C normalmente se modifica qu´ımicamente por un proceso de metilaci´on (adici´on de un grupo metilo (-CH3) a una mol´ecula). Existe una probabilidad bastante alta de que esta metilaci´on de C mute en una T , con la consecuencia de que, en general, los dinucle´otidos CpG son m´as raros en el genoma de lo que podr´ıa esperarse de las probabilidades independientes de C y G. Debido a importantes razones biol´ogicas, el proceso de metilaci´on se suprime en determinadas extensiones cortas del genoma, como alrededor de los promotores o regiones de inicio de muchos genes. En estas regiones, observamos m´as dinucle´otidos CpG que en ning´ un otro lugar, y de hecho m´as nucle´otidos C y G en general. Tales regiones se denominan islas CpG y, en contraste, el resto del genoma es el oc´eano. Normalmente est´an compuestas por cientos de bases de longitud. As´ı, las islas CpG conforman aproximadamente un 40 % de los genes promotores en mam´ıferos. La “p” en CpG representa que est´an enlazados por un fosfato. La definici´on formal de una isla CpG fue dada por Gardiner-Garden y Fommer ([21]), y se expresa como una regi´on con al menos 200 pares de bases, con un porcentaje de CG mayor del 50, y con un promedio del ratio CpG “observado/esperado” mayor de 0.6. Este ratio se calcula dividiendo la proporci´on de dinucle´otidos CpG en la regi´on, entre lo esperado cuando se asume independencia en una distribuci´on multinomial. La f´ormula usada es O/E =

CpG/N C/N × G/N

donde N es el n´ umero de pares de bases en la secuencia considerada. Podemos observar una secuencia de dinucle´otidos pero no sabemos a qu´e tipo de regi´on pertenece cada fragmento. Por ejemplo, en esta secuencia AACAT | {z A} CGTCCG | {z } AT | ACAT {z A}, una cuesti´on relevante es: dado un fragmento de la secuencia gen´omica, ¿c´omo podemos decidir si proviene o no de una isla CpG? Parece que el fragmento se˜ nalado en negrita, X = CGTCCG, s´ı podr´ıa provenir de una isla CpG. Otra pregunta que podr´ıamos hacernos, en general, es: dado un largo segmento de secuencia, ¿c´omo podemos encontrar islas CpG en ´el, si hay alguna? La respuesta a esta pregunta la daremos en el siguiente cap´ıtulo, con los modelos de Markov ocultos. Para poder modelizar las islas CpG, debemos saber primero qu´e tipo de peculiaridades presentan ´estas y los oc´eanos. Sabemos que 1. Hay m´as Ces y Ges en las islas (y m´as Aes y T es en los oc´eanos). 2. La probabilidad de hallar una G despu´es de un nucle´otido ser´a mayor en una isla (menos en un oc´eano) si en la posici´on actual hay una C que si no la hay.

´ LAS ISLAS CPG 4.6. APLICACION:

39

En consecuencia, y como respuesta a la primera pregunta, un modelo de Markov de orden uno puede capturar estas relaciones de dependencia. Las probabilidades de cada transici´on van a depender de si estamos en una isla CpG o no, por lo tanto, construiremos dos modelos de Markov, uno para cada caso.

Figura 4.1: Modelo de Markov para el problema de las islas CpG. Para poder estimar estas probabilidades, obtenemos de bancos de datos gen´omicos secuencias pertenecientes a islas CpG (grupo “+”) y pertenecientes a oc´eanos (grupo “−”). Las probabilidades de transici´on se estimar´an mediante m´axima veroumero de veces que el nucle´otido j sigue al nucle´otido similitud: si n∗ij representa el n´ i en una secuencia, siendo ∗ ∈ {+, −} e i, j ∈ {A, C, G, T }, las probabilidades de transici´on estimadas son pb+ ij

n+ ij =X , n+ im m

pb− ij

n− ij =X , n− im

(4.9)

m

de la misma forma que las obtenidas en la secci´on 4.1. Las siguientes tablas muestran las estimaciones de las probabilidades en dos secuencias: la primera corresponde a un modelo de isla CpG y es S1 = CGAACAGCCTCGACATGGCGTT; y el segundo, S2 = AAACAGCCTGACATGGTTC, no corresponde a una isla CpG. Evidentemente, la suma por filas en ambas tablas es uno. Isla A+ C+ G+ T+ Suma

A+ 0.20 0.29 0.33 0.00 0.82

C+ G+ T+ 0.40 0.20 0.20 0.14 0.43 0.14 0.33 0.17 0.17 0.33 0.33 0.33 1.21 1.13 0.84

Tabla 4.5: Estimaci´on de las probabilidades de transici´on en el modelo de isla CpG. Supongamos ahora que queremos puntuar una secuencia para decidir si corresponde a una isla CpG o a un oc´eano, como present´abamos en la primera pregunta. Disponemos de dos modelos, como antes indicamos, el modelo “+” de las islas CpG y el modelo “−” de los oc´eanos. La idea subyacente tras el sistema de puntuaciones es

40

´ CAP´ITULO 4. ANALISIS DE LAS SECUENCIAS DEL ADN Oc´eano A− C− G− T− Suma

A− 0.33 0.50 0.25 0.00 1.08

C− G− T− 0.33 0.17 0.17 0.25 0.00 0.25 0.25 0.25 0.25 0.25 0.50 0.25 1.08 0.92 0.92

Tabla 4.6: Estimaci´on de las probabilidades de transici´on en el modelo de oc´eano. Al basarse en una secuencia corta aparece un cero en la transici´on C → G. 1. Si la secuencia pertenece a una isla CpG, tendr´a una probabilidad m´as alta sobre el modelo “+” que sobre el “−”. 2. Si la secuencia no es de una isla CpG, la probabilidad que le asignar´a el modelo “−” ser´a mayor. En lugar de multiplicar las probabilidades, sumaremos los logaritmos de las razones de probabilidades seg´ un cada modelo, y calcularemos el logg-odds ratio de la siguiente forma   n Y ! pb+    xa−1 xa  n n X  X  pb+ P (X|+) x x a=1 a a−1   log βxa−1 xa , S(X) = log = = log  n = P (X|−) pb− xa−1 xa  a=1 Y − a=1 pbxa−1 xa a=1

(4.10) donde xa denota el nucle´otido en la posici´on a de la secuencia X, para a = 1, . . . , n; y n es la longitud de la secuencia dada. La decisi´on de si la secuencia es o no una isla CpG depender´a de que los valores de la suma sean m´as o menos altos. En la siguiente tabla se muestra un ejemplo para i, j ∈ E. − βij = log(p+ A C G T ij /pij ) A -0.51 0.18 0.18 0.18 C -0.56 -0.56 1.00 -0.56 G 0.29 0.29 -0.41 -0.41 T 1.00 0.29 -0.41 0.29

Tabla 4.7: Matriz de puntuaciones basada en los modelos de las tablas anteriores. Al tratarse en secuencias cortas, la transici´on C → G se puntuar´a con un 1. Deber´ıan tomarse m´as valores para mejorar la estimaci´on. Ante estas puntuaciones que podr´ıamos considerar altas (menos valores negativos que positivos), dir´ıamos que la secuencia s´ı corresponde a una isla CpG. En este caso,

´ LAS ISLAS CPG 4.6. APLICACION:

41

dada la secuencia X = CGTCCG, se tiene una puntuaci´on S(X) = 1 + (−0.41) + 0.29 + (−0.56) + 1 = 1.32. Esta puntuaci´on normalizada por la longitud de X resulta 1.32 = 0.22. Concluimos entonces, dado que la puntuaci´on es positiva, que s´ı se trata 6 de una isla CpG, como hab´ıamos intuido.

42

´ CAP´ITULO 4. ANALISIS DE LAS SECUENCIAS DEL ADN

Cap´ıtulo 5 Los modelos de Markov ocultos 5.1.

Introducci´ on

Llegamos ya a nuestro principal objeto de estudio: los Modelos de Markov Ocultos. Utilizaremos toda la informaci´on vista en las dos partes anteriores para comprender mejor la teor´ıa general sobre estos modelos, y sus aplicaciones pr´acticas. Nuestro estudio sobre los modelos de Markov ocultos contendr´a 1. Una definici´on de estos modelos y descripci´on de sus propiedades. 2. El principal algoritmo de estimaci´on de par´ametros del modelo. 3. Una aplicaci´on y un caso pr´actico Por u ´ltimo, introduciremos algunas de las nuevas aplicaciones de los modelos de Markov ocultos en diferentes problemas relacionados con la Biolog´ıa y la Bioinform´atica. Antes de comenzar con la definici´on, veamos una breve introducci´on a su historia y una idea intuitiva de dichos modelos. Podemos decir, de manera general, que los modelos de Markov ocultos (HMMs son las siglas en ingl´es para “Hidden Markov Models”) son modelos matem´aticos que capturan informaci´on oculta en secuencias de s´ımbolos u observaciones (nucle´otidos en nuestro caso). Los modelos de Markov ocultos, o simplemente los modelos de Markov, nos ayudan a responder a las siguientes preguntas: “dada una secuencia de ADN, ¿pertenece a una familia particular?” o “asumiendo que la secuencia pertenece a una determinada familia, ¿qu´e podemos saber de su estructura interna?”. Un ejemplo del segundo tipo de problema podr´ıa ser tratar de identificar determinadas regiones de una secuencia de prote´ına. La mayor´ıa de los escritos sobre los modelos de Markov ocultos pertenecen a la literatura del reconocimiento ortogr´afico ([31]), donde fueron aplicados por primera 43

44

CAP´ITULO 5. LOS MODELOS DE MARKOV OCULTOS

vez a principios de 1970. En este campo, la aplicaci´on de los modelos de Markov ocultos trata b´asicamente de, dada una se˜ nal ortogr´afica ac´ ustica (fonemas), reconocer la secuencia de palabras (ocultas) que han emitido tal se˜ nal, es decir, que se han dicho. Muchos problemas en el an´alisis de secuencias biol´ogicas tienen la misma estructura: basados en una secuencia de s´ımbolos de un determinado alfabeto, encontrar qu´e secuencia representa. Para las prote´ınas, las secuencias consisten en s´ımbolos del “alfabeto” de los 20 amino´acidos; y t´ıpicamente nosotros queremos saber a qu´e familia de prote´ınas pertenece la secuencia dada. Aqu´ı, la secuencia primaria de amino´acidos es an´aloga a la se˜ nal ortogr´afica, y la familia de prote´ınas a la palabra dicha que representa. Nuestro alfabeto, como ya sabemos, consta de las iniciales de las cuatro bases de nucle´otidos A, C, G y T . En 1989, Gary Churchill ([22]) introdujo el uso de los modelos de Markov ocultos para la segmentaci´on de secuencias de ADN. La utilizaci´on de estos modelos por Churchill, le permiti´o segmentar una secuencia de ADN alternando regiones de nucle´otidos similares. Desde entonces, los modelos de Markov ocultos se han aplicado cada vez m´as a la gen´omica, incluyendo la b´ usqueda de genes y la predicci´on de funciones de las prote´ınas, sobre todo debido al trabajo pionero de David Haussler, de la universidad de California, Santa Cruz. Hoy en d´ıa, con los m´etodos de alineamiento, los HMMs se encuentran entre los algoritmos m´as representativos en el a´rea de la Bioinform´atica. La necesidad b´asica de los HMMs surge del hecho de que los datos del genoma est´an inherentemente desordenados. Incluso en regiones con un alto contenido en CG, por ejemplo, puede haber largas extensiones de Aes y T s. Los patrones observados en secuencias de ADN son necesariamente una representaci´on “robusta” del estado oculto del genoma (“gran contenido de CG” podr´ıa ser uno de esos estados). Las cadenas de Markov simples no son lo suficientemente flexibles para capturar muchas propiedades de las secuencias de ADN. Los modelos de Markov ocultos las generalizan en un simple y eficiente esquema. La idea b´asica que se encuentra detr´as de la aplicaci´on de los modelos de Markov ocultos es la de modelizar una secuencia como si hubiese sido indirectamente generada por una cadena de Markov. Cada posici´on en la secuencia tiene un estado desconocido (oculto) de la cadena de Markov, pero podemos observar los s´ımbolos generados de acuerdo a una distribuci´on multinomial que depende de ese estado. En otras palabras, la informaci´on que recibimos sobre la cadena de Markov oculta es indirecta. La secuencia que tratamos de analizar es por lo tanto modelizada siendo el resultado de un proceso doblemente estoc´astico: uno generando una cadena de Markov oculta, y otro transformando esta cadena oculta en una secuencia observable. Este segundo proceso sigue una distribuci´on multinomial: en cada estado, se utiliza un conjunto diferente de par´ametros para producir la secuencia observada. Una de las claves de los HMMs es tomar esta secuencia observada e inferir o dedu´ cir los estados ocultos. Estos pueden representar diferentes tipos de secuencias. Los

´ 5.1. INTRODUCCION

45

modelos de Markov ocultos simples solo tienen dos estados, como “ricos en CG” o “ricos en AT ”, mientras que los modelos ocultos m´as complejos pueden tener m´as estados, como “A”, “C”, “G”, “T ”, “regi´on codificante” o “intr´on”. Los modelos de Markov ocultos tienen dos importantes par´ametros: las probabilidades de transici´on y las probabilidades de emisi´on. El par´ametro de transici´on describe la probabilidad de que la cadena de Markov cambie a lo largo de varios estados ocultos. Estos cambios pueden ocurrir muy a menudo o de forma escasa; ya que la cadena puede hacer la transici´on entre solo dos estados o entre muchos estados. El par´ametro de emisi´on describe las probabilidades de que los s´ımbolos en la secuencia observada se produzcan en cada uno de los diferentes estados. Cada uno de los estados ocultos debe poder producir los mismos s´ımbolos, aunque en diferentes frecuencias. El n´ umero de s´ımbolos emitidos puede variar entre dos o muchos; el n´ umero m´as com´ un en an´alisis de secuencias es 4 (para ADN y ARN) o 20 (para los amino´acidos). Vamos ahora a formalizar la notaci´on para los modelos de Markov ocultos. Necesitamos ahora distinguir la secuencia de estados de la secuencia de s´ımbolos. Llamaremos a la secuencia de estados el camino π. El propio camino sigue una cadena de Markov simple, por lo tanto, la probabilidad de un estado depende solamente del estado anterior. El estado en la posici´on a del camino se denota por πa . La cadena est´a caracterizada por los par´ametros puv = P (πa = v|πa−1 = u), que indican las probabilidades de transici´on entre estados contiguos del camino, para a = 1, . . . , n, donde n es la longitud de la secuencia de observaciones que queremos analizar; y u, v son posibles estados de la cadena de Markov oculta. Para modelizar el comienzo del proceso introducimos un estado inicial I. La probabilidad de transici´on pI,u del estado inicial al estado u puede pensarse como la probabilidad (0) de empezar en el estado u, pu , para cualquier u. Es tambi´en posible modelizar los finales, terminando una secuencia de estados con una transici´on en un estado final F , pu,F . Por convenio, tomaremos la transici´on al estado final como cero. Como hemos desacoplado los s´ımbolos i de los estados u, debemos introducir un nuevo conjunto de par´ametros para el modelo, eu (i). En general, un estado puede producir un s´ımbolo sobre todos los s´ımbolos posibles. Por lo tanto, definimos eu (i) = P (xa = i|πa = u), la probabilidad de que se observe el s´ımbolo i cuando nos encontramos en el esta´ do u. Estas se conocen como las probabilidades de emisi´on. La raz´on por la cual las probabilidades de emisi´on se llaman de esta forma es porque a veces es conveniente pensar los modelos de Markov ocultos como modelos generativos, que generen o emitan secuencias. Puede generarse una secuencia de un HMM como sigue: primero se

46

CAP´ITULO 5. LOS MODELOS DE MARKOV OCULTOS (0)

elige un estado π1 de acuerdo a las probabilidades pπa , a = 1, . . . , n. En este estado, se emite una observaci´on de acuerdo a la distribuci´on eπ1 para ese estado. Despu´es, se elige un nuevo estado π2 en concordancia con las probabilidades de transici´on pπ1 πa , a = 1, . . . , n; y as´ı sucesivamente. De esta forma, se genera una secuencia de observaciones aleatorias. En consecuencia, cuando el HMM trabaja hay, en primer lugar, una secuencia de estados visitados, los cuales denotamos por π1 , π2 , π3 , . . .; y en segundo lugar, una secuencia de s´ımbolos emitidos, denotados por x1 , x2 , x3 , . . .. Denotamos la secuencia entera de los πa por Π y la secuencia de los xa por X, y escribimos la secuencia observada como X = x1 , x2 , . . .; y la secuencia de estados como Π = π1 , π2 , . . .. Supongamos que cada s´ımbolo xa toma un valor del conjunto de observaciones, que es finito de tama˜ no k (en nuestro caso k = 4 para A, C, G y T ); y que cada estado πa toma uno de los valores del conjunto de estados de tama˜ no S. Ahora, es f´acil escribir la probabilidad conjunta de una secuencia observada X y una secuencia de estados Π n Y (0) P (X, Π) = pπ1 eπa (xa )pπa πa+1 , (5.1) a=1

donde exigimos que pπn πn+1 = 0. Sin embargo, esta expresi´on no es u ´til en la pr´actica porque normalmente no conocemos el camino oculto (secuencia de estados oculta). M´as adelante, describiremos un algoritmo para estimar el camino encontrando la secuencia de estados m´as probable u o´ptima para el modelo de Markov oculto, llamado algoritmo de Viterbi. Antes de comenzar con el estudio del algoritmo en cuesti´on, conviene se˜ nalar que existen un gran n´ umero de variantes de los modelos de Markov ocultos que modifican y extienden el modelo b´asico, para satisfacer las necesidades de diferentes aplicaciones. Por ejemplo, podemos a˜ nadir estados silenciosos (es decir, estados que no emiten ning´ un s´ımbolo) al modelo para representar la falta de ciertos s´ımbolos que se espera que est´en presentes en posiciones espec´ıficas. Podemos tambi´en hacer que los estados emitan dos s´ımbolos en lugar de un u ´nico s´ımbolo y, por lo tanto, el resultante HMM genere simult´aneamente dos secuencias de s´ımbolos relacionadas. Otra variante de los modelos de Markov ocultos, adecuada para la b´ usqueda de genes, se denomina modelo de Markov semi-oculto o generalizado. En este caso, cada estado oculto emite una secuencia de observaciones, en lugar de una u ´nica observaci´on o s´ımbolo. Veremos una breve descripci´on de este modelo con la aplicaci´on en b´ usqueda de genes. Por otro lado, es tambi´en posible hacer que las probabilidades de ciertos estados dependan de parte de las emisiones previas y as´ı podr´ıamos describir de forma m´as compleja correlaciones entre s´ımbolos. Veamos un ejemplo cl´asico y muy sencillo de c´omo aplicar un modelo de Markov oculto con dos estados, para introducir una idea b´asica e intuitiva de ´este.

´ 5.1. INTRODUCCION

47

Ejemplo: el casino deshonesto. Imaginemos que en un casino utilizan un dado libre o justo (correcto) la mayor´ıa de las veces, pero otras utilizan un dado cargado. Tenemos entonces un par de dados, uno que es justo o normal, y otro que est´a cargado (es decir, no caemos en todas las caras con igual probabilidad). De esta forma, nuestros estados ser´an: {“dado libre”, “dado cargado” }. Un modelo de Markov decide cu´al de los dos dados juega, y dependiendo del estado del modelo, se aplicar´an las probabilidades de emisi´on para el dado cargado o para el normal. As´ı, generamos una secuencia de s´ımbolos que es el resultado de estar en dos estados diferentes. Especificaremos nuestro par´ametro de transici´on de modo que en cualquiera de los dos estados haya un 90 % posibilidades de quedarse en ese estado, y un 10 % de cambiar de estado (el casino cambia un dado libre por uno cargado), como se muestra en la figura inferior.

Figura 5.1: Modelo de Markov oculto asociado al ejemplo. Las transiciones entre el dado cargado y el dado libre se modelizan con una cadena de Markov, y los lanzamientos como emisiones independientes de un modelo multinomial. Para nuestro dado no cargado o libre, la probabilidad de obtener un n´ umero entre 1 y 6 (ambos incluidos) es la misma, y est´a dada por 0.1667 0.1667 0.1667 0.1667 0.1667 0.1667, donde cada columna describe la probabilidad para cada uno de los seis n´ umeros. Para el dado cargado, las probabilidades de emisi´on son 0.1000 0.1000 0.1000 0.1000 0.1000 0.5000. Aqu´ı, las probabilidades de emitir un n´ umero entre 1 y 5 son iguales, pero es mucho mayor la probabilidad de obtener un 6 (50 %).

48

CAP´ITULO 5. LOS MODELOS DE MARKOV OCULTOS

La secuencia visible producida por tal modelo de Markov oculto podr´ıa ser como la siguiente X = 4553653163363555133362665132141636651666. Si conocemos las propiedades de los dos dados y de la cadena de Markov subyacente u oculta, ¿podemos encontrar la secuencia m´as probable de estados ocultos detr´as de X? En otras palabras, ¿podemos averiguar qu´e dado ha sido utilizado por el casino en cada instante o posici´on de la secuencia de resultados? M´as adelante, con el algoritmo de Viterbi, veremos c´omo responder a estas preguntas, aqu´ı simplemente mostramos la secuencia oculta que genera nuestra secuencia visible Oculto : Π = 1111111111111111111122221111111222222222. V isible : X = 4553653163363555133362665132141636651666. Observemos que el s´ımbolo “6” ocurre con una frecuencia mayor cuando el estado en la secuencia oculta es 2, correspondiente al dado cargado, pero esta dependencia es simplemente probabil´ıstica. En las aplicaciones biol´ogicas, aplicaremos nuestro modelo de Markov oculto en un conjunto de datos donde los estados ocultos son conocidos, incluso cuando no sepamos exactamente cuales son las probabilidades de transici´on y emisi´on. Esto permite calcular las matrices de transici´on y emisi´on de nuestro modelo basadas en los datos y, por lo tanto, una mejor inferencia en los estados ocultos de nuevos datos. Este aspecto especialmente u ´til en la b´ usqueda de genes, donde los modelos se aplican a genomas conocidos para despu´es extenderlos a otros genomas semejantes. Seguimos ahora con el ejemplo de las islas CpG dado en el cap´ıtulo anterior (secci´on 4.6), para estudiar su an´alisis con los modelos de Markov ocultos. Aplicaci´ on: las islas CpG. Recordemos la segunda pregunta planteada relativa al an´alisis de las islas CpG: ¿c´omo podemos encontrarlas, si las hay, en una larga secuencia no conocida? Los modelos de cadenas de Markov ya descritos pueden utilizarse para responder esta cuesti´on, calculando las puntuaciones log-odds de un conjunto o “ventana” de, por ejemplo, 100 nucle´otidos “alrededor” de cada nucle´otido en la secuencia. As´ı, las puntuaciones altas indicar´ıan islas CpG potenciales. Sin embargo, esto no nos vale ya que, de hecho, las islas CpG tienen unos l´ımites definidos y son de longitud variable. Por lo tanto, el problema ser´ıa saber qu´e taman ˜o de ventana hay que usar. Una aproximaci´on m´as satisfactoria ser´ıa construir un modelo simple de la secuencia entera, no s´olo 100 nucle´otidos, que incorpore ambas cadenas de Markov, el modelo “+” de las islas CpG, y el modelo “−” del resto de regiones (oc´eano). Para ello, utilizaremos un modelo de Markov oculto. Para simular en un mismo modelo las “islas” en un “oc´eano”, queremos tener ambos modelos de Markov presentes en uno solo, con una peque˜ na probabilidad de cambiar de un modelo al otro en cada punto de transici´on. Sin embargo, esto

5.2. EL ALGORITMO DE VITERBI

49

introduce las posibles dificultades de tener dos estados con el mismo nombre que correspondan a cada nucle´otido. Resolveremos este problema etiquetando de nuevo los estados. Ahora tenemos los estados A+ , C+ , G+ y T+ en las regiones de islas CpG, y los estados A− , C− , G− y T− correspondientes a las regiones de oc´eano. Cada uno de los estados de ambos modelos puede emitir uno de los cuatro s´ımbolos de nuestro alfabeto A, C, G y T . De esta forma, adem´as de las probabilidades de transici´on entre los dos grupos (isla y oc´eano), hay tambi´en un conjunto de transiciones en cada uno de los grupos, como en los modelos de Markov simples. Las probabilidades de transici´on en este modelo est´an establecidas de modo que los elementos dentro de cada grupo tienen probabilidades de transici´on mayores. En general, es m´as probable cambiar de “+” a “−” que al rev´es (hay m´as regiones de oc´eano que de isla en el genoma), por tanto, si lo recorremos libremente, el modelo pasar´a m´as tiempo en los estados “−” de oc´eano que en los estados de islas. El etiquetado es el paso importante. Las probabilidades de emisi´on en este caso son todas 0 o 1, dependiendo de si el s´ımbolo emitido coincide con el estado (tanto positivo como negativo) o no. Por ejemplo, P (xa = C|πa = C− ) = 1 y P (xa = C|πa = A+ ) = 0. As´ı, siguiendo la ecuaci´on 5.1, podemos calcular la probabilidad de la secuencia CGCG que ha sido emitida por la secuencia de estados (C+ , G− , C− , G+ ) (0)

pC+ × 1 × pC+ ,G− × 1 × pG− ,C− × 1 × pC− ,G+ × 1 × 0. La diferencia esencial entre una cadena de Markov y un modelo de Markov oculto es que para este u ´ltimo no hay una correspondencia entre los estados y los s´ımbolos. Es decir, no es posible, a la larga, saber en qu´e estado est´a el modelo cuando se genera el s´ımbolo xa , s´olo observando xa . En nuestro ejemplo, no hay forma posible de saber, observando s´olo un s´ımbolo C, si fue emitido por un estado C+ o C− . Por ello, la f´ormula 5.1 no es v´alida para nuestro inter´es, y tendremos que estudiar un algoritmo que estime la secuencia de estados ocultos.

5.2.

El algoritmo de Viterbi

A pesar de que no es posible saber en qu´e estado est´a el sistema observando los s´ımbolos correspondientes, s´ı podemos estimar la secuencia de estados oculta. Describimos a continuaci´on el algoritmo de programaci´on din´amica m´as com´ un para poder hallarla, denominado algoritmo de Viterbi. En general, es posible que haya muchas secuencias de estados que pueden dar lugar a una particular secuencia de s´ımbolos. Por ejemplo, en nuestro modelo de islas CpG, las secuencias de estados (C+ , G+ , C+ , G+ ), (C− , G− , C− , G− ) y (C+ , G− , C+ , G− ) podr´ıan generar la secuencia de s´ımbolos CGCG. Sin embargo, lo hacen con diferentes probabilidades. La tercera es el producto de m´ ultiples probabilidades alternando entre los dos modelos (del modelo “+” al “−”, del “−” al “+”, y as´ı sucesivamente), y por lo tanto es mucho m´as peque˜ na que la probabilidad de las dos

50

CAP´ITULO 5. LOS MODELOS DE MARKOV OCULTOS

primeras (las probabilidades de transici´on entre modelos son menores). La segunda es significativamente m´as peque˜ na que la primera, ya que contiene dos transiciones de C a G que son significativamente menos probables en el caso del modelo “−” de oc´eanos, que en el modelo “+” de islas. Aplicaremos el algoritmo de Viterbi para hallar la sucesi´on de estados m´as probable de haberla generado. La estimaci´on de un camino a trav´es de un modelo de Markov oculto nos dir´a qu´e parte de la secuencia se predice por una isla CpG, ya que anteriormente asumimos que cada estado se asignaba al modelo de las islas CpG o al del oc´eano. Si tenemos que elegir s´olo un camino, el que tenga la mayor probabilidad deber´a ser el escogido. De entre todas las posibles secuencias de estados Π, queremos encontrar la secuencia de estados que explique mejor la secuencia de s´ımbolos observada. Dicho de otra forma, queremos encontrar “el mejor alineamiento” entre la secuencia de s´ımbolos y el HMM, por lo que a veces este procedimiento se denomina problema del alineamiento ´optimo o decodificaci´on. Formalmente, dada una secuencia observada X = x1 , x2 , x3 , ..., xR de salidas u observaciones, queremos computar eficientemente una secuencia de estados Π = π1 , π2 , π3 , ..., πR que tenga la mayor probabilidad condicionada dada X (tomamos ahora R en lugar de n como longitud de la secuencia para facilitar la notaci´on). En otras palabras, queremos encontrar el camino o´ptimo Π∗ que satisfaga lo siguiente Π∗ = argmaxΠ P (Π|X).

(5.2)

Encontrar la secuencia de estados ´optima Π∗ comparando todas las S R (S el el tama˜ no del conjunto formado por todos los estados) posibles secuencias de estados es computacionalmente inviable. Sin embargo, podemos utilizar un algoritmo de programaci´on din´amica, llamado algoritmo de Viterbi, para encontrar una de las muchas secuencias de estados que pueden maximizar la probabilidad anterior. El algoritmo se divide en dos partes: primero encuentra el maxΠ P (Π|X), y despu´es realiza un procedimiento “hacia atr´as” para encontrar la secuencia Π que satisfaga este m´aximo. En primer lugar se define, para r y u arbitrarios, δr (u) = maxπ1 ,π2 ,...,πr−1 P (π1 , π2 , ..., πr−1 , πr = u; x1 , x2 , ..., xr ),

(5.3)

(δ1 (u) = P (π1 = u; x1 )). En otras palabras, δr (u) es la m´axima probabilidad de todas las formas de terminar en el estado u en el instante o posici´on r; y haber observado la secuencia x1 , x2 , ..., xr . Entonces maxΠ P (Π, X) = maxu δR (u).

(5.4)

La probabilidad en esta expresi´on es la probabilidad conjunta de Π y X, no una probabilidad condicionada. Nuestro objetivo es encontrar la secuencia Π para la cual se alcance la m´axima probabilidad condicionada. Ya que maxΠ P (Π|X) = maxΠ

P (Π, X) , P (X)

(5.5)

5.2. EL ALGORITMO DE VITERBI

51

y ya que el denominador del lado derecho de la igualdad no depende de Π, se tiene que argmaxΠ P (Π|X) = argmaxΠ

P (Π, X) = argmaxΠ P (Π, X). P (X)

(5.6)

El primer paso es calcular los δr (u) inductivamente en r. Despu´es, con el procedimiento “hacia atr´as”, recuperaremos la secuencia que de el mayor δR (u). El paso de inicializaci´on es δ1 (u) = p(0) 1 ≤ u ≤ S. (5.7) u eu (x1 ), teniendo en cuenta que δ0 (0) = 1 y δr (0) = 0, para todo r > 0. El paso de inducci´on es δr (v) = max1≤u≤S {δr−1 (u)puv }ev (xr ),

2 ≤ r ≤ R, 1 ≤ v ≤ S.

(5.8)

Para obtener la observaci´on con la m´axima probabilidad, recuperamos los πa como sigue. Se define ψR = argmax1≤u≤S δR (u). (5.9) Y sea πR = ψR . Entonces, πR es el estado final en la secuencia de estados requerida. Los restantes πr para r ≤ R − 1 se encuentran recursivamente definiendo primero ψr = argmax1≤u≤S {δr (u)puψr+1 }. (5.10) Y despu´es escribiendo πr = ψr . Si el m´aximo no es u ´nico, tomaremos arbitrariamente un valor de u. El algoritmo de Viterbi encuentra la secuencia de estados ´optima con una eficiencia del orden de O(RS 2 ). Veremos c´omo obtener el camino oculto en un problema de predicci´on de islas CpG utilizando el algoritmo de Viterbi en el siguiente cap´ıtulo, donde analizaremos un ejemplo con R. Existen otros dos algoritmos que resuelven problemas relacionados con un HMM, como el problema de puntuaci´on o evaluaci´on y el problema de estimaci´on. El primero de ellos, consiste en calcular la probabilidad P (X), dada la secuencia observada X = x1 , . . . , xR , basada en un modelo de Markov oculto. Un m´etodo eficiente para calcular dicha probabilidad est´a basado en un algoritmo de programaci´on din´amica, llamado algoritmo “hacia delante” (existe un an´alogo llamado algoritmo “hacia atr´as”). Este algoritmo tiene un tiempo de computaci´on del orden de O(RS 2 ). Si queremos clasificar una nueva prote´ına, podemos construir un HMM para cada familia conocida de prote´ınas y utilizarlo para calcular la probabilidad de la nueva secuencia. Asignaremos la prote´ına a la familia en cuyo modelo se obtenga mayor probabilidad. En cuanto al segundo problema, se refiere a c´omo podemos elegir los

52

CAP´ITULO 5. LOS MODELOS DE MARKOV OCULTOS

par´ametros de un HMM bas´andose en un conjunto de secuencias observadas. No existe una manera o´ptima de estimar estos par´ametros de un n´ umero limitado de secuencias de observaciones finitas, pero s´ı existen formas de localizar los par´ametros del modelo que maximicen localmente la probabilidad de la secuencia observada. Por ejemplo, el algoritmo de Baum-Welch estima y actualiza de forma iterativa el conjunto de par´ametros del modelo basado en un procedimiento “hacia delante-hacia atr´as”. Existen muchos m´etodos para la estimaci´on de par´ametros de un HMM, ya que en el fondo se trata de un problema de optimizaci´on. No profundizaremos en el estudio de ambos problemas y algoritmos ([40], [10]). Ejemplo: el casino deshonesto Veamos c´omo aplicar el algoritmo de Viterbi a nuestro ejemplo del casino deshonesto. Recordemos que la probabilidad de quedarse en cualquiera de los dos estados (dado libre o cargado) es 0.9, y la probabilidad de cambiar de estado es 0.1. Adem´as, la probabilidad de obtener un n´ umero en el dado libre es siempre la misma ( 16 ), mientras que en el dado cargado la probabilidad de que salga un 6 es 5 veces mayor que la del resto. Con estos datos, y teniendo en cuenta el diagrama presentado en la figura 5.1, veamos c´omo aplicar el algoritmo de Viterbi a una secuencia de lanzamientos, como por ejemplo X = 6 2 6, correspondiente a 3 lanzamientos. El paso inicial es 1 0 0. Despu´es, para el s´ımbolo 6 se tiene 1 1 × 16 = 12 2 1 × 12 = 14 = 0.25 2

siguiendo con el siguiente n´ umero obtenido, el 2, se obtiene 1 6 1 10

1 × max{( 12 × 0.9, 41 × 0.1) = 0.0125 1 × max{( 12 × 0.1, 41 × 0.9) = 0.0225

y finalmente aplicamos el algoritmo para el u ´ltimo s´ımbolo, de nuevo el 6. 1 × max{(0.0125 × 0.9, 0.0225 × 0.1) = 0.000375 6 1 × max{(0.0125 × 0.1, 0.0225 × 0.9) = 0.005625 2

Las primeras filas representan los c´alculos para el estado “dado libre” y las segundas para el estado “dado cargado”. As´ı, tomando los mayores valores de δ, que se presentan en negrita, obtenemos que el camino o´ptimo es Π = {222} para el conjunto de s´ımbolos X = 6 2 6 (el estado 2 hace referencia al dado cargado). Se proceder´ıa de forma an´aloga para un conjunto de s´ımbolos mayor, es decir, para m´as lanzamientos de los dados. Ejemplo: reconocimiento de regiones en secuencias de ADN. Otro ejemplo sencillo con tres estados y aplicado al estudio del ADN es el siguiente. Imaginemos

5.2. EL ALGORITMO DE VITERBI

53

que tenemos un problema de reconocimiento de una regi´on o posici´on de empalme 5’ (“sitio donor 5”), en una secuencia gen´etica. Asumiremos que nos dan la secuencia de ADN que comienza con un ex´on, contiene un lugar de empalme 5’ y termina con un intr´on. El problema es identificar donde se produce el cambio del ex´on al intr´on, es decir, donde est´a el lugar de empalme 5’. Sabemos que las secuencias de exones, los lugares de empalme y los intrones tienen diferentes propiedades estad´ısticas. Imaginemos algunas diferencias simples: supongamos que los exones tienen, en promedio, una composici´on de base uniforme (25 % cada base), los intrones son ricos en Aes y T s (es decir, 40 % cada una para A y T , 10 % cada una para C y G), y los lugares de empalme 5’ tienen como nucle´otido mayoritario a G (95 % de Gs y 5 % de Aes). Con esta informaci´on, podemos dibujar el modelo de Markov oculto asociado, como se muestra en la figura inferior.

Figura 5.2: Ejemplo de un modelo de Markov oculto para el reconocimiento de una regi´on de empalme 5’. El modelo tiene tres estados (m´as un estado inicial y otro final), uno para cada uno de los tres niveles que podemos asignar a un nucle´otido: E (ex´on), 5 (empalme 5’) e I (intr´on). As´ı, el conjunto de estados es Π = {E, 5, I} y el conjunto de observaciones son los cuatro nucle´otidos X = {A, C, G, T }. Cada estado tiene sus propias probabilidades de emisi´on (mostradas debajo de los estados), que modelan la composici´on b´asica de los exones, intrones y las Gs en los lugares de empalme. Cada estado tiene tambi´en sus probabilidades de transici´on (flechas): las probabilidades de moverse de ese estado a un nuevo estado, desde el estado inicial y hacia el estado final (en este caso las probabilidades de transici´on de los estados inicial y final son no nulas). Las probabilidades de transici´on describen el orden lineal en el cual esperamos que los estados ocurran: uno o m´as Es, un 5, uno o m´as Ies. Consideremos la secuencia mostrada en la figura, compuesta por 26 nucle´otidos, y el camino de estados correspondiente, donde hay 27 transiciones y 26 emisiones en total. Multiplicando todas las 53 probabilidades juntas (y tomando el logaritmo, ya que se trata de n´ umeros peque˜ nos), podemos calcular la P (X, Π) = −41.22.

54

CAP´ITULO 5. LOS MODELOS DE MARKOV OCULTOS

Existen otros 14 posibles caminos que no tienen probabilidad cero, ya que el lugar de empalme 5’ debe estar en una de las 14 Aes o Gs internas. En la figura se enumeran los seis caminos con las mayores puntuaciones (con G en el lugar de empalme). El mejor tiene el valor de −41.22 para el logaritmo de la probabilidad, lo que indica que la posici´on m´as probable para el empalme es el quinto G. El camino alternativo difiere solo ligeramente de ´este, −41.71 frente a −41.22. Teniendo en cuenta esta leve diferencia, ¿c´omo de seguros estamos de que la quinta G es la elecci´on correcta? Aqu´ı se muestra la ventaja del modelado probabil´ıstico: podemos calcular nuestra fiabilidad directamente. La probabilidad de que un estado emita un amino´acido i es la suma de las probabilidades de todos los caminos que utilizan ese estado para generar tal amino´acido, normalizada por la suma de todos los posibles caminos. En este ejemplo, solo hay un camino de estados en el numerador, y una suma de 14 caminos en el denominador. Obtenemos una probabilidad del 46 % de que la mejor puntuaci´on de la quinta G es correcta, y un 28 % de que la sexta G es correcta (parte inferior de la figura). A ´esto se le denomina decodificaci´on posterior. Para problemas grandes, se utilizan algoritmos de programaci´on din´amica, como el algoritmo forward o “hacia deltante”.

5.3.

Aplicaciones: b´ usqueda de genes

Hoy en d´ıa, se producen ya secuencias de genes con longitudes del orden de millones de bases. Tales secuencias consisten en la colecci´on de genes separados unos de otros por largas extensiones de secuencias no funcionales. Es de vital importancia encontrar donde est´an los genes en la secuencia y, por lo tanto, son muy u ´tiles m´etodos computacionales que identifiquen r´apidamente una gran proporci´on de los genes. El problema implica tomar al mismo tiempo una gran cantidad de diversa informaci´on, y se han estudiado muchas aproximaciones para poder hacerlo. Actualmente, un exitoso y popular buscador de genes para secuencias de ADN humano es el GENSCAN ([9]), que est´a basado en una generalizaci´on de los modelos de Markov ocultos. Expondremos a continuaci´on, de manera breve, un algoritmo similar con el esp´ıritu del GENSCAN para ilustrar los conceptos b´asicos de un HMM en la b´ usqueda de genes humanos. Para incrementar la precisi´on del procedimiento es necesario introducir algunos detalles que no describiremos aqu´ı ([40]). Adem´as de la herramienta GENSCAN, existen muchos m´as sistemas para la predicci´on de genes, tanto en organismos eucariotas como procariotas. En la siguiente direcci´on web se describen algunos de los m´as utilizados y precisos: http://cmgm.stanford.edu/classes/genefind/. Antes de comenzar con la introducci´on al algoritmo en el que se basa la herramienta GENSCAN, veamos un ejemplo intuitivo que nos muestra la idea y objetivo del posterior estudio.

´ 5.3. APLICACIONES: BUSQUEDA DE GENES

55

Ejemplo: modelos de Markov ocultos para genes eucariotas. Los modelos de Markov ocultos pueden utilizarse efectivamente para representar secuencias biol´ogicas, como ya hemos visto. Como un ejemplo simple, consideremos un modelo de Markov oculto que modele genes codificantes en prote´ınas de organismos eucariotas. Es conocido que muchas regiones codificantes a prote´ınas muestran tendencias entre codones. El uso no uniforme de los codones resulta en diferentes s´ımbolos estad´ısticos para diferentes posiciones de los codones. Estas propiedades no se observan en los intrones, ya que no se traducen en amino´acidos. Por lo tanto, es importante incorporar estas caracter´ısticas de los codones a la hora de modelar genes codificantes en prote´ınas y construir un buscador de genes. La figura inferior muestra un ejemplo de un HMM para modelizar genes eucariotas. El modelo de Markov oculto dado, trata de capturar las diferencias estad´ısticas en exones e intrones. El modelo tiene cuatro estados, donde E1 , E2 , y E3 se utilizan para modelar las propiedades estad´ısticas b´asicas en los exones. Cada estado Ea utiliza un conjunto de probabilidades de emisi´on diferentes para reflejar el s´ımbolo en la posici´on a de un cod´on. El estado I se utiliza para modelizar los intrones. Observemos que este HMM puede representar genes con m´ ultiples exones, donde los respectivos exones pueden tener un n´ umero variable de codones; y los intrones pueden tener tambi´en longitudes variables. Este ejemplo muestra que si conocemos la estructura y las caracter´ısticas importantes de la secuencia biol´ogica de inter´es, construir el modelo de Markov oculto es relativamente sencillo y puede hacerse de una manera intuitiva.

Figura 5.3: Un modelo de Markov oculto sencillo para modelizar genes eucariotas. El HMM construido se utiliza ahora para analizar nuevas secuencias observadas. Por ejemplo, supongamos que tenemos una nueva secuencia de ADN X = x1 , . . . , x19 = ATGCGACTGCATAGCACTT. ¿C´omo podemos saber si esta secuencia de ADN es la regi´on codificada de un gen o no? O, si asumimos que X es un gen que codifica en prote´ına, ¿c´omo podemos predecir las posiciones de los exones y los intrones en la secuencia dada? Podemos responder a la primera pregunta calculando la probabilidad observada de X basada en el HMM dado, que modeliza genes codificantes. Si esta probabilidad es alta, implica que esta secuencia de ADN es probablemente la de la regi´on codificada de un gen. Por otra parte, podr´ıamos

56

CAP´ITULO 5. LOS MODELOS DE MARKOV OCULTOS

concluir que X no es probable que sea una regi´on codificante ya que no contiene las propiedades estad´ısticas que se observan t´ıpicamente en los genes codificantes en prote´ınas. La segunda cuesti´on trata de la predicci´on de la estructura interna de la secuencia, que no puede ser directamente observada. Para responder a esta pregunta, debemos primero predecir la secuencia de estados Π en el modelo de Markov oculto que mejor describa X. Una vez que tenemos el mejor Π, es sencillo predecir las localizaciones de los exones e intrones. Por ejemplo, supongamos que la secuencia de estados ´optima es Π, como se muestra en la figura (donde se denomina y). Esto implica que las nueve primeras bases x1 , . . . , x9 pertenecen al primer ex´on, las siguientes cuatro bases x10 , . . . , x13 pertenecen a un intr´on, y las u ´ltimas seis bases x14 , . . . , x19 pertenecen al otro ex´on.

5.3.1.

GENSCAN

El sofware para la b´ usqueda de genes GENSCAN (http://genes.mit.edu/GENSCAN.html) ha sido desarrollado por Chris Burge y Samuel Karlin en la Universidad de Stanford, y est´a actualmente alojado en el departamento de biolog´ıa del MIT (Instituto Tecnol´ogico de Massachusetts). Este programa utiliza un modelo probabil´ıstico complejo sobre la estructura del gen, que se basa en la actual informaci´on biol´ogica sobre las propiedades de las se˜ nales de transcripci´on, traducci´on y de empalme. Adem´as, utiliza varias propiedades estad´ısticas de regiones codificantes y no codificantes. Para tener en cuenta la heterogeneidad del genoma humano que afecta a la estructura y densidad de los genes, GENSCAN deriva diferentes conjuntos de modelos de genes de regiones del genoma con diferente contenido en CG. Es muy r´apido y su precisi´on hace que el GENSCAN sea el m´etodo elegido para el an´alisis inicial de grandes (en el rango de megabases) conjuntos de ADN en genomas eucariotas. GENSCAN ha sido utilizado como la herramienta principal para la predicci´on de genes en el Proyecto Internacional del Genoma Humano. A continuaci´on, exponemos la idea para poder responder a la segunda cuesti´on del ejemplo, y por tanto poder predecir estructuras de genes. Como hemos dicho, el siguiente algoritmo est´a basado en el esp´ıritu del GENSCAN. Modelos de Markov semi-ocultos Un modelo de Markov semi-oculto (semiHMM) o generalizado es un modelo de Markov oculto en el que cada estado puede emitir una secuencia de observaciones, y no solamente una observaci´on o s´ımbolo como en los modelos de Markov ocultos estudiados. Veamos c´omo podemos modelizar esta idea. En un HMM, supongamos que p es la probabilidad de transici´on de cualquier estado a ´el mismo. La probabilidad de que el proceso se quede en un estado para n pasos es p(n−1) (1−p) y, por lo tanto, la cantidad de tiempo que el proceso est´a en este estado sigue una distribuci´on geom´etrica de par´ametro (1−p). Para el modelo de ge-

´ 5.3. APLICACIONES: BUSQUEDA DE GENES

57

Figura 5.4: Diagrama de estados del HMM en el que se basa GENSCAN.

nes que construyamos, es necesario permitir otras distribuciones para esta cantidad. En un HMM semi-oculto (semiHMM) todas las probabilidades de transici´on de un estado a ´el mismo son cero; y cuando el proceso visita un estado, ´este produce no un s´olo s´ımbolo del alfabeto sino una secuencia entera. La longitud de la secuencia puede seguir cualquier distribuci´on, y el modelo que genera la secuencia de esa longitud puede seguir tambi´en cualquier distribuci´on. Las posiciones en la secuencia emitida de un estado no necesitan ser iid. El modelo se formula m´as precisamente como sigue. Cada estado π tiene asociado a ´el mismo una variable aleatoria Lπ (L de “longitud”) cuyo rango es un subconjunto de 0, 1, 2, . . .; y para cada valor observado l de Lπ existe una variable aleatoria Yπ,l cuyo rango consiste de todas las secuencias de longitud l. Cuando se visita el estado π, se determina aleatoriamente una longitud l de la distribuci´on de Lπ . Luego, la distribuci´on de Yπ,l se utiliza para determinar una secuencia de longitud l. Despu´es, se toma una transici´on a un nuevo estado y el proceso se repite, generando otra secuencia. Estas secuencias se concatenan para crear la secuencia final del modelo de Markov semi-oculto. Los algoritmos involucrados en este modelo tienen un orden de magnitud m´as complejo que un HMM regular, ya que dada una secuencia observada no s´olo no conocemos el camino de estados que produce, sino tampoco conocemos la divisi´on de puntos en la secuencia indicando d´onde se hace una transici´on a un nuevo estado.

58

CAP´ITULO 5. LOS MODELOS DE MARKOV OCULTOS

La aplicaci´on en la b´ usqueda de genes requiere una generalizaci´on del algoritmo de Viterbi. Es una generalizaci´on natural, sin embargo, ya que se trabaja generalmente con secuencias muy largas, la generalizaci´on natural no funciona en tiempo razonable. En la pr´actica, deben hacerse m´as suposiciones. Burge ([9]) observ´o que si las longitudes de las largas regiones interg´enicas se toman siguiendo distribuciones geom´etricas, y si estas longitudes generan secuencias en una manera relativamente iid, entonces el algoritmo puede ajustarse de modo que puedan obtenerse tiempos de funcionamiento pr´acticos. Estas suposiciones no son irrazonables en nuestro caso y, por tanto, no deber´ıan afectar en gran medida a la precisi´on de las predicciones. Omitiremos los detalles t´ecnicos que subyacen en este tema. Nuestro objetivo es mostrar la idea principal de c´omo trabaja un modelo de Markov oculto buscador de genes. De esta forma, un par φ es una secuencia de estados π1 , π2 , ..., πn y una secuencia de longitudes l1 , l2 , ..., ln . Dada una secuencia observada X de un semiHMM, el algoritmo de Viterbi encuentra un par o´ptimo φopt tal que la P (φopt |X) ≥ P (φ|X) para todos los pares φ. En otras palabras, φopt es el par m´as probable para dar lugar a la secuencia X. El par o´ptimo da las predicciones de genes, as´ı como la predicciones de la estructura completa del gen. En la pr´actica, hay muchas m´as consideraciones para optimizar el m´etodo. Quiz´as una de las m´as importantes es que muchas de las probabilidades estimadas dependen del contenido CG de la regi´on de ADN buscada. Por ejemplo, regiones con mayor contenido de CG (ilas CpG) tienden a contener una densidad de genes significativamente mayor. El modelo de Burge toma esto y otros factores en cuidadosa consideraci´on, y contin´ ua refin´andose tanto como diferentes tipos de datos haya disponibles. Veamos de forma breve la idea de modelizaci´on de los modelos de Markov semiocultos. El modelo. Cada estado del modelo que construyamos es un modelo es su propio sentido. Es necesario “entrenar” a cada estado para producir secuencias que modelen las partes correspondientes del gen en cuesti´on. Para ello, utilizamos un gran conjunto de datos de entrenamiento que consiste en largas cadenas de ADN, donde las estructuras de los genes han sido completamente caracterizadas. Burge et al. ([9]) recopilaron 2.5 millones de pares de bases (Mb) de ADN humano con 380 genes, 142 genes ex´on individuales, y un total de 1492 exones y 1254 intrones. Adem´as, se incluyeron regiones s´olo codificantes (sin intrones) de 1619 genes humanos. Modelamos, por ejemplo, un gen orientado de 5’ a 3’, con 13 estados de un modelo de Markov semi-oculto, como muestra la figura inferior. La primera fila representa la regi´on interg´enica o no funcional. La segunda fila representa la regi´on promotora. La tercera fila es la 5’UTR. La cuarta fila de los cinco estados representa los intrones y exones. La fila final es la 3’UTR y la se˜ nal Poli(A). La regi´on no funcional entre

´ 5.4. CASO PRACTICO: LAS ISLAS CPG

59

Figura 5.5: Diagrama de estados. genes est´a etiquetada por N en la figura. Primero se describe c´omo se forma cada estado dando las distribuciones de Lπ y de Yπ,l para cada estado π. Despu´es, se establecen las probabilidades de transici´on entre estados basadas de los modelos elegidos para cada regi´on del gen. Para cada uno de los estados tendremos un modelo distinto, dependiendo de sus caracter´ısticas gen´eticas. La figura muestra la complejidad de esta modelizaci´on, y las consecuentes estimaciones de par´ametros de forma o´ptima. No profundizaremos m´as en el estudio de los modelos de Markov semi-ocultos en b´ usqueda de genes por no ser objeto de estudio del presente trabajo ([40]).

5.4.

Caso pr´ actico: las islas CpG

Como aplicaci´on pr´actica a toda la teor´ıa estudiada anteriormente, y siguiendo con los ejemplos ya propuestos sobre este tema, nos planteamos la b´ usqueda de islas CpG en un organismo en particular, el “Danio Rerio” o “pez cebra”. Tomaremos la secuencia de datos del cromosoma 10 de tal organismo, para analizarla y encontrar en ella las posibles islas CpG, utilizando los modelos de Markov ocultos. Realizaremos esta b´ usqueda con R y el programa Bioconductor, de donde podemos descargarnos el genoma entero del pez cebra. Finalmente, una vez encontradas las posibles islas, comprobaremos si ´estas se encuentran en las mismas posiciones, aproximadamente, que las estudiadas en el art´ıculo “Redefining CpG islands using hidden Markov models”, por Rafael A. Irizarry et al. ([24]); donde se utiliza una variante “m´as fina”

60

CAP´ITULO 5. LOS MODELOS DE MARKOV OCULTOS

de los modelos de Markov ocultos, llamados modelos de Markov ocultos jer´arquicos. Estas posiciones ya estudiadas ser´an nuestros datos de entrenamiento. Podemos descarg´arnoslos de la siguiente direcci´on web: http://rafalab.jhsph.edu/CGI/. El paquete de R “makeCGI” realiza la b´ usqueda de regiones CpG autom´aticamente, pero no lo utilizaremos para el presente trabajo pues nos interesa programar los pasos y el algoritmo utilizado. Comenzamos por descargarnos el programa Bioconductor y el genoma del pez cebra (versi´on 6 correspondiente a la secuenciaci´on de su genoma en el a˜ no 2008), incluyendo en R las sentencias > source("http://bioconductor.org/biocLite.R") > available.genomes() # genomas disponibles de diferentes organismos [16] "BSgenome.Drerio.UCSC.danRer5" [17] "BSgenome.Drerio.UCSC.danRer6" [18] "BSgenome.Drerio.UCSC.danRer7" > biocLite("BSgenome.Drerio.UCSC.danRer6") # versi´ on 6, como el art´ ıculo > library(BSgenome.Drerio.UCSC.danRer6) > Drerio Zebrafish genome | | organism: Danio rerio (Zebrafish) | provider: UCSC | provider version: danRer6 | release date: Dec. 2008 | release name: Sanger Institute Zv8 | | single sequences (see ’?seqnames’): | chr1 chr2 chr3 chr4 chr5 chr6 chr7 chr8 | chr9 chr10 chr11 chr12 chr13 chr14 chr15 chr16 | chr17 chr18 chr19 chr20 chr21 chr22 chr23 chr24 | chr25 chrM | | multiple sequences (see ’?mseqnames’): | Zv8_NA Zv8_scaffold upstream1000 upstream2000 | upstream5000 | | (use the ’$’ or ’[[’ operator to access a given sequence) > class(Drerio) # objeto clase BSgenome [1] "BSgenome" attr(,"package") [1] "BSgenome"

´ 5.4. CASO PRACTICO: LAS ISLAS CPG > seqnames(Drerio) # [1] "chr1" "chr2" [7] "chr7" "chr8" [13] "chr13" "chr14" [19] "chr19" "chr20" [25] "chr25" "chrM"

elegiremos el cromosoma "chr3" "chr4" "chr5" "chr9" "chr10" "chr11" "chr15" "chr16" "chr17" "chr21" "chr22" "chr23"

61 10 etiquetado por ‘‘chr10’’ "chr6" "chr12" "chr18" "chr24"

> crom10 crom10 43467561-letter "MaskedDNAString" instance (# for masking) seq: CACACACACACACACACACACACACA...TATTCGGAACATTTTAATGTCAGAT masks: maskedwidth maskedratio active names 1 110500 2.542126e-03 TRUE AGAPS 2 413 9.501338e-06 TRUE AMB 3 22431376 5.160486e-01 FALSE RM 4 1151720 2.649608e-02 FALSE TRF desc 1 assembly gaps 2 intra-contig ambiguities 3 RepeatMasker 4 Tandem Repeats Finder [period length(crom10) # 43.467.561 de nucle´ otidos de longitud [1] 43467561 Una vez tenemos el cromosoma 10 disponible, observamos lo que ser´an nuestros datos de entrenamiento, es decir, las posiciones ya estudiadas de las islas que nos hemos descargado. > posic head(posic)

1 2 3 4 5

ı ¨..chr chr10 chr10 chr10 chr10 chr10

start 10146 11319 25650 34657 45703

end length CpGcount GCcontent pctGC obsExp 10266 121 11 72 0.595 1.034 12498 1180 148 731 0.619 1.312 25762 113 10 56 0.496 1.471 34818 162 8 84 0.519 0.756 45986 284 28 157 0.553 1.290

CAP´ITULO 5. LOS MODELOS DE MARKOV OCULTOS

62 6

chr10 46944 47514

571

38

301 0.527

0.967

> posic10 dim(posic10) # es un data.frame de 3889 filas 8 columnas [1] 3889 8 > start end class(crom10) # objeto de la clase MaskedDNAString [1] "MaskedDNAString" attr(,"package") [1] "Biostrings" lo que significa que existen componentes o posiciones en la secuencia del genoma dentro cromosoma 10 que est´an enmascaradas, es decir, que no se muestran. Si nos fijamos en la frecuencia de las letras de nuestro alfabeto, observamos cu´ales son las letras enmascaradas y en qu´e cantidad se encuentran en el genoma > af10_1 af10_1 A 13792137 W 0 D 0

C 7928689 S 0 B 0

G T 7888654 13747168 Y K 0 0 N 0 0

M 0 V 0 + 0

R 0 H 0

> crom10 class(crom10) # objeto de la clase ‘‘DNAString’’ [1] "DNAString" attr(,"package") [1] "Biostrings" > af10_2 af10_2 # N son 110913

´ 5.4. CASO PRACTICO: LAS ISLAS CPG A 13792137 W 0 D 0

C 7928689 S 0 B 0

G T 7888654 13747168 Y K 0 0 N 110913 0

63 M 0 V 0 + 0

R 0 H 0

Como ya indicamos anteriormente, las posiciones etiquetadas con una “N ” son nucle´otidos del genoma todav´ıa desconocidos o sin determinar. Para realizar el an´alisis, y poder calcular las matrices de transici´on y de emisi´on necesarias que son nuestros par´ametros, tenemos que trabajar con la secuencia desenmascarada. Decido eliminar los “gaps” etiquetados con “N” para calcular las frecuencias del resto de nucle´otidos y dinucle´otidos (secuencias de dos nucle´otidos de longitud), y poder calcular ambas matrices. Con ellas, aplicaremos el algoritmo de Viterbi. Calculo la matriz de transici´on a la que llamar´e “trans”, de dimensi´on 8 × 8 de acuerdo a los 8 estados (A+ , C+ , G+ , T+ , A− , C− , G− , T− ; en este orden), de la siguiente forma: calculo las tablas de contingencia de todos los posibles dinucle´otidos en las regiones de islas y de oc´eano, por lo que obtengo dos matrices 4 × 4. En efecto, utilizando las frecuencias de los dinucle´otidos, puedo conocer las transiciones de un estado o nucle´otido a otro. Dimensiono estas matrices en los bloques diagonales de mi matriz “trans”. Los restantes espacios de la matriz de transici´on miden las probabilidades de transici´on de una isla a un oc´eano (bloque superior) y al rev´es (bloque inferior). Para calcularlas, tomo la primera letra antes de una isla y la primera que est´a despu´es de una isla, para todas las islas de “posic10”; y calculo las frecuencias (o transiciones) con la primera de la isla y la u ´ltima de la isla, respectivamente. Dimensiono el resultado en dos matrices 4 × 4 que a˜ nado de nuevo a mi matriz “trans” en el orden correspondiente. > # Modelo de islas > cadenamas for (i in 1:3889) + cadenamas[i,] # juntamos filas de ‘‘cadenamas’’ > juntomas juntomas > > >

cantidad_n + > >

vmas seq > > > > > >

x