Inferencia bayesiana en filogenia

La inferencia bayesiana en filogenia genera la probabilidad posterior de un parámetro, un árbol filogenético y/o un modelo evolutivo, basada en la probabilidad anterior de ese parámetro y la función de verosimilitud de los datos. La aplicación del análisis bayesiano en la inferencia filogenética presenta varias ventajas en comparación con otros métodos de inferencia, como la fácil interpretación de los resultados, la posibilidad de usar información apriorística​ y algunas ventajas computacionales.​ Se calcula la probabilidad de que nuestro árbol sea correcto condicionada por los datos que tenemos: P (árbol|datos). Lo contrario a otros métodos, que calculan la probabilidad de que nuestros datos se adapten al árbol: P (datos|árbol).

Fundamento teórico

Probabilidad posterior

La estadística bayesiana, a diferencia de la estadística clásica, considera los parámetros de un modelo como variables aleatorias con una distribución estadística y no como constantes desconocidas. Puesto que desconocemos los valores de los parámetros del modelo, desde un punto de vista bayesiano resulta más razonable especificar a priori unas distribuciones para describir los valores que éstos pueden tomar. La Inferencia Bayesiana (IB) de la filogenia se basa en el cálculo de la probabilidad posterior de un árbol, esto es, la probabilidad de que dicho árbol sea correcto dados unos datos y un modelo evolutivo.

La probabilidad posterior se calcula mediante el Teorema de Bayes, según el cual la probabilidad posterior de un árbol P(H|D) es proporcional a su probabilidad previa P(H) multiplicada por su verosimilitud P(D|H) (o probabilidad de los datos dados un árbol y modelo evolutivo). El cálculo de la probabilidad posterior implica evaluar todos los posibles árboles y, para cada árbol, investigar todas las posibles combinaciones de longitudes de ramas y parámetros del modelo evolutivo. Sin embargo, este cálculo resulta prohibitivo computacionalmente. Se debe tener en cuenta que, aunque, si bien es cierto que un modelo más rico en parámetros siempre tiene una mayor verosimilitud que un submodelo anidado, su verosimiltud total no tiene por qué ser superior. La razón es que un modelo más rico en parámetros también tiene un espacio de parámetros más amplio y, por tanto, una menor densidad de probabilidad a priori.

Algoritmo empleado

En el caso de la IB, utilizamos algoritmos MCMC (o Markov Chain Monte Carlo) para obtener una aproximación de la distribución de la probabilidad posterior. Los MCMC son algoritmos de simulación estocástica que evitan el cálculo directo de las probabilidades posteriores y permiten obtener una muestra de la distribución posterior. De manera análoga a cómo opera la búsqueda heurística por MV (máxima verosimilitud), la cadena MCMC de la IB parte de un árbol al azar (combinación al azar de topología, longitudes de rama y parámetros del modelo) y realiza un cierto número de visitas a árboles concretos. De este modo, los algoritmos MCMC modifican ligeramente una de las variables del árbol inicial (ya sea la posición o longitud de una rama en concreto, o el valor de un parámetro del modelo evolutivo) y evalúan el nuevo árbol. Esto es lo que se conoce como una generación. Tras cada generación, el cambio propuesto se acepta o rechaza en función del valor de la razón entre las probabilidades posteriores de los estados actual y anterior. Si esta razón es mayor que uno (el nuevo estado es más probable), el cambio se acepta y la cadena MCMC continúa desde este punto. Si por el contrario, la razón es menor de uno, el cambio puede aceptarse o rechazarse en función del valor que toma la razón entre probabilidades posteriores.

Fases en la búsqueda de árboles filogenéticos mediante Inferencia Bayesiana

Debemos tener en cuenta que la aproximación de la probabilidad posterior sería más exacta cuantas más generaciones se evalúen en la cadena de Markov. Durante las generaciones iniciales de las cadenas MCMC, los árboles suelen tener una probabilidad posterior baja como resultado de haber comenzado a partir de combinaciones aleatorias de topología, longitudes de rama y valores de parámetros. Esto es lo que conocemos como la fase de burnin o calentamiento. Tras esta fase inicial, las cadenas MCMC alcanzan una fase estacionaria donde los árboles muestreados tienen una probabilidad posterior elevada. Para la estimación del árbol final mediante IB es necesario descartar los árboles del burnin y todo el resto de árboles con probabilidad posterior elevada se resumen mediante un árbol consenso.​

Para trabajar correctamente en la práctica desde este enfoque, se deben tener en cuenta algunas consideraciones.

Modelo representativo de la búsqueda topológica en filogenética bayesiana mediante el algoritmo Monte Carlo Cadenas de Markov

• En primer lugar, el mecanismo de modificación de estado de la cadena debe ser diseñado con cuidado. Por un lado, los cambios sugeridos deben ser bastante pequeños a fin de evitar una caída drástica en la probabilidad posterior. Por otro lado, las modificaciones deben ser lo suficientemente grandes como para asegurar que la cadena es capaz de explorar un porción suficiente de espacio de parámetros en un período razonable de tiempo.

• En segundo lugar, como la cadena MCMC comienza en una ubicación aleatoria en el espacio de parámetros, generalmente de muy baja probabilidad posterior, se le debe dar tiempo para llegar a la estacionalidad, un estado en el que la cadena es en realidad muestreo de la probabilidad posterior. En la práctica, aproximadamente las primeras 10.000 generaciones, quedarán excluidas de análisis posteriores.

• En tercer lugar, los árboles que se encuentran cerca unos de otros están altamente correlacionados en su construcción, ya que sólo se diferencian por pequeñas modificaciones. Para evitar este problema de autocorrelación, el conjunto final de los árboles de salida se obtiene, por ejemplo, manteniendo el último árbol cada 1000 muestreos que realiza la cadena MCMC.

Extensión del algoritmo MCMC

Como ocurre con las búsquedas heurísticas mediante el método de máxima verosimilitud, la cadena MCMC puede quedar atrapada en máximos locales, de manera que dos análisis con los mismos criterios de búsqueda pueden presentar resultados diferentes para el mismo conjunto de datos. Para evitarlo, la IB utiliza varias cadenas MCMC simultáneamente y además se recomienda hacer cada análisis de IB por duplicado. Este algoritmo se conoce como Metropolis-coupled Montecarlo vía Cadenas de Markov (MC3). Los programas de IB como MrBayes generalmente implementan cuatro cadenas MCMC, una cadena “fría" y tres “calientes". Las cadenas “calientes" elevan la probabilidad posterior de los árboles visitados a una potencia entre 0 y 1, favoreciendo así el salto entre picos locales de probabilidad posterior. Durante el proceso de búsqueda, las cadenas se intercambian, de modo que en cada momento la cadena con mayor probabilidad posterior se comporta como “fría" y el resto como “calientes". Esto favorece una exploración más exhaustiva del conjunto de árboles posibles, puesto que se facilita el salto entre picos locales de probabilidad posterior elevada y con ello se reduce la posibilidad de que la cadena “fría" quede atrapada en máximos locales. Para el resultado final, sólo se tiene en cuenta los resultados de la cadena “fría", pues incluye a los árboles con mayor probabilidad posterior entre todos los visitados por las cuatro cadenas.

El factor "temperatura de las cadenas" es, por tanto, un parámetro importante a tener en cuenta en la búsqueda filogenética mediante IB. Si aplicamos un exceso de temperatura, entonces las cadenas que se mueven en los paisajes calentados caminarán por todo el espacio y la probabilidad de estar en un pico interesante, cuando se intercambie el estado con la cadena fría, será menor. Por lo tanto, la mayoría de los intercambios serán rechazados y la temperatura no aumentará la eficacia de mezcla con la cadena fría. Por otro lado, si no calentamos lo suficiente, entonces las cadenas serán muy similares, y las cadenas calentadas no intercambiarán información a una velocidad adecuada que la cadena fría.

Para saber si las cadenas MCMC de nuestro análisis de IB han explorado de modo eficiente el conjunto de árboles es necesario estudiar la convergencia de las cadenas “frías" entre los dos análisis de IB. La forma más sencilla es observar sus diferencias medias, que en MrBayes conocemos como Average Standard Deviation of Split Frequencies. Una diferencia entre las desviaciones estándar medias menor de 0.01 suele indicarnos que los dos análisis de IB han convergido. Otra herramienta adecuada de diagnóstico de convergencia es el llamado Potential Scale Reduction Factor (PSRF) originalmente propuesto por Gelman y Rubin (1992). El PSRF compara la varianza entre runs con la varianza dentro de los runs. Si se inician las cadenas en puntos de partida dispersos, la varianza entre runs será inicialmente mayor que la varianza dentro de las mismas. Sin embargo, la convergencia de las cadenas hará que esta varianza se asimile y el PSRF se acercará a 1. Sin embargo, es muy recomendable además estudiar la convergencia de las cadenas MCMC a posteriori con herramientas adecuadas desarrolladas a tal efecto, como AWTY.​

Por último, señalar que algunos programas como MrBayes nos permiten estudiar los resultados mediante medidas de calidad de modelos estadísticos como son: el criterio de información de Akaike, el criterio de información Bayesiano o el factor Bayes. Además, una alternativa a la ejecución de un análisis completo de cada modelo escogiendo luego entre ellos, en función de las probabilidades estimadas y factor Bayes que presentan los mismos, consiste en llevar a cabo un solo análisis bayesiano que explore los modelos en un espacio predefinido de éstos (reversible-jump MCMC). Además, en la práctica, también se pueden emplear varios tipos de modificaciones en el flujo de trabajo. Los métodos de intercambio de ramas como NNI (nearest neighbor interchange), SPR (subtree prune and regraft) y TBR (tree bisection and reattachment) se emplean en el contexto de los métodos de máxima parsimonia y máxima verosimilitud para la búsqueda a través del espacio de topologías de los árboles filogenéticos; a efectos de la construcción de una cadena MCMC se necesita acoplar un método para proponer estas modificaciones y los parámetros del modelo evolutivo.​

Fundamento práctico

Teorema de Bayes

La inferencia bayesiana se basa en la interrelación cuantitativa entre la función de verosimilitud y las distribuciones anteriores y posteriores de probabilidad. Esta interrelación viene dada por el Teorema de Bayes, el cual permite calcular la probabilidad posterior a partir de la verosimilitud y probabilidad anterior de los datos.

$P(H|D)={\frac {P(H)P(D|H)}{P(D)}}$

$H$=Hipótesis o parámetro

$D$=Datos

$P(D|H)$=Verosimilitud de los datos, probabilidad de H (o valor del parámetro)dados D

$P(H)$=Probabilidad anterior, es la probabilidad incondicional de H

$P(D)$=Probabilidad incondicional de D, que puede ser obtenida usando la ley de la probabilidad total, calculando ΣH P(H) P(D|H). Funciona como constante normalizadora,asegurando que la sumatoria de $pP=1$

La probabilidad posterior de un árbol filogenético $(A_{i})$ condicional a un alineamiento múltiple $(B)$ puede ser calculado mediante el Teorema de Bayes:

$P(A_{i}|B)={\frac {P(B|A_{i})P(A_{i})}{\sum _{j=1}^{T(e)}P(B|A_{j})P(A_{j})}}$

Siendo la función de verosimilitud:

$P(B|A_{i})=\int _{l}\int _{s}P(B|A_{i},l,s)P(l,s)\mathrm {d} l\mathrm {d} s$

El sumatorio es sobre todos los $T(e)$ árboles posibles para $e$ especies, así como, las integrales son sobre todas las combinaciones de longitud de rama $(l)$ y parámetros de sustitución $(s)$.

La probabilidad anterior de un árbol suele establecerse como: $P(A_{i})={\frac {1}{T(e)}}$

Montecarlo vía cadenas de Markov

El sumatorio y las integrales requeridas en el análisis bayesiano filogenético no pueden ser evaluados analíticamente por ello se requiere de métodos no determinísticos para aproximar las probabilidades posteriores de los árboles. El método más usado es el de Montecarlo vía Cadenas de Markov (MCMC, siglas en inglés) que se basa en el muestreo de una distribución de probabilidad simulada. Su algoritmo parte de un estado inicial aleatorio de la cadena, y a continuación se propone un nuevo estado generado estocásticamente. Se calcula el cociente (R) de la función de densidad probabilística de ambos estados, usando el algoritmo Metropolis-Hastings. La decisión de aceptar T' o mantener Ti se basa en la relación de las probabilidades a posteriori de los dos árboles, que viene dada por:

$R={\frac {P(A'|B)}{P(A_{i}|B)}}={\frac {P(B|A')P(A')P(B)}{P(B|A_{i})P(A_{i})P(B)}}={\frac {P(B|A')P(A')}{P(B|A_{i})P(A_{i})}}$

Siendo $A_{i}$ el árbol filogenético en una posición determinada de la cadena de Markov, $A'$ la modificación propuesta y $B$ el conjunto de datos alineados.

• Si R ≥1 se acepta el nuevo estado y la cadena es actualizada.
• Si R< 1, se genera un número aleatorio entre 0 y 1. Si es mayor que R se rechaza el nuevo estado.

Por tanto, el MCMC trabaja de la siguiente manera:

• Comienza una cadena markoviana con un árbol ya sea elegido al azar o elegido por el investigador.

• Se propone un nuevo árbol: el proceso de cambio del árbol 1 al 2 debe satisfacer las siguientes condiciones:

1) El mecanismo debe ser estocástico

2) Cada árbol posible debe ser obtenido por aplicaciones repetidas del mismo mecanismo.

3) La cadena debe ser aperiódica.

Este proceso se puede repetir varios miles o millones de veces. La probabilidad posterior se puede aproximar a la proporción de veces que un árbol es visitado por la cadena.

Metropolis-coupled Montecarlo vía cadenas de Markov

El MCMC presenta limitaciones cuando la distribución presenta varios picos separados valles profundos, pudiendo la cadena quedar atascada en algún pico local. El algoritmo Metropolis-coupled Montecarlo vía Cadenas de Markov (MC3, siglas en inglés) permite a la cadena saltar valles profundos. En esta variante del MCMC se corren $n$ cadenas en paralelo(normalmente 4), donde $n-1$ son cadenas calentadas. Las cadenas tienen una distribución $P(A_{i}|B)^{C}$. Cuando C es igual a 0 todos los árboles tiene la misma probabilidad, pudiendo las cadenas visitarlos con total libertad. En la cadena fría C es igual a 1, siendo, por tanto, su distribución igual a la distribución de interés. Por ejemplo, se puede aplicar un calentamiento gradual con forma:

$C={\frac {1}{1+(i-1)\alpha }}$

Cuando todas las cadenas han terminado un ciclo se produce un intercambio entre dos cadenas elegidas al azar. Si el cambio se acepta las dos cadenas son actualizadas. La inferencia se realiza basada únicamente en los estados muestreados por la cadena fría. Mediante este método se permite que la cadena fría (que quizás este bloqueada en un pico local) colonice nuevos picos, si se producen intercambios exitosos con cadenas calentadas, que tienen mayor libertad para explorar todo el espacio topográfico. Por tanto, este algoritmo implica:

• Correr algunas cadenas independientemente.

• La primera cadena que cuenta es la fría el resto se denomina cadenas accesorias (calentadas, con más capacidad de exploración).

• Los cambios de estado se establecen al azar entre dos cadenas distintas.

• Se necesita correr varios análisis independientes para confirmar convergencias.

Inferencia Bayesiana y otros métodos computacionales

En la inferencia filogenética bayesiana el objetivo no es determinar un solo árbol filogenético de máxima probabilidad posterior (el objetivo de la MV), sino más bien calcular una muestra de árboles filogenéticos de acuerdo con la probabilidad posterior de distribución de los mismos. Dicha muestra de árboles se procesa adicionalmente, para generar un único árbol de consenso. Esta integración de árboles otorga la ventaja de incluir la incertidumbre presente en la determinación de la topología, longitudes de rama y parámetros del modelo evolutivo.

Sin embargo, los métodos de IB tienen una conexión fuerte con la MV, ya que la hipótesis preferida es aquella con la mayor probabilidad posterior, y éste es un valor que se calcula en función de la verosimilitud. Puesto que la IB se basa directamente en la verosimilitud, comparte sus propiedades de eficiencia y consistencia.

