Tutoral de uso de PhyML v3 para selección de modelos y búsquedas intensivas de árboles

1 Selección y análisis de modelos paramétricos y empíricos mediante LRTs, AIC y BIC con phyml, programación AWK, Bash y R

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/bin

2.1 Publicaciones clave sobre PhyML

1: Lefort V, Longueville JE, Gascuel O. SMS: Smart Model Selection in PhyML. Mol Biol Evol. 2017 Sep 1;34(9):2422-2424. doi: 10.1093/molbev/msx149. PubMed PMID: 28472384; PubMed Central PMCID: PMC5850602.

2: 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.

3: 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.

4: 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.

# para ver la ayuda
phyml -help

# -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 alineamientos de 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/null

3 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.tgz

y confirma que tienes el tarro comprimido (tgz):

ls -ltr

-rw-rw-r-- 1 vinuesa vinuesa 145254 nov 21 19:36 phyml_tutorial_data.tgz

file 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:

tar -xzf phyml_tutorial_data.tgz

---

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:

1. Calcular un árbol de distancias NJ, bajo el modelo más sencillo (\(JC\)), lo cual es muy rápido
2. 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

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.

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

# 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 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:

1. frecuencia de bases
2. número de clases de tasas de sustitución
3. 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, igual para cada tipo o clase 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 $K2P$ o $HKY85$
- 012345 representa un modelo de 6 tasas, es decir, el GTR.

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 las máximas restricciones.

# 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 una \(shell\) de \(unix\) o \(linux\), podemos fácilmente acceder a este valor sin necesidad de abrir el archivo:

# o usando grep:
grep Log-like primates.phy_phyml_stats.txt 

. 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.

grep Log-like primates.phy_phyml_stats.txt > JC.fit

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.

awk '/Log-like/{print $NF}' primates.phy_phyml_stats.txt

-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.fit

3.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. Modelo F81
# phyml -i primates.phy -d nt -m 000000 -f m -c 1 -u primates_JC.ph -o lr -b 0 &> /dev/null
phyml -i primates.phy -d nt -m F81 -f m -c 1 -u primates_JC.ph -o lr -b 0 &> /dev/null

grep Log- primates.phy_phyml_stats.txt > F81.fit

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.fit

NOTA: 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 cada vez más compleja.

for m in JC69 K80 F81 HKY85; do 
  phyml -i primates.phy -d nt -m $m -f m -u primates_JC.ph -c 1 -o lr -b 0 &> /dev/null
  awk '/Log-like/{print $NF}' primates.phy_phyml_stats.txt > ${m}.fit
done

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.

# 5. Modelo HKY85+G
# phyml -i primates.phy -d nt -m 121121 -f m -t e -c 4 -a 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 -u primates_JC.ph -o lr -b 0 &> /dev/null
grep Log-like primates.phy_phyml_stats.txt > HKY85G.fit

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.

# 5. Modelo HKY85+G
# phyml -i primates.phy -d nt -m 121121 -f m -t e -c 1 -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 1 -v e -u primates_JC.ph -o lr -b 0 &> /dev/null
grep Log-like primates.phy_phyml_stats.txt > HKY85I.fit

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 una proporción de sitios invariantes.

# 6. 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.fit

Nota: 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 resultdos 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

sed '1d' model_fits.tsv | sort -grk2 | column -t

HKY85G   -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 tiene 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 abrumadora.

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

eog LRTs_primates_dotchart.png &

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:

1. 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\).
   
2. El estadístico \(LRT=2*(lnL_1 - lnL_0)\) sigue una distribución de chi-cuadrada:

\[\chi^2 = \sum \frac {(O - E)^2}{E}\]

# 6.3 Entre los dos modelos que parecen empatados, veamos si existen o no diferencias significativas, 
#     corriendo una prueba estadistica formal: 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 ALTAMENTE SIGNIFICATIVO.

# cómputo de estadísticas
pchisq(xsq,1)      # 1;  este es el valor critico de la X2 para el estadistico LRT con 1 grado de libertad; [1] 1
1-(pchisq(xsq,1))  # 0; esta es la p asociada a dicho valor crítico, ABSOLUTAMENTE SIGNIFICATIVA =>
                   # rechazamos el modelo nulo                                                       

q()         # para salir de R
            # decimos que NO quieren salvar sesion

En 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?

phyml -i primates.phy -d nt -m HKY85 -f m -t e -c 4 -a e -o tlr

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 al 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.

# despliega el árbol: 
figtree primates.phy_phyml_tree.txt &

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 de 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.

# GTR+G+I
phyml -i primates.phy -m GTR -a e -f m -v e -u primates_JC.ph -b 0 -o lr 
grep Log-like primates.phy_phyml_stats.txt > GTRGI.fit  

# visualizamos las tasas del modelo GTR+G+I
grep '<->' primates.phy_phyml_stats.txt

# 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

3.3.2 Evaluación de modelos simplificados de 2, 3 y 4 tasas

En base a los resultados arriba mostrados, vamos a evaluar modelos alternativos simplificados.

La simplificación más obvia es igualar ambas tasas de transiciones, ya que son muy parecidas entre sí y alrededor de un órden de magnitud más altas que las transversiones.

Evaluaremos otras opciones de simplificación o reducción de número de 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

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 -t

GTRGI          -5708.92372
TN93G          -5710.55182
HKYG           -5711.93899

Noten que si bien el modelo de más alta parametrización (\(GTR+G+I\)) tiene el puntaje más alto de \(lnL\), la diferencia contra el competidor \(HKY+G\) es pequeña, < 1-2 unidades de \(lnL\).

Usemos el criterio de información de Akaike (AIC) para evaluar estos resultados formalmente.

3.3.3 Cálculo de los \(AICi\) (de cada modelo) y \(\Delta_i\) 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 dataframe model_AICs, con los vectores "model" y "AIC"
# vector model con los nombres de los modelos
model <- c("AIC_GTRGI", "AIC_TN93T", "AIC_HKYG")

AIC_GTRGI <- -2*(-5708.92372) + 2*(10+nbranches)
AIC_TN93G <- -2*(-5710.55182) + 2*(7+nbranches)
AIC_HKYG <- -2*(-5711.93898) + 2*(6+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)
AICmin

Explorando el dataframe \(model\_AICs\) ordenado por valores crecientes de \(AICi\), vemos que el \(AICmin\) corresponde al modelo AIC_5ratesG

model_AICs[order(model_AICs$AIC),]
      model      AIC
2 AIC_TN93G 11471.10
3  AIC_HKYG 11471.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 2.74380

3.3.4 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.77\).
• El parámetro de \(pInv\) es perfectamente prescindible, y un modelo de 2 o tres 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\ =\ 2.74\)), indicando que está claramente sobre-parametrizado.

En mi experiencia, esta estrategia es buena por dos motivos:

1. es muy rápida (posiblemente la más rápida)
2. 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, habría que escribir un \(script\) para automatizarla.

3.3.5 Tarea de ajuste de modelos paramétricos mediante criterios de información AIC y BIC

1. Evalúa mediante AICs las otras simplificaciones en las clases de transversiones propuestas en el esquema de tasas del modelo \(GTR+I+G\).
2. 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.
3. Como veremos en la siguiente sección, el criterio de información de Akaike \(AIC\) tiende a seleccionar modelos más parametrizados que el citerio 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

Recuerden 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 \(AIC\) o \(BIC\) para evaluar su mérito o ajuste relativo.

En el siguiente ejercicio evaluaremos primeramente una gama de 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 con un \(bucle\ for\)

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)
done

4.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.csv

JTT+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?

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\)

eog dotchart_sorted_model_scores_for_primates_21_AA.png &

Figura 4. Puntaciones de máxima verosimilitud para los modelos evualuados 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 5

4.3.1 Despliegue y edición del árbol con figtree

# Despliegue y edición del árbol con figtree 
figtree primates_21_AA.phy_phyml_tree.txt &

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

cat primates_21_AA.phy_phyml_stats.txt

...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 por \(BIC\) con phyml_DNAmodelFinder.sh y phyml_protModelFinder.sh

Si bien en el ejemplo anterior hemos intentado optimizar nuestro 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.

Los pasos seguidos son fácilmente automatizables mediante un \(script\) en \(Bash\) u otro lenguaje. 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 facilidad.

Vamos a cerrar este tutorial corriendo los scripts phyml_protModelFinder.sh y phyml_DNAmodelFinder.sh, que pueden copiar a su directorio ~/bin con los siguientes comandos:

[ ! -d ~/bin ] && mkdir ~/bin

cp /home/vinuesa/cursos/LCG-ENS/scripts/phyml_protModelFinder.sh ~/bin
cp /home/vinuesa/cursos/LCG-ENS/scripts/phyml_DNAmodelFinder.sh ~/bin

o descargar directamente desde el repositorio GitHub de los TIB, usando \(wget\) desde la línea de comandos:

