library(phangorn)Loading required package: apemtdna <- read.phyDat('../../data/mtDNA.fasta', format = 'fasta')Reconstrucción filogenética
PRÁCTICA 9
Introducción
Aunque existen muy buenos programas para inferir filogenias fuera del entorno de R (es decir, para el terminal de Bash), utilizar las herramientas en R tiene la ventaja de poder representar los árboles, manipularlos y extraer información de ellos con toda la flexibilidad que permite un lenguaje de programación. Otros programas como RAxML o IQtree producirían los resultados en forma de archivos de texto con formato Newick o NEXUS. En R, los árboles serán objetos de clase phylo, que pueden también ser exportados como archivos de texto con formato Newick con la función write.tree() del paquete ape. Este paquete será cargado automáticamente cuando carguemos el paquete phangorn, que proporciona las principales funciones para inferencia filogenética.
Objetivos
- Se busca obtener una filogenia calibrada de hominidos y primates a partir de secuencias mitocondriales.
Datos
Descargamos un alineamiento en formato FASTA de 42 cromosomas mitocondriales completos de: 1 Gorila, 1 Homo heidelbergensis (de hace 400.000 años), 4 humanos denisovanos (de hace entre 30.000 y 110.000 años), 19 humanos neandertales (de hace entre 30.000 y 180.000 años), 9 humanos anatómicamente modernos (actuales).
El archivo descargado recibe el nombre de mtDNA.fasta y se ha guardado en la carpeta “data”.
A continuación ejecutamos el siguiente bloque para cargar el paquete phangorn (que previamente hemos instalado):
Tras ejecutar el bloque anterior, se cargan los datos en el objeto de clase phyDat llamado mtdna. Con la función glance() observamos los datos o con image() podemos generar una imagen del alineamiento.
glance(mtdna) nseq nchar unique_site_pattern parsimony_informative_sites const_sites
1 42 16603 974 2047 13482image(mtdna)
Ejercicio 1
¿Puedes confirmar el número de secuencias? ¿Qué longitud tiene el alineamiento? ¿Qué función usarías para extraer los nombres de las secuencias?
Con la función glance() que ha sido usada anteriormente, descubrimos que hay 42 secuencias (nseq) y que la longitud del alineamiento es de 16603 caracters (nchar). Para extraer los nombres de las secuencias utilizamos la función names().
names(mtdna) [1] "G.gorilla_D38114" "P.troglodytes_JF727212"
[3] "P.troglodytes_JN191207" "P.troglodytes_JF727168"
[5] "P.troglodytes_HM068590" "P.paniscus_KX211934"
[7] "P.paniscus_KX211928" "P.paniscus_GU189658"
[9] "P.paniscus_GU189664" "H.heidelbergensis_NC_023100"
[11] "H.s.denisova_NC_013993" "H.s.denisova_FR695060"
[13] "H.s.denisova_KX663333" "H.s.denisova_KT780370"
[15] "H.s.neandertal_KX198085" "H.s.neandertal_KX198086"
[17] "H.s.neandertal_MG025538" "H.s.neandertal_NC_011137"
[19] "H.s.neandertal_FM865409" "H.s.neandertal_FM865407"
[21] "H.s.neandertal_KX198088" "H.s.neandertal_KX198083"
[23] "H.s.neandertal_MG025540" "H.s.neandertal_MG025536"
[25] "H.s.neandertal_MG025537" "H.s.neandertal_KX198087"
[27] "H.s.neandertal_KX198082" "H.s.neandertal_FM865408"
[29] "H.s.neandertal_KU131206" "H.s.neandertal_KF982693"
[31] "H.s.neandertal_FM865411" "H.s.neandertal_KY751400"
[33] "H.s.neandertal_MK033602" "H.s.modern_EU092846"
[35] "H.s.modern_EU092665" "H.s.modern_JQ045086"
[37] "H.s.modern_DQ779927" "H.s.modern_DQ779932"
[39] "H.s.modern_AP008640" "H.s.modern_AY495317"
[41] "H.s.modern_J01415" "H.s.modern_FJ236981" Árboles basados en distancias
Los árboles NJ y UPGMA pueden proporcionar una idea rápida de cómo podría ser la relación filogenética entre las secuencias. Como métodos basados en distancias que son, necesitamos calcular primero las distancias. Podemos hacerlo en el mismo paso en el que generamos el árbol: nj(dist.ml(mtdna)). O bien, podemos hacerlo en dos pasos: guardar la matriz de distancias en un objeto de clase dist, y después generar los árboles.
La función dist.ml() (de phangorn) solo ofrece un par de modelos de evolución de nucleótidos para estimar distancias: F81 y JC69 (ofrece 17 modelos para proteínas). La función dist.dna() de ape (que no se aplica sobre objetos phyDat, sino DNAbin) ofrece algunas opciones más. Pero en el fondo, los métodos NJ y UPGMA son sólo heurísticos rápidos, con los que no hace falta usar distancias de alta precisión.
distancias <- dist.ml(mtdna, model = 'F81')
mt_nj <- NJ(distancias)
mt_upgma <- upgma(distancias)Ejercicio 2
A) ¿De qué clase son los objetos mt_nj y mt_upgma?
Para conocer la clase de los objetos utilizo class(). En este caso se observa que ambos son de clase “phylo”, es decir, son objetos en los que se guarda la información sobre la estructura del árbol filogénetico.
class(mt_nj)[1] "phylo"class(mt_upgma)[1] "phylo"B) Consulta la estructura de estos objetos.
Para observar la estructura de ambos objetos utilizo str(). Esta función nos proporciona información sobre la estructura interna de un objeto de R.
(str(mt_nj))List of 4
$ edge : int [1:81, 1:2] 47 47 54 54 53 53 52 52 51 51 ...
$ edge.length: num [1:81] 6.03e-05 5.20e-08 -4.37e-08 4.37e-08 8.51e-08 ...
$ tip.label : chr [1:42] "G.gorilla_D38114" "P.troglodytes_JF727212" "P.troglodytes_JN191207" "P.troglodytes_JF727168" ...
$ Nnode : int 40
- attr(*, "class")= chr "phylo"
- attr(*, "order")= chr "postorder"NULL(str(mt_upgma))List of 4
$ edge : int [1:82, 1:2] 63 63 49 49 48 48 59 59 55 55 ...
$ edge.length: num [1:82] 0.00127 0.00127 0.001422 0.000151 0.011063 ...
$ tip.label : chr [1:42] "G.gorilla_D38114" "P.troglodytes_JF727212" "P.troglodytes_JN191207" "P.troglodytes_JF727168" ...
$ Nnode : int 41
- attr(*, "class")= chr "phylo"
- attr(*, "order")= chr "postorder"NULLC) Represéntalos gráficamente (plot(mt_nj)) y consulta la ayuda de la función plot.phylo() para saber qué opciones de representación gráfica tienes.
(plot(mt_nj))
$type
[1] "phylogram"
$use.edge.length
[1] TRUE
$node.pos
[1] 1
$node.depth
[1] 1
$show.tip.label
[1] TRUE
$show.node.label
[1] FALSE
$font
[1] 3
$cex
[1] 1
$adj
[1] 0
$srt
[1] 0
$no.margin
[1] FALSE
$label.offset
[1] 0
$x.lim
[1] 0.0000000 0.1535281
$y.lim
[1] 1 42
$direction
[1] "rightwards"
$tip.color
[1] "black"
$Ntip
[1] 42
$Nnode
[1] 40
$root.time
NULL
$align.tip.label
[1] FALSEplot.phylo(mt_nj)
D) Para qué sirven las funciones add.scale.bar() y axisPhylo()?
Para contestar esta pregunta he consultado help(add.scale.bar) y help(axisPhylo).
add.scale.bar: Esta función añade una barra horizontal que indica la escala de las longitudes de las ramas a la representación gráfica de un árbol filogenético en el dispositivo gráfico actual.
axisPhylo: Esta función añade un eje escalado al lado de un gráfico filogenético.
Enraizar
Escogiendo el nodo
El árbol mt_nj no está enraizado, se comprueba de la siguiente manera:
is.rooted(mt_nj)[1] FALSEAl ejecutarlo nos dice [FALSE], por lo que se confirma que no esta enraizado. Para enraizarlo podemos usar la función root() pero primero hay que indicarr el outgroup o nodo interno en el que se quiere situar la raíz. No obstante, primero para saber identificar los nodos internos invocamos nodelabels() tras haber usado plot(mt_nj). Por ejemplo, el nodo 75 parece ser el ancestro entre gorilas y chimpancés, y por tanto también es ancestro de los humanos.
Por tanto:
mt_nj_rooted <- root(mt_nj, node = 75)
mt_nj_rooted
Phylogenetic tree with 42 tips and 40 internal nodes.
Tip labels:
G.gorilla_D38114, P.troglodytes_JF727212, P.troglodytes_JN191207, P.troglodytes_JF727168, P.troglodytes_HM068590, P.paniscus_KX211934, ...
Unrooted; includes branch length(s).En este caso obtenemos también un árbol filogenético con 42 puntas y 40 puntos internos.
Con un outgroup
La numeración de los nodos internos siempre sigue a la de las hojas (o puntas) del árbol. El orden es caprichoso y después de enraizar, lo que era el 75 en mt_nj ahora será el 43 en mt_nj_rooted. Para garantizar que el enraizado no depende de la numeración, es más seguro usar el nombre del outgroup:
mt_nj_rooted <- root(mt_nj, outgroup = 'G.gorilla_D38114')
mt_nj_rooted
Phylogenetic tree with 42 tips and 40 internal nodes.
Tip labels:
G.gorilla_D38114, P.troglodytes_JF727212, P.troglodytes_JN191207, P.troglodytes_JF727168, P.troglodytes_HM068590, P.paniscus_KX211934, ...
Unrooted; includes branch length(s).En este caso también da un árbol filogenético con 42 puntos y 40 nodos internos.
En un punto medio
Al representar el árbol enraizado (plot(mt_nj_rooted)) verás que no todas las hojas del árbol están a la misma altura. El punto exacto en el que se sitúa la raíz en una rama es arbitrario y root() no se esfuerza mucho en intentar que el árbol parezca ultramétrico. La otra opción es la función midpoint(), que escoge la posición de la raíz para que el árbol esté más equilibrado:
# Esto sobreescribe el objeto anterior:
#|fig-width: 12
#|fig-height: 12
mt_nj_rooted <- midpoint(mt_nj)
plot(mt_nj_rooted, align.tip.label = TRUE, main = 'NJ, enraizado en punto medio')
nodelabels(frame = 'none')
add.scale.bar()
Con datación de hojas
El enraizado en punto medio asume que todas las hojas son contemporáneas, como sucede en muchos casos. Pero, en nuestros datos tenemos secuencias contemporáneas y también arcaicas. Conociendo estas edades, podemos aplicar un método de enraizado llamado root-to-tip regression, mediante la función rtt() de ape.
Lo primero es crear un vector de fechas (o edades) para todas las hojas del árbol. Podemos expressar las edades de las secuencias arcaicas como el número (negativo) de años que hace que dejaron de evolucionar. Para saber el orden en que debemos especificar las edades, podemos usar mt_nj$tip.label.
fechas <- c(0, 0, 0, 0, 0, 0, 0, 0, 0,
-400000,
-39000, -50000, -110000, -110000,
-39820, -40096, -38515, -38310, -39000, -40000, -41210,
-42430, -42540, -43230, -43780, -44290, -44770, -49000,
-50000, -82752, -65000, -122287, -110450,
0, 0, 0, 0, 0, 0, 0, 0, 0)
# Observa que hay que asignar el resultado de setNames() a "fechas"
fechas <- setNames(fechas, mt_nj$tip.label)
fechas G.gorilla_D38114 P.troglodytes_JF727212
0 0
P.troglodytes_JN191207 P.troglodytes_JF727168
0 0
P.troglodytes_HM068590 P.paniscus_KX211934
0 0
P.paniscus_KX211928 P.paniscus_GU189658
0 0
P.paniscus_GU189664 H.heidelbergensis_NC_023100
0 -400000
H.s.denisova_NC_013993 H.s.denisova_FR695060
-39000 -50000
H.s.denisova_KX663333 H.s.denisova_KT780370
-110000 -110000
H.s.neandertal_KX198085 H.s.neandertal_KX198086
-39820 -40096
H.s.neandertal_MG025538 H.s.neandertal_NC_011137
-38515 -38310
H.s.neandertal_FM865409 H.s.neandertal_FM865407
-39000 -40000
H.s.neandertal_KX198088 H.s.neandertal_KX198083
-41210 -42430
H.s.neandertal_MG025540 H.s.neandertal_MG025536
-42540 -43230
H.s.neandertal_MG025537 H.s.neandertal_KX198087
-43780 -44290
H.s.neandertal_KX198082 H.s.neandertal_FM865408
-44770 -49000
H.s.neandertal_KU131206 H.s.neandertal_KF982693
-50000 -82752
H.s.neandertal_FM865411 H.s.neandertal_KY751400
-65000 -122287
H.s.neandertal_MK033602 H.s.modern_EU092846
-110450 0
H.s.modern_EU092665 H.s.modern_JQ045086
0 0
H.s.modern_DQ779927 H.s.modern_DQ779932
0 0
H.s.modern_AP008640 H.s.modern_AY495317
0 0
H.s.modern_J01415 H.s.modern_FJ236981
0 0 A continuación se escoge la posición de la raíz que optimiza la regresión entre distancia de raíz a puntas y edades de las puntas:
mt_nj_rtt <- rtt(mt_nj, tip.dates = fechas, objective = 'rms')
mt_nj_rtt
Phylogenetic tree with 42 tips and 41 internal nodes.
Tip labels:
G.gorilla_D38114, P.troglodytes_JF727212, P.troglodytes_JN191207, P.troglodytes_JF727168, P.troglodytes_HM068590, P.paniscus_KX211934, ...
Rooted; includes branch length(s).Ejercicio 3
A) ¿Se corresponden bien las distancias evolutivas estimadas con las dataciones?
Primero representamos el árbol filogenético teniendo en cuenta las dataciones:
#|fig-width: 12
#|fig-height: 7
plot(mt_nj_rtt, align.tip.label = TRUE, show.tip.label=TRUE)
add.scale.bar()
CONCLUSIÓN: Las muestras arqueológicas deben de haber dejado de acumular mutaciones hace tiempo. Homo heldeibergensis es el que mejor se ve, su rama se queda más atrás que otros. Lo que no es normal que los denisovanos tengan una rama tan larga. La función rtt busca la mejor correlación entre longitudes de rama y el tiempo que ha pasado desde que dejaron de existir. Este árbol nos dice que no se aplicaría a la perfección el reloj molecular. Podemos desconfiar de que la datación sea correcta, podemos haber asignado a los denisovanos una edad más antigua de la que parece (igual nos hemos equivocado al transcribir los datos, la media no es representativa de la edad real).
B) ¿Qué secuencias crees que se ajustan menos a su supuesta edad?
Eliminamos la secuencia del gorila del dataset porque es más antigua que las otras. Queremos un outgroup más cercano que P.troglodytes, pero como no hay nada más cercanos a los humanos nos quedamos con ellos. También hay pares de secuencias de neandertales que son idénticos entre ellos, de cada uno de esos grupos nos quedamos con una sola copia.
C) ¿Qué razones puede haber para que la edad de una secuencia no se corresponda con la cantidad de cambios que ha acumulado?
Las muestras arqueologicas pueden tener errores - el método estadístico para corregir los errores de secuenciación y degradación del DNA puede que haya sido mejor en algunas muestras que en otras. Esto puede producir una tendencia a alargar las ramas que pueden parecer mutaciones pero son efectos de la degradación.
Distancias
Al crear anres el objeto distancias, se puede usar para conocer la cantidad estimada de cambios entre los cromosomas mitocondriales humanos y de neandertales, por ejemplo. El objeto de clase dist es incómodo, pero podemos obtener una simple matriz fácilmente:
distMat <- as.matrix(distancias)Aprovechando que los nombres de las secuencias tienen una estructura coherente, podemos usar startsWith() para seleccionar a los humanos modernos o a los neandertales, etc.:
# Esto son vectores lógicos que podemos usar después para seleccionar filas
# y columnas de la matriz de distancias.
modernos <- startsWith(colnames(distMat), 'H.s.modern')
neandertales <- startsWith(colnames(distMat), 'H.s.neandertal')
denisovanos <- startsWith(colnames(distMat), 'H.s.denisova')
chimpances <- startsWith(colnames(distMat), 'P.troglodytes')
bonobos <- startsWith(colnames(distMat), 'P.paniscus')La distancia estimada media entre humanos y neandertales se puede obtener como:
mean(distMat[modernos, neandertales])[1] 0.01214282En este caso la distancia estimada es de 0.01214282.
Ejercicio 4
A) Calcula las distancias medias entre todos los grupos.
Como previamente hemos calculado la distancia para cada uno de los grupos, ahora calculamos la media de ambos (de dos en dos). Para ello utilizo la función mean().
mean(distMat[modernos, denisovanos])[1] 0.02250501mean(distMat[modernos, chimpances])[1] 0.09104153mean(distMat[modernos, bonobos])[1] 0.09239182mean(distMat[modernos, neandertales])[1] 0.01214282mean(distMat[neandertales, bonobos])[1] 0.09089352mean(distMat[neandertales, denisovanos])[1] 0.02201398mean(distMat[neandertales, chimpances])[1] 0.09041225mean(distMat[denisovanos, bonobos])[1] 0.09297814mean(distMat[denisovanos, chimpances])[1] 0.09129525mean(distMat[chimpances, bonobos])[1] 0.04059035B) ¿Cuántas veces mayor es la distancia entre humanos y chimpancés que entre humanos y neandertales?
Para calcular esto lo hago en modo de proporción, en este caso, divido la distancia media entre humanos y chimpancés por la distancia media entre humanos y neandertales. Vemos que es 7,5 veces (aprox.) más grande la distancia de los primeros entre los segundos.
mean(distMat[modernos, chimpances])/mean(distMat[modernos, neandertales])[1] 7.497561C) ¿Cómo calcularías la distancia media entre dos secuencias de humanos modernos?
Siguiendo la misma lógica que en los apartados anteriores:
mean(distMat[modernos, modernos])[1] 0.003331521Este resultado es bajo pero más alto de lo que encontraríamos en cromosomas nucleares.
Otra alternativa presentada en clase es:
d <- distMat[modernos,modernos]
mean(d[lower.tri(d)])[1] 0.003747961El resultado es similar al obtenido previamente. Esta en el mismo orden de magnitud.
Por otra parte, atendiendo a la solución presentada en el guion de prácticas del Aula Virtual:
# La submatriz siguiente contiene ceros en la diagonal. Si queremos excluir
# del cálculo de la media las distancias entre cada secuencia consigo misma,
# habrá que tenerlo en cuenta:
distMat[modernos, modernos] H.s.modern_EU092846 H.s.modern_EU092665 H.s.modern_JQ045086
H.s.modern_EU092846 0.000000000 0.004115877 0.005454272
H.s.modern_EU092665 0.004115877 0.000000000 0.005393179
H.s.modern_JQ045086 0.005454272 0.005393179 0.000000000
H.s.modern_DQ779927 0.005392590 0.005270729 0.004724566
H.s.modern_DQ779932 0.005513658 0.005391783 0.004906390
H.s.modern_AP008640 0.005517618 0.005152317 0.004787395
H.s.modern_AY495317 0.005635614 0.005148861 0.004784781
H.s.modern_J01415 0.005453806 0.005210372 0.004785411
H.s.modern_FJ236981 0.005331436 0.004966433 0.004481205
H.s.modern_DQ779927 H.s.modern_DQ779932 H.s.modern_AP008640
H.s.modern_EU092846 0.005392590 0.005513658 0.005517618
H.s.modern_EU092665 0.005270729 0.005391783 0.005152317
H.s.modern_JQ045086 0.004724566 0.004906390 0.004787395
H.s.modern_DQ779927 0.000000000 0.001087163 0.002359394
H.s.modern_DQ779932 0.001087163 0.000000000 0.002298600
H.s.modern_AP008640 0.002359394 0.002298600 0.000000000
H.s.modern_AY495317 0.002478978 0.002660427 0.001692960
H.s.modern_J01415 0.002418729 0.002600192 0.001693182
H.s.modern_FJ236981 0.002055212 0.002236605 0.001269395
H.s.modern_AY495317 H.s.modern_J01415 H.s.modern_FJ236981
H.s.modern_EU092846 0.005635614 0.0054538063 0.0053314359
H.s.modern_EU092665 0.005148861 0.0052103722 0.0049664328
H.s.modern_JQ045086 0.004784781 0.0047854113 0.0044812052
H.s.modern_DQ779927 0.002478978 0.0024187289 0.0020552122
H.s.modern_DQ779932 0.002660427 0.0026001920 0.0022366046
H.s.modern_AP008640 0.001692960 0.0016931824 0.0012693954
H.s.modern_AY495317 0.000000000 0.0010873068 0.0009058850
H.s.modern_J01415 0.001087307 0.0000000000 0.0006642889
H.s.modern_FJ236981 0.000905885 0.0006642889 0.0000000000# Por ejemplo:
z <- distMat[modernos, modernos]
z[z == 0] <- NA
mean(z, na.rm = TRUE)[1] 0.003747961Las distancias estimadas no son exactas. Cuanto mayor es la distancia, mayor es el error en la estimación. Además el modelo disponible podia no ser el adecuado. Pero nos ayudan a estimar el tiempo transcurrido desde que dos linajes se separaron, si asumimos el reloj molecular. Por ejemplo, estas distancias sugieren que los humanos modernos divergieron de los neandertales hace una cantidad de tiempo aproximadamente igual al 13.3% del tiempo que nos separa de los chimpancés, que podrían ser más de 850000 años si del chimpancé nos separan 6.5 millones de años.
En los árboles las distancias se representan mediante la suma de las longitudes que separan las hojas. Las podemos obtener con la función cophenetic():
distCofen <- cophenetic(mt_nj_rooted)Ejercicio 5
¿Cuánto se parecen las distancias estimadas a las distancias cofenéticas?
Para valorar esto representamos gráficamente de la siguiente manera:
# Aquí usamos las matricees como vectores:
plot(distMat, distCofen, xlab = 'Distancias F81', ylab = 'Dist. cofenéticas')
abline(a = 0, b = 1, col = 'red', lty = 2)
Limpieza o filtrado del alineamiento
Las distancias crudas, como número de diferencias, nos pueden informar de la presencia de pares de secuencias idénticas. Podríamos dejarlas, pero algunos de los métodos que usaremos a continuación pueden generar errores, por ejemplo cuando dos secuencias idénticas se supone que pertenecen a épocas diferentes.
crudas <- dist.dna(as.DNAbin(mtdna), model = 'raw')
crudas <- as.matrix(crudas)
identicas <- which(crudas == 0 & lower.tri(crudas), arr.ind = TRUE)
# posiciones de la semimatriz inferior donde las distancias son 0:
identicas row col
H.s.neandertal_KX198086 16 15
H.s.neandertal_MG025540 23 15
H.s.neandertal_MG025540 23 16
H.s.neandertal_KX198088 21 20
H.s.neandertal_KX198083 22 20
H.s.neandertal_KX198082 27 20
H.s.neandertal_KX198083 22 21
H.s.neandertal_KX198082 27 21
H.s.neandertal_KX198082 27 22repetidas <- unique(c(identicas[, 'row'], identicas[, 'col']))
heatmap(crudas[repetidas, repetidas], cexRow = 0.7, cexCol = 0.7)
Así pues, las secuencias de Neandertales KX198085, KX198086 y MG025540 (todas ellas con una edad estimada de unos 40000 años) son idénticas entre ellas; y las secuencias, también Neandertales FM865407, KX198082, KX198088 y KX198083 son idénticas entre ellas. Bueno, en realidad, la comparación ignora las posiciones en las que una secuencia contiene ambigüedades (Ns). Dejaremos sólo una de cada grupo: MG025540 y FM865407, que son las que menos ambigüedades tienen en la secuencia. Como no conocemos una función que elimine secuencias de un objeto phyDat, necesitamos especificar aquellas con las que nos queremos quedar. Podríamos usar la función subset() o bien la función [, tal como la aplicaríamos a una matriz. También existe una función unique(), pero si la usamos perdemos control sobre qué representantes conservamos de cada grupo de secuencias idénticas.
# La exclamación ("!") invierte el valor lógico de lo que la sigue.
# La función "grepl()" produce un vector lógico con "TRUE" donde los nombres
# de las secuencias contienen alguna de esas palabras.
mtdna37 <- mtdna[! grepl('(KX198085|KX198086|KX198082|KX198088|KX198083)',
names(mtdna)), ]
mtdna3737 sequences with 16603 character and 974 different site patterns.
The states are a c g t Hay 37 secuencias con 16603 carácteres y 974 site patterns. Los estados son: a c g t.
Ejercicio 6
¿Serías capaz de crear otro conjunto de datos, mtdna28, que sólo contenga secuencias del género Homo?
mtdna28 <- mtdna37[grepl('^H.', names(mtdna37)), ]
mtdna2828 sequences with 16603 character and 974 different site patterns.
The states are a c g t Identifica 28 secuencias.
Por otra parte, la secuencia del gorila no nos hace falta: es demasiado divergente y puede haber evolucionado a un ritmo diferente, por tanto:
mtdna36 <- mtdna37[! startsWith(names(mtdna37), 'G.gorilla'),]Escoger el modelo de evolución molecular
Las secuencias mitocondriales proporcionan suficiente información para determinar con confianza la topología del árbol que relaciona los grupos principales. Sólo algunas regiones de árbol, como la relación exacta entre los miembros de un mismo grupo, presentan dificultad. Para estimar mejor la topología y las longitudes de las ramas, debemos usar el método de máxima verosimilitud, para el cual sí es fundamental escoger el modelo evolutivo más adecuado. Para ello usamos la funcion modelTest(), que devuelve un objeto de clase modelTest, que también puede usarse como un simple marco de datos.
modelos <- modelTest(mtdna36, control = pml.control(trace = 0))
head(modelos) Model df logLik AIC AICw AICc AICcw BIC
1 JC 69 -43394.16 86926.33 0 86926.91 0 87458.83
2 JC+I 70 -42831.03 85802.05 0 85802.66 0 86342.27
3 JC+G(4) 70 -42827.52 85795.04 0 85795.64 0 86335.25
4 JC+G(4)+I 71 -42807.54 85757.08 0 85757.70 0 86305.01
5 F81 72 -42567.47 85278.94 0 85279.57 0 85834.58
6 F81+I 73 -41987.61 84121.22 0 84121.87 0 84684.58Lo que hace modelTest() es calcular la verosimilitud de cada modelo sobre un mismo árbol. Al hacerlo, debe optimizar los parámetros de cada modelo, pero no optimiza el árbol. Sería una pérdida de tiempo no guardar los parámetros optimizados, que son un excelente punto de partida para la búsqueda del árbol más verosímil. Podemos extraer cualquiera de los modelos comparados con la función as.pml().
fit <- as.pml(modelos, 'GTR+G(4)+I')
fitmodel: GTR+G(4)+I
loglikelihood: -39275.11
unconstrained loglikelihood: -45918.42
Proportion of invariant sites: 0.6583584
Model of rate heterogeneity: Discrete gamma model
Number of rate categories: 4
Shape parameter: 0.9797904
Rate Proportion
1 0.0000000 0.6583584
2 0.3890159 0.0854104
3 1.3774432 0.0854104
4 2.9170011 0.0854104
5 7.0247157 0.0854104
Rate matrix:
a c g t
a 0.000000 1.9682572 63.0009475 1.212889
c 1.968257 0.0000000 0.6641116 46.671218
g 63.000948 0.6641116 0.0000000 1.000000
t 1.212889 46.6712184 1.0000000 0.000000
Base frequencies:
a c g t
0.3094827 0.3112863 0.1305929 0.2486381 Aquí se demuestra que modelos no es un simple marco de datos, sino que tiene escondida mucha información, incluyendo de hecho el árbol y el alineamiento. También podemos usar as.pml() para extraer el modelo más adecuado a los datos según alguno de los criterios disponibles (BIC, AIC, AICc). El motivo por el que el modelo más adecuado no es necesariamente el más verosímil es que la verosimilitud siempre será mayor para un modelo más complejo (con más parámetros), porque tiene más grados de libertad para ajustarse a los datos. Los criterios de selección penalizan los modelos con más parámetros para hacer la comparación más justa. El peligro de usar un modelo con demasiados parámetros es el de sobreajustar el modelo, es decir, tomar el ruído natural de los datos como algo que el modelo también debe explicar.
fit <- as.pml(modelos, 'BIC')
fitmodel: TIM2+G(4)+I
loglikelihood: -39278.26
unconstrained loglikelihood: -45918.42
Proportion of invariant sites: 0.6579267
Model of rate heterogeneity: Discrete gamma model
Number of rate categories: 4
Shape parameter: 0.9785258
Rate Proportion
1 0.0000000 0.65792675
2 0.3877809 0.08551831
3 1.3745614 0.08551831
4 2.9126763 0.08551831
5 7.0183832 0.08551831
Rate matrix:
a c g t
a 0.000000 2.033328 78.09853 2.033328
c 2.033328 0.000000 1.00000 57.925290
g 78.098532 1.000000 0.00000 1.000000
t 2.033328 57.925290 1.00000 0.000000
Base frequencies:
a c g t
0.3094827 0.3112863 0.1305929 0.2486381 Ejercicio 7
¿Cuál es el mejor modelo para secuencias sólo del género Homo?
Para hacer esto he seguido el mismo planteamiento que en el caso anterior.
modelos2 <- modelTest(mtdna28, control = pml.control(trace = 0))
fit2 <- as.pml(modelos2, 'BIC')
fit2 model: TrN+G(4)+I
loglikelihood: -28823.51
unconstrained loglikelihood: -45918.42
Proportion of invariant sites: 0.8149294
Model of rate heterogeneity: Discrete gamma model
Number of rate categories: 4
Shape parameter: 1.003691
Rate Proportion
1 0.0000000 0.81492938
2 0.7439846 0.04626766
3 2.5820334 0.04626766
4 5.4066329 0.04626766
5 12.8807197 0.04626766
Rate matrix:
a c g t
a 0.000 1.00000 53.837 1.00000
c 1.000 0.00000 1.000 32.09634
g 53.837 1.00000 0.000 1.00000
t 1.000 32.09634 1.000 0.00000
Base frequencies:
a c g t
0.3088407 0.3122559 0.1312366 0.2476667 fit2 <- as.pml(modelos2, 'BIC')
fit2model: TrN+G(4)+I
loglikelihood: -28823.51
unconstrained loglikelihood: -45918.42
Proportion of invariant sites: 0.8149294
Model of rate heterogeneity: Discrete gamma model
Number of rate categories: 4
Shape parameter: 1.003691
Rate Proportion
1 0.0000000 0.81492938
2 0.7439846 0.04626766
3 2.5820334 0.04626766
4 5.4066329 0.04626766
5 12.8807197 0.04626766
Rate matrix:
a c g t
a 0.000 1.00000 53.837 1.00000
c 1.000 0.00000 1.000 32.09634
g 53.837 1.00000 0.000 1.00000
t 1.000 32.09634 1.000 0.00000
Base frequencies:
a c g t
0.3088407 0.3122559 0.1312366 0.2476667 El mejor modelo en este caso es TrN+G(4)+I
Máxima verosimilitud
Tras seleccionar el modelo, optimizamos la topología del árbol, longitudes de la ramas y otros parámetros del modelo. Usaremos la función pml_bb(), que internamente usa las funciones pml() y optim.pml(). Podemos modular el comportamiento de pml_bb() (cuánta información muestra del proceso, cuántas iteraciones usa, cuándo decide que no puede mejorar el modelo…) con los parámetros definidos mediante pml.control(). A la función pml_bb() podemos proporcionarle el alineamiento (objeto de clase phyDat) y el nombre del modelo de evolución molecular, o bien un objeto de clase modelTest, de donde toma automáticamente el alineamiento y el modelo con un valor de BIC (Bayesian Information Criterion) más bajo.
mt_ml <- pml_bb(modelos, control = pml.control(trace = 0))El objeto mt_ml, de clase pml, contiene el árbol (mt_ml$tree), además de los parámetros optimizados del modelo. Antes de ver el árbol, echémosle un vistazo al modelo ajustado, para ver cuánto se aparta del que habíamos ajustado con modelTest().
mt_mlmodel: TIM2+G(4)+I
loglikelihood: -39277.23
unconstrained loglikelihood: -45918.42
Proportion of invariant sites: 0.6575751
Model of rate heterogeneity: Discrete gamma model
Number of rate categories: 4
Shape parameter: 0.9767073
Rate Proportion
1 0.000000 0.65757508
2 0.386313 0.08560623
3 1.371503 0.08560623
4 2.908756 0.08560623
5 7.014820 0.08560623
Rate matrix:
a c g t
a 0.000000 2.032831 77.89999 2.032831
c 2.032831 0.000000 1.00000 57.982110
g 77.899986 1.000000 0.00000 1.000000
t 2.032831 57.982110 1.00000 0.000000
Base frequencies:
a c g t
0.3094827 0.3112863 0.1305929 0.2486381 El modelo recomendado es TIM2+G(4)+I
Por defecto, pml_bb() usa un método de reordenación del árbol denominado “stochastic”, genera un árbol no enraizado y 1000 pseudo-réplicas de bootstrap. Al pedirle la representación gráfica, añade automáticamente los valores de soporte de bootstrap y lo enraiza en el punto medio:
plot(mt_ml, align.tip.label = TRUE, main = 'Árbol más verosímil')
Ejercicio 8
A) El árbol no está enraizado. Como las secuencias no son todas de la misma época, utiliza rtt() para situar la raíz en el punto más compatible con la hipótesis de reloj molecular.
Para ello he seguido el mismo procedimiento que en el que se ha realizado anteriormente.
fechas2 <- c(0, 0, 0, 0, -400000, 0, 0, 0, 0, -39000, -50000, -110000, -110000, -122287, -43230, -82752, -50000, -43780, -39000, -38310, -38515, -42540, -40000, -44290, -49000, -65000, -110450, 0, 0, 0, 0, 0, 0, 0, 0, 0)
fechas2 <- setNames(fechas2, mt_ml$tip.label)
fechas2 [1] 0 0 0 0 -400000 0 0 0 0
[10] -39000 -50000 -110000 -110000 -122287 -43230 -82752 -50000 -43780
[19] -39000 -38310 -38515 -42540 -40000 -44290 -49000 -65000 -110450
[28] 0 0 0 0 0 0 0 0 0mt_ml_rtt <- rtt(mt_ml$tree, tip.dates = fechas2, objective = 'rms')
plot(mt_ml_rtt, align.tip.label = TRUE)
B) ¿Podrías obtener el árbol más verosímil para el subconjunto de secuencias del género Homo?
mt_ml2 <- pml_bb(modelos2,
control = pml.control(trace = 0))
plot(mt_ml2, align.tip.label = TRUE, main = 'Árbol más verosímil')
Con datación de puntas (tip dating)
El método de máxima verosimilitud genera árboles no enraizados por defecto, como consecuencia de la reversibilidad temporal. En esos árboles las longitudes de ramas són el número esperado de sustituciones por sitio, sin ninguna referencia explícita al tiempo. Pero cuando no todas las secuencias son isócronas, podemos usar la información de sus edades para constreñir la búsqueda del árbol más verosímil. Al usar las edades de las secuencias, estamos asumiendo un reloj molecular: una relación constante entre el tiempo y el número esperado de sustituciones acumuladas. Entonces, la posición de la raíz afecta la verosimilitud y el resultado es un árbol enraizado.
# Habíamos creado el vector "fechas", pero incluye al gorila.
fechas36 <- fechas[names(mtdna36)]
mt_ml_clock <- pml_bb(modelos,
method = 'tipdated',
tip.dates = fechas36, rearrangement = 'NNI',
control = pml.control(trace = 0))
mt_ml_clockmodel: TIM2+G(4)+I
loglikelihood: -41234.09
unconstrained loglikelihood: -45918.42
Proportion of invariant sites: 0.7646238
Model of rate heterogeneity: Discrete gamma model
Number of rate categories: 4
Shape parameter: 1.092121
Rate Proportion
1 0.000000e+00 0.76462377
2 7.627412e-09 0.05884406
3 2.474842e-08 0.05884406
4 4.989794e-08 0.05884406
5 1.146144e-07 0.05884406
Rate: 1.15857e-08
Rate matrix:
a c g t
a 0.000000 2.319011 57.89519 2.319011
c 2.319011 0.000000 1.00000 47.953244
g 57.895191 1.000000 0.00000 1.000000
t 2.319011 47.953244 1.00000 0.000000
Base frequencies:
a c g t
0.3094827 0.3112863 0.1305929 0.2486381
Rate: 1.15857e-08 En este caso el modelo TIM2+G(4)+I.
plot(mt_ml_clock, align.tip.label = TRUE)
Ejercicio 9
A)¿Cuál es la tasa de mutación global estimada?
A partir de los datos obtenidos anteriores, concluyo que la tasa de mutación es de 1.15857e-08.
B) ¿En qué unidades está expresada?
Corresponde al número de mutaciones por sitio, por lo que estas son las unidades.
C) ¿Cuánto tiempo dice el árbol que divergimos de los chimpancés? (Puedes usar la función cophenetic(), puesto que en el árbol datado las longitudes de las ramas están en años).
He consultado la ayuda de la función cophenetic () y luego he procedido tal y como se ha hecho anteriormente en esta práctica. Esta función calcula la distancia entre dos ramas del árbol (como se ha visto antes, en los árboles las distancias se representan mediante la suma de las longitudes que separan las hojas). En mi caso, he visualizado en la consola de R studio la distancia entre humanos-chimpancés que es de 1804164,5 años.
distCofen2 <- cophenetic(mt_ml_clock$tree)Estimación de edades
Las edades de las secuencias aportan información sobre el reloj molecular en la medida en que se echen de menos las mutaciones que dejaron de acumularse en las secuencias arcaicas. Esta información puede ser muy útil en filogenias de virus de ARN muestreados a lo largo de los años. Pero en una filogenia de millones de años, la poca información de las puntas no es suficiente para datar con exactitud las divergencias antiguas. Hay otras estrategias. Podríamos fijar una tasa de mutación conocida y estimar sólo los otros parámetros. Alternativamente, podemos dar la fecha de algún nodo interno para ayudar a calibrar el reloj. Por ejemplo, supondremos que el ancestro común más reciente entre humanos y chimpancés vivió hace 6.5 millones de años.
Esta información se la proporcionamos a la función estimate.dates() mediante un vector de fechas o edades de las secuencias que incluya no sólo las hojas del árbol sino también los nodos internos, o al menos aquellos de los que conocemos la edad. El vector debe tener tantos componentes como nodos internos y hojas tiene el árbol, con valores NA en las edades desconocidas. Para saber en qué orden están numerados los nodos, podemos ejecutar lo siguiente:
mt_ml_rtt <- rtt(mt_ml$tree, fechas36)
plot(mt_ml_rtt, align.tip.label = TRUE)
nodelabels()
Como suele suceder, la raíz (en los árboles enraizados) es el nodo que sigue a las hojas, es decir el 37, porque los 36 primeros números son las puntas.
fechasNodos <- c(fechas36, -6500000, rep(NA, 34))
fechasEstimadas <- estimate.dates(mt_ml_rtt, fechasNodos)
tasaEstimada <- estimate.mu(mt_ml_rtt, fechasNodos)
# Creamos una copia del árbol...
mt_ml_años <- mt_ml_rtt
# ...y actualizamos las longitudes de sus ramas como la resta # entre las fechas de los nodos que conectan:
mt_ml_años$edge.length <- fechasEstimadas[mt_ml_años$edge[,2]] -
fechasEstimadas[mt_ml_años$edge[,1]]
plot(mt_ml_años)
axisPhylo()
tasaEstimada[1] 1.314216e-08La tasa estimada es de 1.31428e-08.
Al forzar que la raíz esté a unos 6.5 millones de años, la tasa de mutación estimada disminuye. Es habitual que las tasas parezcan diferentes a diferentes profundidades del árbol, o entre linajes diferentes.