Tutoral de uso de PhyML v3 para selección de modelos y búsquedas intensivas de árboles
Pablo Vinuesa, CCG-UNAM, Cuernavaca, México. GitHub repos; @pvinmex
2024-11-23
- 1 Selección y análisis de modelos paramétricos y empíricos mediante LRTs, AIC y BIC con \(PhyML\) desde la línea de comandos de \(Linux\)
- 2 NOTAS DE INSTALACION DE PHYML DESDE EL REPOSITORIO GitHub
- 3 SELECCION DE MODELOS PARAMÉTRICOS DE
DNA USANDO PHYML - Prácticas
- 3.1 Descarga de secuencias para la práctica
- 3.2 Pruebas de razones de
verosimilitud (\(LRTs\)) entre pares de
modelos anidados
- 3.2.1 Estima del árbol NJ de referencia, sobre el que estimaremos los parámetros de todos los modelos en competición
- 3.2.2 Optimización de parámetros y
estimas de \(lnL\) para una gama de
modelos anidados a evaluar
- 3.2.2.1 Jukes y Cantor de 1969 (\(JC\) o \(JC69\))
- 3.2.2.2 El modelo de Kimura de 1980 (\(K80\)) o Kimura de 2 parámetros (\(K2P\))
- 3.2.2.3 Modelo de Felsenstein de 1981 (\(F81\))
- 3.2.2.4 Modelo de Hasegawa, Kishino y Yano de 1985 (\(HKY85\))
- 3.2.2.5 Modelo de Hasegawa, Kishino y Yano de 1985 + distribución Gamma (\(HKY85+G\))
- 3.2.2.6 Modelo de Hasegawa, Kishino y Yano de 1985 + proporción de sitios invariantes (\(HKY85+I\))
- 3.2.2.7 Modelo de Hasegawa, Kishino y Yano de 1985 + distribución Gamma + proporción de sitios invariantes (\(HKY85+G+I\))
- 3.2.3 Análisis gráfico y estadístico de resultados usando herramientas de \(Linux\) y \(R\)
- 3.2.4 Estima de la filogenia bajo el mejor modelo encontrado: HKY85+G
- 3.2.5 Despliegue y edición del árbol con figtree
- 3.3 Estrategia alternativa de
selección de modelos de DNA: reducción de parámetros del modelo GTR+G+I,
y evaluación de los modelos resultantes mediante \(AIC\)
- 3.3.1 Ajuste del modelo \(GTR+G+I\) a primates.phy
- 3.3.2 Exploración del archivo stats para identificar posibles simplificaciones al modelo GTR+G+I
- 3.3.3 Evaluación de modelos simplificados de 2, 3 y 4 tasas
- 3.3.4 Cálculo de los \(AICi\) y \(\Delta_i\) (de cada modelo) con respecto al
\(AICmin\) en
R - 3.3.5 Conclusiones
- 3.3.6 Tarea de ajuste de modelos paramétricos mediante criterios de información AIC y BIC
- 4 Cálculo de \(AICs\) y \(BICs\) para selección modelos de
aminoácidos
- 4.1 Evaluación rápida de un subconjunto de matrices empíricas con un \(bucle\ for\) y un \(hash\)
- 4.2 Análisis gráfico y numérico de resultados de \(AIC\) usando herramientas de filtrado de \(Linux\) y \(R\)
- 4.3 Búsqueda intensiva bajo el mejor modelo, usando rearreglos tipo BEST y 5 árboles de semilla
- 5 Selección automatizada de modelos de
DNA y proteína con phyml_DNAmodelFinder.sh y
phyml_protModelFinder.sh
- 5.1 Ejemplo de invocación del script phyml_protModelFinder.sh
- 5.2 El criterio bayesiano de
información (\(BIC\))
- 5.2.1 Cálculo de distancias (\({\Delta}AIC_i\) o \({\Delta}BIC_i\)) y ponderaciones de Akaike (\(AIC_{wi}\)) o bayesianas (\(BIC_{wi}\)) para cada modelo en competición
- 5.2.2 Cálculo de ponderaciones de cada modelo y su peso acumulado en base a un conjunto de valores ordenados de (\(\Delta BIC\))
- 5.2.3 Código de AWK para calcular los pesos y pesos acumulados de cada modelo
- 5.2.4 Tarea: selección de modelos con BIC
- 5.3 Selección “inteligente” de subconjuntos de modelos de sustitución para DNA con phyml_DNAmodelFinder.sh
- 6 Comparación del rendimiento (benchmarking) de phyml_DNAmodelFinder con jModelTest2
- 7 Cita para este tutorial y los scripts phyml_DNAmodelFinder y phyml_protModelFinder
1 Selección y análisis de modelos paramétricos y empíricos mediante LRTs, AIC y BIC con \(PhyML\) desde la línea de comandos de \(Linux\)
1.1 Presentación
Estas prácticas fueron diseñadas originalmente para los Talleres Internacionales de Bionformática, TIB2022. Este tutorial está basado en el material del repositorio GitHub - TIB-filoinfo.
La versión actual, es una adaptación para el Tópico selecto de: Introducción al Biocómputo en Sistemas Linux para Genómica, impartido para alumn@s de la Facultad de Ciencias - UNAM que cursan el N3 del Taller de Ciencias Genómicas - de Moléculas a Ecosistemas (semestre 2023-1), impartido en el Centro de Ciencias Genómicas - UNAM.
El material fue extendido para mi curso en la Licenciatura en Ciencias Genómicas, LCG-UNAM, incluyendo los \(scripts\) phyml_DNAmodelFinder.sh y phyml_ProtModelFinder.sh que acompañan a este tutoral.
1.2 Objetivos
En esta práctica aprenderemos a usar \(phyml\) para seleccionar modelos de sustitución, paramétricos (DNA) y empíricos (prot), usando pruebas de cocientes de verosimilitudes (\(LTRs\)) y métodos basados en criterios de información como los de Akaike (\(AIC\)) y bayesiano (\(BIC\)). Usaremos scripts sencillos de \(R\) para despliegue gráfico de resultados y análisis numérico de valores de \(lnL\) para los modelos en competición.
A lo largo de las prácticas iremos optimizando los comandos usados, agregando \(blucles\), \(arreglos\), \(hashes\) entre otros, para hacer más eficiente el análisis, culminando con la presentación de los \(scripts\) phyml_DNAmodelFinder y phyml_protModelFinder, con los que automatizaremos el proceso. Ambos programas útiles para la investigación en filogenética.
Además, aprenderemos a hacer búsquedas intensivas bajo el mejor modelo seleccionado, usando múltiples árboles de semilla, así como el despliegue y edición de los mejores árboles encontrados con \(FigTree\).
Al concluir el tutorial habrán adquirido una sólida formación en selección de modelos evolutivos para filogenética, y aprendido o reforzado aspectos fundamentales de biocómputo desde la línea de comandos en un entorno de \(Linux\).
1.3 Licencia
El material del Taller, TIB-filoinfo lo distribuyo públicamente a través de este repositorio GitHub bajo la Licencia No Comercial Creative Commons 4.0
This
work is licensed under a
Creative
Commons Attribution-NonCommercial 4.0
2 NOTAS DE INSTALACION DE PHYML DESDE EL REPOSITORIO GitHub
Para correr las prácticas que siguen a continuación es necesario que tengas instalada una versión reciente de \(phyml\), idealmente v3.3.20220408.
1. Descarga la ultima version estable desde el [repositorio GitHub de phyml](https://github.com/stephaneguindon/phyml/releases)
wget -c https://github.com/stephaneguindon/phyml/archive/refs/tags/v3.3.20220408.tar.gz
2. desempaqueta y ve al direcotorio
tar -xvzf v3.3.20220408.tar.gz
cd phyml-v3.3.20220408.tar.gz
3. Para compilarlo, ejecuta:
./autogen.sh
./configure
make
4. El binario aparecerá en el subdirecorio src
cd src
./phyml --help
5. Sugiero copiar phyml a un directorio en el $PATH
cp $HOME/bin
# o para instalarlo a nivel del sistema, con poderes de administrador
sudo cp phyml /usr/local/bin2.1 Publicaciones clave sobre PhyML
1: Criscuolo A. morePhyML: improving the phylogenetic tree space exploration with PhyML 3. Mol Phylogenet Evol. 2011 Dec;61(3):944-8. doi: 10.1016/j.ympev.2011.08.029. Epub 2011 Sep 8. PubMed PMID: 21925283.
2: Guindon S, Dufayard JF, Lefort V, Anisimova M, Hordijk W, Gascuel O. New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0. Syst Biol. 2010 May;59(3):307-21. doi: 10.1093/sysbio/syq010. Epub 2010 Mar 29. PubMed PMID: 20525638.
3: Guindon S, Delsuc F, Dufayard JF, Gascuel O. Estimating maximum likelihood phylogenies with PhyML. Methods Mol Biol. 2009;537:113-37. doi: 10.1007/978-1-59745-251-9_6. PubMed PMID: 19378142.
2.2 Opciones más usadas de PhyML v3
Notas:
- Phyml requiere alineamientos múltiples de DNA o proteína en formato phylip.
- Puedes usar el \(script\) convert_aln_format_batch_bp.pl para interconvertir alineamientos múltiples entre diferentes formatos, incluyendo FASTA, PHYLIP, NEXUS, entre muchos otros.
# -i input_file # -d nt|aa # -b -4 usa SH-like branch support values (-b 100 corre 100 repl de bootstrap) # -m PARA NT: (modelo de sustitucion HKY85 (default) | JC69 | K80 | F81 | F84 | TN93 | GTR o definido con modificadores de tasa 012345 (ac ag at cg ct gt) # -m PARA AA LG (default) WAG | JTT | MtREV | Dayhoff | DCMut | RtREV | CpREV | VT | Blosum62 | MtMam | MtArt | HIVw | HIVb | custom # -f frec. de nt [e empirica| m estimada por ML] # -t ts/tv # -v prop inv [0,1] o estimarla con e # -c no. clases para discretizar la funcion gamma # -a e (estima gamma) -a digito (fija su valor) # -s (define la busqueda o search) NNI SPR BEST (mejor de ambas) # --rand_start (inicia busquedas de arboles de adic. sec. aleatorizada) # --n_rand_starts no. (no. de arboles semilla) # -o params [tlr|tl|lr|l|r|n] permite optimizar combinaciones de t=topología, l=longitud de ramas, r parámetros de tasas (rates) de modelos; n=ninguno # -u usertree permite pasarle a phyml un árbol precalculado por el usuario
2.3 Ejemplos de invocación de PhyML v3 para secuencias de DNA (A) y proteína (B)
# A: Ejemplo de corrida para NT; recureda que los modificatodres de tasa indicados con -m se aplican en este orden: (ac ag at cg ct gt)
phyml -i primates.phy -m HKY85 -c 4 -a e -o tlr -s BEST --rand_start --n_rand_starts 5 &> /dev/null
# B: Ejemplo de corrida para AA
phyml -i 10_GDP_eucar_prokar_cluO_prof2prof.faa -d aa -m WAG -c 4 -a e -v e -s BEST &> /dev/null3 SELECCION DE MODELOS PARAMÉTRICOS DE DNA USANDO PHYML - Prácticas
Iniciaremos estas prácticas enfocándonos primeramente en el uso de \(phyml\), y en conceptos básicos de la selección de modelos paramétricos. Empezaremos por las pruebas de cocientes o razones de verosimilitudes (\(LRTs\)).
Una vez familiarizados con el uso de \(phyml\) y los conceptos básicos de contraste de hipótesis anidadas mediante \(LRTs\), nos enfocaremos a eficientizar el análisis computacional, introduciendo estructuras de programación básicas como bucles y estructuras de datos.
Esta progresión culmina con la presentación de los \(scripts\) phyml_DNAmodelFinder y phyml_protModelFinder, que automatizan el proceso de selección de modelos de manera eficiente, haciendo un riguroso y exhaustivo análisis de los modelos en competición basado en el criterio de información bayesiano \(BIC\).
3.1 Descarga de secuencias para la práctica
Para iniciar las prácticas, descarga las secuencias desde el repositorio GitHub del TIB22.
Lo haremos llamando a \(wget\) desde la línea de comandos.
Haz un \(cd\) al directorio de trabajo, y desde éste ejecuta el siguiente comando:
wget -c https://raw.githubusercontent.com/vinuesa/TIB-filoinfo/master/docs/sesion6_maxima_verosimilitud/phyml_tutorial_data.tgzy confirma que tienes el tarro comprimido (tgz):
-rw-rw-r-- 1 vinuesa vinuesa 145254 nov 21 19:36 phyml_tutorial_data.tgz
phyml_tutorial_data.tgz: gzip compressed data, last modified: Mon Jul 29 01:29:04 2019, from Unix, original size modulo 2^32 40960
el cual puedes desempacar con el siguiente comando:
3.2 Pruebas de razones de verosimilitud (\(LRTs\)) entre pares de modelos anidados
Las pruebas de razones o cocientes de verosimilitudes son una manera natural de evaluar parejas de modelos en competición, en los que comparamos el ajuste relativo de un modelo nulo (\(lnL_{H_0}\)) y uno alternativo (\(lnL_{H_1}\)), el cual relaja una o más constricciones en parámetros impuestos por \(H_0\). Es decir, compararemos modelos anidados, donde:
LRT=2*(\(lnL_{H_1}\) - \(lnL_{H_0}\))
El estadístico \(LRT\) se distribuye como una \(\chi^2\) con \(q\) grados de libertad, donde \(q\)=diferencia en parámetros libres a estimar entre \(H_1\) y \(H_0\).
Para acelerar el procedimiento de seleccion de modelos, típicamemente seguiremos los siguientes pasos:
- Calcular un árbol de distancias NJ, bajo el modelo más sencillo (\(JC\)), lo cual es muy rápido
- Este árbol lo usaremos para calcular los valores de los parámetros de los diferentes modelos a ajustar, y determinar la puntuación de verosimilitud de los modelos evaluados
Es decir, recordando la definición de la función de máxima verosimilitud:
lnL = Pr(D|H), que para modelos filogenéticos
podemos reescribir como:
lnL = Pr(D|\(\tau\nu\rho\)), donde:
- D=Datos (alineamiento)
- \(\tau\)=topología
- \(\nu\)=conjunto de longitudes de rama
- \(\rho\)= parámetros del modelo de sustitución
Veamos cómo hacerlo en el siguiente ejemplo.
3.2.1 Estima del árbol NJ de referencia, sobre el que estimaremos los parámetros de todos los modelos en competición
El objetivo de calcular inicialmente un un árbol NJ (método
algorítmico de agrupación muy rápido) es fijar la topología
(\(\tau\)), que es el
parámetro más costoso computacionalmente de calcular, pidiendo a \(phyml\) que optimice sólo los parámetros de
longitud de rama (\(\nu\)) y del modelo
de sustitución (\(\rho\)) a evaluar
(opción -o lr), calculando el puntaje (\(score\)) de verosimilitud
correspondiente.
# Visualizemos el alineamiento para confirmar que es tal, y que está en formato phylip (requerido por phyml)
head primates.phy && tail primates.phy
# Estima de árbol BioNJ con distancias JC69 (una sola tasa de sustitución e igual frecuencia de bases)
# phyml -i primates.phy -d nt -m 000000 -f 0.25,0.25,0.25,0.25 -b 0 -c 1 -o n &> /dev/null
phyml -i primates.phy -d nt -m JC69 -b 0 -c 1 -o n &> /dev/null
# Renombremos el archivo
mv primates.phy_phyml_tree.txt primates_JC.ph
# Podemos visualizar el arbol con:
figtree primates_JC.ph &3.2.2 Optimización de parámetros y estimas de \(lnL\) para una gama de modelos anidados a evaluar
Después de calcular el árbol NJ-JC, iniciamos el proceso de seleccion de modelos.
Existen 203 modelos en la familia \(GTR+I+G\), los cuales se distinguen por su grado de parametrización con respecto a:
- frecuencia de bases
- número de clases de tasas de sustitución
- modelado de la heterogeneidad de tasas de sustutución entre sitios
El objetivo de la selección de modelos es encontrar uno con un nivel óptimo de parámetros para nuestros datos. Es decir, un modelo con suficientes parámetros para modelar con realismo las características más salientes de nuestros datos (alineamiento), haciendo uso de tres clases de parámetros listados en el párrafo entarior.
La idea es agregar los parámetros necesarios para que el modelo presente un ajuste óptimo, pero sin incurrir en sobreparametrización.
Para entender el proceso, vamos a evaluar una serie de pares de modelos anidados de parametrización creciente, concretamente:
- JC vs K2P
- JC vs F81
- JC vs HKY
- HKY vs HKY+G
- HKY+G vs HKY+G+I
Recordemos cómo especificar la codificación de tasas de sustitución
relativas en phyml, haciendo uso de los siguientes códigos
numéricos:
tipo de sust. : Tv Ti Tv Tv Ti Tv
clase de tasas: 0 1 2 3 4 5
sustituciones: ac ag at cg ct gt
---------------------------------
A == G Ti (purinas)
|| X || Tv
C == T Ti (pirimidinas)
donde:
- 000000 ó 111111 representa una sola tasa de sustitución, como en los modelos JC69 o F81
- 121121 ó 010010 representa una clase de tasas para las transiciones (Ti)
y otra clase de tasas para las transversiones (Tv), como en K2$ o HKY85
- 012345 representa un modelo de 6 tasas, es decir, el GTR (General Time Reversible).
3.2.2.1 Jukes y Cantor de 1969 (\(JC\) o \(JC69\))
Comencemos la evaluacion de modelos. Empezaremos arbitrariamente por el modelo más sencillo, el JC.
- Parametrización de JC: igual frecuencia de bases (\(\pi_a = \pi_c = \pi_g = \pi_c = 0.25\) ) y tasa de Ti = Tv (\(\alpha = \beta\)).
- Este es el modelo nulo por excelencia, es decir, el que impone restricción total a los parámetros de tasa y frecuencia.
# 1. Modelo JC
phyml -i primates.phy -d nt -m JC69 -c 1 -b 0 -u primates_JC.ph -o lr &> /dev/null
# Anotemos el valor de Log-likelihood para este modelo
cat primates.phy_phyml_stats.txt...truncado . Sequence filename: primates.phy . Data set: #1 . Initial tree: user tree (primates_JC.ph) . Model of nucleotides substitution: JC69 . Number of taxa: 12 . Log-likelihood: -6424.20245 <=== ESTE ES EL PUNTAJE (lnL) QUE NECESITAMOS . Unconstrained log-likelihood: -4184.81795 . Composite log-likelihood: -14938.70804 . Parsimony: 1153 . Tree size: 1.43328 . Discrete gamma model: Yes - Number of classes: 1 - Gamma shape parameter: 1.000 - Relative rate in class 1: 1.00000 [freq=1.000000] . Nucleotides frequencies: - f(A)= 0.25000 - f(C)= 0.25000 - f(G)= 0.25000 - f(T)= 0.25000 ...truncado
Trabajando desde un \(shell\) de \(Unix\) o \(linux\), podemos acceder a este valor sin necesidad de abrir el archivo:
. Log-likelihood: -6424.20245
Guardaremos el valor de \(lnL\) para cada uno de los modelos a evaluar en un archivo adecuadamente nombrado para posterior referencia.
También podríamos haber usado \(awk\), que nos permite encontrar la línea que contiene ‘Log-like’, y además imprimir el último campo de la misma, que contiene el valor que buscamos.
-6424.20245
3.2.2.2 El modelo de Kimura de 1980 (\(K80\)) o Kimura de 2 parámetros (\(K2P\))
- Parametrizacaión de K80: asume igual frecuencia de bases (\(\pi_a = \pi_c = \pi_g = \pi_c = 0.25\) ), pero permite acomodar Ti \(\ne\) Tv (\(\alpha \ne \beta\)).
- Con este modelo evaluamos si existe un sesgo transicional en los datos
# 2. Modelo K80; veamos lo que tarda al usar el "user-tree" vs. tener que estimarlo
#phyml -i primates.phy -d nt -m 010010 -f 0.25,0.25,0.25,0.25 -t e -c 1 -u primates_JC.ph -o lr -b 0 &> /dev/null
phyml -i primates.phy -d nt -m K80 -t e -c 1 -u primates_JC.ph -o lr -b 0 &> /dev/null
grep Log- primates.phy_phyml_stats.txt > K80.fit3.2.2.3 Modelo de Felsenstein de 1981 (\(F81\))
- Parametrización: acomoda diferentes frecuencias de bases (\(\pi_a \ne \pi_c \ne \pi_g \ne \pi_c\)) pero asume que la tasa de Ti = Tv (\(\alpha = \beta\)).
- Con este modelo evaluamos si existe sesgo composicional,
evaluando frecuencias bajo máxima verosimilitud
-f m.
3.2.2.4 Modelo de Hasegawa, Kishino y Yano de 1985 (\(HKY85\))
- Parametrización: acomoda diferentes frecuencias de bases (\(\pi_a \ne \pi_c \ne \pi_g \ne \pi_c\)) y tasa de Ti \(\ne\) Tv (\(\alpha \ne \beta\)).
- Con este modelo evaluamos si existe sesgo composicional de bases y sesgo transicional en tasas de sustitución.
# 4. Modelo HKY85
# phyml -i primates.phy -d nt -m 010010 -f m -t e -c 1 -u primates_JC.ph -o lr -b 0 &> /dev/null
phyml -i primates.phy -d nt -m HKY85 -f m -t e -c 1 -o lr -b 0 &> /dev/null
grep Log-like primates.phy_phyml_stats.txt > HKY85.fitNOTA: Como veremos en ejemplos posteriores, podíamos haber optimizado
el proceso corriendo un bucle for como el mostrado abajo.
Por ahora nos concetramos en entender cómo llamar a \(phyml\) con modelos de parametrización
creciente.
3.2.2.5 Modelo de Hasegawa, Kishino y Yano de 1985 + distribución Gamma (\(HKY85+G\))
- Parametrización: (\(\pi_a \ne \pi_c \ne \pi_g \ne \pi_c\)), Ti \(\ne\) Tv (\(\alpha \ne \beta\)) y distribución Gamma (\(+\Gamma\)) para acomodar heterogeneidad de tasas entre sitios o columnas del alineamiento
- Con este modelo podemos evaluar si existe sesgo composicional de bases, sesgo transicional en tasas de sustitución y heterogeneidad de tasas entre sitios.
3.2.2.6 Modelo de Hasegawa, Kishino y Yano de 1985 + proporción de sitios invariantes (\(HKY85+I\))
- Parametrización: (\(\pi_a \ne \pi_c \ne \pi_g \ne \pi_c\)), Ti \(\ne\) Tv (\(\alpha \ne \beta\)) y proporción de sitios invariantes (\(+I\)) para acomodar heterogeneidad de tasas entre sitios o columnas del alineamiento
- Con este modelo podemos evaluar si existe sesgo composicional de bases, sesgo transicional en tasas de sustitución y heterogeneidad de tasas entre sitios, considerando sólo una proporción de sitios invariantes.
3.2.2.7 Modelo de Hasegawa, Kishino y Yano de 1985 + distribución Gamma + proporción de sitios invariantes (\(HKY85+G+I\))
- Parametrización: (\(\pi_a \ne \pi_c \ne \pi_g \ne \pi_c\)), Ti \(\ne\) Tv (\(\alpha \ne \beta\)), distribución Gamma (\(+\Gamma\)) y proporción de sitios invariantes (\(pInv\)) para modelar heterogeneidad de tasas entre sitios
- Con este modelo podemos evaluar si existe sesgo composicional de bases, sesgo transicional en tasas de sustitución, heterogeneidad de tasas entre sitios, asumiendo además una proporción de sitios invariantes.
# 7. Modelo HKY85+G+I
# phyml -i primates.phy -d nt -m 121121 -f m -t e -c 4 -a e -v e -u primates_JC.ph -o lr -b 0 &> /dev/null
phyml -i primates.phy -d nt -m HKY85 -f m -t e -c 4 -a e -v e -u primates_JC.ph -o lr -b 0 &> /dev/null
grep Log-like primates.phy_phyml_stats.txt > HKY85GI.fitNota: en el ejercicio de selección de matrices de sustitución para proteínas veremos cómo optimizar este tedioso proceso, haciendo uso de un \(bucle\ for\) y un \(hash\) para guardar los valores de \(lnL\) directamente en una estructura de datos, sin escribir archivos intermedios, como hicimos en los ejemplos anteriores.
3.2.3 Análisis gráfico y estadístico de resultados usando herramientas de \(Linux\) y \(R\)
Una vez evaluados los diferentes modelos en competición, haremos un simple análisis gráfico y estadístico de los resultados.
- Parseo de valores \(lnL\) de los archivos *.fit, y generación de archivos tsv y csv
# generemos un archivo separado por tabuladores con su cabecera
awk 'BEGIN{OFS="\t"; print "model\tlnL"}{print FILENAME , $NF}' *.fit | sed 's/\.fit//' > model_fits.tsv
# limpiemos un poco el directorio de trabajo
rm *fit Visualizemos los resultados obtenidos, ordenados de mayor a menor
\(lnL\) con sort -grk2
# sort -g, --general-numeric-sort: compare according to general numerical value
sed '1d' model_fits.tsv | sort -grk2 | column -tHKY85G -5711.93899 <== mejor (mayor) lnL HKY85GI -5711.94038 HKY85I -5768.00412 HKY85 -5981.72029 K80 -6142.42909 F81 -6284.99565 JC69 -6424.20245
Los resultados arriba mostrados indican que el mejor modelo es \(HKY+G\). Noten que en este caso existe una interacción negativa entre los parámetros \(\Gamma\) e \(I\) en el modelo \(HKY+G+I\), ya que este último tiene un puntaje de \(lnL\) ligeramente más bajo que el del modelo \(HKY+G\), a pesar de contar con un parámetro más.
Esta interacción negativa entre \(\Gamma\) y \(I\) la he observado con cierta frecuencia, por lo que suelo ignorar \(I\) para modelar la heterogeneidad de tasas entre sitios, a no ser que la evidencia en favor de usar \(I\) sea contundente.
3.2.3.1 Graficado con R de valores de verosimilitud de los modelos en competición, ordenados de \(lnL\)
# 6. Evaluacion del merito relativo de los modelos HKY85+G vs HKY85+G+I
# mediante la prueba de cocientes de verosimilitud (LRT), en R
# carga R
R
# leamos los datos en una estructura tipo dataframe
# system("ls")
list.files(pattern="\\.tsv$")
dfr <- read.csv(file="model_fits.tsv", header = TRUE, sep="\t")
# ver la estructura del data frame
str(dfr)
# imprimelo a pantalla
head(dfr)
# 6.2 veamos una grafica tipo "dotchart". Que puedes concluir?
dfrord <- dfr[order(dfr$lnL),]
dfrord
dotchart(dfrord$lnL, labels = dfrord$model, main = "sorted model scores", xlab="lnL score")
# guardemos salida como archivo png
png(file="LRTs_primates_dotchart.png")
dotchart(dfrord$lnL, labels = dfrord$model, main = "sorted model scores", xlab="lnL score")
dev.off()- visualizemos la gráfica
Figura 1. Dotchart mostrando resultados ordenados de puntuación de \(lnL\) para los modelos indicados, estimados sobre una filogenia \(NJ-JC\), y el alineamiento primates.phy.
La gráfica revela claramente la importancia los parámetros de frecuencia de bases, ti/tv y heterogeneidad de tasas entre sitios modelada con una distribución gamma (\(\Gamma\)), reflejados en los incrementos altamente significativos de \(lnL\) de los modelos de complejidad creciente analizados.
3.2.3.2 Evaluación formal de ajuste relativo entre pares de modelos anidados con grado creciente de parametrización, mediante pruebas de razones de verosimilitud (\(LRTs\))
Examinando la figura 1, vemos que los modelos \(JC69\) y \(F81\) tienen puntajes de verosimilitud muy distintos. La pregunta es: ¿Son significativamente distintos?. Vamos a hacer una prueba de cocientes de verosimilitudes o \(LRT\) para evaluar estadísticamente si existe un sesgo significativo en frecuencia de baes (\(-f\ e\)), que justifique estimar tres parámetros libres adicionales (por \(F81\)).
Recordemos que:
- Las pruebas de razones de verosimilitud (\(LRTs\)) sólo permiten comparar
parejas de modelos anidados, en los que el modelo nulo
(\(L_{H_0}\)) está anidado dentro del
alternativo (\(L_{H_1}\)), es
decir, \(L_0\) es un caso particular de
\(L_1\).
- El estadístico \(LRT=2*(lnL_1 - lnL_0)\) sigue una distribución \(\chi^2\) (chi-cuadrada):
\[\chi^2 = \sum \frac {(O - E)^2}{E}\]
Para ganar una apreciación de las magnitudes de LRTs que son altamente significativas, y las que no, hagamos los siguientes contrastes de hipótesis:
- JC vs. F81
# Dado que los modelos JC69 y F81 estan anidados,
# podemos correr un simple LRT (prueba de razones de verosimilitudes): LRT = 2(lnL1 - lnL 0)
dfr
# LRT = 2( lnL1 - lnL0 )
# 2 veces la diferencia de lnL entre el alternativo (F81) - modelo nulo (JC69)
xsq <- 2*(-6284.99565 - -6424.20245)
xsq # [1] 278.4136 # Este es un valor gigantesco de la ji-cuadrada, y por tanto esperamos que sea ALTAMENTE SIGNIFICATIVO.
# cómputo del valor crítico de la distibución ji-cuadrada para p=0.05 y p=0.01, con 3 grados de libertad:
qchisq(0.95,3) # 7.814728
qchisq(0.99,3) # 11.34487
# cómputo de estadísticas
pchisq(xsq,3,lower.tail=FALSE) # 4.66758e-60; esta es la probabilidad en la distribución de ji-cuadrada a la derecha de xsq (del LRT) con 3 grado2 de libertad;
# como es << 0.05, es un valor altamente significativo, e implica que podemos rechazar la hipótesis o modelo nulo!
1-(pchisq(xsq,1)) # 0; esto es equivalente a lo escrito arriba, es una p ABSOLUTAMENTE SIGNIFICATIVA => rechazamos el modelo nulo
q() # para salir de R
# decimos que NO quieren salvar sesionEn general, cuando observamos una diferencia de \(lnL\) de decenas de unidades entre modelos anidados, la diferencia será significativa, como podrás comprobar en el siguiente ejercicio:
- Tarea: repite el ejercicio de evaluación de
cocientes de similitud para las siguientes parejas de modelos anidados:
- K80 vs. HKY85 (¿qué evaluamos con este contraste de hipótesis, y los dos que siguen?)
- HKY85+I vs. KY85+I+G
3.2.4 Estima de la filogenia bajo el mejor modelo encontrado: HKY85+G
Estimemos la filogenia de máxima verosimilitud bajo el modelo seleccionado: \(HKY+G\)
Nótese que usamos sólo la parameterización del modelo
seleccionado (-m HKY85 -c 4 -a e), pero no los
varlores de los parámetros, ya que fueron optimizados sobre el árbol
\(NJ-JC\).
Pediremos a phyml que optimice todos los parámetros,
incluyendo ahora la topología, longitudes de ramas y
parámetros del modelo de sustitución con la opción -o
tlr.
- Pregunta: ¿qué tipo de rearreglos de rama usa \(phyml\) en la llamada mostrada seguidamente?
Una vez terminada la búsqueda, veamos las estimas de máxima verosimilitud de los parámetros del modelo de sustitución HKG+G.
# veamos las estimas de máxima verosimilitud de los parámetros del modelo de sustitución HKG+G
cat primates.phy_phyml_stats.txt...truncado . Sequence filename: primates.phy . Data set: #1 . Initial tree: BioNJ . Model of nucleotides substitution: HKY85 . Number of taxa: 12 . Log-likelihood: -5711.93881 . Unconstrained log-likelihood: -4184.81795 . Composite log-likelihood: -14275.02293 . Parsimony: 1153 . Tree size: 3.38341 . Discrete gamma model: Yes - Number of classes: 4 - Gamma shape parameter: 0.363 - Relative rate in class 1: 0.01179 [freq=0.250000] - Relative rate in class 2: 0.15381 [freq=0.250000] - Relative rate in class 3: 0.68822 [freq=0.250000] - Relative rate in class 4: 3.14618 [freq=0.250000] . Transition/transversion ratio: 12.726587 . Nucleotides frequencies: - f(A)= 0.36372 - f(C)= 0.31987 - f(G)= 0.08153 - f(T)= 0.23488 ...truncado
Compara estos parámetros con los obtenidos durante el proceso de selección de modelos usando el árbol \(NJ-JC\). ¿Notan diferencias?
Los resultados de estimas de máxima verosimilitud de los parámetros del modelo \(HKY+G\) arriba mostrados (del archivo primates.phy_phyml_stats.txt) son fáciles de entender, pero:
- Pregunta: ¿qué significan los valores “Relative rate in class 1 - 4”?
3.2.5 Despliegue y edición del árbol con figtree
Podemos desplegar y editar los árboles con figtree del Prof. Andrew Rambaut.
Figura 2. Mejor árbol de máxima verosimilitud encontrado para primates.phy bajo el modelo \(HKY+G\), seleccionado mediante LRTs
3.3 Estrategia alternativa de selección de modelos de DNA: reducción de parámetros del modelo GTR+G+I, y evaluación de los modelos resultantes mediante \(AIC\)
En el ejercicio anterior evaluamos una serie de modelos anidados, iniciando arbitrariamente por el más sencillo (JC69), progresando hasta el modelo \(HKY85+I+G\).
Pero podríamos haber construido y analizado otras jerarquías de modelos anidados, lo cual, posiblemente, nos hubiera llevado a seleccionar otros modelos. Esta es una limitación de las pruebas de cocientes de verosimilitudes.
Una alternativa es ajustar el modelo más rico en parámetros \(GTR+G+I\) (10 parámetros libres: 3 de frecuencias de bases + 5 de tasas + 1 de gamma + 1 de propInv) a los datos, y examinar la matriz de transición resultante para simplificarla (reducir parámetros) donde parezca obvio, por ser tasas muy similares (< ~20-30 % de diferencia), y evaluar los ajustes de estos modelos usando criterios de información.
3.3.1 Ajuste del modelo \(GTR+G+I\) a primates.phy
Igual que en ejemplos anteriores, usamos el árbol NJ-JC para estimar sobre éste los parámetros del modelo de sustitución y evaluar el ajuste (\(lnL\)) resultante, sin necesidad de invertir tiempo de cómputo en reestimar la topología y longitudes de rama.
Con \(grep\) filtramos la línea con el ‘Log-likelihood’ y la guardamos en un archivo adecuadamente nombrado.
3.3.2 Exploración del archivo stats para identificar posibles simplificaciones al modelo GTR+G+I
Analicemos los valores obtenidos para las tres clases de parámetros de modelos para secuencias de DNA: frecuencias, tasas, y de heterogeneidad de tasas entre sitios
# 1. Analiza las frecuencias del modelo GTR+G+I
grep '\- f(' primates.phy_phyml_stats.txt
# 2. Analiza las tasas del modelo GTR+G+I
grep '<->' primates.phy_phyml_stats.txt
# 3. Analiza Prop. invariantes y gamma
grep 'Gamma shape parameter' primates.phy_phyml_stats.txt
grep 'Proportion of invariant' primates.phy_phyml_stats.txt# Frecuencias de nucleótidos: - f(A)= 0.35497 - f(C)= 0.32231 - f(G)= 0.08008 - f(T)= 0.24264 # GTRIG 2 3 3 4 rates << clases de tasas para modelos simplificados de 5 y 4 tasas A <-> C 3.60550 0 0 0 0 A <-> G 43.95008 1 1 1 1 A <-> T 3.05954 0 0 0 0 C <-> G 2.37102 0 0 0 2 C <-> T 34.89961 1 2 1 1 G <-> T 1.00000 0 0 2 3 - Gamma shape parameter: 0.369 . Proportion of invariant: 0.000
En base a los resultados arriba mostrados, concluimos que:
- Las frecuencias de bases son muy heterogéneas, por lo que usaremos modelos que las acomoden
- Las Ti >> Tv, es decir, existe un muy marcado sesgo transicional, que debe ser contemplado por el modelo
- La proporción de sitios invariantes es despreciable si usamos la distribución \(\Gamma\) para modelar la heterogeneidad de tasas entre sitios
En resumen, usaremos modelos que acomodan diferentes frecuencias de bases, separan al menos tasas entre transiciones y transversiones, y usan sólo la distribución \(\Gamma\) para modelar la heterogeneidad de tasas entre sitios.
En la siguiente sección evaluaremos algunas opciones de simplificación de número de tasas (2 ó 3).
3.3.3 Evaluación de modelos simplificados de 2, 3 y 4 tasas
La simplificación más obvia es igualar ambas tasas de transiciones (010010), ya que son parecidas entre sí y alrededor de un órden de magnitud más altas que las de transversiones.
Evaluaremos también si las tasas de trasiciones entre purinas y entre pirimididnas son significativamente diferentes (010020) entre sí, lo que equivale a usar un modelo TN93 (Tamura-Nei 93)
Habría otras opciones de simplificación, incluyendo modelos de 4 tasas, como se indica en las columnas numéricas mostradas arriba, junto a las tasas del \(GTR+I+G\)
Copia el siguiente bloque completo en la terminal para ejecutar todos los comandos secuencialmente.
# Evaluación del modelo alternativo (simplificado) de 2 tasas (010010) + G # NOTA: es el HKY85+G
phyml -i primates.phy -m HKY85 -f m -a e -u primates_JC.ph -b 0 -o lr &> /dev/null
grep Log- primates.phy_phyml_stats.txt > model_HKY85G.fit
# Evaluación del modelo alternativo (simplificado) de 3 tasas (010020) + G # NOTA: es el TN93+G
phyml -i primates.phy -m TN93 -f m -a e -u primates_JC.ph -b 0 -o lr &> /dev/null
grep Log- primates.phy_phyml_stats.txt > model_TN93G.fit
# Visualicemos los ajustes, ordenados por valores decrecientes de lnL:
for f in *fit; do echo -ne "$f\t"; cat $f; done | awk '{print $1,$4}' \
| sed 's/\.fit//' | sort -k2 | column -tGTRGI -5708.92372 TN93G -5710.55182 HKYG -5711.93899
Nota que si bien el modelo de más alta parametrización (\(GTR+G+I\)) tiene el puntaje de verosimilitud (\(lnL\)) más alto, la diferencia contra el competidor \(HKY+G\) es pequeña, de sólo ~3 unidades de \(lnL\).
Usemos el criterio de información de Akaike (AIC) para evaluar estos resultados formalmente.
3.3.4 Cálculo de los
\(AICi\) y \(\Delta_i\) (de cada modelo) con respecto al
\(AICmin\) en R
Recordemos la fórmula del criterio de información de Akaike para un modelo particular: \(AIC_i=-2*(lnL_i)+(2N_i)\)
donde:
- \(lnL_i=log-Likelihood\ del\ modelo\ i\)
- \(N_i=número\ de\ parámetros\ libres\ del\ modelo\ i\),
Noten que \(N_i\) debería considerar además el número de ramas libres a estimar (\(2n-3\); \(n = num\_secuencias\)), por lo que:
- \(N_i=número\ de\ parámetros\ libres\ del\ modelo\ i + (2n - 3)\)
R
# cálculo del número de ramas libres, donde n=12 taxa
n <- 12 # número de secuencias o nodos terminales
nbranches <- 2*(n-3) # número de ramas libres a estimar para n secuencias
# construyamos el data frame model_AICs, con los vectores "model" y "AIC"
# vector model con los nombres de los modelos
model <- c("AIC_GTRGI", "AIC_TN93G", "AIC_HKYG")
# Parámetros libres del GTR+G+I = 10 (frecuencias=3, tasas=5, pInv=1, Gamma=1)
# Parámetros libres del TN+G = 6 (frecuencias=3, tasas=2, Gamma=1)
# Parámetros libres del HKY+G = 6 (frecuencias=3, tasas=1, Gamma=1)
AIC_GTRGI <- -2*(-5708.92372) + 2*(10+nbranches)
AIC_TN93G <- -2*(-5710.55182) + 2*(6+nbranches)
AIC_HKYG <- -2*(-5711.93898) + 2*(5+nbranches)
# vector AIC, con los AICi de cada modelo, en el mismo orden.
AIC <- c(AIC_GTRGI, AIC_TN93G, AIC_HKYG)
# el dataframe
model_AICs <- data.frame(model, AIC)
model_AICs
# mejor ordenado por valor de AIC
model_AICs[order(model_AICs$AIC),]
# y guardamos en AIC_min
AICmin <- min(model_AICs$AIC)
AICminExplorando el dataframe \(model\_AICs\) ordenado por valores crecientes de \(AICi\), vemos que el \(AICmin\) corresponde al modelo AIC_TN93G
model_AICs[order(model_AICs$AIC),]
model AIC
2 AIC_TN93G 11469.10
3 AIC_HKYG 11469.88
1 AIC_GTRGI 11473.85
Construyamos un nuevo \(data\ fame\) \(model\_deltas\) con los vectores \(models\) y \(deltas\), con los valores de \(\Delta_i\) para cada modelo, con respecto al \(AIC\_min\)
# compute deltas
delta_GTRGI <- AIC_GTRGI - AICmin
delta_HKYG <- AIC_HKYG - AICmin
delta_TN93G <- AIC_TN93G - AICmin
# build vectors
deltas <- c(delta_GTRGI, delta_HKYG, delta_TN93G)
models <- c("GTRGI","HKYG", "TN93G")
# build and print model_deltas data frame
model_deltas <- data.frame(models, deltas)
model_deltas[order(model_deltas$deltas),]models deltas 3 TN93G 0.00000 2 HKYG 0.77432 1 GTRGI 4.74380
3.3.5 Conclusiones
- El mejor modelo (AICmin) evaluado es TN93G (0100020), pero el modelo HKYG (0100010) es un serio competidor, ya que está a una distancia \(\Delta_i\ =\ 0.774\).
- El parámetro de \(pInv\) es perfectamente prescindible, y un modelo de 2 ó 3 tasas parece suficiente para obtener un buen ajuste.
- Queda claro que, si bien el modelo \(GTR+G+I\) presenta el valor más alto de \(lnL\), al ponderar los parámetros libres a estimar, resulta ser el peor modelo de los evaluados por AIC (\(\Delta_i\ =\ 4.744\)), indicando que está claramente sobre-parametrizado.
En mi experiencia, esta estrategia es buena por dos motivos:
- es muy rápida (posiblemente la más rápida)
- encuentra muy buenos modelos, ya que los métodos tradicionales,
incluyendo estrategias automáticas, como las implementadas en \(jModelTest2\), no evalúan por defecto a los
203 modelos posibles de la familia \(GTR+G+I\) (aunque es posible con la opción
-s 203).
Para una implementación eficiente de la estrategia arriba descrita, sería deseable escribir un \(script\) para automatizarla.
3.3.6 Tarea de ajuste de modelos paramétricos mediante criterios de información AIC y BIC
- Evalúa mediante AICs las otras simplificaciones en las clases de transversiones propuestas en el esquema de tasas del modelo \(GTR+I+G\).
- Compara los puntajes de \(lnL\) de todos estos modelos entre sí mediante análisis de \(\Delta_i\), incluyendo los modelos \(HKY+G\) y \(TN93G\) evaluados previamente.
- Como veremos en la siguiente sección, el criterio de información de Akaike \(AIC\) tiende a seleccionar modelos más parametrizados que el criterio de información bayesiano \(BIC\) que presentaremos en la siguiente sección. Repite el análisis usando \(BIC\).
4 Cálculo de \(AICs\) y \(BICs\) para selección modelos de aminoácidos
Recuerda que para la inferencia de filogenias de proteínas bajo máxima verosimilitud se usan matrices empíricas de sustitución, como la \(BLOSUM62\) o \(JTT\) … Dado que estas matrices NO están anidadas entre sí, necesitamos usar criterios de información como el de Akaike (\(AIC\)) o bayesiano (\(BIC\)) para evaluar su mérito o ajuste relativo.
En el siguiente ejercicio evaluaremos diversas matrices empíricas, ajustando o no frecuencia de los residuos (\(-f\ e\)), y/o distribución gamma (\(+G\)) para modelado de la heterogeneidad de tasas entre los sitios o columnas del alineamiento.
Como hicimos en ejemplos anteriores, para minimizar el tiempo de
cómputo, estimaremos un un árbol NJ-LG, el cual pasaremos a \(phyml\) como filogenia pre-definida por el
usuario (-u usertree).
Copia todo el bloque en la terminal, a partir de
infile=primates_21_AA.phy, para ejecutar todas las líneas
secuencialmente.
4.1 Evaluación rápida de un subconjunto de matrices empíricas con un \(bucle\ for\) y un \(hash\)
En este ejercicio haremos uso más eficiente de las herramientas que nos brindan el \(Shell\) y \(linux\), como primeros pasos hacia la automatización del proceso.
# Evaluemos un subconjunto de los modelos empiricos (matrices de sustitución) implementados en phyml para proteínas
# - Amino-acid based models : LG (default) | WAG | JTT | MtREV | Dayhoff | DCMut | RtREV | CpREV | VT | AB
# Blosum62 | MtMam | MtArt | HIVw | HIVb | custom
# primero seleccionemos la mejor matriz de base y además ajustando con ajuste empírico de frecuencias y
# una distribución gamma para modelar heterogeneidad de tasas entre sitios
# NOTA: copia todo el bloque, hasta pasado el done, y pégalo en tu consola
# guardemos el archivo de entrada en una variable
infile=primates_21_AA.phy
# primero estimemos un árbol de distancias NJ-LG, de cómputo rápido
phyml -i "$infile" -d aa -m LG -c 1 -b 0 -o n &> /dev/null
# y renombrémoslo, para su uso como ´usertree'
mv "${infile}"_phyml_tree.txt primates_LG.ph
# verifica que se generó el archivo;
# si no fuera el caso, termina el proceso con mensaje explícito de error
[[ ! -s primates_LG.ph ]] && { echo "ERROR: primates_LG.ph not found"; exit 1 ; }
# declaramos el hash model_fits para guardar los valores de lnL asociados a cada modelo
declare -A model_scores
# Corramos las matrices derivadas de genes nucleares en un bucle for para un procesado más eficiente
# A las matrices de base le agregamos sucesivamente parámetros adicionales: -a e; -f e; -f m -a e
# Noten el uso de awk para procesar directamente los valores de lnL de los archivos *_phyml_stats.txt,
# los cuales se guardan en el hash model_scores, usando como llave el nombre del modelo :)
for mat in AB BLOSUM62 DAYHOFF JTT LG VT WAG; do
echo " >>> running: phyml -i ${infile} -d aa -m $mat -u primates_LG.ph -c 1 -v 0 -o lr &> /dev/null"
phyml -i ${infile} -d aa -m $mat -f m -u primates_LG.ph -c 1 -o lr &> /dev/null
model_scores["$mat"]=$(awk '/Log-like/{print $NF}' ${infile}_phyml_stats.txt)
echo " >>> running: phyml -i ${infile} -d aa -m $mat -c 4 -a e -u primates_LG.ph -o lr &> /dev/null"
phyml -i ${infile} -d aa -m $mat -f m -c 4 -a e -u primates_LG.ph -o lr &> /dev/null
model_scores["${mat}+G"]=$(awk '/Log-like/{print $NF}' ${infile}_phyml_stats.txt)
# noten el uso de '-f e' para estimar las frequencias empíricas de los aminoácidos del alineamiento
echo " >>> running: phyml -i ${infile} -d aa -m $mat -f m -c 1 -u primates_LG.ph -o lr &> /dev/null"
phyml -i ${infile} -d aa -m $mat -f e -c 1 -u primates_LG.ph -o lr &> /dev/null
model_scores["${mat}+f"]=$(awk '/Log-like/{print $NF}' ${infile}_phyml_stats.txt)
echo " >>> running: phyml -i ${infile} -d aa -m $mat -f e -a e -u primates_LG.ph -o lr &> /dev/null"
phyml -i ${infile} -d aa -m $mat -f e -c 4 -a e -u primates_LG.ph -o lr &> /dev/null
model_scores["${mat}+f+G"]=$(awk '/Log-like/{print $NF}' ${infile}_phyml_stats.txt)
done4.2 Análisis gráfico y numérico de resultados de \(AIC\) usando herramientas de filtrado de \(Linux\) y \(R\)
Primero escribiremos en un archivo los valores de \(lnL\) para cada modelo evaluado guardados
en nuestro hash model_scores, el cual leeremos en \(R\) para generar las gráficas de
puntuaciones de \(lnL\), y calcular los
\(AIC_i\) de cada modelo.
# vamos a recorrer el hash model_scores para imprimir los resultados en forma tabular
# y ordenada por valores decrecientes de lnL
echo "# writing ${infile}_sorted_model_fits.tsv, sorted by lnL"
for m in "${!model_scores[@]}"; do
echo -e "$m\t${model_scores[$m]}"
done | sort -bk2 > "${infile}"_sorted_model_fits.tsv
# si se generó el archivo, y no está vacío, lo imprimimos a STDOUT
[[ -s "${infile}"_sorted_model_fits.tsv ]] && cat "${infile}"_sorted_model_fits.tsv
# es trivial generar también un archivo con campos delimitados por comas (csv)
sed 's/\t/,/' "${infile}"_sorted_model_fits.tsv > "${infile}"_sorted_model_fits.csv
# visualizemos la tabla
column -t -s $',' "${infile}"_sorted_model_fits.csvJTT+f+G -4778.48904 AB+f+G -4814.44420 JTT+G -4815.80148 JTT+f -4865.61699 WAG+f+G -4870.83180 VT+f+G -4880.25893 JTT -4904.05775 AB+G -4910.17657 WAG+G -4910.24192 LG+f+G -4913.60180 VT+G -4916.81421 DAYHOFF+f+G -4920.31326 AB+f -4933.84484 LG+G -4946.47152 WAG+f -4961.23284 VT+f -4968.11339 DAYHOFF+G -4969.87622 BLOSUM62+f+G -4976.30876 WAG -4998.94774 BLOSUM62+G -5002.18584 VT -5004.20867 LG+f -5008.95318 DAYHOFF+f -5021.06263 AB -5040.59363 LG -5042.35342 DAYHOFF -5069.90036 BLOSUM62+f -5071.35214 BLOSUM62 -5097.14921
- ¿Qué podemos concluir de esta salida con respecto a las matrices y parámetros adicionales usados?
4.2.1 Análisis gráfico con R de modelos ordenados por valor de \(lnL\)
# carga R
R
# 6.1 leamos los datos en una estructura tipo dataframe
list.files(pattern='\\.csv')
dfr <- read.csv(file="primates_21_AA.phy_sorted_model_fits.csv", header = FALSE)
# ver la estructura del data frame
str(dfr)
# imprimelo a pantalla
dfr
# aniadimos un header (nombramos las columnas con las variables)
names(dfr) <- c("model", "lnL")
head(dfr)
# veamos un dotchart()
dfrord <- dfr[order(dfr$lnL),]
dotchart(dfrord$lnL, labels = dfrord$model, main = "sorted model scores", xlab="lnL score")
# guardémoslo como pdf
pdf(file="dotchart_sorted_model_scores_for_primates_21_AA.pdf")
dotchart(dfrord$lnL, labels = dfrord$model, main = "sorted model scores", xlab="lnL score")
dev.off()
# guardémoslo como png
png(file="dotchart_sorted_model_scores_for_primates_21_AA.png")
dotchart(dfrord$lnL, labels = dfrord$model, main = "sorted model scores", xlab="lnL score")
dev.off()- Veamos la gráfica resultante de modelos ordenados por sus valores de \(lnL\)
Figura 4. Puntaciones de máxima verosimilitud para los modelos evaluados usando un árbol \(NJ-LG\) con el alineamiento \(primates\_21\_AA.phy\)
# veamos los resultados ordenados por valor decreciente de lnL
dfr[order(dfr$lnL, dfr$model, decreasing = TRUE), ]
# ... truncado
model lnL
14 JTT+f+G -4778.489
2 AB+f+G -4814.444
16 JTT+G -4815.801
13 JTT+f -4865.617
26 WAG+f+G -4870.832
22 VT+f+G -4880.259
15 JTT -4904.058
4 AB+G -4910.177
28 WAG+G -4910.242
# ... truncado
- Pregungas
- ¿Qué podemos concluir de estos resultados?
- ¿Usarían \(JTT+f+G\) como el mejor modelo?
4.2.2 Cálculo formal de \(AIC_{min}\) y \(\Delta_{AIC_{i}}\) para los modelos en competición
# Computo de AICmin y deltaAICi
# recuerda: AICi = 2Ki - 2(lnLi); donde Ki = numero parametros libres del modelo i; lnLi = lnL del modelo i
# número de ramas libres a estimar = 2(n-3)
# G.L:
# - 19 para los 20 AAs (estimo 19 frecuencias y la 20ava es 1-(sum(19 freqs estimadas)))
# - 1 para el parametro alpha de la distriucion gamma
# - 39 ramas libres a estimar (2n-3)
AIC_JTTfg <- 2*(19+1+39) -2*(-4778.48904)
AIC_JTTG <- 2*(1+39) -2*(-4815.80148)
AIC_JTTf <- 2*(19+39) -2*(-4865.61699)
cat("AIC_JTTfg: ", AIC_JTTfg, "\n", "AIC_JTTG: ", AIC_JTTG, "\n", "AIC_JTTf: ", AIC_JTTf, "\n")AIC_JTTfg: 9674.978 <<< AIC_min AIC_JTTG: 9711.603 AIC_JTTf: 9847.234
# compute and print the deltaAICs for the 2 top models
AICmin <- AIC_JTTfg
deltaJTTG <- AIC_JTTG - AICmin
deltaJTTf <- AIC_JTTf - AICmin
cat ("deltaAIC", "\t", "value", "\n", "deltaJTTG", "\t", deltaJTTG, "\n", "deltaJTTf", "\t", deltaJTTf, "\n")# deltaAIC value # deltaJTTG 36.62 # dista >> 10 unidades de AIC del AICmin, por tanto, son modelos que podemos descartar # deltaJTTf 172.26
4.2.3 CONCLUSIÓN
El análisis de \(AIC\) claramente muestra que por mucho el mejor modelo de los evaluados es JTT+f+G.
Noten que aunque con \(-f\ e\) se tuvieron que estimar 19 parámetros libres, el incremento del ajuste fue mayor que la penalización consiguiente \(2Ni=38+39\).
4.3 Búsqueda intensiva bajo el mejor modelo, usando rearreglos tipo BEST y 5 árboles de semilla
# finalmente, hagamos una busqueda intensiva bajo este modelo, usando el algoritmo BEST de rearreglo de ramas (evalúa NNI y SPR en cada iteración),
# con 5 arboles de semilla y -o tlr
phyml -i primates_21_AA.phy -d aa -m JTT -c 4 -a e -f m -s BEST --rand_start --n_rand_starts 54.3.1 Despliegue y edición del árbol con figtree
Figura 5. Mejor árbol de máxima verosimilitud encontrado para primates_21_AA.phy bajo el modelo \(JTT+G+f\)
- Veamos los parámetros relevantes
...truncado . Sequence filename: primates_21_AA.phy . Data set: #1 . Initial tree: BioNJ . Model of amino acids substitution: JTT . Number of taxa: 21 . Log-likelihood: -4775.99427 . Unconstrained log-likelihood: -2434.66382 . Composite log-likelihood: -30488.21256 . Parsimony: 560 . Tree size: 1.29999 . Discrete gamma model: Yes - Number of classes: 4 - Gamma shape parameter: 0.718 - Relative rate in class 1: 0.07800 [freq=0.250000] - Relative rate in class 2: 0.37181 [freq=0.250000] - Relative rate in class 3: 0.93129 [freq=0.250000] - Relative rate in class 4: 2.61890 [freq=0.250000] . Run ID: none . Random seed: 1659906799 . Subtree patterns aliasing: no . Version: 3.3.3:3.3.20190909-1 . Time used: 0h0m14s (14 seconds) ...truncado
5 Selección automatizada de modelos de DNA y proteína con phyml_DNAmodelFinder.sh y phyml_protModelFinder.sh
Si bien en el ejemplo anterior hemos optimizado significativamente el análisis corriendo un \(bucle\ for\) y guardando los valores de \(lnL\) para cada modelo en el \(hash\) \(model\_scores\), este trabajo sigue siendo tedioso, con mucha intervención manual, lo cual, además de lento, nos expone a cometer errores.
No obstante, los pasos seguidos son fácilmente automatizables mediante un \(script\) en \(Bash\) u otro lenguaje interpretado.
Además del \(bucle\ for\) y el \(hash\), vamos a poder optimizar el proceso si escribimos \(funciones\) que encapsulen el trabajo de las diferentes operaciones que tenemos que hacer de manera repetida.
Como discutimos en el tema de funciones, éstas son fundamentales para modularizar el código en unidades lógicas autocontenidas, lo cual es esencial para escribir programas complejos y poderlos mantener y desarrollar con mayor facilidad.
Vamos a cerrar este tutorial corriendo los \(scripts\)
phyml_protModelFinder.sh y
phyml_DNAmodelFinder.sh, que puedes descargar a tu
directorio ~/bin desde el repositorio GitHub de los
TIB, usando una llamada a \(wget\)
con el siguiente código:
# si no existe el directorio ~/bin, entonces genéralo
[ ! -d ~/bin ] && mkdir ~/bin
# cámbiate al directorio ~/bin y llama a wget para decargar los scripts
cd ~/bin
wget -c https://raw.githubusercontent.com/vinuesa/TIB-filoinfo/master/phyml_protModelFinder.sh
wget -c https://raw.githubusercontent.com/vinuesa/TIB-filoinfo/master/phyml_DNAmodelFinder.sh
# haz ejecutables los scripts de shell (*.sh) descargados
chmod 755 *sh
# regresa al directorio previo (de trabajo)
cd -Estos \(scripts\), además de hacer las llamadas a \(phyml\), evalúan los resultados haciendo uso de los criterios de información de Akaike (\(AIC\)) y bayesiano (\(BIC\)).
Veremos que el \(AIC\) usa un criterio más laxo que el \(BIC\), tendiendo a seleccionar modelos más complejos que el segundo.
Si te interesa la programación en Bash, te invito a que explores en detalle el código de ambos \(scripts\).
NOTA: Los scripts phyml_protModelFinder.sh, y en particular phyml_DNAmodelFinder.sh, hacen uso de sintaxis moderna de Bash, por lo que requieren una versión reciente (>= 5.1.0) del intérprete de comandos \(bash\) para poder hacer uso de toda la funcionalidad que ofrecen.
Puedes ver la versión de \(bash\) de
tu sistema con el siguiente comando bash –version.
Aquí encuentras la versión actual del manual oficial de Bash.
5.1 Ejemplo de invocación del script phyml_protModelFinder.sh
Vamos a correr el \(phyml\_protModelFinder.sh\) para automatizar la selección de modelos para proteínas y la inferencia de la filogenia bajo el mejor modelo encontrado.
phyml_protModelFinder.sh v0.9_2023-11-17 requires two arguments provided on the command line: phyml_protModelFinder.sh# model sets to choose from: 1 -> nuclear genes (AB BLOSUM62 DAYHOFF DCMut JTT LG VT WAG) 2 -> organellar genes (CpREV MTMAM MtREV MtArt) 3 -> nuclear and organellar (1 + 2) 4 -> retroviral genes (HIVw HIVb RtREV) 5 -> nuclear + retroviral genes (1 + 4) 6 -> all (1+2+3+4+5) 7 -> test (JTT LG) EXAMPLE: phyml_protModelFinder.sh primates_21_AA.phy 5 10 AIM: phyml_protModelFinder.sh v0.9_2023-11-17 will evaluate the fit of the the seleced model set, combined or not with +G and/or +f, computing AICi, BICi, deltaBIC, BICw and inferring the ML tree under the BIC-selected model PROCEDURE - Models are fitted using a fixed NJ-LG tree, optimizing branch lenghts and rates to calculate their AICi, BICi, delta_BIC and BICw - The best model is selected by BIC - SPR searches can be launched starting from multiple random trees - Default single seed tree searches use a BioNJ with BEST moves SOURCE: the latest version of the program is available from GitHub at: https://github.com/vinuesa/TIB-filoinfo/blob/master/phyml_protModelFinder.sh LICENSE: GPL v3.0. See https://github.com/vinuesa/get_phylomarkers/blob/master/LICENSE CITATION: Vinuesa 2023. Fast and easy selection of protein models for phylogenetics with phyml_protModelFinder. Available from GitHub at https://github.com/vinuesa/TIB-filoinfo/.
Necesitamos pasarle dos parámetros a phyml_protModelFinder.sh: i) el archivo \(phylip\) con las secuencias alineadas, y ii) elegir un set de matrices a evaluar.
Hay un tercer parámetro \(<int>\) opcional, que nos permite indicarle al programa que queremos hacer la búsqueda final de árboles a partir de \(n\) árboles de semilla de adición secuencial aleatorizada, bajo el mejor modelo seleccionado.
Probemos con el set (#5) de modelos para genes nucleares y retrovirales.
Abajo vemos el STDOUT del \(script\).
Antes de iniciar, el \(script\) comprueba que estamos corriendo una versión de \(bash >= 4.4\), ya que hace uso de idiomas modernos de \(Bash\), incluyendo \(hashes\). De hecho algunas funcionalidades requieren \(bash\ >=\ 4.8\). Si usas MacOS, probablemente debas hacer una actualización, disponible vía \(MacPorts\), \(Homebrew\) o \(Fink\).
Seguidamente comprueba los argumentos pasados al \(script\), determinado si existe el archivo de entrada y si está en formato phylip. Esto lo hace llamando a la función \(check\_is\_phylip\).
El primer cálculo consiste en evaluar las frecuencias empríricas de los aminoácidos con la función \(compute\_AA\_freq\_in\_phylip()\), que se muestra seguidamente:
function compute_AA_freq_in_phylip()
{
local phylip
phylip=$1
awk '
BEGIN{print "idx\tAA\tobs_freq"}
{
# ignore first row and column
if( NR > 1 && NF > 1){
# remove empty spaces
gsub(/[ ]+/," ")
l=length($0)
for(i=1; i<=l; i++){
c = substr($0, i, 1)
# count only standard amino acids
if (c ~ /[ARNDCQEGHILKMFPSTWYV]/){
ccounts[c]++
letters++
}
}
}
}
# print relative frequency of each residue
END {
for (c in ccounts){
aa++
printf "%i\t%s\t%.4f\n", aa, c, (ccounts[c] / letters )
}
}' "$phylip"
}Puedes copiar la función a la terminal y llamarla desde la línea de comandos:
Seguidamente construye un árbol NJ-LG que será usado, como en los ejemplos anteriores, para fijar la topología en las llamadas a \(phyml\) para optimizar parámetros y calcular el ajuste de \(lnL\)
Se entra al bucle principal del \(script\), que hace llamadas a \(phyml\) con cada uno de los modelos a evaluar, llenando diversos \(hashes\) con los resultados para su despliegue y posterior análisis estadístico.
Una vez estimados los parámetros, y calculados los puntajes de verosimilitud correspondientes, se realizan llamadas a funciones para calcular los \(AIC_i\), \(AIC_c\), \(BIC_i\), \(\Delta_{BIC}\), \(BIC_w\) y \(BIC_{cumW}\) asociados a cada modelo, se guardan los resultados en un \(hash\), y se imprime finalmente la tabla de resultados, ordenada por valores crecientes de BIC.
Por último se usa el mejor modelo (seleccionado por \(BIC\)), para estimar la filogenia. Por defecto ésta se estima a partir de un árbol semilla \(BioNJ\) y haciendo rearreglos tipo \(BEST\) (por defecto). Si se le pasa al \(script\) un íntegro como tercer parámetro, generará ese número de árboles de semilla de adición secuencial aleatorizada, a partir de los cuales hace búsquedas con rearreglos de tipo \(SPR\).
Es importante observar que el criterio bayesiano de información (BIC) es más conservador que los LRTs (no aplicables a proteínas) y AICs, los cuales tienden a elegir modelos sobreparametrizados. Como ejemplo de ésto, vean en la tabla que se muestra abajo que BIC selecciona a JTT+G, mientras que AIC selecciona JTT+f+G.
$ phyml_protModelFinder.sh primates_21_AA.phy 5
=========================================================================================
# Bash version 5.1 OK
# phyml_protModelFinder.sh v0.9_2023-11-17 running on alisio. Run started on: 2023-11-19 at 18:35:44
# running with phyml v.3.3_2023
/home/vinuesa/anaconda3/bin/mpirun
/usr/bin/phyml-mpi
/usr/bin/awk
/home/vinuesa/anaconda3/bin/bc
/usr/bin/sed
/home/vinuesa/anaconda3/bin/perl
/usr/local/bin/phyml
# Run check_dependencies() ... looks good: all required binaries are in place.
# infile:primates_21_AA.phy; model_set:5; mpi_OK:1; seed trees: 1
=========================================================================================
18:35:44 # 1. Computing sequence stats for primates_21_AA.phy:
- number of sequences: 21
- number of sites: 500
- number of branches: 39
- observed amino acid frequencies:
idx AA obs_freq
1 N 0.0423
2 A 0.0442
3 P 0.0439
4 C 0.0396
5 D 0.0399
6 Q 0.0656
7 R 0.0545
8 E 0.0939
9 F 0.0381
10 S 0.0713
11 T 0.0424
12 G 0.0414
13 H 0.0287
14 V 0.0763
15 I 0.0491
16 W 0.0135
17 K 0.0666
18 Y 0.0285
19 L 0.0961
20 M 0.0243
--------------------------------------------------------------------------------
18:35:44 1. Computing NJ-LG tree for input file primates_21_AA.phy with 21 sequences
--------------------------------------------------------------------------------
2. running in a for loop to combine all base models in model_set 5=>nuclear_and_viral,
with (or not) +G and or +f, and compute the model lnL, after optimizing branch lengths and rates
--------------------------------------------------------------------------------
18:35:44 # running: phyml -i primates_21_AA.phy -d aa -m AB -u primates_21_AA.phy_LG-NJ.nwk -c 1 -v 0 -o lr
18:35:44 # running: phyml -i primates_21_AA.phy -d aa -m AB -f m -c 4 -a e -u primates_21_AA.phy_LG-NJ.nwk -o lr
18:35:45 # running: phyml -i primates_21_AA.phy -d aa -m AB -f e -c 1 -u primates_21_AA.phy_LG-NJ.nwk -o lr
18:35:45 # running: phyml -i primates_21_AA.phy -d aa -m AB -u primates_21_AA.phy_LG-NJ.nwk -f e -a e -o lr
18:35:46 # running: phyml -i primates_21_AA.phy -d aa -m BLOSUM62 -u primates_21_AA.phy_LG-NJ.nwk -c 1 -v 0 -o lr
18:35:46 # running: phyml -i primates_21_AA.phy -d aa -m BLOSUM62 -f m -c 4 -a e -u primates_21_AA.phy_LG-NJ.nwk -o lr
18:35:47 # running: phyml -i primates_21_AA.phy -d aa -m BLOSUM62 -f e -c 1 -u primates_21_AA.phy_LG-NJ.nwk -o lr
18:35:47 # running: phyml -i primates_21_AA.phy -d aa -m BLOSUM62 -u primates_21_AA.phy_LG-NJ.nwk -f e -a e -o lr
18:35:48 # running: phyml -i primates_21_AA.phy -d aa -m DAYHOFF -u primates_21_AA.phy_LG-NJ.nwk -c 1 -v 0 -o lr
18:35:48 # running: phyml -i primates_21_AA.phy -d aa -m DAYHOFF -f m -c 4 -a e -u primates_21_AA.phy_LG-NJ.nwk -o lr
18:35:49 # running: phyml -i primates_21_AA.phy -d aa -m DAYHOFF -f e -c 1 -u primates_21_AA.phy_LG-NJ.nwk -o lr
18:35:49 # running: phyml -i primates_21_AA.phy -d aa -m DAYHOFF -u primates_21_AA.phy_LG-NJ.nwk -f e -a e -o lr
18:35:50 # running: phyml -i primates_21_AA.phy -d aa -m DCMut -u primates_21_AA.phy_LG-NJ.nwk -c 1 -v 0 -o lr
18:35:50 # running: phyml -i primates_21_AA.phy -d aa -m DCMut -f m -c 4 -a e -u primates_21_AA.phy_LG-NJ.nwk -o lr
18:35:51 # running: phyml -i primates_21_AA.phy -d aa -m DCMut -f e -c 1 -u primates_21_AA.phy_LG-NJ.nwk -o lr
18:35:51 # running: phyml -i primates_21_AA.phy -d aa -m DCMut -u primates_21_AA.phy_LG-NJ.nwk -f e -a e -o lr
18:35:52 # running: phyml -i primates_21_AA.phy -d aa -m JTT -u primates_21_AA.phy_LG-NJ.nwk -c 1 -v 0 -o lr
18:35:53 # running: phyml -i primates_21_AA.phy -d aa -m JTT -f m -c 4 -a e -u primates_21_AA.phy_LG-NJ.nwk -o lr
18:35:53 # running: phyml -i primates_21_AA.phy -d aa -m JTT -f e -c 1 -u primates_21_AA.phy_LG-NJ.nwk -o lr
18:35:54 # running: phyml -i primates_21_AA.phy -d aa -m JTT -u primates_21_AA.phy_LG-NJ.nwk -f e -a e -o lr
18:35:54 # running: phyml -i primates_21_AA.phy -d aa -m LG -u primates_21_AA.phy_LG-NJ.nwk -c 1 -v 0 -o lr
18:35:55 # running: phyml -i primates_21_AA.phy -d aa -m LG -f m -c 4 -a e -u primates_21_AA.phy_LG-NJ.nwk -o lr
18:35:56 # running: phyml -i primates_21_AA.phy -d aa -m LG -f e -c 1 -u primates_21_AA.phy_LG-NJ.nwk -o lr
18:35:56 # running: phyml -i primates_21_AA.phy -d aa -m LG -u primates_21_AA.phy_LG-NJ.nwk -f e -a e -o lr
18:35:57 # running: phyml -i primates_21_AA.phy -d aa -m VT -u primates_21_AA.phy_LG-NJ.nwk -c 1 -v 0 -o lr
18:35:57 # running: phyml -i primates_21_AA.phy -d aa -m VT -f m -c 4 -a e -u primates_21_AA.phy_LG-NJ.nwk -o lr
18:35:58 # running: phyml -i primates_21_AA.phy -d aa -m VT -f e -c 1 -u primates_21_AA.phy_LG-NJ.nwk -o lr
18:35:58 # running: phyml -i primates_21_AA.phy -d aa -m VT -u primates_21_AA.phy_LG-NJ.nwk -f e -a e -o lr
18:35:59 # running: phyml -i primates_21_AA.phy -d aa -m WAG -u primates_21_AA.phy_LG-NJ.nwk -c 1 -v 0 -o lr
18:35:59 # running: phyml -i primates_21_AA.phy -d aa -m WAG -f m -c 4 -a e -u primates_21_AA.phy_LG-NJ.nwk -o lr
18:36:00 # running: phyml -i primates_21_AA.phy -d aa -m WAG -f e -c 1 -u primates_21_AA.phy_LG-NJ.nwk -o lr
18:36:00 # running: phyml -i primates_21_AA.phy -d aa -m WAG -u primates_21_AA.phy_LG-NJ.nwk -f e -a e -o lr
18:36:01 # running: phyml -i primates_21_AA.phy -d aa -m HIVw -u primates_21_AA.phy_LG-NJ.nwk -c 1 -v 0 -o lr
18:36:01 # running: phyml -i primates_21_AA.phy -d aa -m HIVw -f m -c 4 -a e -u primates_21_AA.phy_LG-NJ.nwk -o lr
18:36:02 # running: phyml -i primates_21_AA.phy -d aa -m HIVw -f e -c 1 -u primates_21_AA.phy_LG-NJ.nwk -o lr
18:36:02 # running: phyml -i primates_21_AA.phy -d aa -m HIVw -u primates_21_AA.phy_LG-NJ.nwk -f e -a e -o lr
18:36:03 # running: phyml -i primates_21_AA.phy -d aa -m HIVb -u primates_21_AA.phy_LG-NJ.nwk -c 1 -v 0 -o lr
18:36:03 # running: phyml -i primates_21_AA.phy -d aa -m HIVb -f m -c 4 -a e -u primates_21_AA.phy_LG-NJ.nwk -o lr
18:36:04 # running: phyml -i primates_21_AA.phy -d aa -m HIVb -f e -c 1 -u primates_21_AA.phy_LG-NJ.nwk -o lr
18:36:04 # running: phyml -i primates_21_AA.phy -d aa -m HIVb -u primates_21_AA.phy_LG-NJ.nwk -f e -a e -o lr
18:36:05 # running: phyml -i primates_21_AA.phy -d aa -m RtREV -u primates_21_AA.phy_LG-NJ.nwk -c 1 -v 0 -o lr
18:36:05 # running: phyml -i primates_21_AA.phy -d aa -m RtREV -f m -c 4 -a e -u primates_21_AA.phy_LG-NJ.nwk -o lr
18:36:06 # running: phyml -i primates_21_AA.phy -d aa -m RtREV -f e -c 1 -u primates_21_AA.phy_LG-NJ.nwk -o lr
18:36:06 # running: phyml -i primates_21_AA.phy -d aa -m RtREV -u primates_21_AA.phy_LG-NJ.nwk -f e -a e -o lr
18:36:07 # writing primates_21_AA.phy_sorted_model_set_5_fits.tsv, sorted by BIC
--------------------------------------------------------------------------------------------------
model K sites/K lnL AIC AICc BIC deltaBIC BICw BICcumW
HIVb+G 40 12.50000 -4785.64708 9651.29416 9658.44013 9819.87848 0 0.97 0.97
HIVw+f+G 59 8.47000 -4729.97721 9577.95442 9594.04533 9826.61630 6.73782 0.03 1.00
HIVb+f+G 59 8.47000 -4745.24175 9608.48350 9624.57441 9857.14538 37.26690 0.00 1.00
HIVw+G 40 12.50000 -4811.89672 9703.79344 9710.93941 9872.37776 52.49928 0.00 1.00
JTT+G 40 12.50000 -4815.80148 9711.60296 9718.74893 9880.18728 60.30880 0.00 1.00
HIVw 39 12.82000 -4839.90549 9757.81098 9764.59359 9922.18070 102.30222 0.00 1.00
JTT+f+G 59 8.47000 -4778.48904 9674.97808 9691.06899 9923.63996 103.76148 0.00 1.00
HIVb 39 12.82000 -4845.99268 9769.98536 9776.76797 9934.35508 114.47660 0.00 1.00
JTT 39 12.82000 -4865.61699 9809.23398 9816.01659 9973.60370 153.72522 0.00 1.00
AB+f+G 59 8.47000 -4814.44420 9746.88840 9762.97931 9995.55028 175.67180 0.00 1.00
HIVw+f 58 8.62000 -4839.90549 9795.81098 9811.33025 10040.25825 220.37977 0.00 1.00
HIVb+f 58 8.62000 -4845.99268 9807.98536 9823.50463 10052.43263 232.55415 0.00 1.00
AB+G 40 12.50000 -4910.22076 9900.44152 9907.58749 10069.02584 249.14736 0.00 1.00
WAG+G 40 12.50000 -4910.24193 9900.48386 9907.62983 10069.06818 249.18970 0.00 1.00
VT+G 40 12.50000 -4916.81420 9913.62840 9920.77437 10082.21272 262.33424 0.00 1.00
JTT+f 58 8.62000 -4865.61699 9847.23398 9862.75325 10091.68125 271.80277 0.00 1.00
WAG+f+G 59 8.47000 -4870.83180 9859.66360 9875.75451 10108.32548 288.44700 0.00 1.00
AB 39 12.82000 -4933.84484 9945.68968 9952.47229 10110.05940 290.18092 0.00 1.00
VT+f+G 59 8.47000 -4880.25893 9878.51786 9894.60877 10127.17974 307.30126 0.00 1.00
LG+G 40 12.50000 -4946.47152 9972.94304 9980.08901 10141.52736 321.64888 0.00 1.00
WAG 39 12.82000 -4961.23284 10000.46568 10007.24829 10164.83540 344.95692 0.00 1.00
VT 39 12.82000 -4968.11339 10014.22678 10021.00939 10178.59650 358.71802 0.00 1.00
DAYHOFF+G 40 12.50000 -4969.87622 10019.75244 10026.89841 10188.33676 368.45828 0.00 1.00
DCMut+G 40 12.50000 -4970.16584 10020.33168 10027.47765 10188.91600 369.03752 0.00 1.00
LG+f+G 59 8.47000 -4913.60179 9945.20358 9961.29449 10193.86546 373.98698 0.00 1.00
DAYHOFF+f+G 59 8.47000 -4920.31326 9958.62652 9974.71743 10207.28840 387.40992 0.00 1.00
DCMut+f+G 59 8.47000 -4920.53388 9959.06776 9975.15867 10207.72964 387.85116 0.00 1.00
AB+f 58 8.62000 -4933.84484 9983.68968 9999.20895 10228.13695 408.25847 0.00 1.00
BLOSUM62+G 40 12.50000 -5002.18584 10084.37168 10091.51765 10252.95600 433.07752 0.00 1.00
LG 39 12.82000 -5008.95318 10095.90636 10102.68897 10260.27608 440.39760 0.00 1.00
WAG+f 58 8.62000 -4961.23284 10038.46568 10053.98495 10282.91295 463.03447 0.00 1.00
DAYHOFF 39 12.82000 -5021.06263 10120.12526 10126.90787 10284.49498 464.61650 0.00 1.00
DCMut 39 12.82000 -5021.17208 10120.34416 10127.12677 10284.71388 464.83540 0.00 1.00
VT+f 58 8.62000 -4968.11339 10052.22678 10067.74605 10296.67405 476.79557 0.00 1.00
RtREV+G 40 12.50000 -5027.95301 10135.90602 10143.05199 10304.49034 484.61186 0.00 1.00
BLOSUM62+f+G 59 8.47000 -4976.30875 10070.61750 10086.70841 10319.27938 499.40090 0.00 1.00
RtREV+f+G 59 8.47000 -4984.38124 10086.76248 10102.85339 10335.42436 515.54588 0.00 1.00
LG+f 58 8.62000 -5008.95318 10133.90636 10149.42563 10378.35363 558.47515 0.00 1.00
BLOSUM62 39 12.82000 -5071.35214 10220.70428 10227.48689 10385.07400 565.19552 0.00 1.00
DAYHOFF+f 58 8.62000 -5021.06263 10158.12526 10173.64453 10402.57253 582.69405 0.00 1.00
DCMut+f 58 8.62000 -5021.17208 10158.34416 10173.86343 10402.79143 582.91295 0.00 1.00
RtREV 39 12.82000 -5084.78701 10247.57402 10254.35663 10411.94374 592.06526 0.00 1.00
BLOSUM62+f 58 8.62000 -5071.35214 10258.70428 10274.22355 10503.15155 683.27307 0.00 1.00
RtREV+f 58 8.62000 -5084.78701 10285.57402 10301.09329 10530.02129 710.14281 0.00 1.00
--------------------------------------------------------------------------------------------------
* NOTE 1: when sites/K < 40, the AICc is recommended over AIC.
* NOTE 2: The best model is selected by BIC, because AIC is biased, favoring parameter-rich models.
==================================================================================================
... will estimate the ML tree under best-fitting model HIVb+G selected by BIC
18:36:09 # running: phyml -i primates_21_AA.phy -d aa -m HIVb -f m -c 4 -a e -o tlr -s BEST
# Your results:
- primates_21_AA.phy_HIVb+G_BESTmoves_phyml_stats.txt
- primates_21_AA.phy_HIVb+G_BESTmoves_phyml_tree.txt
--------------------------------------------------------------------------------------------------
Elapsed time: 0 days, 00 hr, 00 min, 33 sec
Done!
========================================================================================
If you use phyml_protModelFinder.sh v.0.9_2023-11-17 for your research,
I would appreciate that you:
1. Cite the code in your work as:
Pablo Vinuesa. phyml_protModelFinder.sh v.0.9_2023-11-17
https://github.com/vinuesa/TIB-filoinfo/blob/master/phyml_protModelFinder.sh
2. Give it a like on the https://github.com/vinuesa/TIB-filoinfo/ repo
Thanks!
Podemos desplegar el árbol directamente en la terminal usando herramientas del paquete \(newick\_utilities\).
+-+ Gorilla
|
+---+ 0.999542Human
| +-+ 0.988925
+-------------+ 0.999985Chimp
| |
| |+-------------------------+ Gibbon
| ++ 0.731672
| +----------+ Orangutan
|
+--------------------------------------------------------+ 1.000000 +--------+ Patas
| | |
| | +---+ 1.000000AGM cDNA
| | | | +--+ 1.000000
| | | | | +-+ Tant cDNA
| | | +--+ 0.999995
| | | |-+ Baboon
| +----------------------+ 1.000000.977694
| | +----+ Rhes cDNA
| |
=| |+--+ Colobus
| ++ 0.776704
| +-+ DLangur
|
| +-----------------------+ Owl
| |
| +-+ 0.550063 +---------------+ Tamarin
| | +-------------+ 1.000000
| | +---------+ PMarmoset
| |
+--------------------------------------------------------+ 0.419386----------------------------------+ Squirrel
| |
|+-+ 0.633907+-------------+ Titi
|| +---------+ 1.000000
++ 1.000000 +----------------+ Saki
|
| +---------------------------------+ Howler
| |
+---------------+ 0.999170-------------+ Woolly
+-----+
+----------+ Spider
|----------------|-----------------|----------------|----------------|-----------------|----------------|----
0 0.05 0.1 0.15 0.2 0.25 0.3
substitutions/site
Figura 6. Filogenia de máxima verosimilitud bajo el modelo \(HIVb+G\) seleccionado automáticamente por \(phyml\_protModelFinder.sh\) en base a valores de \(BIC\), y desplegada directamente en la terminal con programas del paquete \(newick\_utilities\).
Nota cómo llamamos a \(nw\_reroot\)
para hacer un enraizamiento en el punto medio, pasando
el resultado a \(nw\_display\),
haciendo uso de lo que se conoce como sustitución de
proceso comando2 <(comando1), que pasa la
salida de comando1 como si fuera un archivo que lee
comando2.
Gracias a la sustitución de proceso evitamos tener que escribir un archivo intermedio con la filogenia enraizada, para ser desplegada por \(nw\_display\)
Abajo te muestro un ejemplo más avanzado del uso de sustitución de proceso que hacen los \(scripts\) \(phyml\_DNAmodelFinder\) y \(phyml\_protModelFinder\) para añadir las columnas \(\Delta_{BIC}\), \(BIC_w\) y \(BIC_{cumW}\) a la tabla de resultados.
# 6.4 paste the BIC_deltas_a & BICw_a values as a new column to
# "${infile}"_sorted_model_set_"${model_set}"_fits.tsv
paste "${infile}"_sorted_model_set_"${model_set}"_fits.tsv \
<(for i in "${BIC_deltas_a[@]}"; do echo "$i"; done) \
<(for i in "${BICw_a[@]}"; do echo "$i"; done) \
<(for i in "${BICcumW_a[@]}"; do echo "$i"; done) > t
[[ -s t ]] && mv t "${infile}"_sorted_model_set_"${model_set}"_fits.tsv
5.2 El criterio bayesiano de información (\(BIC\))
Tal vez notaste que en nuestro caso el mejor modelo seleccionado por \(BIC\) es el \(HIVb+G\), mientras que \(AIC\) selecciona \(HIVw+f+G\). Es más, vemos que \(HIVw+f+G\) tiene un \(\Delta_{BIC}\) > 6, por lo que es un modelo claramente peor que el \(HIVb+G\). Es tan grande la distancia \(\Delta_{BIC}\) entre ambos modelos, que \(HIVw+f+G\) tiene un \(BIC_w\ =\ 0.03\).
El criterio de información de Akaike (\(AIC\)) es sabido que tiende a favorecer modelos con más parámetros, y con frecuencia selecciona modelos sobreparametrizados.
Por ello los \(scripts\) \(phyml\_protModelFinder.sh\) y \(phyml\_DNAmodelFinder.sh\) usan el \(BIC\) para la selección de modelos.
El \(BIC\) para un modelo se define de la siguiente manera:
\(BIC_i=-2lnL_i + K_i * log(N)\)
donde:
- \(lnL_i=log-Likelihood\ del\ modelo\ i\)
- \(K=número\ de\ parámetros\ libres\ del\ modelo\ i\)
- \(N=número\ de\ sitios\ del\ alineamiento\),
Noten que \(K_i\) debería considerar además el número de ramas libres a estimar (\(2n-3\); \(n = num\_secuencias\)), por lo que:
- \(K_i=número\ de\ parámetros\ libres\ del\ modelo\ i + (2n - 3)\)
Al igual que con los \(AIC\), el modelo con menor \(BIC\) es el preferido, el cual equivale al modelo de máxima probabilidad posterior.
- ¿A qué se debe esta discrepancia entre \(AIC\) y \(BIC\)?
Dado que para muchos alineamientos \(log(N) > 2\), el \(BIC\) suele seleccionar modelos con menos parámetros que \(AIC\).
5.2.1 Cálculo de distancias (\({\Delta}AIC_i\) o \({\Delta}BIC_i\)) y ponderaciones de Akaike (\(AIC_{wi}\)) o bayesianas (\(BIC_{wi}\)) para cada modelo en competición
Con el \(AIC\) y el \(BIC\) comparamos el ajuste de todos los modelos en competición simultáneamente, independientemente de si están o no anidados entre sí.
Por ello para cada \(modelo_i\) podemos:
- Ordenar a los modelos en competición de menor (mejor) a mayor (peor) AIC o BIC. El mejor modelo tendrá el del \(AIC_{min}\) o \(BIC_{min}\)
- Calcular la diferencia en unidades de \(AIC\) o \(BIC\) de cada modelo con respecto al \(AIC_{min}\) o \(BIC_{min}\), es decir, su \({\Delta}AIC_i\) o \({\Delta}BIC_i\)
Recuerda que, de acuerdo con Burnham and Anderson (2002). Model Selection and Multimodel Inference, Springer, p.70:
- Modelos a una distancia de 0-2 unidades de \(AIC\) o \(BIC\) del \(AIC_{min}\) o \(BIC_{min}\) tienen un soporte sustancial
- Modelos a una distancia de 4-7 unidades de \(AIC\) o \(BIC\) del \(AIC_{min}\) o \(BIC_{min}\) alternativas considerablemente peores
- Modelos a una distancia de > 10 unidades de \(AIC\) o \(BIC\) del \(AIC_{min}\) o \(BIC_{min}\) los podemos despreciar, ya que carecen de soporte
Esto queda claramente reflejado en los pesos o ponderaciones de Akaike y bayesianas para cada modelo, su \({\Delta}AIC_i\) o \({\Delta}BIC_i\), como te muestro seguidamente:
5.2.2 Cálculo de ponderaciones de cada modelo y su peso acumulado en base a un conjunto de valores ordenados de (\(\Delta BIC\))
Cálculo de pesos de cada modelo \(BIC_{wi}\): El peso de cada modelo viene dado por la fórmula: \[ w_i = \frac{e^{-\Delta BIC_i / 2}}{\sum_{j=1}^{N} e^{-\Delta BIC_j / 2}} \] donde \(w_i\) es el peso del modelo \(i\), y el denominador es la suma de pesos de todos los modelos evaluados, normalizando las ponderaciones para que sumen 1.
Cálculo de pesos acumulados: El peso acumulado de cada modelo el la sumatoria consecutiva de los pesos de una lista ordenada de modelos.
5.2.2.1 Ejemplo detallado, paso a paso de cálculo de \(BIC_{wi}\) a partir de una lista ordenada de \(\Delta BIC\)
Dada la siguiente lista ordenada de valores de \(\Delta BIC\): \[ \Delta BIC = [0, 4.02575, 6.80293, 10.82696, 26.16869, 27.94213] \]
Calcula \(e^{-\Delta BIC_i / 2}\) para cada \(\Delta BIC\):
- Para \(\Delta BIC_0 = 0\): \(e^{0 / 2} = 1.00000\)
- Para \(\Delta BIC_1 = 4.02575\): \(e^{-4.02575 / 2} = 0.18648\)
- Para \(\Delta BIC_2 = 6.80293\): \(e^{-6.80293 / 2} = 0.03305\)
- Para \(\Delta BIC_3 = 10.82696\): \(e^{-10.82696 / 2} = 0.00447\)
- Para \(\Delta BIC_4 = 26.16869\): \(e^{-26.16869 / 2} = 0.0000076\)
- Para \(\Delta BIC_5 = 27.94213\): \(e^{-27.94213 / 2} = 0.0000053\)
Calcula la sumatoria de estos pesos: \[ \text{Sumatoria de pesos} = 1.00000 + 0.18648 + 0.03305 + 0.00447 + 0.0000076 + 0.0000053 \approx 1.22402 \]
Calcula el peso normalizado de cada modelo:
- \(w_0 = \frac{1.00000}{1.22402} \approx 0.81719\)
- \(w_1 = \frac{0.18648}{1.22402} \approx 0.15236\)
- \(w_2 = \frac{0.03305}{1.22402} \approx 0.02700\)
- \(w_3 = \frac{0.00447}{1.22402} \approx 0.00365\)
- \(w_4 = \frac{0.0000076}{1.22402} \approx 0.0000062\)
- \(w_5 = \frac{0.0000053}{1.22402} \approx 0.0000043\)
Calcula los pesos o ponderaciones acumuladas:
- Peso acumulado para el modelo 0: \(0.81719\)
- Peso acumulado para el modelo 1: \(0.81719 + 0.15236 = 0.96955\)
- Peso acumulado para el modelo 2: \(0.96955 + 0.02700 = 0.99655\)
- Peso acumulado para el modelo 3: \(0.99655 + 0.00365 = 1.00020\) (approx 1.000 due to rounding)
- Peso acumulado para modelos 4 y 5 queda en \(1.000\).
A partir de estos datos podemos determinar que con los modelos 0 y 1 alcanzamos el intervalo de credibilidad del 95%.
5.2.3 Código de AWK para calcular los pesos y pesos acumulados de cada modelo
Aquí te dejo un pequeño fragmento de código \(AWK\) para el cálculo de \(BIC_{wi}\) y \(BIC_{Cum_{wi}}\), en base a una lista ordenada de valores de \(\Delta BIC\), similar al que usan los \(scripts\) \(phyml\_protModelFinder.sh\) y \(phyml\_DNAmodelFinder.sh\)
BEGIN {
# Delta BIC values
delta_bics[0] = 0
delta_bics[1] = 4.02575
delta_bics[2] = 6.80293
delta_bics[3] = 10.82696
delta_bics[4] = 26.16869
delta_bics[5] = 27.94213
n = 6 # Number of models
# Step 1: Calculate exp(-Delta BIC / 2) for each model
sum_exp = 0
for (i = 0; i < n; i++) {
exp_values[i] = exp(-delta_bics[i] / 2)
sum_exp += exp_values[i]
}
# Step 2: Calculate weights and cumulative weights
cumulative_weight = 0
printf "%-10s %-10s %-10s %-10s\n", "Model", "Delta BIC", "Weight", "Cum. Weight"
for (i = 0; i < n; i++) {
weight = exp_values[i] / sum_exp
cumulative_weight += weight
printf "%-10d %-10.5f %-10.5f %-10.5f\n", i+1, delta_bics[i], weight, cumulative_weight
}
}5.2.3.1 Explicación de cada función
exp(-delta_bics[i] / 2): calcula el exponencial de \(-\Delta BIC / 2\) para cada modelo, requerido para derivar su peso.sum_exp += exp_values[i]: Esta línea acumula la suma de todos los valores \(e^{-\Delta BIC / 2}\) para la normalización.weight = exp_values[i] / sum_exp: El peso de cada modelo es calculado dividiendo su valor exponencial por la sumatoria total.cumulative_weight += weight: El peso acumulado secuencialmente de la lista ordenada (decreciente) de pesos para cada modelo.
5.2.4 Tarea: selección de modelos con BIC
- Calcula manualmente los valores de \(BICi\), \(BIC_{MIN}\) y \(\Delta_{BIC}\) para los modelo \(HIVb+G\) y \(HIVb+f+G\) del alineamiento de primates_21_AA.phy, usando la fórmula mostrada en la sección de \(BIC\) y los datos que imprime \(phyml\_protModelFinder.sh\) al inicio de su salida.
- Vuelve a correr el script, pero evaluando todos los modelos
phyml_protModelFinder.sh primates_21_AA.phy 6¿Hay alguna diferencia? Comenta el resultado. - Si no lo has hecho aún, repite el ejercicio de selección de modelos paramétricos (de DNA) con las secuencias primates.phy, usando \(BIC\) para evaluar los modelos simplificados a partir de los parámetros estiamdos por GTR+I+G.
5.3 Selección “inteligente” de subconjuntos de modelos de sustitución para DNA con phyml_DNAmodelFinder.sh
El \(script\) phyml_DNAmodelFinder.sh tiene una estructura similar a phyml_protModelFinder.sh, pero es algo más complejo, ya que el usuario puede pedirle que evalúe el conjunto de modelos estándar (JC69, K80, F81, HKY85, TN93 y GTR) más un conjunto amplio de modelos transicionales (TIM) y transversionales (TVM), con 3, 4 y 5 clases de tasas, los cuales son elegidos después de comprobar si existen o no sesgos en frecuencia de bases y transicional significativos.
NOTA: \(phyml\_DNAmodelFinder.sh\) hace uso de
idiomas más avanzados y recientes del lenguaje Bash que \(phyml\_protModelFinder.sh\), ya que además
de arreglos y hashes, usa
referencias a los mismos para pasarlos a funciones. Por
tanto la funcionalidad completa del script sólo está disponible si usas
una versión de \(bash >= 5.1\).
Veamos la ayuda de \(phyml\_DNAmodelFinder.sh\):
phyml_DNAmodelFinder.sh v1.2.0_2023-11-18 requires two arguments provided on the command line, the third one being optional: phyml_DNAmodelFinder.sh# model sets to choose from: 1 -> standard models (JC69 K80 F81 HKY85 TN93 GTR) 2 -> standard + 64 extended_ef_models, OR standard + 62 extended_uf_models NOTE: phyml_DNAmodelFinder.sh automatically chooses the proper extended set (ef|uf) to evaluate, based on delta_BIC evaluation of compositional bias (JC69 vs. F81) 3 -> minimal test set (JC69 F81 HKY85 TN93) EXAMPLE: phyml_DNAmodelFinder.sh primates.phy 2 AIM: phyml_DNAmodelFinder.sh v1.2.0_2023-11-18 will evaluate the fit of the selected model set, combined or not with +G and/or +f, computing AICi, BICi, deltaBIC and BICw, inferring the ML tree under the BIC-selected model. PROCEDURE: - Models are fitted using a fixed BioNJ-JC tree, optimizing branch lengths and rates to calculate their AICi, BICi, delta_BIC and BICw - Only relevant models among the extended set are evaluated, based on delta_BIC comparisons between JC69-F81, to decide if ef|uf models should be evaluated, and comparisons between KHY85-TN93 to determine if models with two Ti rates should be evaluated - pInv is automatically excluded in the extended model set, if the delta_BICi_HKY+G is =< 2 when compared with the BIC_HKY+G+I - The best model is selected by BIC - SPR searches can be launched starting from multiple random trees - Default single seed tree searches use a BioNJ tree with BEST moves SOURCE: the latest version of the program is available on GitHub: https://github.com/vinuesa/TIB-filoinfo LICENSE: GPL v3.0. See https://github.com/vinuesa/get_phylomarkers/blob/master/LICENSE CITATION: Vinuesa P. (2023). Efficient two-step selection of DNA models with phyml_DNAmodelFinder. https://github.com/vinuesa/TIB-filoinfo
Al elegir el set 2 de modelos, en vez de evaluar a todos los modelos estándar+TIM+TVM, phyml_DNAmodelFinder.sh hace unas evaluaciones de valores de \(\Delta_{BIC}\) entre los siguientes pares de modelos:
- *JC69 vs. F81* para determinar si existe un **sesgo composicional** significativo - *HKY85 vs. TN93* para evaluar si existe un **sesgo transicional entre purinas y pirimidinas** significativo - *HKY85+G vs. HKY85+I+G* para calcular si amerita o no incluir una **proporción de sitios invariantes**
En función de estos análisis, efectuados sobre los modelos del conjunto estándar, el \(script\) decide qué conjuntos de modelos adicionales evaluar:
- TIM con con uf|ef
- TVM con con uf|ef
- añadir o no +I
5.3.1 Ejemplo de corrida phyml_DNAmodelFinder.sh en modo de selección de set extendido de modelos
Al igual que phyml_protModelFinder.sh, phyml_DNAmodelFinder.sh primero verifica:
- el ambiente en el que corre (verión de \(Bash\) y dependencias)
- que las secuencias están en formato \(phylip\)
- calcula unas estadísticas descriptivas de las secuencias
========================================================================================= # Bash version 5.1 OK # running with phyml v.3.3_2023 # phyml_DNAmodelFinder.sh v1.2.0_2023-11-18 running on alisio. Run started on: 2023-11-19 at 18:54:14 # workding directory: /home/vinuesa/Cursos/LCG/ENS_y_PANGENOMICA/tema12_maxima_verosimilitud/data /home/vinuesa/anaconda3/bin/mpirun /usr/bin/phyml-mpi /usr/bin/awk /home/vinuesa/anaconda3/bin/bc /usr/bin/sed /home/vinuesa/anaconda3/bin/perl /usr/local/bin/phyml # Run check_dependencies() ... looks good: all required binaries are in place. # infile:primates.phy; model_set:2 => extended_models; seed trees: 1; delta_BIC_cutoff=2 ========================================================================================= 18:54:14 # 1. Computing sequence stats for primates.phy: - number of sequences: 12 - number of sites: 898 - number of branches: 21 - observed nucleotide frequencies: idx NT obs_freq 1 A 0.3241 2 C 0.3040 3 T 0.2664 4 G 0.1056 --------------------------------------------------------------------------------
Queda claro del análisis de frecuencias empríricas que existe un sesgo composicional importante en los datos, por lo que phyml_DNAmodelFinder.sh debería favorecer modelos que lo consideren.
En los siguientes pasos:
- Calcula un árbol NJ-JC que servirá de topología de referencia
- Ajusta los modelos de base, haciendo uso del árbol NJ-JC
- Calucula los \(AIC\), \(AICc\) y \(BIC\) para estos modelos de base e imprime a pantalla los valores de \(BIC_i\) para los modelos arriba enlistados, en base a los cuales toma la decisión de qué tipo de modelos adicionales evaluar.
18:54:14 1. Computing NJ-JC tree for input file primates.phy with 12 sequences
--------------------------------------------------------------------------------
2.1. running in a for loop to combine all base models in model_set 2=>extended_models,
with (or without) +G and or +I, and compute the model lnL, after optimizing branch lengths and rates
18:54:14 # running: phyml -i primates.phy -d nt -m JC69 -f m -u primates.phy_JC-NJ.nwk -c 1 -v 0 -o lr
JC_BICi: 12991.20847
18:54:15 # running: phyml -i primates.phy -d nt -m JC69 -f m -c 4 -a e -u primates.phy_JC-NJ.nwk -o lr
18:54:15 # running: phyml -i primates.phy -d nt -m JC69 -f m -v e -c 1 -u primates.phy_JC-NJ.nwk -o lr
18:54:15 # running: phyml -i primates.phy -d nt -m JC69 -f m -u primates.phy_JC-NJ.nwk -v e -a e -o lr
18:54:15 # running: phyml -i primates.phy -d nt -m K80 -f m -u primates.phy_JC-NJ.nwk -c 1 -v 0 -o lr
18:54:15 # running: phyml -i primates.phy -d nt -m K80 -f m -c 4 -a e -u primates.phy_JC-NJ.nwk -o lr
18:54:15 # running: phyml -i primates.phy -d nt -m K80 -f m -v e -c 1 -u primates.phy_JC-NJ.nwk -o lr
18:54:15 # running: phyml -i primates.phy -d nt -m K80 -f m -u primates.phy_JC-NJ.nwk -v e -a e -o lr
18:54:16 # running: phyml -i primates.phy -d nt -m F81 -f m -u primates.phy_JC-NJ.nwk -c 1 -v 0 -o lr
F81_BICi: 12733.19534
18:54:16 # running: phyml -i primates.phy -d nt -m F81 -f m -c 4 -a e -u primates.phy_JC-NJ.nwk -o lr
18:54:16 # running: phyml -i primates.phy -d nt -m F81 -f m -v e -c 1 -u primates.phy_JC-NJ.nwk -o lr
18:54:16 # running: phyml -i primates.phy -d nt -m F81 -f m -u primates.phy_JC-NJ.nwk -v e -a e -o lr
18:54:16 # running: phyml -i primates.phy -d nt -m HKY85 -f m -u primates.phy_JC-NJ.nwk -c 1 -v 0 -o lr
HKY85_BICi: 12133.44475
18:54:16 # running: phyml -i primates.phy -d nt -m HKY85 -f m -c 4 -a e -u primates.phy_JC-NJ.nwk -o lr
HKY85+G_BICi: 11475.87798
18:54:17 # running: phyml -i primates.phy -d nt -m HKY85 -f m -v e -c 1 -u primates.phy_JC-NJ.nwk -o lr
18:54:17 # running: phyml -i primates.phy -d nt -m HKY85 -f m -u primates.phy_JC-NJ.nwk -v e -a e -o lr
HKY85+I+G_BICi: 11477.88060
18:54:18 # running: phyml -i primates.phy -d nt -m TN93 -f m -u primates.phy_JC-NJ.nwk -c 1 -v 0 -o lr
TN93_BICi: 12134.51460
18:54:18 # running: phyml -i primates.phy -d nt -m TN93 -f m -c 4 -a e -u primates.phy_JC-NJ.nwk -o lr
18:54:18 # running: phyml -i primates.phy -d nt -m TN93 -f m -v e -c 1 -u primates.phy_JC-NJ.nwk -o lr
18:54:18 # running: phyml -i primates.phy -d nt -m TN93 -f m -u primates.phy_JC-NJ.nwk -v e -a e -o lr
--------------------------------------------------------------------------------
18:54:20 # Starting evaluation and automatic selection of extended model set
# setting compositional_heterogeneity flag to: 1
--------------------------------------------------------------------------------
18:54:20 # ... evaluation and automatic selection of extended model set
# setting transitional_heterogeneity flag to: 0
--------------------------------------------------------------------------------
18:54:20 # ... evaluation and automatic selection of extended model set
# setting use_pInv flag to: 0
# will evaluate models with unequal frequencies
--------------------------------------------------------------------------------
Como ven arriba, el \(script\) detecta:
- un sesgo composicional significativo,
compositional_heterogeneity_flag=1 - No considera significativa la diferencia de tasas de
transiciones entre purinas y pirimidinas
transitional_heterogeneity_flag=0 - Tampoco considera que el parámetro \(pInv\) incremente significativamente el
ajuste, por lo que
use_pInv_flag=0
En el siguiente ciclo, el programa evaluará sólo modelos \(TVM*\) del set extendido, con estima empírica de la frecuencia de las bases, sin \(pInv\)
18:54:20 2.2. running a loop to combine the extended models in model_set 2=>extended_models,
with (or not) +G and or +I, and compute the model lnL, after optimizing branch lengths and rates
skipping TN|TIM|SYM matrix 012230
skipping TN|TIM|SYM matrix 010023
skipping TN|TIM|SYM matrix 010232
skipping TN|TIM|SYM matrix 012232
skipping TN|TIM|SYM matrix 012332
skipping TN|TIM|SYM matrix 012342
skipping TN|TIM|SYM matrix 012343
skipping TN|TIM|SYM matrix 012340
skipping TN|TIM|SYM matrix 010021
skipping TN|TIM|SYM matrix 010022
skipping TN|TIM|SYM matrix 010120
skipping TN|TIM|SYM matrix 011123
skipping TN|TIM|SYM matrix 012223
skipping TN|TIM|SYM matrix 012222
skipping TN|TIM|SYM matrix 000120
skipping TN|TIM|SYM matrix 000121
skipping TN|TIM|SYM matrix 001021
skipping TN|TIM|SYM matrix 012234
skipping TN|TIM|SYM matrix 010231
skipping TN|TIM|SYM matrix 011020
skipping TN|TIM|SYM matrix 011230
skipping TN|TIM|SYM matrix 012130
skipping TN|TIM|SYM matrix 010121
skipping TN|TIM|SYM matrix 010122
skipping TN|TIM|SYM matrix 010123
skipping TN|TIM|SYM matrix 001020
skipping TN|TIM|SYM matrix 000123
skipping TN|TIM|SYM matrix 010203
skipping TN|TIM|SYM matrix 010223
skipping TN|TIM|SYM matrix 010230
skipping TN|TIM|SYM matrix 012032
skipping TN|TIM|SYM matrix 010233
skipping TN|TIM|SYM matrix 001234
skipping TN|TIM|SYM matrix 010234
skipping TN|TIM|SYM matrix 011234
skipping TN|TIM|SYM matrix 012134
skipping TN|TIM|SYM matrix 012304
skipping TN|TIM|SYM matrix 012324
skipping TN|TIM|SYM matrix 012334
skipping TN|TIM|SYM matrix 012341
skipping TN|TIM|SYM matrix 012344
18:54:20 # running: phyml -i primates.phy -d nt -m 012210 -f m -u primates.phy_JC-NJ.nwk -c 1 -v 0 -o lr
18:54:20 # running: phyml -i primates.phy -d nt -m 012210 -f m -c 4 -a e -u primates.phy_JC-NJ.nwk -o lr
skipping 012210+I and 012210+I+G
18:54:20 # running: phyml -i primates.phy -d nt -m 012314 -f m -u primates.phy_JC-NJ.nwk -c 1 -v 0 -o lr
18:54:20 # running: phyml -i primates.phy -d nt -m 012314 -f m -c 4 -a e -u primates.phy_JC-NJ.nwk -o lr
skipping 012314+I and 012314+I+G
18:54:21 # running: phyml -i primates.phy -d nt -m 012310 -f m -u primates.phy_JC-NJ.nwk -c 1 -v 0 -o lr
18:54:21 # running: phyml -i primates.phy -d nt -m 012310 -f m -c 4 -a e -u primates.phy_JC-NJ.nwk -o lr
skipping 012310+I and 012310+I+G
18:54:21 # running: phyml -i primates.phy -d nt -m 010213 -f m -u primates.phy_JC-NJ.nwk -c 1 -v 0 -o lr
18:54:21 # running: phyml -i primates.phy -d nt -m 010213 -f m -c 4 -a e -u primates.phy_JC-NJ.nwk -o lr
skipping 010213+I and 010213+I+G
18:54:22 # running: phyml -i primates.phy -d nt -m 012213 -f m -u primates.phy_JC-NJ.nwk -c 1 -v 0 -o lr
18:54:22 # running: phyml -i primates.phy -d nt -m 012213 -f m -c 4 -a e -u primates.phy_JC-NJ.nwk -o lr
skipping 012213+I and 012213+I+G
18:54:22 # running: phyml -i primates.phy -d nt -m 012013 -f m -u primates.phy_JC-NJ.nwk -c 1 -v 0 -o lr
18:54:22 # running: phyml -i primates.phy -d nt -m 012013 -f m -c 4 -a e -u primates.phy_JC-NJ.nwk -o lr
skipping 012013+I and 012013+I+G
18:54:23 # running: phyml -i primates.phy -d nt -m 010012 -f m -u primates.phy_JC-NJ.nwk -c 1 -v 0 -o lr
18:54:23 # running: phyml -i primates.phy -d nt -m 010012 -f m -c 4 -a e -u primates.phy_JC-NJ.nwk -o lr
skipping 010012+I and 010012+I+G
18:54:23 # running: phyml -i primates.phy -d nt -m 012012 -f m -u primates.phy_JC-NJ.nwk -c 1 -v 0 -o lr
18:54:23 # running: phyml -i primates.phy -d nt -m 012012 -f m -c 4 -a e -u primates.phy_JC-NJ.nwk -o lr
skipping 012012+I and 012012+I+G
18:54:24 # running: phyml -i primates.phy -d nt -m 010212 -f m -u primates.phy_JC-NJ.nwk -c 1 -v 0 -o lr
18:54:24 # running: phyml -i primates.phy -d nt -m 010212 -f m -c 4 -a e -u primates.phy_JC-NJ.nwk -o lr
skipping 010212+I and 010212+I+G
18:54:24 # running: phyml -i primates.phy -d nt -m 012313 -f m -u primates.phy_JC-NJ.nwk -c 1 -v 0 -o lr
18:54:24 # running: phyml -i primates.phy -d nt -m 012313 -f m -c 4 -a e -u primates.phy_JC-NJ.nwk -o lr
skipping 012313+I and 012313+I+G
18:54:24 # running: phyml -i primates.phy -d nt -m 010011 -f m -u primates.phy_JC-NJ.nwk -c 1 -v 0 -o lr
18:54:24 # running: phyml -i primates.phy -d nt -m 010011 -f m -c 4 -a e -u primates.phy_JC-NJ.nwk -o lr
skipping 010011+I and 010011+I+G
18:54:25 # running: phyml -i primates.phy -d nt -m 012212 -f m -u primates.phy_JC-NJ.nwk -c 1 -v 0 -o lr
18:54:25 # running: phyml -i primates.phy -d nt -m 012212 -f m -c 4 -a e -u primates.phy_JC-NJ.nwk -o lr
skipping 012212+I and 012212+I+G
18:54:25 # running: phyml -i primates.phy -d nt -m 010210 -f m -u primates.phy_JC-NJ.nwk -c 1 -v 0 -o lr
18:54:25 # running: phyml -i primates.phy -d nt -m 010210 -f m -c 4 -a e -u primates.phy_JC-NJ.nwk -o lr
skipping 010210+I and 010210+I+G
18:54:25 # running: phyml -i primates.phy -d nt -m 011010 -f m -u primates.phy_JC-NJ.nwk -c 1 -v 0 -o lr
18:54:26 # running: phyml -i primates.phy -d nt -m 011010 -f m -c 4 -a e -u primates.phy_JC-NJ.nwk -o lr
skipping 011010+I and 011010+I+G
18:54:26 # running: phyml -i primates.phy -d nt -m 001101 -f m -u primates.phy_JC-NJ.nwk -c 1 -v 0 -o lr
18:54:26 # running: phyml -i primates.phy -d nt -m 001101 -f m -c 4 -a e -u primates.phy_JC-NJ.nwk -o lr
skipping 001101+I and 001101+I+G
18:54:26 # running: phyml -i primates.phy -d nt -m 000100 -f m -u primates.phy_JC-NJ.nwk -c 1 -v 0 -o lr
18:54:26 # running: phyml -i primates.phy -d nt -m 000100 -f m -c 4 -a e -u primates.phy_JC-NJ.nwk -o lr
skipping 000100+I and 000100+I+G
18:54:26 # running: phyml -i primates.phy -d nt -m 000101 -f m -u primates.phy_JC-NJ.nwk -c 1 -v 0 -o lr
18:54:26 # running: phyml -i primates.phy -d nt -m 000101 -f m -c 4 -a e -u primates.phy_JC-NJ.nwk -o lr
skipping 000101+I and 000101+I+G
18:54:26 # running: phyml -i primates.phy -d nt -m 010110 -f m -u primates.phy_JC-NJ.nwk -c 1 -v 0 -o lr
18:54:27 # running: phyml -i primates.phy -d nt -m 010110 -f m -c 4 -a e -u primates.phy_JC-NJ.nwk -o lr
skipping 010110+I and 010110+I+G
18:54:27 # running: phyml -i primates.phy -d nt -m 000102 -f m -u primates.phy_JC-NJ.nwk -c 1 -v 0 -o lr
18:54:27 # running: phyml -i primates.phy -d nt -m 000102 -f m -c 4 -a e -u primates.phy_JC-NJ.nwk -o lr
skipping 000102+I and 000102+I+G
18:54:27 # running: phyml -i primates.phy -d nt -m 010112 -f m -u primates.phy_JC-NJ.nwk -c 1 -v 0 -o lr
18:54:27 # running: phyml -i primates.phy -d nt -m 010112 -f m -c 4 -a e -u primates.phy_JC-NJ.nwk -o lr
skipping 010112+I and 010112+I+G
18:54:27 # running: phyml -i primates.phy -d nt -m 010211 -f m -u primates.phy_JC-NJ.nwk -c 1 -v 0 -o lr
18:54:27 # running: phyml -i primates.phy -d nt -m 010211 -f m -c 4 -a e -u primates.phy_JC-NJ.nwk -o lr
skipping 010211+I and 010211+I+G
--------------------------------------------------------------------------------
Como puede verse en la salida del programa arriba mostrada, ajusta modelos \(TVM*+G\), estimando frecuencia de bases \(-f\ e\), e ignorando los modelos \(Tn|TIM|SYM\), con dos tasas de transiciones.
Una vez evaluados los dos conjuntos de modelos (base + selección de modelos extendidos), calcula diversos estadísticos de cirterios de información, tanto de Akaike como bayesiano, ordenando los resultados por valores crecientes de \(BIC\), e imprimiendo también los valores de \(\Delta_{BIC}\), \(BIC_{w_i}\) y \(BICcumW\).
18:54:28 # writing primates.phy_sorted_model_set_2_fits.tsv, sorted by BIC; extmodels+F true -------------------------------------------------------------------------------- model K sites/K lnL AIC AICc BIC deltaBIC BICw BICcumW HKY85+G 26 34.53000 -5711.93899 11475.87798 11477.48992 11600.68240 0 0.85 0.85 TN93+G 27 33.25000 -5710.55176 11475.10352 11476.84145 11604.70811 4.02571 0.11 0.97 HKY85+I+G 27 33.25000 -5711.94030 11477.88060 11479.61853 11607.48519 6.80279 0.03 1.00 TN93+I+G 28 32.07000 -5710.55229 11477.10458 11478.97339 11611.50934 10.82694 0.00 1.00 010212+G+F 27 33.25000 -5721.62318 11497.24636 11498.98429 11626.85095 26.16855 0.00 1.00 012212+G+F 27 33.25000 -5722.51008 11499.02016 11500.75809 11628.62475 27.94235 0.00 1.00 010012+G+F 27 33.25000 -5722.63081 11499.26162 11500.99955 11628.86621 28.18381 0.00 1.00 012313+G+F 28 32.07000 -5719.92282 11495.84564 11497.71445 11630.25040 29.56800 0.00 1.00 012213+G+F 28 32.07000 -5720.26837 11496.53674 11498.40555 11630.94150 30.25910 0.00 1.00 010213+G+F 28 32.07000 -5720.86864 11497.73728 11499.60609 11632.14204 31.45964 0.00 1.00 012013+G+F 28 32.07000 -5722.34471 11500.68942 11502.55823 11635.09418 34.41178 0.00 1.00 010210+G+F 27 33.25000 -5725.84777 11505.69554 11507.43347 11635.30013 34.61773 0.00 1.00 012012+G+F 27 33.25000 -5726.25000 11506.50000 11508.23793 11636.10459 35.42219 0.00 1.00 012314+G+F 29 30.96000 -5719.45652 11496.91304 11498.91765 11636.11797 35.43557 0.00 1.00 012210+G+F 27 33.25000 -5726.65816 11507.31632 11509.05425 11636.92091 36.23851 0.00 1.00 012310+G+F 28 32.07000 -5725.20152 11506.40304 11508.27185 11640.80780 40.12540 0.00 1.00 HKY85+I 26 34.53000 -5768.00399 11588.00798 11589.61992 11712.81240 112.13000 0.00 1.00 TN93+I 27 33.25000 -5766.10222 11586.20444 11587.94237 11715.80903 115.12663 0.00 1.00 010211+G+F 27 33.25000 -5804.50691 11663.01382 11664.75175 11792.61841 191.93601 0.00 1.00 010011+G+F 26 34.53000 -5810.83201 11673.66402 11675.27596 11798.46844 197.78604 0.00 1.00 010112+G+F 27 33.25000 -5816.57851 11687.15702 11688.89495 11816.76161 216.07921 0.00 1.00 010110+G+F 26 34.53000 -5823.31701 11698.63402 11700.24596 11823.43844 222.75604 0.00 1.00 001101+G+F 26 34.53000 -5913.85413 11879.70826 11881.32020 12004.51268 403.83028 0.00 1.00 011010+G+F 26 34.53000 -5917.56175 11887.12350 11888.73544 12011.92792 411.24552 0.00 1.00 K80+G 23 39.04000 -5950.87340 11947.74680 11949.00996 12058.15071 457.46831 0.00 1.00 K80+I+G 24 37.41000 -5950.57605 11949.15210 11950.52667 12064.35618 463.67378 0.00 1.00 012313+F 27 33.25000 -5950.31559 11954.63118 11956.36911 12084.23577 483.55337 0.00 1.00 012213+F 27 33.25000 -5952.11058 11958.22116 11959.95909 12087.82575 487.14335 0.00 1.00 012314+F 28 32.07000 -5949.09947 11954.19894 11956.06775 12088.60370 487.92130 0.00 1.00 010212+F 26 34.53000 -5959.18463 11970.36926 11971.98120 12095.17368 494.49128 0.00 1.00 012212+F 26 34.53000 -5959.65943 11971.31886 11972.93080 12096.12328 495.44088 0.00 1.00 010213+F 27 33.25000 -5957.06064 11968.12128 11969.85921 12097.72587 497.04347 0.00 1.00 K80+I 23 39.04000 -5974.90514 11995.81028 11997.07344 12106.21419 505.53179 0.00 1.00 010012+F 26 34.53000 -5964.74421 11981.48842 11983.10036 12106.29284 505.61044 0.00 1.00 012013+F 27 33.25000 -5962.12925 11978.25850 11979.99643 12107.86309 507.18069 0.00 1.00 012012+F 26 34.53000 -5974.64960 12001.29920 12002.91114 12126.10362 525.42122 0.00 1.00 010210+F 26 34.53000 -5976.57458 12005.14916 12006.76110 12129.95358 529.27118 0.00 1.00 012310+F 27 33.25000 -5973.19079 12000.38158 12002.11951 12129.98617 529.30377 0.00 1.00 012210+F 26 34.53000 -5977.70123 12007.40246 12009.01440 12132.20688 531.52448 0.00 1.00 HKY85 25 35.92000 -5981.72025 12013.44050 12014.93133 12133.44475 532.76235 0.00 1.00 TN93 26 34.53000 -5978.85509 12009.71018 12011.32212 12134.51460 533.83220 0.00 1.00 000101+G+F 26 34.53000 -6028.02955 12108.05910 12109.67104 12232.86352 632.18112 0.00 1.00 000102+G+F 27 33.25000 -6027.15138 12108.30276 12110.04069 12237.90735 637.22495 0.00 1.00 010112+F 26 34.53000 -6065.15408 12182.30816 12183.92010 12307.11258 706.43018 0.00 1.00 000100+G+F 26 34.53000 -6079.11681 12210.23362 12211.84556 12335.03804 734.35564 0.00 1.00 010211+F 26 34.53000 -6080.81607 12213.63214 12215.24408 12338.43656 737.75416 0.00 1.00 001101+F 25 35.92000 -6089.44631 12228.89262 12230.38345 12348.89687 748.21447 0.00 1.00 010110+F 25 35.92000 -6090.63982 12231.27964 12232.77047 12351.28389 750.60149 0.00 1.00 010011+F 25 35.92000 -6094.74661 12239.49322 12240.98405 12359.49747 758.81507 0.00 1.00 F81+G 25 35.92000 -6106.18486 12262.36972 12263.86055 12382.37397 781.69157 0.00 1.00 F81+I+G 26 34.53000 -6103.95758 12259.91516 12261.52710 12384.71958 784.03718 0.00 1.00 F81+I 25 35.92000 -6117.43632 12284.87264 12286.36347 12404.87689 804.19449 0.00 1.00 K80 22 40.81000 -6142.42909 12328.85818 12330.01475 12434.46192 833.77952 0.00 1.00 011010+F 25 35.92000 -6167.29688 12384.59376 12386.08459 12504.59801 903.91561 0.00 1.00 000101+F 25 35.92000 -6198.63537 12447.27074 12448.76157 12567.27499 966.59259 0.00 1.00 000102+F 26 34.53000 -6197.84929 12447.69858 12449.31052 12572.50300 971.82060 0.00 1.00 000100+F 25 35.92000 -6259.08395 12568.16790 12569.65873 12688.17215 1087.48975 0.00 1.00 JC69+G 22 40.81000 -6272.46906 12588.93812 12590.09469 12694.54186 1093.85946 0.00 1.00 JC69+I+G 23 39.04000 -6269.14954 12584.29908 12585.56224 12694.70299 1094.02059 0.00 1.00 JC69+I 22 40.81000 -6277.78094 12599.56188 12600.71845 12705.16562 1104.48322 0.00 1.00 F81 24 37.41000 -6284.99563 12617.99126 12619.36583 12733.19534 1132.51294 0.00 1.00 JC69 21 42.76000 -6424.20245 12890.40490 12891.45969 12991.20847 1390.52607 0.00 1.00 -------------------------------------------------------------------------------------------------- * NOTE 1: when sites/K < 40, the AICc is recommended over AIC * NOTE 2: The best model is selected by BIC, because AIC is biased, favoring parameter-rich models
En la columna de \(BIC_{w_i}\) se aprecia que:
- El modelo de dos tasas HKY85+G (010010+G) es el mejor, con un peso \(BIC_w\ =\ 0.85\).
- El segundo mejor (\(TN93+G\)) está a una \(\Delta_{BIC}\ =\ 4.026\) con respecto al \(BIC_{min}\), con un peso relativo \(BIC_w\ =\ 0.11\).
- Recuerden que el modelo seleccionado por \(AIC\) fue el \(TN93G\), demostrando que tiende a seleccionar modelos más parametrizados que el \(BIC\).
Finalmente, se infiere la filogenia bajo el mejor modelo seleccionado por BIC, en este caso \(HKY85+G\)
==================================================================================================
# Will estimate the ML tree under best-fitting model HKY85+G selected by BIC
--------------------------------------------------------------------------------------------------
18:54:31 # running: phyml -i primates.phy -d nt -m HKY85 -c 4 -a e -o tlr -s BEST
# Your results:
- primates.phy_HKY85+G_BESTmoves_phyml_stats.txt
- primates.phy_HKY85+G_BESTmoves_phyml_tree.txt
--------------------------------------------------------------------------------------------------
Elapsed time: 0 days, 00 hr, 00 min, 18 sec
Done!
========================================================================================
If you use phyml_DNAmodelFinder.sh v.1.2.0_2023-11-18 for your research,
I would appreciate that you:
1. Cite the code in your work as:
Pablo Vinuesa. phyml_DNAmodelFinder.sh v.1.2.0_2023-11-18
https://github.com/vinuesa/TIB-filoinfo/blob/master/phyml_DNAmodelFinder.sh
2. Give it a like on the https://github.com/vinuesa/TIB-filoinfo/ repo
Thanks!
Podemos desplegar el árbol directamente en la línea de comandos con \(nw\_display\), enraizando previamente en Tarsius syrichta y Lemur catta con \(nw\_reroot\), como se muestra abajo:
+----------------+ Tarsius sy
|
=| +-----------------------+ Lemur catt
| |
| | +-----------------------------+ Saimiri sc
+----------------+ |
| | +---+ M. sylvanu
| | |
+------------------+ 1.000000---------------+ 1.000000 fascicu
| | | |
| | +-+ 0.987818ulatta
| | +--+ 0.999992
+------+ 0.992860 ++ Macaca fus
|
| +----------+ Hylobates
| |
+--------+ 1.000000--+ Pongo
| |
+--+ 0.994098 Gorilla
| |
+----+ 1.000000n
+-+ 0.999171
+--+ Homo sapie
|---------------|----------------|---------------|---------------|--
0 0.25 0.5 0.75 1
substitutions/site
Figura 7. Filogenia de máxima verosimilitud bajo el model \(HKY85+G\) seleccionado automáticamente por phyml_DNAmodelFinder mediante \(BIC\) y desplegada con nw_display.
6 Comparación del rendimiento (benchmarking) de phyml_DNAmodelFinder con jModelTest2
Para finalizar, vamos a comparar los resultados y rendimiento computacional de phyml_DNAmodelFinder con respecto a jModelTesst2, que puede ser considerado el software de referencia para selección de modelos de DNA.
Las corridas de ambos programas se realizaron sobre la misma máquina con las siguientes características:
- Linux 5.15.0-56-generic, arch: amd64, bits: 64, numcores: 12 - Intel Corporation 8th Gen Core Processor Host Bridge/DRAM Registers (rev 07)
6.1 Rendimiento de phyml_DNAmodelFinder
- Llamada a phyml_DNAmodelFinder v1.2.0_2023-11-18 corriendo phyml v.3.3_2023 (single-threaded) sobre un solo núcleo.
real 0m17.044s user 0m15.052s sys 0m2.221s
con los siguientes modelos seleccionados por \(BIC\) conformando el intervalo de credibilidad del 100%
model K sites/K lnL AIC AICc BIC deltaBIC BICw BICcumW HKY85+G 26 34.53000 -5711.93896 11475.87792 11477.48986 11600.68234 0 0.85 0.85 TN93+G 27 33.25000 -5710.55184 11475.10368 11476.84161 11604.70827 4.02593 0.11 0.97 HKY85+I+G 27 33.25000 -5711.94038 11477.88076 11479.61869 11607.48535 6.80301 0.03 1.00 TN93+I+G 28 32.07000 -5710.55238 11477.10476 11478.97357 11611.50952 10.82718 0.00 1.00 010212+G+F 27 33.25000 -5721.62321 11497.24642 11498.98435 11626.85101 26.16867 0.00 1.00 012212+G+F 27 33.25000 -5722.51024 11499.02048 11500.75841 11628.62507 27.94273 0.00 1.00 ... TRUNCADO
Nota: todos los modelos arriba enlistados como códigos incorporan la estima de frecuencia de bases, como indica ‘la decoración’ \(+F\), (compositional_heterogeneity flag to: 1).
6.2 Rendimiento de jModelTest2
La versión de \(jmodeltest2\) empleada es la que viene pre-empacada para Ubuntu 22.04.1 LTS
- Llamada a jModelTest2 v2.1.10 v20160303 corriendo phyml-mpi PhyML_3.0_linux64 v20130103 sobre 12 cores asignados automáticamente
real 1m30.086s user 13m27.152s sys 0m1.791s
con los siguientes modelos seleccionados por \(BIC\) conformando el intervalo de credibilidad del 100%
* BIC MODEL SELECTION : Selection uncertainty Model -lnL K BIC delta weight cumWeight ------------------------------------------------------------------------- HKY+G 5711.93890 26 11600.682222 0.000000 0.708136 0.708136 010012+G+F 5710.19575 27 11603.996092 3.313870 0.135058 0.843194 TrN+G 5710.55146 27 11604.707512 4.025290 0.094632 0.937825 HKY+I+G 5711.93893 27 11607.482452 6.800230 0.023630 0.961455 010023+G+F 5709.27633 28 11608.957422 8.275200 0.011303 0.972758 012232+G+F 5709.82603 28 11610.056822 9.374600 0.006523 0.979281 012213+G+F 5709.89075 28 11610.186262 9.504040 0.006114 0.985395 010012+I+G+F 5710.19559 28 11610.795942 10.113720 0.004508 0.989903 TIM1+G 5710.43657 28 11611.277902 10.595680 0.003542 0.993445 TrN+I+G 5710.55181 28 11611.508382 10.826160 0.003157 0.996602 ... truncado, basado en redondeo de 0.996602 ~= 1
6.3 Benchmark entre \(phyml\_DNAmodelFinder.sh\) y \(jModelTest2\)
Un análisis comparativo, basado en 75 alineamientos con 7-30 secuencias entre 132 y 1234 nucleótidos, resultó en los siguientes resultados:
- jmodeltest2_12_threads-time.out
real 57m38.413s user 526m10.365s sys 2m39.542s
- phyml_DNAmodelFinder_v08_phyml22_PARALLEL_BICcutoff2_fm-time.out
real 4m45.518s user 45m9.525s sys 0m55.629s
En promedio, ambos programas encuentran modelos muy parecidos para los 75 alineamientos, como se deduce de:
\(\sum_{n=1}^{75}BIC_{phyml\_DNAmodelFinder} - \sum_{n=1}^{75}BICcum_{jmodeltest2} = -0.411463\),
resultado que favorece mínimamente a \(phyml\_DNAmodelFinder\).
Que ambos programas seleccionan modelos idénticos o muy parecidos se deduce del hecho de que los 75 modelos seleccionados por cada uno de los programas tiene el mismo promedio de 43 parámetros:
(\(\overline{K}_{phyml\_DNAmodelFinder} = \overline{K}_{jModelTest2} = 43\)).
6.4 Conclusión
El \(script\) \(phyml\_DNAmodelFinder\) es una herramienta muy ligera, eficiente y competitiva para encontrar modelos de sustitución de nucleótidos de entre los 203 modelos de la familia \(GTR+G+I\), haciendo uso exclusivamente de \(phyml\), \(Bash\) y \(awk\), junto a algunas otras herramientas estándar de \(Linux\), como \(bc\), \(sed\) y \(sort\).
Es de hecho una herramienta más rápida que \(jModelTest2\) (Vinuesa, benchmark sin publicar), seleccionando el mismo modelo o modelos muy parecidos.
\(phyml\_DNAmodelFinder.sh\) usa exclusivamente el criterio bayesiano de información (\(BIC\)), lo cual minimiza el riesgo de sobre-parametrización.
6.5 Impacto de \(modelTest\) y su sucesor \(jModelTest2\) en el campo de la filogenética
\(modelTest\), y especialmente su sucesor \(jModelTest2\), sin duda alguna representan un hito en la popularización del uso riguroso de selección de modelos para filogenética.
\(jModelTest2\) es una herramienta muy versátil y poderosa, que implementa cinco estrategias para la selección de modelos:
- LRTs jerárquicos(\(hLRT\)) y dinámicos (\(dLRT\))
- criterios de información de Akaike (\(AIC\)) y bayesianos (\(BIC\))
- un método de decisión teórica (\(DT\)).
Además calcula ponderaciones o pesos de Akaike o bayesianos de los parámetros de todos los modelos evaluados, y hace promediado de modelos y de sus parámetros, ¡una joya!
El trabajo del Dr. David Posada y su grupo sobre selección de modelos filogenéticos ha tenido un enorme impacto en el campo, y recomiendo enfáticamente que lean sus trabajos. ¡Gracias David!.
No obstante, algunos autores son también críticos con respecto al uso generalizado de \(modeltest\). En palabras del Prof. Ziheng Yang, la aplicación masiva de \(jModelTest2\) a datos biológicos reales ha contribuido al uso extendido de los modelos más complejos, incluyendo los “modelos patológicos” \(+I+\Gamma\), como describe en su influyente libro (Yang Z. 2014. Molecular Evolution - a statistical approach, pp. 149), disponible en línea.
7 Cita para este tutorial y los scripts phyml_DNAmodelFinder y phyml_protModelFinder
Si este material te es de utilidad, agradeceré que lo cites en tu tesis o publicaciones, y/o le des un “like” en el repositorio GitHub de los TIB.
- Citas: