Árboles filogenéticos
Martín Guacaneme
Mayo de 2021
Nota: El contenido de este cuaderno proviene de clases del profesor Luis Fernando García, del departamento de Biología de la Universidad Nacional de Colombia.
1 Librerías
La librería phytools es usada para visualizar árboles.
library(ape)
library(phangorn)
library(phytools)## Loading required package: maps2 Matriz de distancias
Vamos a utilizar secuencias del virus de la influenza para diferentes especies. El archivo con dichas secuencias alineadas está en formato fasta y se llama chars2.txt. Para empezar volvemos el archivo en un formato compatible con la librería phangorn.
options(width = 300)
influenza = phangorn::read.phyDat(file = "chars2.txt", format = "fasta", type = "DNA")
influenza## 21 sequences with 1698 character and 217 different site patterns.
## The states are a c g tHay 21 secuencias con 1698 caracteres y 217 columnas en las que hay diferentes posiciones. Los estados son adenina timina, etc. Ahora lo que se hace es crear una matriz de distancia para poder crear árboles de distancia o parsimonia a través de ella. En este caso la clase debe ser DNAbin, por eso transformamos influenza en este tipo de clase. La función dist.dna (método o función de objeto o clase) calcula una matriz de distancias por pares de las secuencias de ADN y a partir de una instancia de objecto de clase GNAbin utilizando un modelo de evolución del ADN (por ejemplo K80, Kimura80).
matrizdist = ape::as.DNAbin(influenza)
matrizdist = ape::dist.dna(matrizdist)
matrizdist## IV_PATO PATOVIETNAM HUMANO TAIPOLLO VIETNAM MALLARD HONGKONG II_HONGKONG II_POLLO GANSO PATO POLLO III_POLLO II_PATO GANSO_SHANTOU FALCON GORRION II_GORRION III_GORRION CISNE
## PATOVIETNAM 0.024763883
## HUMANO 0.022720037 0.019713335
## TAIPOLLO 0.019746607 0.018724219 0.020719371
## VIETNAM 0.013740529 0.012736124 0.012727193 0.009785048
## MALLARD 0.017721210 0.016707313 0.016695867 0.013740529 0.007805741
## HONGKONG 0.053711759 0.059037274 0.058993731 0.051618174 0.049492772 0.047342246
## II_HONGKONG 0.047376367 0.054781278 0.052606317 0.045302310 0.043201467 0.041080376 0.005843030
## II_POLLO 0.047376367 0.054781278 0.052606317 0.045302310 0.043201467 0.041080376 0.005843030 0.000000000
## GANSO 0.049456414 0.056884418 0.054701859 0.047376367 0.045268882 0.043140541 0.007800249 0.001940571 0.001940571
## PATO 0.049456414 0.056884418 0.054701859 0.047376367 0.045268882 0.043140541 0.007800249 0.001940571 0.001940571 0.000000000
## POLLO 0.062284039 0.065501415 0.063276104 0.058003631 0.055853090 0.060117188 0.068734673 0.062236986 0.062236986 0.064363061 0.064363061
## III_POLLO 0.037034561 0.035969066 0.033878426 0.032924472 0.026776854 0.029817941 0.050480916 0.044187661 0.044187661 0.046257026 0.046257026 0.050446514
## II_PATO 0.042289991 0.041205597 0.041172026 0.038139166 0.031921469 0.034990591 0.053711759 0.047376367 0.047376367 0.049456414 0.049456414 0.056884418 0.022704000
## GANSO_SHANTOU 0.040193552 0.039115729 0.039084341 0.036057997 0.029865968 0.032924472 0.051579559 0.045268882 0.045268882 0.047342246 0.047342246 0.055811457 0.018683878 0.013729173
## FALCON 0.040193552 0.039115729 0.039084341 0.033985455 0.029865968 0.032924472 0.049456414 0.043170243 0.043170243 0.045236980 0.045236980 0.055811457 0.018683878 0.015721207 0.011745043
## GORRION 0.048433450 0.047342246 0.045206605 0.044250785 0.037987655 0.042124218 0.054664481 0.048332346 0.048332346 0.050413658 0.050413658 0.047221117 0.030743211 0.040038709 0.035893608 0.035893608
## II_GORRION 0.035997216 0.034990591 0.036977610 0.033956464 0.027796222 0.031869226 0.039987351 0.033813845 0.033813845 0.035850823 0.035850823 0.047278605 0.026705125 0.031822925 0.027737100 0.027737100 0.021659674
## III_GORRION 0.037005336 0.038014954 0.035917259 0.034962019 0.028795446 0.032873081 0.045236980 0.038999208 0.038999208 0.041052855 0.041052855 0.044104426 0.025678748 0.034885267 0.030778467 0.030778467 0.020666275 0.014695820
## CISNE 0.039115729 0.046429698 0.042185990 0.039148619 0.034990591 0.039115729 0.051618174 0.045302310 0.045302310 0.047376367 0.047376367 0.042154346 0.033855400 0.039054456 0.036977610 0.036977610 0.038973850 0.028774697 0.031822925
## IV_GORRION 0.046429698 0.049570090 0.047412018 0.044356868 0.040159571 0.044319991 0.052725564 0.046392114 0.046392114 0.048470219 0.048470219 0.034909354 0.032849619 0.042219147 0.038043758 0.038043758 0.023660581 0.031822925 0.028737648 0.0167345503 Árboles de distancias
Estos árboles son árboles de similaridad, no representan inferencia evolutiva.
3.1 Método UPGMA
Con la matriz de distancia perdemos los caracteres a favor de las diferencias en caracteres entre especies. Cuando aparece un valor de 0, significa que no hay diferencia al nivel de los caracteres (nucléotidos) entre las sequencias de dos especies. Ahora creamos un árbol con el método de grupo de pares no ponderados con media aritmética (UPGMA) usando la matriz de distancia.
arbolUPGMA = phangorn::upgma(matrizdist)
plot(arbolUPGMA)
Si la longitud de dos ramas es indéntica, significa que las secuancias también son indénticas y en la matriz de distancia la diferencia es de 0. Por ejemplo, los virus de influenza de pato y ganso son idénticos.
3.2 Método NJ
Ahora hacemos un árbol con el método de unión de vecinos (NJ) usando la misma matriz de distancias:
arbolNJ = nj(matrizdist)
plot(arbolNJ)
Este último árbol puede ser distinto del árbol creado con el método UPGMA. Para personalizar los árboles podemos agregar argumentos a parámetros como cex, para el tamaño de la letra, edge.color, para el grosos de las ramas, etc. También se puede escoger entre diferentes visualizaciones de árbol como filograma, cladograma, radial y demás.
plot(arbolUPGMA, type= "p", cex=0.8, edge.width=2, edge.color="red", font=3)
plot(arbolUPGMA, type= "c", cex=0.8, edge.width=2, edge.color="blue", font=3)
plot(arbolUPGMA, type= "p", label.offset=0.0005, edge.lty=1,
node.pos=2, cex=0.8, edge.width=2, edge.color="black", font=3)
4 Personalización del árbol con phytools
Además de plot podemos graficar árboles con el método plotTree del paquete phytools, el cual es compatible con ape y con phangorn.
plotTree(arbolNJ)
A este también le podemos hacer modificaciones.
plotTree(arbolNJ, ftype="b", fsize=0.8, offset=1, color="red", lwd=2)
En los árboles, sin cambiar la topología, se puede cambiar el orden en que los grupos son visualizados. Por ejemplo, se pueden ordena las puntas de manera alfabética (en la medida de lo posible), o con los grupos más derivados hacia uno de los lados del árbo. Para escalarizar hacia la derecha:
plotTree(ladderize(arbolNJ))
Para escalizar hacia la izquierda:
plotTree(ladderize(arbolNJ, right=FALSE))
Fialmente, a notar que para guardar un árbol podemos usar el comando write.tree(tree, file = "file_name.nex"). Este archivo puede ser leído usando read.tree(file = file_name.nex).
5 Enraizar
Hasta ahora los árboles que hemos construido no están enraízados. Para enraizarlos podemos usar la función root del paquete ape. Por ejemplo, para poner por raíz el virus de la influenza presente en humano pasamos por argumento el nombre de la secuencia que corresponde a esta secuencia en el parámetro outgroup.
arbolNJraiz <-ape::root(arbolNJ, outgroup = "HUMANO", r = TRUE)
plot(arbolNJraiz)
Se puede hacer lo mismo con el árbol creado a partir del método UPGMA.
arbolUPGMAraiz <-ape::root(arbolUPGMA, outgroup = "HUMANO", r=TRUE)
plot(arbolUPGMAraiz)
Además podemos visualizar los dos árboles al tiempo con los siguientes comandos:
layout(matrix(c(1,2)), height=c(10,10))
par(mar=c(1,1,1,1))
plot(arbolUPGMAraiz, label.offset=0.0005, main="ARBOL UPGMA", cex=0.4)
plot(arbolNJraiz, label.offset=0.0005, main="ARBOL NJ", cex=0.4)
6 Árboles de parsimonia
La parsimonia busca disminuir el número de pasos que explican un árbol evolutivo contando el número de cambios de cada uno de los caracteres. En los métodos de distancia (UPGMA y NJ) sólo se llega a un único árbol, mientras que en parsimonia se evalúan múltiples árboles. Se suman todos los pasos reqerridos para explicar un árbol obteniendo un número que se compara entre otros árboles con diferentes pasos. El árbol con menor número de cambios es el ideal.
Hay que tener en cuenta que la parsimonia no usa todos los caracteres. Son eliminados aquellos que son constantes en todos los taxones, como los caracteres 1, 6 y 8 en la matrix de abajo; o aquellos que son varaibles pero no informativos, como 2, 3 y 4. Un caracter es informativo si tiene al menos dos estados de caracter y por lo menos dos estos estados ocurren co una frecuencia mínima de dos.
| Taxa / Caracter | 1 | 2 | 2 | 5 | 5 | 6 | 6 | 8 | 9 |
|---|---|---|---|---|---|---|---|---|---|
| A | A | A | G | A | G | T | G | C | A |
| B | A | G | C | C | G | T | G | C | G |
| C | A | G | A | T | A | T | C | C | A |
| D | A | G | A | G | A | T | C | C | G |
Para estimar árboles de máxima parsiomonia existen varias posbilidades, la más sencilla es partir de árboles de distancia. Se utiliza un árbol de inicio obtenido por distancia y se cuenta su número de pasos. Podemos estimar, por ejemplo, el número de pasos del árbol arbolUPGMAraiz.
phangorn::parsimony(arbolUPGMAraiz, influenza)## [1] 349El árbol arbolUPGMAraiz tiene 349 pasos, algo importante es que aunque esté con raíz o no el número de pasos debe ser el mismo. Probémoslo usando el árbol sin raíz:
phangorn::parsimony(arbolUPGMA, influenza)## [1] 349Con el método optim.parsimony se obtienen el árbol con mejor parsimonia. Este método permite encontrar árboles bajo máxima parsimonia usando árboles de distancia de inicio.
mejorUPGMA = phangorn::optim.parsimony(arbolUPGMAraiz, influenza)## Final p-score 324 after 4 nni operationsAhora hagámoslo con el árbol de NJ:
mejorNJ = phangorn::optim.parsimony(arbolNJraiz, influenza)## Final p-score 324 after 0 nni operationsHacer parsimonia con un número de taxones muy alto se convierte en una labor computacional ardua. Más allá de 10 taxones hay un número de combinaciones difíciles de obtener computacionalmente y se deben hacer búsquedas heurísticas, esto significa que no se evalúan todas las posibilidades de árboles.
Los métodos heurísticos buscan algunos árboles seleccionando un árbol inicial y cortando ramas y poniéndolas en otro lugar del árbol (branch swapping). Se puede interpretar como una muestra representativa de una problación para hacer inferencias a partir de ella. Hay varias formas de búsquedas heurísticas como:
Intercambio de vecino más cercano (Nearest Neighbor Interchange, NNI): Consiste en intercambir ramas adyacentes hacia una rama interna.
Podado de subárbol y reinjerto (Subtree pruning and regrafting, SPR): Se corta un pedazo de rama y se inserta en otra rama de manera azarosa.
Bisección de árbol y reconexión (Tree bisection and reconnection, TBR): Se corta el árbol en la mitad y una porción del árbol que queda de coloca en alguna parte que queda de las ramas del otro árbol.
Los anteriores métodos heurísticos son lentos computacionalente, a veces si crear cientos o miles de árboles con una parsimonia idéntica. Usalmente se debe, entonces, poner una restricción de búsqueda, por ejemplo un límite de árboles con los que trabajar.
Otra estrategia para hacer el proceso de búsqueda de árbol con mayor parsimonia es con el algoritmo de búsqueda pratchet. La cual tiene los siguientes pasos:
- Generar un árbol de inicion con algún nivel de intercambio de ramas
- Seleccionar al azar un sunconjunto de caracteres (nucleótidos) y darles más peso, esto quiere decir que se usan con mayor frecuencia dichos caraceres para hacer los análisis. La cantida de caracteres que se seleccionan es establecida por el usuario, típicamente entre 5 y 25% de los caracteres.
- Realizar itercambio de ramas (branch swapping).
- Iterar pasos 1 a 3 entre 50 y 200 veces.
El algoritmo anterior es computacionalmente más amigable y rápido. Probémoslo usando la función pratchet:
virus_parsimonia = phangorn::pratchet(influenza, all = TRUE)Se genera un árbol, se pesan caracteres por defecto y se itera un número determiando de veces. El algoritmo ha encontrado múltiples ocasiones que el árbol más corto es de 324 pasos ero sólo 4 árboles con igual longitud y número de pasos aunque presenten diferente topología.
virus_parsimonia## 4 phylogenetic treesPara poderlos comparar es necesario enraizarlos. en este caso los enreizamos con el virus de la influenza del humano, de que sabemos que debe de estar menos emparentado con los virus de las aves.
virus_parsimoniaR = ape::root(phy = virus_parsimonia, outgroup = "HUMANO")
plot(virus_parsimoniaR, cex = 0.6)
Para escoger el sólo un árbol con igual parsimonia hay que un árbol de consenso. Una posibilidad es un árbol de consenso estricto, sólo los grupos monofiléticos presentes en todos los árboles son incluidos; si un grupo monofilétco no está presente en todos los árboles colapsa formando una politomía en el árbol consenso. Este tipo de consenso es muy riguroso y suele tener muchas politomías, de modo que se puede hacer un consenso con una regla de mayoría. Por ejemplo, si un grupo está presente en una mayoría del 51% de árboles, entonces el árbol consenso lo incluye, si sólo está presente en el 50%, el consenso no lo usa.
Para hacer un árbol de consenso estricto podemos usar el método ape con parámetro p de 1:
estrictode100 = ape::consensus(virus_parsimoniaR, p = 1)
plot(estrictode100, cex = 0.6)
Para un árbol menos estricto podemos cambiar el valor del parámetro p:
estricto50 = ape::consensus(virus_parsimoniaR, p = 0.5)
plot(estricto50, cex = .6)
El resultado es un árbol con menos politomías, más resuelto.
7 Bootstrap
Para dar soporte a los árboles se puede hacer una serie de seudoréplicas con remplazamiento (bootstrapping) de una matriz. Los cambios entre seudoréplicas consisten en el uso diferencial de los caracteres. Se crean árboles en los que un mismo caracter se repite u otro no se usa.
Por ejemplo, si la matriz original es la siguiente, con los caracteres 1, 2, 3, 4, 5, 6, 7 y 8:
| Taxa / Caracter | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 |
|---|---|---|---|---|---|---|---|---|
| A | A | A | G | C | C | T | A | G |
| B | A | G | C | C | T | A | G | C |
| C | G | C | C | T | G | A | C | T |
| D | G | G | T | C | T | A | G | C |
| Outgroup | G | T | C | G | T | T | A | C |
Una posible seudoréplica es una en la que se usen sólo los caracteres 1, 2, 5, 6, y 8 repiéndolos de modo que hay la misma dimensión de la matriz original.
| Taxa / Caracter | 1 | 2 | 2 | 5 | 5 | 6 | 6 | 8 |
|---|---|---|---|---|---|---|---|---|
| A | A | A | A | C | C | T | T | G |
| B | A | G | G | T | T | A | A | C |
| C | G | C | C | G | G | A | A | T |
| D | G | G | G | T | T | A | A | C |
| Outgroup | G | T | T | T | T | T | T | C |
Con cada réplica se hace un árbol consenso; la veces que un grupo se repita en el conjunto de réplica es el valor de soporte del nodo. Esta técnica, ampliamente utilizada, proporciona evaluaciones de «confianza» para cada clado de un árbol observado, basándose en la proporción de árboles bootstrap que muestran ese mismo clado. Normalmente un nodo con un soporte mayor a 80% es un buen soporte.
arbolesbootstrap = phangorn::bootstrap.phyDat(influenza, FUN = pratchet, bs = 10)En este caso usamos la función pratchet, si usáramos otra se demoraría más. También usamos un número de réplicas irrisosio pues este es un ejercicio demostrativo y no es necesrio gatar mucho tiempo en él. La rutina anterior genera entonces 10 árboles seudoréplicas.
plot(arbolesbootstrap, cex = .6)
Ahora bien, generamos un consenso; en este caso con un consenso al 60%:
estricto100 = ape::consensus(arbolesbootstrap, p = 0.6)
plot(estricto100, cex = .6)
Lo cual genera un árbol con bastante parafilia.
8 Modelos probabilísticos
Métodos de reconstrucción filogenética basados en modelos probabilísticos. Los caracteres en este caso no son sólo caracteres si no que también representn dinámicas evolutivas. A diferencia de la parsimonia, con estos modelos todos los caracteres son útiles, ncluso los que son constantes entre todos los taxones. Además estos cumplen con un modelo evolutivo.
8.1 Árboles de máxima verosimilitud
Se calcula la verosimilitud (probabilidad de obtener un dato según un modelo) de un árbol de acuerdo a un alineamiento de secuencias usando un modelo de evolución de nucleótidos. Este modelo se presenta como una matriz de substitución que muestra las probabilidades de cambio entre los cuatro nucleótidos (por lo tanto esta matriz será 4x4). También se tienen en cuenta la frecuencia de de ocurrencia de las bases. Los pasos para determinar la verosimilitud de un árbol son:
- Se genera un árbol enraizado de inicio de cualqueir tipo, puede ser al azar.
- Se calcula la verosimilitud de casa sitio (de bases) usando un árbol.
- Se calcula la verosimilitud total del árbol
Por nodo, se consideran todos los escenario posibles que pudieron dar origen a los estados de caracter observados para dicho sitio. La verosimilitud de un sitio es, entonces, la suma de verosimilitudes de la reconstrución de este sitio (en forma de árbol) dado un modelo de evolución o sustitución —los cuales se mencionana abajo—.
La verosimilitud total del árbol, \(L\), a diferencia del cálculo por sitio en que se suma, es el producto de las verosimilitudes por sitio.
\[ L = L_1 * L_2 * ... * L_n = \prod_{i = 1}^{n}L_i \] Al final, se suele trabajar con el logaritmo nepetiano de las verosimilitudes para hacer los valores más grandes y manejables:
\[
ln(L) = ln(L_1) * ln(L_2) * ... * ln(L_n) = \prod_{i = 1}^{n}ln(L_i)
\] Por lo anterior, la verosimilitud es computacionalmente costosa. Para el ejemplo usamos de nuevo nuestro objeto influenza de clase phy y creamos un árbol al azar de 21 ramas con rtree como punto de partida.
arbolazar = ape::rtree(n = 21, tip.label = names(influenza))
plot(arbolazar, cex = .5)
En seguida lo enraizamos por el virus de la influenza humana para poderlo visualizar mejor. Además los «escalerizamos» hacia la derecha y le agregamos escala; aquí la longitud de la rama sí es significativa, indica cantidad de cambio en cuanto a sustituciones nucleotídicas (pues se usan nucleótidos en el alineamiento múltiple de las secuencias, chars2.txt).
arbolazarR = ape::root(phy = arbolazar, outgroup = "HUMANO")
plot(ladderize(arbolazarR), cex = .5); add.scale.bar()
A partir del arbol de arriba se peude iniciar la búsqueda del mejor árbol por máxima verosimilitud. Lo primero que se hace es carcular la verosimilitud del árbol dadas las secuencias. Con pml (Phylogenetic maximum likelihood), podemos computar tal verosimilitud.
ajustado = phangorn::pml(arbolazarR, influenza)
ajustado##
## loglikelihood: -32234.36
##
## unconstrained loglikelihood: -5180.792
##
## Rate matrix:
## a c g t
## a 0 1 1 1
## c 1 0 1 1
## g 1 1 0 1
## t 1 1 1 0
##
## Base frequencies:
## 0.25 0.25 0.25 0.25La información que tiene el objeto ajustado nos reporta la verosimilitud del árbol al azar que habíamos creado, que es -5180.792. También reporta un modeo de substitución general, el cual tal vez no se ajuste bien a los datos. Lo que hay que hacer es encontrar un árbol que optimice la verosimilitud usando un modelo de sustitución; para esto vamos a usar el método phangorn::optim.pml, el cual computa la verosimilitud de un árbol filogenético dado un alineamiento múltiple de secuencias y un modelo de evolución de DNA. Toma como argumentos un objeto de clase pml, el tipo de modelo que se queire usar así como el tiempo de rearreglo para los árboles.
ajustadoconGTR = phangorn::optim.pml(object = ajustado, model = "GTR", rearrangement = "ratchet")Para ver el árbol oculto usamos $tree. También lo enraizamos.
ajustadoconGTR$tree
ajustadoconGTRraíz = ape::root(ajustadoconGTR$tree, outgroup = "HUMANO")
plot(ladderize(ajustadoconGTRraíz), cex = .5); add.scale.bar()
El árbol anterior fue generado por el modelo evolutivo GTR. Pero se pueden usar diferentes modelos. Por ejemplo, con HKY y F81.
ajustadoconJC = phangorn::optim.pml(object = ajustado, model = "JC", rearrangement = "ratchet")
ajustadoconHKY = phangorn::optim.pml(object = ajustado, model = "HKY", rearrangement = "ratchet")
ajustadoconF81 = phangorn::optim.pml(object = ajustado, model = "F81", rearrangement = "ratchet")8.1.1 Modelo de sustitución (o de evolución del DNA)
Representaciones del proceso de evolución. Tienen en cuenta la frecuencia de los nucléotidos y la tasa de transición o transversión. Usualmente son más comunes las transiciones porque suelen producir mutaciones sinónimas, lo cual hace que no sean tan sensibles a la selección. También se tiene en cuenta en estos modelos es la heterogeneidad de las secuencias (\(G\)), esto es una medida para la variabilidad de las bases. Otro aspecto relevante en los modelos de sustitución o de evolución son los mismos árboles filogenéticos, de los que se tiene en cuenta las longitud de las ramas y su topología.
El modelo más sencillo es el de Jukes-Cantor, según el cual la frecuencia de las bases en idéntica y las tasas de substitución tampoco varían. Por ello, sólo cuenta con un prámetro, \(\lambda\).
\[ Q = \begin{pmatrix} * & \lambda & \lambda \\ \lambda & * & \lambda \\ \lambda & \lambda & * \end{pmatrix} \]
Por otro lado, mientras Kimura 2 tiene los supuestos de que hay igual frecuencia de bases y que las tasas de substitución son diferentes, F81 supone iguales tasas de sustitución pero diferentes frecuencuencias de bases. Normalmente, a mayor número de parámetros, mayor complejidad.
Hay muchos otros modelos. Se precisa, etnocnes, decidir cuál modelo se ajusta mejor con los datos de trabajo en un filogenía por máxima verosimilitud. Sin embargo, se ha demostrado que el modelo de sustitucion usado no tiene mucho efecto en la filogenia resultante.
Existen programas para el modelo de sustitución que mejor se ajusta a los datos, entre estos están Jmodeltest, PAUP. Para determinar el mejor modelo se compra la verosimilutd de un conjunto de dados usando diferentes modelos, esto mdiante pruebas dde hipóteisis (Likelihood Ratio Test, LRT). El modelo con mejor verosimilitud (valor más bajo, \(L\)) es el que mejor se ajusta a los datos.
«Esencialmente, todos los modelos están mal, pero algunos son útiles» —George E. P. Box
En R una manera de comparar los diferentes modelos es mediante phangorn::AIC.
AIC(ajustadoconGTR, ajustadoconJC, ajustadoconHKY, ajustadoconF81)## df AIC
## ajustadoconGTR 47 8869.358
## ajustadoconJC 39 9284.922
## ajustadoconHKY 43 8870.705
## ajustadoconF81 42 9182.652La primera columna corresponde a los grados libertad. Según el criterio anterior, el mejor modelo que se ajusta con los datos es GTR (con el AIC más bajo) y optimiza la verosimilitud. Pero para evaluar más modelos, podemos usar modelTest de phangorn, el cual necesita un árbol de inicio.
influenzaModelTest = phangorn::modelTest(influenza, ajustadoconGTRraíz)## [1] "JC+I"
## [1] "JC+G"
## [1] "JC+G+I"
## [1] "F81+I"
## [1] "F81+G"
## [1] "F81+G+I"
## [1] "K80+I"
## [1] "K80+G"
## [1] "K80+G+I"
## [1] "HKY+I"
## [1] "HKY+G"
## [1] "HKY+G+I"
## [1] "SYM+I"
## [1] "SYM+G"
## [1] "SYM+G+I"
## [1] "GTR+I"
## [1] "GTR+G"
## [1] "GTR+G+I"influenzaModelTest## Model df logLik AIC AICw AICc AICcw BIC
## 1 JC 39 -4603.461 9284.922 1.480290e-112 9286.804 2.280895e-112 9496.973
## 2 JC+I 40 -4555.867 9191.733 2.547111e-92 9193.713 3.737592e-92 9409.222
## 3 JC+G 40 -4554.713 9189.425 8.077474e-92 9191.405 1.185276e-91 9406.913
## 4 JC+G+I 41 -4553.518 9189.036 9.812024e-92 9191.116 1.369425e-91 9411.962
## 5 F81 42 -4549.326 9182.652 2.388330e-90 9184.834 3.166345e-90 9411.014
## 6 F81+I 43 -4501.910 9089.819 3.438913e-70 9092.107 4.325312e-70 9323.619
## 7 F81+G 43 -4500.739 9087.478 1.108578e-69 9089.766 1.394320e-69 9321.278
## 8 F81+G+I 44 -4499.574 9087.148 1.307396e-69 9089.544 1.558054e-69 9326.385
## 9 K80 40 -4444.749 8969.498 4.612162e-44 8971.477 6.767818e-44 9186.986
## 10 K80+I 41 -4395.506 8873.013 4.124354e-23 8875.092 5.756196e-23 9095.938
## 11 K80+G 41 -4394.356 8870.712 1.302867e-22 8872.792 1.818359e-22 9093.638
## 12 K80+G+I 42 -4392.978 8869.957 1.900660e-22 8872.139 2.519814e-22 9098.320
## 13 HKY 43 -4392.352 8870.705 1.307864e-22 8872.992 1.644973e-22 9104.504
## 14 HKY+I 44 -4344.060 8776.121 4.520840e-02 8778.516 5.387589e-02 9015.358
## 15 HKY+G 44 -4342.978 8773.955 1.334790e-01 8776.351 1.590700e-01 9013.192
## 16 HKY+G+I 45 -4341.876 8773.752 1.477445e-01 8776.258 1.666140e-01 9018.426
## 17 SYM 44 -4439.958 8967.917 1.016726e-43 8970.312 1.211656e-43 9207.154
## 18 SYM+I 45 -4390.044 8870.088 1.779809e-22 8872.594 2.007122e-22 9114.763
## 19 SYM+G 45 -4388.867 8867.735 5.773291e-22 8870.241 6.510642e-22 9112.409
## 20 SYM+G+I 46 -4387.367 8866.734 9.524618e-22 8869.353 1.015121e-21 9116.845
## 21 GTR 47 -4387.679 8869.358 2.564353e-22 8872.093 2.579655e-22 9124.907
## 22 GTR+I 48 -4339.464 8774.929 8.203164e-02 8777.782 7.778965e-02 9035.915
## 23 GTR+G 48 -4338.303 8772.607 2.619839e-01 8775.459 2.484363e-01 9033.593
## 24 GTR+G+I 49 -4337.074 8772.148 3.295527e-01 8775.121 2.942142e-01 9038.571El modelo que mejor se ajusta, según ModelTest es GTR+G+I, el cual tienen en cuanta la heterogneidad. Una vez se sabe el mejor modelo, se puede buscar el mejor árbol, pero como el modelo GTR+G+I no existe como opción hay que expecificarlo con condiciones. De modo que incluya los sitios invariantes y tenga en cuenta la heterogeneidad.
mejorarbol = phangorn::optim.pml(
object = ajustadoconGTR,
model = "GTR",
optInv = TRUE,
optGamma = TRUE,
rearrangement = "ratchet")## only one rate class, ignored optGamma## optimize edge weights: -4387.679 --> -4387.679
## optimize base frequencies: -4387.679 --> -4387.679
## optimize rate matrix: -4387.679 --> -4387.679
## optimize invariant sites: -4387.679 --> -4339.704
## optimize edge weights: -4339.704 --> -4339.55
## optimize topology: -4339.55 --> -4339.55
## 0
## [1] "Ratchet iteration 1 , best pscore\n so far: -4339.4835468905"
## [1] "Ratchet iteration 2 , best pscore\n so far: -4339.4835468905"
## [1] "Ratchet iteration 3 , best pscore\n so far: -4339.4835468905"
## [1] "Ratchet iteration 4 , best pscore\n so far: -4339.4835468905"
## [1] "Ratchet iteration 5 , best pscore\n so far: -4339.4835468905"
## [1] "Ratchet iteration 6 , best pscore\n so far: -4339.4835468905"
## [1] "Ratchet iteration 7 , best pscore\n so far: -4339.4835468905"
## [1] "Ratchet iteration 8 , best pscore\n so far: -4339.4835468905"
## [1] "Ratchet iteration 9 , best pscore\n so far: -4339.4835468905"
## [1] "Ratchet iteration 10 , best pscore\n so far: -4339.4835468905"
## [1] "Ratchet iteration 11 , best pscore\n so far: -4339.4835468905"
## [1] "Ratchet iteration 12 , best pscore\n so far: -4339.4835468905"
## [1] "Ratchet iteration 13 , best pscore\n so far: -4339.4835468905"
## [1] "Ratchet iteration 14 , best pscore\n so far: -4339.4835468905"
## [1] "Ratchet iteration 15 , best pscore\n so far: -4339.4835468905"
## [1] "Ratchet iteration 16 , best pscore\n so far: -4339.4835468905"
## [1] "Ratchet iteration 17 , best pscore\n so far: -4339.4835468905"
## [1] "Ratchet iteration 18 , best pscore\n so far: -4339.4835468905"
## [1] "Ratchet iteration 19 , best pscore\n so far: -4339.4835468905"
## [1] "Ratchet iteration 20 , best pscore\n so far: -4339.4835468905"
## [1] "Ratchet iteration 21 , best pscore\n so far: -4339.4835468905"
## optimize base frequencies: -4339.484 --> -4339.467
## optimize rate matrix: -4339.467 --> -4339.39
## optimize invariant sites: -4339.39 --> -4339.389
## optimize edge weights: -4339.389 --> -4339.389
## optimize topology: -4339.389 --> -4339.389
## 0
## optimize base frequencies: -4339.389 --> -4339.386
## optimize rate matrix: -4339.386 --> -4339.386
## optimize invariant sites: -4339.386 --> -4339.386
## optimize edge weights: -4339.386 --> -4339.386
## optimize base frequencies: -4339.386 --> -4339.386
## optimize rate matrix: -4339.386 --> -4339.386
## optimize invariant sites: -4339.386 --> -4339.386
## optimize edge weights: -4339.386 --> -4339.386mejorarbol##
## loglikelihood: -4339.386
##
## unconstrained loglikelihood: -5180.792
## Proportion of invariant sites: 0.6974688
##
## Rate matrix:
## a c g t
## a 0.0000000 1.055574 6.89374 0.8678635
## c 1.0555742 0.000000 0.29449 8.9680753
## g 6.8937398 0.294490 0.00000 1.0000000
## t 0.8678635 8.968075 1.00000 0.0000000
##
## Base frequencies:
## 0.346711 0.1844641 0.2266299 0.242195La verosimilitud ahora es de -4339.386, también se reporta que el 70% de las secuencias son regioes invariantes. Ahora veamos el árbol.
mejorarbolR = ape::root(mejorarbol$tree, outgroup = "HUMANO")
plot(ladderize(mejorarbolR), cex = 0.5); add.scale.bar()
Este sería nuestro mejor mejor árbol bajo el modelo de máxima verosimilitud.
8.2 Árboles bayesianos
Para esto no usamos R.