Repositorio Institucional de Documentos

Centro Polit´ ecnico Superior Universidad de Zaragoza Trabajo Fin de M´ aster Curso 2010-2011 M´ aster en Ingenier´ ıa de Sistemas e Inform´ atica An´alisis filogen´etico mediante l´ogica temporal y model checking Jos´e Ignacio Requeno Director: Jos´e Manuel Colom Departamento de Inform´atica de Ingenier´ıa de Sistemas, CPS 20 de Noviembre del 2010 An´alisis filogen´etico mediante l´ogica temporal y model checking Resumen El an´alisis filogen´etico es la rama de la bioinform´atica que se encarga de estudiar la clasificaci´on de un conjunto de especies distintas en funci´on de su relaci´on de proximidad evolutiva. La metodolog´ıa de trabajo que sigue actualmente es muy “ad-hoc”; es decir, cuando especificamos una nueva propiedad biol´ogica o a˜nadimos otra especie al modelo necesitamos reimplementar el algoritmo de verificaci´on para que considere los nuevos datos. Para solucionarlo, nuestra aproximaci´on trata de sistematizar el m´etodo de trabajo en el an´alisis filogen´etico, adem´as de reaprovechar conceptos y m´etodos consolidados en otras ´areas de an´alisis de sistemas. Para ello proponemos el uso de “model checking”, una t´ecnica de verificaci´on autom´atica madura en el mundo de la verificaci´on. Nuestra aproximaci´on se sustenta sobre tres pilares principales. Primero, el modelado de la din´amica evolutiva es an´alogamente equivalente al modelado de un sistema de transiciones, donde las caracter´ısticas biol´ogicas de las especies (ADN, fenotipo, etc) definen las variables de los estados; y la mutaci´on o salto entre especies se corresponde con la ejecuci´on de una transici´on. Segundo, el uso de l´ogicas temporales sistematiza el proceso de formulaci´on y aporta formalismo a la especificaci´on de propiedades filogen´eticas. Por ´ultimo, la verificaci´on de estas f´ormulas se realiza de manera autom´atica por cualquiera de las herramientas inform´aticas existentes, devolviendo un contraejemplo en el caso de resultar falsas. Adem´as de la aportaci´on te´orica, hemos evaluado emp´ıricamente las prestaciones del sistema sobre una herramienta software de model checking real: SMV. Esta tecnolog´ıa ha demostrado ser viable cuando trabaja aisladamente con genes o fragmentos de ADNmt, ya que su tama˜no no es excesivamente grande. Sin embargo, su coste temporal escala cuadr´aticamente con la longitud de las secuencias, mientras que el consumo de memoria crece linealmente hasta alcanzar m´as de 2,5 GBytes en el caso peor. Por estos motivos, en estudios filogen´eticos futuros sobre el ADNmt completo (o el ADN nucle- ar, que es m´as grande) ser´a necesario paralelizar y componer la soluci´on a partir de los resultados parciales. Tambi´en se ha propuesto una estructura de datos basada en los diagramas de decisi´on para manejar cadenas de ADNmt de forma m´as eficiente en memoria a c´omo se realiza en la actualidad. ´ Indice general 1. Introducci´on 7 1.1. Contexto.............................. 7 1.2. Motivaci´on............................. 8 1.3. Objetivos ............................. 8 1.4. Organizaci´on del documento . . . . . . . . . . . . . . . . . . . 10 2. Uniendo mundos: Filogen´etica y Model Checking 12 2.1. Contexto.............................. 12 2.2. Evoluci´on como un Sistema de Transiciones . . . . . . . . . . 12 2.3. L´ogica Temporal como un Lenguaje de Especificaci´on . . . . . 16 2.4. Model Checking como un Entorno de Inferencia . . . . . . . . 17 3. Modelado de propiedades estructurales 19 3.1. Contexto.............................. 19 3.2. Propiedades del ´ Arbol ...................... 19 3.3. Propiedades de Secuencia . . . . . . . . . . . . . . . . . . . . 23 4. Implementaci´on sobre un paquete real 25 4.1. Traducci´onaSMV ........................ 25 4.2. Rendimiento............................ 27 5. Representaci´on compacta de secuencias de ADN mediante diagramas de decisi´on 34 5.1. Contexto.............................. 34 5.2. Estadodelarte .......................... 35 5.3. Herramientas y m´etodos . . . . . . . . . . . . . . . . . . . . . 36 5.4. Heur´ısticas............................. 37 5.4.1. Reordenaci´on inicial . . . . . . . . . . . . . . . . . . . 37 5.4.2. Reglas de reducci´on . . . . . . . . . . . . . . . . . . . . 40 5.4.3. Alternativas........................ 42 2 ´ INDICE GENERAL 3 5.5. Complejidad m´ınima . . . . . . . . . . . . . . . . . . . . . . . 44 6. Conclusiones 49 6.1. Conclusiones............................ 49 6.2. Trabajofuturo .......................... 50 Bibliograf´ıa 53 ´ Indice de figuras 2.1. Traducci´on de un ´arbol a una estructura de Kripke . . . . . . 15 2.2. Evaluaci´on de operadores l´ogicos temporales . . . . . . . . . . 17 3.1. Distintos tipos de grupos filogen´eticos . . . . . . . . . . . . . . 22 4.1. ADN mitocondrial . . . . . . . . . . . . . . . . . . . . . . . . 28 4.2. Gr´afica de crecimiento temporal lineal respecto del n´umero de muestras.............................. 31 4.3. Gr´afica de crecimiento temporal cuadr´atico respecto de la longitud de las secuencias . . . . . . . . . . . . . . . . . . . . . . 31 5.1. Representaci´on en ROBDD de una secuencia (SeqBDD) . . . . 36 5.2. Set Decision Diagram . . . . . . . . . . . . . . . . . . . . . . . 40 5.3. Fusi´ondenodos.......................... 41 5.4. Uni´ondecaminos......................... 41 5.5. Permutaci´on local . . . . . . . . . . . . . . . . . . . . . . . . . 42 5.6. ncadenas completamente distintas (|y1|=0).......... 45 5.7. ncadenas completamente iguales (|y1|=k)........... 45 5.8. Esquema de la estructura final (|ym|=0)............ 46 5.9. Esquema de la estructura final (|yn|>0)............ 48 4 ´ Indice de tablas 4.1. Coste temporal de creaci´on de la estructura de Kripke y almacenamiento de las secuencias . . . . . . . . . . . . . . . . . . . 29 4.2. Coste en memoria para la creaci´on de la estructura de Kripke y almacenamiento de las secuencias . . . . . . . . . . . . . . . 30 4.3. Coste temporal de verificaci´on de propiedades . . . . . . . . . 32 5.1. Secuencias de ADN . . . . . . . . . . . . . . . . . . . . . . . . 38 5.2. Reordenaci´on inicial . . . . . . . . . . . . . . . . . . . . . . . 39 5.3. Estructuraideal.......................... 42 5.4. Reordenaci´on alternativa . . . . . . . . . . . . . . . . . . . . . 43 5.5. Primer paso de la inducci´on . . . . . . . . . . . . . . . . . . . 44 5.6. Paso de inducci´on con (m≤n) ................. 47 5.7. Paso de inducci´on con (m > n) ................. 47 5 Agradecimientos Es de bien nacido ser agradecido. En primer lugar me gustar´ıa agradecer a Elvira y Jos´e Manuel por su sabidur´ıa y consejos, sin los cuales este trabajo nunca hubiera llegado a buen puerto. Tambi´en, c´omo no, a mi fiel compa˜nero de armas Roberto, que me apadrin´o desde el primer momento en mi retorno a la universidad por el sendero de la bioinform´atica. Recordar´e adem´as siempre en esta etapa a Goyo, por su visi´on num´erica de la vida; y a Jorge, que con sus caf´es, erizos y conversaciones distendidas consigui´o amenizar mis momentos de descanso. Querr´ıa hacer una menci´on especial a mis padres, que me apoyaron en todo momento y sin los cuales no estar´ıa aqu´ı ahora. Por ´ultimo, mencionar expl´ıcitamente al I3A por su beca de iniciaci´on a la investigaci´on y a la DGA por su beca doctoral [17030 G/5423/480072/91002], que con su financiaci´on me permitieron continuar con mi la labor investigadora. 6 Cap´ıtulo 1 Introducci´on 1.1. Contexto La patolog´ıa mitocondrial, provocada por la presencia de mutaciones en el ADN mitocondrial (en adelante, ADNmt), es una de las causas de enfermedades gen´eticas raras m´as prevalente. El ADNmt tiene importantes caracter´ısticas diferenciales respecto al ADN nuclear. En primer lugar, su dimensi´on m´as reducida y una mayor abundancia de copias por c´elula lo hacen m´as manejable. Adem´as, su herencia es estrictamente matrilineal. Esto implica que no hay recombinaci´on gen´etica, simplificando considerablemente la trazabilidad de las mutaciones y los problemas derivados de la reordenaci´on de genes por la interacci´on cromos´omica. Por ´ultimo, su alta tasa de mutaci´on, mucho m´as elevada que la del ADN nuclear, est´a debida principalmente a su localizaci´on en un entorno reactivo como es el de la mitocondria, org´anulo encargado de la respiraci´on aer´obica celular. Gracias a estas tres peculiaridades, hacen al ADNmt objeto habitual de estudio en el ´ambito bioqu´ımico. En especial, este trabajo ha tomado el ADNmt como punto de partida dado el inter´es que suscita, por todos los comentarios anteriores, entre miembros de la Facultad de Veterinaria (con quienes se ha colaborado en el desarrollo de este trabajo). Una de las herramientas con las que cuentan los cient´ıficos para estudiar el ADNmt es el an´alisis filogen´etico. El an´alisis filogen´etico es una rama de la bioinform´atica que se ocupa de la aplicaci´on de modelos y m´etodos computacionales para la construcci´on de un ´arbol (o grafo) que plasme las relaciones evolutivas de un conjunto de genes, individuos o especies (tambi´en llamados taxones). La sucesiva incorporaci´on de m´etodos formales para caracterizar conceptos claves, como la distancia evolutiva, ha proporcionado resultados valiosos. 7 CAP´ ITULO 1. INTRODUCCI ´ ON 8 La representaci´on num´erica de relaciones entre conjuntos de poblaciones o individuos ha sido explorada por m´etodos arborescentes basados en distancias (neigbour joining), y por m´etodos basados en minimizar el n´umero de cambios de caracter´ısticas biol´ogicas (parsimonia, verosimilitud, inferencia Bayesiana) [15]. Desde un punto de vista biol´ogico, el ´arbol (o potenciales ´arboles) evolutivo que obtenemos como resultado tiene una aplicaci´on directa para inferir propiedades biol´ogicas interesantes: identificar mutaciones pat´ogenas causantes de graves enfermedades (pues estar´an en hojas terminales que no se han propagado hasta el presente), deducir el hist´orico de las migraciones poblacionales, estimar el momento de divergencia entre dos especies o incluso encontrar regiones codificantes de gran importancia. 1.2. Motivaci´on Para averiguar propiedades biol´ogicas como las comentadas anteriormente es importante detectar patrones caracter´ısticos en las secuencias de ADN (propiedades como la conservaci´on o covariaci´on de nucle´otidos y amino´acidos [27]), o analizar la organizaci´on en grupos clad´ısticos (familias de especies). Sin embargo, actualmente el an´alisis filogen´etico carece de cierta flexibilidad en el sentido de que el modelo evolutivo, la especificaci´on de propiedades y su proceso de verificaci´on no est´an completamente desacoplados. Esto conduce a una forma r´ıgida de trabajo, donde la comprobaci´on de una nueva propiedad, la modificaci´on del modelo o la optimizaci´on del proceso de verificaci´on deriva en el redise˜no desde cero de todo el sistema software. Adem´as, la inherente naturaleza temporal de los datos filogen´eticos sugiere la posibilidad de introducir nuevos m´etodos formales novedosos, capaces de mejorar esta falta de flexibilidad en los modelos convencionales, para incorporar caracter´ısticas como: predicci´on futura, embeber informaci´on del pasado y combinaci´on de reglas evolutivas a varios niveles de abstracci´on. 1.3. Objetivos La b´usqueda de un entorno de trabajo flexible, que desacople a´un mejor el modelado y especificaci´on de la tarea de verificaci´on, ha sido nuestro principal objetivo y el eje central del trabajo realizado. De esta forma, los investigadores podr´an enfocar sus esfuerzos exclusivamente al modelado de la evoluci´on y la especificaci´on de propiedades, mientras que los algoritmos gen´ericos de verificaci´on se encargar´an del resto. CAP´ ITULO 2. UNIENDO MUNDOS:FILOGEN´ ETICA Y MODEL CHECKING15 En este punto, podemos definir una estructura de branching-time para los ´arboles filogen´eticos, y que servir´a de base para interpretar las f´ormulas de l´ogica temporal que emplearemos para expresar propiedades del ´arbol. La f´ormula de estado m´as sencilla determina si el estado actual est´a asociado con la secuencia s=σ1σ2. . . σl∈Σl(o con una subcadena de la secuencia, o con un conjunto de secuencias). Las secuencias pueden manipularse simb´olicamente (como si fueran conjuntos de datos) mediante una agregaci´on de sus partes. En t´erminos l´ogicos, quedar´ıa representado como: s≡ l ^ i=1 s[i] = σi(2.1) Definici´on 3 (Filogenia Branching-time) Un ´arbol (por Def. 1) P= (T, r, D) est´a inequ´ıvocamente definido por una estructura de Kripke an´aloga MP= (V, {r}, RP, LP), donde: RPes la relaci´on de transici´on compuesta por el conjunto de arcos del ´arbol (dirigidos desde la ra´ız r) m´as los autobucles en las hojas: RP=E∪ {(v, v) : @(v, w)∈E∧v, w ∈V}, y LPes una funci´on est´andar de etiquetado definida por APP, donde un estado vmapeado como D(v) = σ1σ2. . . σlsatisface la familia de propiedades s[i] = σi,1≤i≤l, adem´as de otras necesarias para mantener la unicidad l´ogica de cada estado (por ejemplo, un identificador num´erico propio). Figura 2.1: Traducci´on de un ´arbol a una estructura de Kripke La figura 2.1 ilustra este proceso de traducci´on desde un ´arbol (filogen´etico) etiquetado (Def. 1) a una estructura de Kripke (Def. 2). Recalcamos que pueden considerarse f´acilmente otros modelos de evoluci´on (redes evolutivas, etc) ajustando las definiciones Defs. 1–3 seg´un se requiera. CAP´ ITULO 2. UNIENDO MUNDOS:FILOGEN´ ETICA Y MODEL CHECKING16 2.3. L´ogica Temporal como un Lenguaje de Especificaci´on La l´ogica temporal es un sistema formal que permite la representaci´on y manipulaci´on de proposiciones l´ogicas cualificadas en t´erminos de tiempo [21]. En el contexto de los sistemas de transici´on, se usan para definir propiedades en secuencias de transiciones entre estados de un sistema previamente definido seg´un la abstracci´on correspondiente (en nuestro caso, un tipo espec´ıfico de estructura de Kripke). Por ejemplo, las propiedades pueden indicar si desde un determinado punto es posible alcanzar un estado en particular, o si una propiedad siempre se cumple. Las l´ogicas temporales se clasifican en dos grandes familias dependiendo de c´omo tratan las secuencias de eventos: mientras que las l´ogicas linear-time engloban las l´ogicas que trabajan ´unicamente con caminos individuales π= s0s1s2. . ., las l´ogicas branching-time tienen en cuenta el conjunto de posibles opciones desde un estado (y por tanto, razonando sobre el total del ´arbol de computaci´on). Computational Tree Logic (CTL) es el m´aximo exponente de este ´ultimo tipo de familias y el que ha sido ampliamente adoptado por la comunidad de model checking [13]. Dado que las filogenias representan procesos evolutivos que por naturaleza son esencialmente bifurcaciones, las l´ogicas de tipo branching-time en general, y CTL en particular, son adecuadas para su correcta descripci´on. CTL reinterpreta los cuantificadores de la l´ogica de primer orden como cuantificadores de camino. Esos cuantificadores expresan el cumplimiento de una propiedad a lo largo de todos los caminos de computaci´on (A), o al menos a lo largo de uno de sus caminos (E). Estos dos cuantificadores deben estar inmediatamente cualificados por alguno de los cinco operadores temporales, de los cuales tres expresan la satisfacci´on de una propiedad eventualmente en cualquier momento del futuro (F), para todos los instantes (G), o en el estado siguiente (X). Los otros dos, son construcciones condicionales en donde la propiedad precedente se cumple hasta que el consecuente se haga verdad (U); o la consecuente se satisface hasta el momento en el que la precedente lo haga (si lo hace, y adem´as ambas se cumplen a la vez en el estado de transici´on) (R). Una gram´atica y sem´antica completa de las f´ormulas CTL puede definirse como un subconjunto representativo de operadores l´ogicos de la manera siguiente; (Fig. 2.2 ilustra la sem´antica de todas ellas como un conjunto.) Definici´on 4 (´ Arbol Filogen´etico L´ogico) Una f´ormula arbitraria de l´ogica temporal φse define a partir de la siguiente gram´atica m´ınima, donde CAP´ ITULO 2. UNIENDO MUNDOS:FILOGEN´ ETICA Y MODEL CHECKING17 p∈AP: φ::= p| ¬φ|φ∨φ|EX (φ)|EG (φ)|E[φUφ] (2.2) Las f´ormulas se chequean sobre una estructura Mconsiderando todos los caminos πque se originan desde un determinado estado inicial s0.M, s0φ significa que s0satisface la propiedad φ. La sem´antica para la verificaci´on de f´ormulas bien formadas es la siguiente (con π=s0s1s2. . .): M, s0p⇔p∈L(s0), M, s0¬φ⇔M, s02φ, M, s0φ∨ψ⇔M, s0φoM, s0ψ, M, s0EX (φ)⇔ ∃π:M, s1φ, M, s0EG (φ)⇔ ∃π:M, siφ, ∀i∈N, y M, s0E[φUψ]⇔ ∃π, i ∈N:M, siψy M, sjφ, 0≤j < i. Figura 2.2: Evaluaci´on de operadores l´ogicos temporales Una f´ormula CTL φrepresenta una propiedad que debe verificarse en determinados estados del ´arbol de computaci´on. En este contexto, un sistema Msatisface φsi y s´olo si todos y cada uno de los estados iniciales cumplen: Vs0∈S0M, s0φ. Como licencia particular y por temas de legibilidad, adaptaremos la notaci´on de los operadores condicionales para reescribirlos como E[φUψ]≡ {φ}EU {ψ}. Una l´ogica as´ı definida permite la expresi´on formal de propiedades gen´ericas relacionadas con la evoluci´on biol´ogica de secuencias, adem´as de su eventual verificaci´on autom´atica. 2.4. Model Checking como un Entorno de Inferencia El principio que opera bajo el model checking es simple: la ejecuci´on de software de verificaci´on en un ordenador para comprobar la correcci´on de CAP´ ITULO 2. UNIENDO MUNDOS:FILOGEN´ ETICA Y MODEL CHECKING18 un sistema. Para lograrlo, como entrada el software necesita un modelo del sistema junto con la especificaci´on de sus requerimientos, ambos provistos por el usuario. En el supuesto de que alguna de las especificaciones propuestas no se cumpla, el software devuelve mediante contraejemplos los escenarios en los que se infringe la propiedad. Sin embargo, los algoritmos b´asicos no pueden gestionar por s´ı solos la explosi´on del espacio de estados, por lo que se necesita una manipulaci´on simb´olica de los sistemas y las f´ormulas l´ogicas [12]. Y lo que es a´un m´as importante, el uso de las t´ecnicas de model checking es completamente transparente respecto del sistema que estamos verificando, dado que son independientes del dominio de aplicaci´on y de su interpretaci´on. Esto significa que los filogenetistas pueden redirigir los esfuerzos que antes dedicaban a temas de implementaci´on para enfocarlos al modelado (antes de ejecutar el model checking, estableciendo las propiedades deseables u obligatorias del sistema) y el an´alisis de los resultados (despu´es del model checking, observando el ´exito o fracaso del proceso y estudiando los par´ametros devueltos). Cap´ıtulo 3 Modelado de propiedades estructurales 3.1. Contexto Este cap´ıtulo est´a dedicado a ilustrar la metodolog´ıa a seguir para modelizar propiedades biol´ogicas filogen´eticas empleando el entorno l´ogico aqu´ı presentado. Generalmente, una propiedad no trivial puede descomponerse en otras m´as sencillas (pero igualmente expresivas) y despu´es sintetizar el resultado a partir de estas. Los beneficios que aporta semejante descomposici´on l´ogica se resumen en dos puntos. En primer lugar, se simplifica la formulaci´on y favorece la legibilidad. Y en segundo lugar, esas propiedades pueden reutilizarse modularmente para construir otras m´as complejas (bastan peque˜nos ajustes locales para variar la sem´antica de una f´ormula, o incluso producir otras propiedades completamente distintas). Existen dos clases generales de consultas: sobre propiedades globales o asociadas al ´arbol, en donde lo que se inspecciona es la propia estructura de la filogenia; o sobre propiedades locales oasociadas a la secuencia, donde la composici´on de caracter´ısticas de las secuencias es el eje central (y que probablemente estar´an complementadas por restricciones auxiliares). En resumen, en este cap´ıtulo nos dedicaremos a modelar formalmente ejemplos de ambos tipos de propiedades. 3.2. Propiedades del ´ Arbol Muchas de las propiedades t´ıpicas que se estudian sobre el ´arbol son de naturaleza clad´ıstica. Por ejemplo, una de las consultas m´as frecuentes pregunta si, dado un conjunto determinado de organismos S, constituyen un 19 CAP´ ITULO 3. MODELADO DE PROPIEDADES ESTRUCTURALES 20 grupo monofil´etico oclado bajo la filogenia particular de estudio. En otras palabras, ¿existe un sub´arbol que contiene exactamente esos organismos como sus hojas?. Formalmente, buscamos un nodo en alguna parte del ´arbol (EF) que funcione como ra´ız del grupo y, entonces: a) todo lo que debe estar dentro del sub´arbol, lo est´a; b) todo lo que debe estar fuera, tambi´en lo est´a fuera. Resumiendo, no hay ning´un intruso o espurio que contamine el sub´arbol, y s´olo los miembros de Sest´an en ´el. clado (S)≡EF (in (S)∧out (S)) (3.1) La regla de inclusi´on (in(S)) afirma que para cada secuencia del conjunto S(aplicando un and l´ogico global V) hay un camino que lo encuentra como hoja. Por su parte, la regla de exclusi´on (out(S)) demanda que todos los caminos finalicen en una hoja del mismo conjunto. En casos de computaciones infinitas como este, las hojas finales se localizan mediante los patrones F◦ AG (s) que definen los autobucles de los estados terminales: por construcci´on es el ´unico lugar a partir de donde la computaci´on de un sub´arbol permanece uniforme (es siempre igual). in (S)≡^ s∈S EF ◦AG (s) (3.2) out (S)≡AF ◦AG _ s∈S s!(3.3) Tal y como queda patente, las propiedades individuales por separado tienen la misma carga sem´antica. Aqu´ı, in (S) se satisface por todos los clados que contienen a s, y out (S) por todos los subclados estrictos. Es m´as, si la propiedad clad´ıstica se extiende para englobar tanto a las secuencias ancestrales como a las hojas, la estructura de la f´ormula principal permanece inalterada: s´olo las reglas de inclusi´on y exclusi´on necesitar´an un refinamien- to extra, que comentamos a continuaci´on, para que las secuencias objetivo puedan localizarse en cualquier lugar del sub´arbol y sin contaminaci´on de intrusos externos. in0(S)≡^ s∈S EF (s) (3.4) out0(S)≡AG _ s∈S s!(3.5) Una propiedad m´as excitante todav´ıa es averiguar, dado un ´arbol filogen´etico y una partici´on de sus hojas (que nace de una partici´on seg´un la CAP´ ITULO 3. MODELADO DE PROPIEDADES ESTRUCTURALES 21 clasificaci´on de las secuencias), si esa partici´on constituye una clasificaci´on en haplogrupos del ´arbol. Los haplogrupos son agregaciones de haplotipos relacionados que se identifican con polimorfismos caracter´ısticos comunes. Por tanto, definen poblaciones gen´eticas, que adem´as pueden marcarse geogr´aficamente tambi´en. Este estudio est´a enfocado a las regiones no recombinables del genoma, en especial del ADN mitocondrial (de donde surgi´o la notaci´on clad´ıstica original para haplogrupos [29]); y la mayor parte del cromosoma Y (para el que se utiliza esta notaci´on [32]). Esencialmente, un haplogrupo junto con un conjunto de poblaciones (haplogrupos hijo) que han nacido de ´el con el paso del tiempo, deben formar un clado. En otras palabras, un haplogrupo es un clado anidado: todos sus miembros ocupan las hojas del ´arbol, excepto posiblemente un n´umero de sub-sub´arboles que carecen completamente de miembros (y tienen a su vez una estructura anidada de clado). Una filogenia tiene una clasificaci´on v´alida si cada parte tiene una estructura de haplogrupo. clasificador (S1, S2, . . . , Sh)≡ h ^ i=1 haplogrupo (Si) (3.6) La comprobaci´on es trivial si se conocen las relaciones entre padres e hijos del haplogrupo, simb´olicamente definidas por una funci´on hijos (S) (que no una f´ormula). haplogrupo0(S)≡clado (S∪hijos (S)) (3.7) Sin embargo, normalmente nos interesa permitir flexibilidad en la ubicaci´on de los haplogrupos de estudio, el refinamiento de haplogrupos de grano grueso y la exploraci´on de hip´otesis alternativas. La ecuaci´on (3.7) puede extenderse para determinar si un conjunto de haplogrupos hijo existen, pero a efectos pr´acticos es suficiente con comprobar la calidad local de los haplogrupos para cada parte por separado, sin recurrir a ninguna informaci´on adicional (m´as all´a de la composici´on de las subf´ormulas). Formalmente, el haplogrupo es una relajaci´on de clado (S) a lo largo de su estructura local. Mientras que la regla de inclusi´on se preserva, al igual que en la b´usqueda de la ra´ız del haplogrupo, la regla de exclusi´on se reemplaza por una propiedad de clado anidada. Bajo esta situaci´on, todos los caminos eventualmente alcanzan un punto (AU) donde, o bien encuentran un miembro del subclado, o bien la ra´ız de un clado distinto, pero siempre siguiendo una traza a lo largo de nodos pertenecientes al haplogrupo en cuesti´on (para ello, se debe proporcionar de antemano una funci´on de pertenencia hipara cada Si). haplogrupo (S, h)≡EF (in (S)∧nested (S, h)) (3.8) CAP´ ITULO 3. MODELADO DE PROPIEDADES ESTRUCTURALES 22 nested (S, h)≡ {h}AU {out (S)∨nesting (S)}(3.9) nesting (S)≡AF ◦AG ¬_ s∈S s!(3.10) Fij´emonos bien que nesting (S) es lo opuesto de out (S). Obviamente, los haplogrupos terminales (esto es, clados o monofilias) son v´alidos para esta f´ormula. En t´erminos de clad´ıstica, la estructura de haplogrupos locales se corresponde con el concepto de grupo polifil´etico. Es precisamente por esto ´ultimo que los hace indistinguibles de los grupos parafil´eticos (basados solamente en el contenido de las hojas) por lo que la informaci´on de miembros ancestrales se torna necesaria (mirar el cap´ıtulo 3.3 e imagen 3.1). Figura 3.1: Distintos tipos de grupos filogen´eticos Por ´ultimo, la asunci´on de que miembros ancestrales del haplogrupo satisfacen las propiedades definidas a priori es razonable (es decir, no aparecen casos especiales). Al igual que antes, la incorporaci´on de datos ancestrales deriva en una familia relacionada de propiedades que permiten una evaluaci´on m´as comprensiva del proceso de evoluci´on. CAP´ ITULO 3. MODELADO DE PROPIEDADES ESTRUCTURALES 23 3.3. Propiedades de Secuencia En general, las propiedades asociadas a secuencias est´an basadas en f´ormulas de estado (es decir, aquellas que se eval´uan dentro de un nodo sin necesidad de recurrir a operadores temporales) compuestas junto con patrones temporales simples, principalmente para extender su estudio aplic´andolo a la filogenia completa o para explorar sus alrededores. A ese estilo de f´ormulas de estado las llamaremos patrones (p). Ofrecen un potente formalismo de descripci´on para formular restricciones generales sin las limitaciones de una aproximaci´on “ad-hoc”. A menudo, esas propiedades se usan no necesariamente para prohibir patrones de mutaci´on, sino como par´ametros de calidad (consulta y alerta de se˜nales inusuales posiblemente debidas a un comportamiento an´omalo) y marcarlas para analizar en detalle posteriormente. Un primer grupo de patrones lo constituyen las restricciones que estudian la correcci´on a nivel global, esto es, que supuestamente deben cumplirse a lo largo de toda la filogenia. Pueden clasificarse de la siguiente manera: La conservaci´on, que se modela como una restricci´on en los s´ımbolos que pueden aparecer en una posici´on dada de la secuencia. Com´unmente, el patr´on se codifica como un vector unidimensional de booleanos que clasifica cada s´ımbolo como permisible o no permisible. Sin embargo, entre otras extensiones, es posible definir familias de elementos compatibles (no necesariamente ligados a una determinada posici´on), y que puedan conmutarse entre s´ı para obtener secuencias “equivalentes”. La covariaci´on, que impone una relaci´on de dependencia entre dos (o m´as) posiciones en una secuencia. De forma general, se representa como una matriz bidimensional de booleanos en la que, para cada s´ımbolo de la primera columna, se marca el conjunto de s´ımbolos admisibles en la segunda columna. T´ıpicamente, para que la propiedad sea significativa, las asociaciones deben producirse entre s´ımbolos dispersos. Una combinaci´on de ambas. Un patr´on global as´ı definido es f´acilmente verificable con s´olo extenderlo a todo el ´arbol de computaci´on. global (p)≡AG (p) (3.11) Las excepciones que surjan a las propiedades mencionadas anteriormente pueden indicar que esas mutaciones son potencialmente peligrosas (o al menos CAP´ ITULO 3. MODELADO DE PROPIEDADES ESTRUCTURALES 24 sospechosas), lo que es de gran inter´es en los estudios aplicados a filogen´etica [27]. Adem´as, las (potenciales) mutaciones que se conocen a priori pueden modelarse expl´ıcitamente como patrones y evaluarse en las posiciones adecuadas de la filogenia. Especialmente si esas mutaciones afectan a funciones metab´olicas importantes, es previsible que sean pat´ogenas y limiten o impidan la reproducci´on del organismo, confin´andolo a las hojas terminales (o cerca de ellas). Aunque para algunas mutaciones est´a sistem´aticamente prohibido que aparezcan como patrones globales, para otras mutaciones pat´ogenas se les est´a permitido aparecer bajo a ciertas condiciones. Especialmente, es necesario que si aparece un patr´on de mutaci´on peligroso, ´este no tenga sucesores (es decir, sea una hoja en la filogenia); o para proveerlo de cierta flexibilidad, desaparezca al cabo de un m´aximo de kpasos o sucesores (AXk). terminal (p)≡AG (p→hoja) (3.12) terminal (p, k)≡AG p→AXk(hoja)(3.13) En este caso, las hojas (autobucles en la estructura de Kripke) deben detectarse sin ninguna referencia a una secuencia particular. Para conseguirlo f´acilmente es suficiente con comparar los vectores de estado (los valores de AP) propios del estado objetivo con los de todos sus sucesores. hoja ≡^ p∈AP p↔AX (p) (3.14) Este ´ultimo ejemplo representa las propiedades que desarrollan exploraciones condicionales de la filogenia. La verificaci´on de una l´ınea evolutiva espec´ıfica ser´ıa el paso siguiente, donde los patrones servir´ıan para definir conjuntos de estados relevantes para una f´ormula, aunque su estudio est´a fuera del ´ambito de alcance de este trabajo. Merece la pena rese˜nar que el chequeo de patrones en la clasificaci´on de haplogrupos cae dentro de esta categor´ıa. CAP´ ITULO 4. IMPLEMENTACI ´ ON SOBRE UN PAQUETE REAL 31 Figura 4.2: Gr´afica de crecimiento temporal lineal respecto del n´umero de muestras Figura 4.3: Gr´afica de crecimiento temporal cuadr´atico respecto de la longitud de las secuencias CAP´ ITULO 4. IMPLEMENTACI ´ ON SOBRE UN PAQUETE REAL 32 A continuaci´on, la siguiente prueba (cuyos resultados se pueden ver en la tabla 4.3) difiere de la anterior en que aqu´ı se considera el tiempo consumido para la verificaci´on de propiedades. Los resultados representan los segundos extra que hay que a˜nadir a la estructura inicial para validar 190 ecuaciones sencillas en l´ogica CTL. La verificaci´on de propiedades es realmente dependiente del m´etodo de exploraci´on escogido por el model checker (profundidad o anchura) y del primer estado donde la f´ormula CTL se torna falsa (esto es, rompiendo el proceso de verificaci´on y saliendo con un contraejemplo bajo el brazo). Por tanto, las ecuaciones de l´ogica temporal han sido generadas autom´aticamente, con el fin de que el resultado de su evaluaci´on sea lo m´as diverso posible. De esta forma, podremos considerar los resultados de la tabla 4.3 como una estimaci´on del coste medio introducido por el proceso de verificaci´on. NoMuestras Longitud Secuencias 500 1000 1500 2000 98 (ND4L) 8,75 15,04 26,08 38,57 115 (ND3) 4,48 14,11 25,25 36,89 174 (ND6) 8,92 14,65 25,91 35,47 318 (ND1) 6,33 23,14 34,65 46,58 347 (ND2) 8,89 14,88 29,13 28,74 459 (ND4) 9,28 8,66 27,81 50,05 603 (ND5) 13,28 20,37 29,05 38,88 Tabla 4.3: Coste temporal de verificaci´on de propiedades De esas 190 ecuaciones creadas, la mitad corresponden a f´ormulas del estilo: SPEC EF AG ((adn.secuencia[i] = A1) || (adn.secuencia[i] =-1)); La f´ormula es una variante de la definici´on de conservaci´on: busca un subgrupo de secuencias para el cual la posici´on iest´a completamente conservada en todos ellos y se corresponde con un amino´acido“A”. En particular, el valor “-1”, que en nuestro caso es el valor comod´ın por defecto de las secuencias, hace de guarda para forzar una b´usqueda m´as profunda y as´ı conseguir una estimaci´on m´as realista. El resto de f´ormulas sirven para ejemplificar la covariaci´on y tienen como patr´on: SPEC AF ((adn.secuencia[i] = B1)-> (EX (adn.secuencia[i+1] = V1))); CAP´ ITULO 4. IMPLEMENTACI ´ ON SOBRE UN PAQUETE REAL 33 La ecuaci´on equivale a decir: para todo nodo (especie) que en el futuro tenga un amino´acido “B” en la posici´on ide la cadena, alg´un descendiente directo suyo tendr´a un amino´acido en la posici´on i+ 1 que ser´a “V”. Ambos tipos de ecuaciones se comprueban para todos los valores de secuencia desde la posici´on i= 1 hasta la i= 95. Los benchmark han sido ejecutados sobre una m´aquina con procesador Intel Core 2 Duo E6750 a 2,66 GHz, 8 GB de RAM, Linux Debian con kernel 2.6.32 como sistema operativo, Cadence SMV versi´on 10.11.02p1 y el monitor de memoria memmon (del paquete UTILIB 4.1 [6]) para medir autom´aticamente el consumo de memoria m´aximo. La implementaci´on del software SMV es secuencial, as´ı que ´unicamente se ha utilizado uno de los dos cores del procesador. Para concluir este punto, desear´ıamos resumir dos aspectos importantes. En primer lugar, el sistema propuesto funciona aceptablemente r´apido para analizar genes aislados y con una carga notable de muestras, lo cual demuestra la viabilidad del entorno para estudios concretos. Sin embargo, el entorno escala mal con la longitud de las cadenas. De hecho, se ha comprobado que el sistema es incapaz de manejar cadenas que se mueven entorno a las decenas de miles de caracteres y unos pocos cientos de muestras (u ´ordenes de magnitud similares). Eso significa que, en el caso de desear verificar propiedades complejas y que engloben a la totalidad del ADNmt, ser´a necesario descomponerlas en varias proposiciones l´ogicas separadas m´as sencillas que afecten solamente a segmentos de tama˜no m´as reducido. Posteriormente habr´a que reconstruir la soluci´on final a partir de resultados parciales. Esta idea de descomposici´on autom´atica de proposiciones l´ogicas es innovadora porque hasta ahora los art´ıculos de investigaci´on publicados se enfocaban en reducir el tama˜no de la estructura de Kripke o la longitud de la computaci´on: por ejemplo, en nuestro contexto, dividir el total de muestras en subconjuntos m´as peque˜nos y verificar las propiedades sobre cada subconjunto por separado [9]. El principal problema latente es que al verificar cada subpropiedad independientemente, en la reconstrucci´on de la soluci´on final debemos comprobar que los caminos escogidos sean iguales. Cap´ıtulo 5 Representaci´on compacta de secuencias de ADN mediante diagramas de decisi´on 5.1. Contexto La representaci´on plana cl´asica de una secuencia de ADN, entendida como una sucesi´on lineal de caracteres, es inviable: incluso para cadenas de tama˜no relativamente peque˜no como el ADN mitocondrial (en comparaci´on con el ADN nuclear), el tratamiento de sus 16Kcaracteres se complica a la hora de procesar y almacenar conjuntos elevados de muestras. Esto queda especialmente patente en el coste de almacenamiento en memoria y en el coste de ejecuci´on de las herramientas de model checking (cap´ıtulo 4.2). Como soluci´on, proponemos el uso de diagramas de decisi´on para representar secuencias de texto. Los diagramas de decisi´on son una estructura de datos equivalente a los grafos ac´ıclicos dirigidos. Normalmente, los model checkers trabajan con una versi´on binaria, los diagramas de decisi´on binarios (BDD), para manejar grandes conjuntos de datos. En particular, utilizan los reduced ordered binary decision diagrams (ROBDD) [11], que son un caso especial donde se impone una relaci´on de orden entre las variables, para representar las estructuras de Kripke. Como principal ventaja, los ROBDD destacan por obtener la representaci´on can´onica m´ınima para un conjunto de elementos de acuerdo con una relaci´on de orden predefinida: los diagramas de decisi´on permiten solapar caminos comunes del grafo para consumir menos espacio, lo cual es bastante ventajoso para representar cadenas de texto en estudios poblacionales, donde m´as del 90 % del ADN de una especie est´a completamente conserva- 34 CAP´ ITULO 5. REPRESENTACI ´ ON COMPACTA DE SECUENCIAS 35 do en todos sus individuos. Adem´as, facilita la comparaci´on entre conjuntos aparentemente distintos gracias a la reducci´on a su forma can´onica. Sin embargo, encontrar la funci´on de ordenaci´on ´optima para los datos de entrada es un problema NP-completo. Por ´ultimo, los diagramas de decisi´on tambi´en colaboran impl´ıcitamente en el an´alisis de conservaci´on y covariaci´on de las secuencias. Por ejemplo, un an´alisis topol´ogico de los caminos cr´ıticos del grafo nos proporciona los nodos (posiciones de la cadena) m´as conservados. A su vez, gracias a la reducci´on del espacio de estados, los algoritmos actuales de covariaci´on se ver´an beneficiados: mejorar´an el rendimiento y potencialmente reducir´an el ruido de fondo que interfiere en los falsos positivos. Salvo que se indique lo contrario, a lo largo del presente cap´ıtulo usaremos indistintamente“secuencia”,“cadena (de texto)”o “palabra”como sin´onimos. Aunque el estudio est´a enfocado hacia el ADN, tambi´en es extensible a otro tipo de cadenas, como por ejemplo prote´ınas, palabras de diccionario, etc. 5.2. Estado del arte La aplicaci´on de diagramas de decisi´on a conjuntos de secuencias todav´ıa no ha sido ampliamente explorado, dado que existen otras estructuras de datos basadas en diccionario (como los prefix tree oradix tree [17]) que funcionan bastante bien y no necesitan buscar la funci´on de ordenaci´on ´optima para los datos de entrada. De hecho, el uso de diagramas de decisi´on binarios para representar conjuntos de secuencias es relativamente moderno [20]. Inicialmente, la propuesta m´as sencilla es definir una variable booleana por posici´on de la secuencia y car´acter del alfabeto, para indicar la presencia o ausencia de dicho car´acter en esa posici´on determinada de la cadena. En particular, como sus variables son binarias, se necesita un m´aximo de dlog2|Σ|e “bits” por posici´on para alfabetos con una cardinalidad superior a dos elementos, siendo Σ nuestro alfabeto. En [20] a esta aproximaci´on la llaman seqBDD. Tenemos un ejemplo de ella en la imagen 5.1. En ese art´ıculo, despu´es de aplicar ciertas reducciones, demuestran que el tama˜no de su diagrama de decisi´on es menor que el de un prefix tree equivalente, debido principalmente a la fusi´on de caminos comunes a lo largo del grafo y no ´unicamente de los prefijos compartidos. CAP´ ITULO 5. REPRESENTACI ´ ON COMPACTA DE SECUENCIAS 36 Figura 5.1: Representaci´on en ROBDD de una secuencia (SeqBDD) Sin embargo, su representaci´on naif es cuestionable y a´un podr´ıa mejorarse m´as porque no aprovecha toda la compactaci´on disponible. Los nodos, al ser binarios, sobrecargan el sistema cuando representan el alfabeto del ADN: el alfabeto del ADN (y en general, cualquier otro) tiene m´as de dos caracteres. Aumentar el rango de valores que admiten las variables de los diagramas de decisi´on, para pasar de binarias a multivaluadas, ser´ıa un primer paso suficiente. Para ello, el punto 5.3 estudia las distintas opciones de diagramas de decisi´on dise˜nados para valores multivaluados, con el fin de analizar c´omo de bien encajan en nuestro contexto. Junto con una serie de heur´ısticas expuestas en el punto 5.4, el apartado 5.5 demuestra que el coste de almacenamiento en memoria se ve reducido. 5.3. Herramientas y m´etodos Los data decision diagrams [14] y los multiple-valued decision diagrams [23, 24] son una ampliaci´on de los ROBDD para variables multivaluadas. Los set decision diagrams (SDD) [33, 18] son un caso especial de los anteriores, donde se admite etiquetar las aristas del diagrama con conjuntos de datos en vez de con una variable simple (la imagen 5.2 es un ejemplo de SDD). Esto sirve, sobretodo, para simplificar la notaci´on en determinados casos cuando los arcos de salida de todas las variables asociadas a la posici´on ide la cadena van a parar al mismo nodo siguiente (imagen 5.4). Por tanto, los SDD suponen el caso m´as general y la mejor opci´on para conseguir adaptar la estructura de datos a las cadenas de texto que proponemos. Como breve inciso, en la notaci´on a lo largo de este texto se han mantenido los arcos al CAP´ ITULO 5. REPRESENTACI ´ ON COMPACTA DE SECUENCIAS 37 terminal cero (no pertenencia al conjunto) por tradici´on con la representaci´on cl´asica de los diagramas de decisi´on, aunque es razonable omitirlos en futuras revisiones para aumentar la legibilidad, tal y como hacen los zero binary decision diagrams [25]. Las heur´ısticas aqu´ı propuestas est´an reunidas en dos grandes grupos. Por un lado, est´a la parte considerada como “tratamiento previo” de las secuencias, que hace referencia a la fase de reordenaci´on (punto 5.4.1) y c´omo organizar las cadenas para sacar el m´aximo provecho a priori. Existen otras alternativas y consideraciones suplementarias que explicamos en el punto 5.4.3. Por el otro, est´an las reglas de reducci´on (punto 5.4.2), que se aplican localmente una vez que hemos extra´ıdo toda la compactaci´on inicial. Dado que aplicamos las reducciones sobre una estructura tipo seqBDD (que primeramente traducimos a un diagrama de decisi´on multivaluado SDD, inherentemente menor), con |seqBDD|≤|prefix tree|por [20], y como las simplificaciones son an´alogas a las que convierten un prefix tree en un radix tree (por ejemplo, colapsando los nodos que tienen un s´olo hijo), todo esto implica que |SDD|≤|radix tree|. Es decir, conseguimos una cota superior m´as peque˜na para los diagramas de decisi´on multivaluados, lo que supone una mejora en el ahorro de memoria respecto implementaciones anteriores. 5.4. Heur´ısticas 5.4.1. Reordenaci´on inicial Supongamos que tenemos un pu˜nado de secuencias como las de la tabla 5.1. Idealmente, est´an alineadas correctamente de forma que las zonas comunes aparecen por bloques. En color rojo resaltan las columnas que poseen todos sus valores id´enticos, y en verde aquellas con el elemento m´as frecuente. Si seleccionamos el mayor n´umero posible de columnas que tengan todos sus componentes iguales, las colocamos al comienzo de la cadena y las fusionamos en un ´unico nodo del grafo, reduciremos la dimensi´on total considerablemente. CAP´ ITULO 5. REPRESENTACI ´ ON COMPACTA DE SECUENCIAS 38 1 2 3 4 5 s1A A CGC s2A A AGC s3A A TGC s4A A TGA s5A A TGT Tabla 5.1: Secuencias de ADN Gracias a esta reordenaci´on (figura 5.2), la implementaci´on sobre el diagrama de decisi´on disminuir´a potencialmente su n´umero total de nodos respecto a un seqBDD (incluso menos que un prefix tree): las sentencias de ADN que comparten rcolumnas iguales almacenar´an su cabecera en un mismo nodo, en vez de en rnodos independientes. Adicionalmente, retrasamos la bifurcaci´on de las ramas, evitando que los caminos diverjan demasiado pronto. Esto nos permitir´a aprovechar mejor la estructura de los SDD. Esta filosof´ıa de intercambio a priori de columnas, junto con un par de reducciones que expondremos en el cap´ıtulo 5.4.2, es un buen avance. La recuperaci´on posterior del contenido desde la estructura de datos ser´a sencilla porque sabemos la permutaci´on escogida de los datos en todo momento. No obstante, encontrar la mejor reordenaci´on global para un ´arbol de decisi´on es una tarea NP-completa. La ordenaci´on inicial se limita a resolver el problema del conjunto de subsecuencias com´un m´as largo [34], LCS por sus siglas en ingl´es. Existe un algoritmo de programaci´on din´amica para dos cadenas de texto cuyo tiempo de ejecuci´on es de orden polin´omico con la longitud de los datos de entrada. La extensi´on del problema a un conjunto de mcadenas es NP-completo, pero al estar preprocesado el multialineamiento como en el ejemplo, se simplifica enormemente hasta hacerlo tratable tambi´en en tiempo polinomial: mientras recorremos las secuencias, leeremos los bloques que vayan apareciendo (O(m∗ k), siendo mel n´umero de secuencias y ksu longitud). El objetivo final es simplificar el grafo mediante reordenaci´on, para ahorrar espacio en memoria y evitar registrar linealmente las mcadenas. En todo momento es indispensable registrar la permutaci´on elegida que transforma si=xi1xi2. . . xiken σ(si) = s0 i=x0 i1x0 i2. . . x0 ik. Nos permitir´a revertir posteriormente el efecto y definir la funci´on caracter´ıstica de las secuencias. Es evidente que la sobrecarga introducida por guardar el patr´on de permutaci´on es tanto menor cuanto mayor n´umero de secuencias haya. σ=xi1xi2. . . xik x0 i1x0 i2. . . x0 ik(5.1) CAP´ ITULO 5. REPRESENTACI ´ ON COMPACTA DE SECUENCIAS 39 La permutaci´on es ´unica y se aplica a todo el conjunto total del alineamiento. Es decir, por sencillez, no admitimos de momento reordenaciones locales para subconjuntos de cadenas. El ejemplo anterior quedar´a intercambiado de la siguiente forma: 1 2 4 5 3 s0 1AAGCC s0 2AAGCA s0 3AAGC T s0 4AAGAT s0 5AAGTT Tabla 5.2: Reordenaci´on inicial El conjunto de cadenas de texto S={s1, s2, . . . , sm}(definidas como “estados finales” o “alcanzables” de nuestro sistema) viene definido por la funci´on caracter´ıstica del alineamiento. Esto es, a cada cadena le asignamos una funci´on caracter´ıstica que la representa, y el conjunto Ses la uni´on l´ogica de todas las funciones caracter´ısticas de sus cadenas: χRS =PsiS χ(si). La funci´on caracter´ıstica de una cadena si=xi1xi2. . . xikes el valor resultante de la permutaci´on σ(): χ(si) = x0 i1x0 i2. . . x0 ik(5.2) Como bot´on de muestra, la ecuaci´on 5.3 refleja el tratamiento de las funciones. El operador l´ogico and (·) simboliza la concatenaci´on de subcadenas y el operador l´ogico or (+) la elecci´on. χRS =χ(s1) + χ(s2) + χ(s3) + χ(s4) + χ(s5) (5.3) =AAGCC +AAGCA +AAGCT +AAGAT +AAGTT =AAG ·[C·(C+A+T)+(A+T)T] Finalmente, en el menor grafo SDD posible habr´a tantas variables como subcadenas comunes. Gr´aficamente, el ´arbol resultado es el de la imagen 5.2. AAG est´a en la ra´ız. En el primer nivel de anidamiento, C·(C+A+T) y (A+T)·Tson las dos subramas principales. C·(C+A+T) est´a repartido en dos nodos, siendo Cel primero. Dada su condici´on de caracteres terminales, (C+A+T) quedar´an unidas en un mismo conjunto. {C, A, T}es el arco del ´ultimo nodo. An´alogamente se construye (A+T)·T. CAP´ ITULO 5. REPRESENTACI ´ ON COMPACTA DE SECUENCIAS 40 Figura 5.2: Set Decision Diagram 5.4.2. Reglas de reducci´on En esta secci´on definimos algunas reglas de simplificaci´on que permitir´an transformar un seqBDD en una estructura m´as compacta. Si aplicamos dichas reglas sobre un seqBDD m´ınimo (previamente reducido mediante sus propias heur´ısticas [20]), demostraremos impl´ıcitamente que nuestra nueva estructura es m´as peque˜na. La representaci´on de los seqBDD es binaria. Esto significa que cada car´acter del alfabeto tiene dos posibles valores en la estructura: cero o uno (pertenencia o no pertenencia a la secuencia). Por consiguiente, habr´a un m´aximo de nnodos por posici´on de la cadena, siendo |Σ|=n. El paso inicial consiste en ampliar el rango de valores de los diagramas de decisi´on utilizando los SDD expuestos anteriormente, ya conseguimos un gran avance. Mediante una funci´on sobreyectiva fi:Bn→Ntraducimos las nvariables booleanas de un seqBDD en la posici´on ide la cadena de texto en una ´unica variable num´erica. As´ı, de esta manera dividimos como m´aximo por nla dimensi´on inicial. La diferencia entre un radix tree y un prefix tree est´a en que, en el primero, los nodos con un solo hijo se fusionan con su descendiente para crear un nodo m´as grueso y as´ı economizar memoria por la sobrecarga de la estructura (gesti´on de punteros, etc). Por tanto, podemos aplicar aqu´ı mismo tambi´en la idea de fusi´on de nodos. Si las etiquetas aiyai+1 de la imagen 5.3 son un car´acter (o una subsecuencia de caracteres), pueden colapsarse en un ´unico nuevo arco conformado por la subcadena obtenida de su concatenaci´on. La fusi´on de nodos es la reducci´on m´as prioritaria, dado que interviene en la compactaci´on de caminos comunes del grafo. CAP´ ITULO 5. REPRESENTACI ´ ON COMPACTA DE SECUENCIAS 47 de O(m∗k) a O(m∗(k−α+ 1)) porque s´olo necesitamos almacenar una cabecera. ←k→ ←m−1→ ↑sim(m) sim−1(m) m... ↓si1(2) y1y2. . . ym Tabla 5.6: Paso de inducci´on con (m≤n) Por contra, si m>n, entonces |yn|>0. En consecuencia, todas las columnas del alineamiento tendr´an al menos un car´acter aiΣ tal que su soporte sea ∀j[1, k], sop(ai, j)>1. Esto significa que habr´a elementos duplicados del alfabeto en las columnas: habr´a bloques aislados del esquema 5.7 con soporte m´aximo sop(ai, j) = (m−n) para un car´acter. ←k→ ←n→ ↑sim(n+ 1) (m−n). . . ↓. . . ↑sin(n+1) n... ↓si1(2) y1y2. . . yn+1 Tabla 5.7: Paso de inducci´on con (m>n) Definimos una cota superior para averiguar el tama˜no m´aximo de este caso, ya que encontrar una aproximaci´on exacta es dif´ıcil por la distribuci´on aleatoria de los subbloques. Imaginemos por un momento que dichos subbloques no existen y ´unicamente compartimos caracteres comunes hasta la cabecera yn. Ser´a necesario un ´ultimo nodo adicional, tal y como se aprecia en la imagen 5.9. El arco final recoger´a Sn={sin[n+ 1..k], . . . , sim[n+ 1..k]}. El n´umero total de nodos es n+ 1 y el de arcos 2n. CAP´ ITULO 5. REPRESENTACI ´ ON COMPACTA DE SECUENCIAS 48 Figura 5.9: Esquema de la estructura final (|yn|>0) La longitud media m´ınima de la cabecera en la escalera inicial es de (n+1) 2. En el tramo restante de las (m−n) primeras secuencias, su longitud media es de nelementos. Simplificando, el coste es de α=n m(n+1) 2+ (m−n). En total, la ocupaci´on en memoria pasar´a de O(m∗k) a O(m∗(k−α+ 1)) por el mismo motivo que en el caso anterior. L´ogicamente, aumentar el n´umero de muestras mno penaliza de igual manera que la longitud kde la secuencia, dado que αdecrece proporcionalmente m´as r´apido con el tama˜no del conjunto de datos. En resumen, se puede decir que este apartado ha demostrado que, incluso en el caso peor, el conjunto de muestras mantiene su tama˜no tras un procesamiento en tiempo polin´omico. Mientras tanto, en media, la reducci´on habitual es proporcional al porcentaje de similitud absoluto de las secuencias: cuanto mayor sea el n´umero de muestras my menor sea su longitud k, mejor se comportar´a la compactaci´on. Cap´ıtulo 6 Conclusiones 6.1. Conclusiones El prop´osito de este trabajo ha sido enlazar conceptos entre dos mundos aparentemente alejados como es la filogen´etica y la verificaci´on formal. Esquem´aticamente, se ha hecho: 1. La interpretaci´on de un ´arbol filogen´etico tradicional como un sistema de transiciones, cuyos estados corresponden a caracter´ısticas compartidas por poblaciones de individuos definidas por un patr´on gen´etico com´un (herencia de ADN). Cuando miembros de esas poblaciones mutan y cambian alguna posici´on de su ADN (“evolucionan”), una transici´on se dispara entre estados del sistema. Entonces, se ha demostrado que existen transformaciones directas que asimilan un ´arbol filogen´etico a una estructura de Kripke equivalente sobre la que podr´an usarse t´ecnicas de model checking. 2. La definici´on de una l´ogica temporal adecuada para la definici´on de propiedades filogen´eticas, sem´antica en la que se basan las estructuras de Kripke derivadas de los ´arboles filogen´eticos. En el contexto de la filogen´etica, donde la evoluci´on se representa como una sucesi´on de eventos causales a lo largo de una filogenia, los operadores temporales con sem´antica branching-time son indispensables. Las propiedades temporales preguntan entonces sobre relaciones de parentesco (eventos pasados o futuros) en el ´arbol. 3. El uso de t´ecnicas est´andar de model checking para la verificaci´on autom´atica de propiedades filogen´eticas descritas en f´ormulas de l´ogica temporal sobre ´arboles filogen´eticos (descritos como estructuras de Kripke). Actualmente, el model checking es una t´ecnica madura con 49 CAP´ ITULO 6. CONCLUSIONES 50 una buena base te´orica y soportado por multitud de herramientas software. A nuestro parecer, esta aportaci´on es novedosa por la aplicaci´on de un entorno de chequeo de modelos a la filogen´etica. De ah´ı se desprenden una serie de conclusiones interesantes que son, a saber: 1. El estudio de propiedades filogen´eticas bajo diferentes modelos de evoluci´on. Se ha tomado el ´arbol filogen´etico por defecto como modelo est´andar a lo largo del trabajo; aunque f´acilmente pueden considerarse otros entornos alternativos (redes filogen´eticas, etc) que incluyan expl´ıcitamente eventos de recombinaci´on y otras transformaciones complejas. 2. El entorno formal de razonamiento, derivado de la l´ogica temporal escogida, ayuda a comprender la estructura subyacente de una propiedad y su relaci´on con el resto de propiedades. Entre otras cosas, esto permite descomponer la estructura subyacente de una propiedad en otras subproposiciones m´as sencillas, posiblemente verificadas y/o almacenadas anteriormente. De aqu´ı se puede lograr una mejora importante en el rendimiento del proceso de model checking. 3. La utilidad de los resultados de verificaci´on. Mientras que la verificaci´on de una propiedad verdadera no reporta gran inter´es, el fallo en el cumplimiento de una f´ormula involucra la generaci´on de contraejemplos con los estados conflictivos (o cadena de eventos) que los causan. Los filogenetistas pueden utilizar esos resultados para refinar sus propiedades o descubrir otras nuevas. 4. Por ´ultimo, pero no por ello menos importante, la importancia de software en el an´alisis filogen´etico. Gracias a que los algoritmos de verificaci´on son de prop´osito general, no est´an sujetos a grandes cambios. Ahora los esfuerzos que antes se dedicaban a los detalles de implementaci´on, ahora los investigadores podr´an dedicarlos a la especificaci´on de propiedades e interpretaci´on de resultados. 6.2. Trabajo futuro A pesar de las conclusiones positivas que hemos extra´ıdo de este trabajo, todav´ıa quedan algunos aspectos que mejorar. Por un lado, est´a el tema de la escalabilidad. Aunque se ha demostrado que el rendimiento es aceptable para analizar genes individuales, est´a la asignatura pendiente de expandir CAP´ ITULO 6. CONCLUSIONES 51 el sistema a propiedades complejas que involucren relaciones entre distintos trozos del ADNmt completo. Adem´as, puesto que el coste temporal se incrementa cuadr´aticamente con la longitud de la cadena, y el consumo de memoria se dispara para conjuntos de datos de tama˜no medio, el siguiente paso natural es paralelizar el entorno en segmentos de ADN m´as peque˜nos que puedan analizarse en un tiempo razonable. El rebanado de propiedades en subproposiciones m´as sencillas es el que m´as nos interesa (frente a la descomposici´on tradicional en subestructuras de Kripke o computaciones m´as sencillas y cortas). Nuestro objetivo m´as inmediato consiste en estudiar c´omo se componen los resultados parciales para recuperar la soluci´on global de una f´ormula CTL inicial, ya que adem´as de verificar cada subpropiedad independientemente hemos de comprobar que existe al menos un camino com´un a todas ellas. Por otro lado, queda la implementaci´on de la estructura de datos propuesta para mejorar la eficiencia en el almacenamiento de conjuntos de secuencias de ADN. Aunque se ha estudiado su viabilidad te´oricamente, y se han localizado las herramientas software con las que construir el sistema final, todav´ıa queda por delante la fase de implementaci´on, integraci´on y experimentaci´on dentro un model checker con pruebas de datos reales. Finalmente, como ´ultimo aspecto, se deja la puerta abierta para la definici´on y uso de l´ogicas m´as expresivas (por ejemplo, con la inclusi´on del “pasado”); o la extensi´on del conjunto de propiedades at´omicas para enunciar propiedades cuantitativas en filogen´etica. Nomenclatura ADNmt ADN Mitocondrial, p´agina 7 AP Atomic Propositions, p´agina 13 BDD Binary Decision Diagram, p´agina 34 CTL Computational Tree Logic, p´agina 16 LCS Longest Common Subsequence, p´agina 38 ROBDD Reduced Ordered Binary Decision Diagram, p´agina 34 SDD Set Decision Diagram, p´agina 36 seqBDD Sequence Binary Decision Diagram, p´agina 35 52 Bibliograf´ıa [1] Bioperl. http://www.bioperl.org. [2] Clustal. http://www.clustal.org/. [3] GenBank. www.ncbi.nlm.nih.gov/genbank/. [4] Cadence SMV Model Checker. http://www.kenmcmil.com/, 2001. [5] SMV Model Checker. http://www.cs.cmu.edu/~modelcheck/smv. html, 2001. [6] UTILIB. https://software.sandia.gov/trac/utilib, 2010. [7] Andr´e Arnold. Syst`emes de transitions finis et s´emantique des processus communicants. Masson, Paris, 1992. [8] Christel Baier and Joost-Pieter Katoen. Principles of model checking. The MIT Press, Cambridge, MA, 2008. [9] Sergey Berezin, S´ergio Campos, and Edmund Clarke. Compositional reasoning in model checking. In Willem-Paul de Roever, Hans Langmaack, and Amir Pnueli, editors, Compositionality: The Significant Difference, volume 1536 of Lecture Notes in Computer Science, pages 81– 102. Springer Berlin / Heidelberg, 1998. [10] Roberto Blanco, Gregorio de Miguel Casado, Jos´e Ignacio Requeno, and Jos´e Manuel Colom. Temporal logics for phylogenetic analysis via model checking. In IEEE International Conference on Bioinformatics & Biomedicine 2010, December 2010. [11] Randal E. Bryant. Symbolic manipulation of boolean functions using a graphical representation. In IEEE 22nd Design Automation Conference, 1985. 53 BIBLIOGRAF´ IA 54 [12] Randal E. Bryant. Graph-based algorithms for Boolean function manipulation. IEEE T. Comput., C-35:677–691, Aug. 1986. [13] Edmund M. Clarke and E. Allen Emerson. Design and synthesis of synchronization skeletons using branching time temporal logic. In Proc. Workshop on Logics of Programs, pages 52–71, Yorktown Heights, NY, May 1981. [14] Jean-Michel Couvreur, Emmanuelle Encrenaz, Emmanuel Paviot-Adet, Denis Poitrenaud, and Pierre-Andr´e Wacrenier. Data decision diagrams for petri nets analysis. Technical report, LIP6, 2002. [15] Joseph Felsenstein. Inferring phylogenies. Sinauer, Sunderland, MA, 2003. [16] Orna Grumberg and Helmut Veith, editors. 25 years of model checking: history, achievements, perspectives. Springer, Berlin, 2008. [17] Donald E. Knuth. Art of Computer Programming, Volume 3: Sorting and Searching, chapter 6.3. Addison-Wesley, 1998. [18] Alexandre Hamez; Yann Thierry-Mieg; Fabrice Kordon. Hierarchical set decision diagrams and automatic saturation. Technical report, LIP6, 2008. [19] Christopher James Langmead and Sumit Kumar Jha. Predicting protein folding kinetics via temporal logic model checking. In Proc. WABI 2007, pages 252–264, Philadelphia, PA, Sep. 2007. [20] Elsa Loekito, James Bailey, and Jian Pei. A binary decision diagram based approach for mining frequent subsequences. Knowledge and Information Systems, pages 1–34, September 2009. [21] Zohar Manna and Amir Pnueli. The temporal logic of reactive and concurrent systems: specification. Springer, Berlin, 1991. [22] Kenneth L. McMillan. Symbolic Model Checking. Kluwer Academic Publishers, Norwell, MA, USA, 1993. [23] D. Michael Miller and Rolf Drechsler. Implementing a multiple-valued decision diagram package. Technical report, Department of Computer Science, University of Victoria, Canada, 1998. BIBLIOGRAF´ IA 55 [24] D. Michael Miller and Rolf Drechsler. On the construction of multiplevalued decision diagrams. Technical report, Department of Computer Science, University of Victoria, Canada, 2002. [25] Shin-ichi Minato. Zero-suppressed bdds for set manipulation in combinatorial problems. In Design Automation, 1993. 30th Conference on, pages 272 – 277, June 1993. [26] Pedro T. Monteiro, Delphine Ropers, Radu Mateescu, Ana T. Freitas, and Hidde de Jong. Temporal logic patterns for querying dynamic models of cellular interaction networks. Bioinformatics, 24:i227–i233, 15 Aug. 2008. [27] Julio Montoya, Ester L´opez-Gallardo, Carmen D´ıez-S´anchez, Manuel J. L´opez-P´erez, and Eduardo Ruiz-Pesini. 20 years of human mtDNA pathologic point mutations: carefully reading the pathogenicity criteria. Biochim. Biophys. Acta, 1787:476–483, May 2009. [28] Michael Arnold; Enno Ohlebusch. Algorithmica, chapter Linear Time Algorithms for Generalizations of the Longest Common Substring Problem, pages 1–13. Springer New York, 2009. [29] Martin B. Richards, Vincent A. Macaulay, Hans-J¨ urgen Bandelt, and Bryan C. Sykes. Phylogeography of mitochondrial DNA in western Europe. Ann. Hum. Genet., 62:241–260, May 1998. [30] Aur´elien Rizk, Gr´egory Batt, Fran¸cois Fages, and Sylvain Soliman. On a continuous degree of satisfaction of temporal logic formulae with applications to systems biology. In Proc. CMSB 2008, pages 251–268, Rostock, Oct. 2008. [31] Richard Rudell. Dynamic variable ordering for ordered binary decision diagrams. In ICCAD ’93: Proceedings of the 1993 IEEE/ACM international conference on Computer-aided design, pages 42–47, Los Alamitos, CA, USA, 1993. IEEE Computer Society Press. [32] The Y Chromosome Consortium. A nomenclature system for the tree of human Y-chromosomal binary haplogroups. Genome Res., 12:339–348, Feb. 2002. [33] Jean-Michel Couvreur; Yann Thierry-Mieg. Hierarchical decision diagrams to exploit model structure. Technical report, LIP6, 2005. BIBLIOGRAF´ IA 56 [34] V. G. Timkovskii. Complexity of common subsequence and supersequence problems and related problems. Cybernetics and Systems Analysis, 25(5):565–580, 1989.