# si no existe el directorio ~/bin, entonces generalo
[ ! -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úa 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.

Te invito a que explores en detalle el código de ambos.

NOTA: Los scripts phyml_protModelFinder.sh y en particular phyml_DNAmodelFinder.sh hacen uso de sintaxis moderna de Bash, por lo que requiren 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.

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 correspondiente.

# despliegue de la ayuda
phyml_protModelFinder.sh

phyml_protModelFinder.sh v0.6 requires two arguments provided on the command line; the third one is optional

phyml_protModelFinder.sh [input phylip file of aligned PROTEIN sequences] [model sets:1-7] [num random seed trees; default:1]
 
   # 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.6 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

LICENSE: GPL v3.0. See https://github.com/vinuesa/get_phylomarkers/blob/master/LICENSE

Vemos que necesitamos pasarle dos parámetros: el archivo phylip con las secuencias alineadas y elejir 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 de modelos para genes nucleares

phyml_protModelFinder.sh primates_21_AA.phy 5

Abajo vemos el STDOUT del \(script\).

• Antes de inciar, comprueba que estamos corriendo una versión de \(bash >= 4.4\), ya que el \(script\) hace uso de idiomas modernos de \(bash\), incluyendo hashes. De hecho algunas funcionalidades requieren \(bash\ >=\ 4.8\). Si usas MacOS, probablemente debas hacer un upgrade, 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:

compute_AA_freq_in_phylip primtates_21_AA.phy

• 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 del modelo 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 calculan los \(AIC_i\), \(AIC_c\), \(BIC_i\), \(\Delta_{BIC}\), \(BIC_w\) y \(BIC_{cumW}\) asociados a cada modelo, haciendo llamadas a funciones, 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, bien sea a partir de un árbol semilla \(BioNJ\) con 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\).

nw_display -w 120 <(nw_reroot primates_21_AA.phy_HIVb+G_BESTmoves_phyml_tree.txt)

                                                                            +-+ 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 herramientas del paquete \(newick\_utilities\)

Noten 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\)

Les muestro abajo 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.1.1 El criterio bayesiano de información (\(BIC\))

Noten que el mejor modelo seleccionado por \(BIC\) caso es \(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\).

¿A qué se debe esta disparidad entre \(AIC\) y \(BIC\)? 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.

Dado que para muchos alineamientos \(log(n) > 2\), el \(BIC\) suele seleccionar modelos con menos parámetros que \(AIC\).

Con \(BIC\) podemos comparar todos los modelos en competición, independientemente de si están o no anidados.

5.1.2 Tarea: selección de modelos con BIC

1. 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. 
2. Vuelve a correr el script, pero evaluando todos los modelos. $phyml\_protModelFinder.sh$ ¿Hay alguna diferencia? Comenta el resultado.
3. 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$.

---

5.2 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 existe o no un sesgo transicional importante.

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 v0.7_2022-12-01 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 + 24 extended_ef_models OR  standard + 23 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 (K80 F81 HKY85 GTR)

EXAMPLE: phyml_DNAmodelFinder.sh primates.phy 2
 
AIM:  phyml_DNAmodelFinder.sh v0.7_2022-12-01 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-JC tree, optimizing branch lenghts and rates 
      to calculate their AICi, BICi, delta_BIC and BICw
 - Only relevant matrices 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 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. (2022). 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.2.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

phyml_DNAmodelFinder.sh primates.phy 2

=========================================================================================
# 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:

1. Calcula un árbol NJ-JC que servirá de topología de referencia
2. Ajusta los modelos de base, haciendo uso del árbol NJ-JC
3. 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:

1. un sesgo composicional significativo, compositional_heterogeneity_flag=1
2. No considera significativa la diferencia de tasas de transiciones entre purinas y pirimidinas transitional_heterogeneity_flag=0
3. 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:

1. El modelo de dos tasas HKY85+G (010010+G) es el mejor, con un peso \(BIC_w\ =\ 0.85\).
2. 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\).
3. 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:

nw_display <(nw_reroot primates.phy_HKY85+G_BESTmoves_phyml_tree.txt Tarsius_sy)

 +----------------+ 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.

5.2.2 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 excluisvamente de \(phyml\), \(Bash\) y \(awk\), junto a algunas otras herramientas 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 sobreparametrización.

---

6 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.

• Cita: Vinuesa P. (2023). Efficient two-step selection of DNA models with \(phyml\_DNAmodelFinder\). https://github.com/vinuesa/TIB-filoinfo