Algunos investigadores prefieren usar IB sobre MV debido a que el primer método utiliza “atajos” para los cálculos al emplear el algoritmo Markov Chain Monte Carlo, el cual permite realizar búsquedas a través de un número menor de árboles según sus valores de probabilidades posteriores. Esto permite que la IB demande menos poder computacional y sea más rápida que MV. Además, no es tan necesario en IB el empleo de métodos posteriores de análisis de bootstrap para comprobar la robustez del árbol.

Por otro lado, la IB nos permite introducir cierto conocimiento previo en el análisis de datos, ya que resulta necesario especificar la distribución previa de varios parámetros​ (aunque, en general, el conocimiento de la mayoría de los parámetros es escaso). De hecho, las áreas de investigación en IB aplicada a filogenias indagan en cómo asegurar la convergencia de cadenas y también cómo escoger las distribuciones a priori apropiadas.

Sin embargo, ningún método de análisis filogenético se encuentra exento de críticas:

• La máxima parsimonia se ve afectada por el fenómeno conocido como atracción de ramas largas, y su concepto teórico no resulta aplicable en muchos casos.

• El método de máxima verosimilitud suele dar lugar a fenómenos de “repulsión” de grupos hermanos cuando estos se ubican en ramas largas de los árboles.​

• La MV y la IB son poco fiables cuando las tasas de evolución de ADN no son homogéneas en el tiempo ni entre linajes.​

Debido a la gran dificultad que se presenta a la hora de escoger el método a implementar en una reconstrucción filogenética, un procedimiento conciliatorio consiste en emplear los tres métodos. Si las topologías de los árboles obtenidos, usando los diferentes métodos, son concordantes, la hipótesis filogenética resultante es considerada robusta

Software más utilizado

• MrBayes
• BEAST

Véase también

• Filogenia
• Análisis moleculares de ADN
• Filogenética computacional

Bibliografía

• Hall, Barrie G. "Pgylogenetic Trees Made Easy: A How-To Manual". Second Edition. Publishers Sunderland, Massachusetts. (2004)
• Holder, M. and Lewis, P. O. “Phylogeny estimation: traditional and bayesian approaches," Nature Reviews Genetics, vol. 4, pp. 275-284. (2003)
• Huelsenbeck, John P., Fredrik Ronquist, Rasmus Nielsen, and Jonathan P. Bollback. “Bayesian Inference of Phylogeny and Its Impact on Evolutionary Biology.” Science 294, no. 5550 (December 14, 2001): 2310–2314. doi:10.1126/science.1065889.
• Huson, David H.; Rupp Regula and Scornavacca Celine."Phylogenetic Networks: Concepts, Algorithms and Applications." Cambridge University Press. (2011)
• Lemey, Philippe; Salemi, Marco and Vandamme, Anne-Mieke. "The Phylogenetic Handbook: A Practical Approach to Phylogenetic Analysis and Hypothesis Testing". Cambridge University Press. (2009)
• Nylander J. A.; Wilgenbusch, J. C.; Warren, D. L. and Swofford, D. L. “Awty (are we there yet?): A system for graphical exploration of mcmc convergence in bayesian phylogenetics,". Bioinformatics, vol. 24, no. 4, pp. 581-583. (2008)
• Peña, C. "Métodos de inferencia filogenética". Rev. peru. biol. 18(2): 265 - 267. (2011)
• Ronquist, F.; van der Mark P. and Huelsenbeck, J. P. “Bayesian phylogenetic analysis using MrBayes”. Second edition, ed. New York: Cambridge University Press, pp. 210-266. (2009)
• Yang, Ziheng and Rannala, Bruce. “Molecular Phylogenetics: Principles and Practice.” Nature Reviews Genetics 13, no. 5 (May 2012): 303–314. doi:10.1038/nrg3186.
• Yang, Z. and Rannala, B.. "Bayesian phylogenetic inference using DNA sequences: A Markov chain Monte Carlo method." Molecular Biology and Evolution, 14, 717–724.(1997)

• Datos: Q4874472