## Instala R y Rstudio (Posit Cloud)
Por favor instala en tu computadora
- Version gratuita para estudiantes
o si quieres trabajar online solo abre una cuenta de (tambien gratuita para estudiantes)
## Correr codigo en Rmarkdown
En cada seccion vas a encontrar texto y codigo. El codigo esta agregado entre entre tres acentos al reves y empieza con las llaves {r echo=TRUE}. Junto a esa linea vas a ver una flecha verde. Esa flecha ejecuta el codigo. Los archivos .Rmd son una herramienta muy util para reproducir el codigo. Son una combinacion de texto con codigo ejecutable y nos permiten crear un documento de manera eficiente donde todos los datos se pueden manejar.
## Lee el arbol y la filogenia
Los datos que vas a leer son datos de la polinizacion de Thalictrum con una filogenia de 103 especies. Vamos a generar el modelo de Markov y su estimacion utilizando los paquetes mas populares de R .
Por favor installa los siguientes paquetes (ejecuta con la flecha verde)
```{r echo=TRUE}
#include=FALSE
#install.packages("ape")
#install.packages("geiger")
#install.packages("corHMM")
```
Y cargalos para poder empezar a analizar datos
```{r echo=TRUE}
#Libraries required
library(ape)
library(geiger)
library(phytools)
library(corHMM)
```
## Lectura de los datos
Lo primero es saber donde tenemos guardamos nuestros archivos. Yo por ejemplo los tengo en este folder `~Dropbox/Teaching/Workshops/sureste2024`. La funcion `setwd()` nos ayuda a marcar el folder donde estaremos trabajando.
```{r}
## Leer los datos
setwd("~/Dropbox/Teaching/Workshops/sureste2024")
pollination_tree<-read.tree("poliniza_arbol.tre")
pollination_data<-read.csv("poliniza_datos.csv", header=TRUE)
head(pollination_data)
## Graficar un arbol
plotTree(pollination_tree,fsize=0.4,ftype="i")
## Graficar un arbol circular
#plotTree()
```
El paquete `phytools` requiere que los datos sean leidos como lista
```{r}
fmode<-setNames(pollination_data$pollination,pollination_data$name)
fmode
```
Ahora grafiquemos los datos con el arbol filogenetico
```{r}
plotTree(pollination_tree,fsize=0.4,ftype="i",type="fan")
cols<-setNames(c("blue","red"),sort(unique(fmode)))
tiplabels(pie=to.matrix(fmode,sort(unique(fmode))),piecol=cols,cex=0.5)
add.simmap.legend(colors=cols,prompt=FALSE,x=0.9*par()$usr[1],
y=-max(nodeHeights(pollination_tree)),fsize=0.8)
```
## Utilizando los modelos de Markov
El modelo mas sencillo asume que las tasas de transicion entre 0 y 1 son iguales. Vamos a ajustarlo en estos datos
```{r}
# la variable fitER se refiere a "fit equal rates model" ajustar el model con tasas iguales para dos estados
fitER<-ace(fmode,pollination_tree,model="ER",type="discrete",C=TRUE)
fitER
fitER$rates
```
### Reconstruccion ancestral en los nodos
Ya que tenemos un modelo ajustado con estimadores de maxima verosimilitud, podemos reconstruir el pasado utilizando verosimilitudes marginales en cada nodo. A este grafico se le llama arbol con reconstruccion ancestral
```{r}
plotTree(pollination_tree,type="fan",fsize=0.2,ftype="i")
nodelabels(node=1:pollination_tree$Nnode+Ntip(pollination_tree),
pie=fitER$lik.anc,piecol=cols,cex=0.5)
tiplabels(pie=to.matrix(fmode,sort(unique(fmode))),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=0.9*par()$usr[1],
y=-max(nodeHeights(pollination_tree)),fsize=0.5)
```
Tambien podemos hacer simulaciones condicionales para saber lo que sucedio en las ramas. A este grafico se le llama mapa estocastico de reconstruccion ancestral
```{r}
mtree<-make.simmap(pollination_tree,fmode,model="ER")
plotSimmap(mtree,cols,type="fan",fsize=0.3,ftype="i")
add.simmap.legend(colors=cols,prompt=FALSE,x=0.9*par()$usr[1],
y=-max(nodeHeights(pollination_tree)),fsize=0.8)
```
## Modelos de Markov con mas variabilidad
En macroevolucion no siempre podemos asumir que las transiciones entre dos estados (tasas evolutivas) ocurren al mismo paso. Ajustar un modelo mucho mas libre es en general una mejor opcion
```{r}
mk2<-corHMM(pollination_tree, pollination_data, rate.cat=1)
plotMKmodel(mk2)
```
Como se ve el mapa estocastico?
```{r}
phy= mk2$phy
data=mk2$data
model= mk2$solution
model[is.na(model)]=0
diag(model)= -rowSums(model)
simmap<-makeSimmap (tree=phy, data=data, model=model, rate.cat=1, nSim=1, nCores=1)
plotSimmap(simmap[[1]], fsize=0.2, type="fan")
```
### Modelos con estados escondidos
En evolucion, a traves de millones de años debemos pensar en la posibilidad de que las tasas de transicion ocurren a differentes pasos. Para esto utilizamos modelos de Markov de estados escondidos.
```{r}
mk2_hidden<-corHMM(pollination_tree, pollination_data, rate.cat=2,model="ARD")
plotMKmodel(mk2_hidden)
```
Como se ve el mapa estocastico?
```{r}
phy= mk2_hidden$phy
data=mk2_hidden$data
model= mk2_hidden$solution
model[is.na(model)]=0
diag(model)= -rowSums(model)
simmap<-makeSimmap (tree=phy, data=data, model=model, rate.cat=2, nSim=1, nCores=1)
plotSimmap(simmap[[1]], fsize=0.2, type="fan")
```