• ¿Cómo cambió el sindrome de polinización a traves de la evolución?

• ¿En que momento surgió el canto de las hembras en aves?

• ¿Que capacidades cognitivas estuvieron presentes en el ancestro común de los primates?

Estimación del estado ancestral de un rasgo a partir de su expresion en las especies existentes

• "Paleontología estadística" (Pagel 1997)
• Escencial para el pensamiento evolutivo
• Permite determinar la historia evolutiva de rasgos

Permite inferir:

• Cuántas veces evolucionó un rasgo
• Qué condiciones propiciaron su evolución
• En qué linajes evolucionó el rasgo

• No es una construcción del árbol filogenético
• Necesita:
  • un árbol previamente construido
  • datos sobre el valor de un rasgo en las especies actuales

Limitaciones

• Solo puede reconstruir valores que estén dentro del rango actual del rasgo
• Asume patrones evolutivos en los datos (i.e. depende de modelos)
• Es afectada por la exclusión de especies (filogenias incompletas)

para rasgos discretos (i.e. sexo, syndrome de polinización)

• Parsimonia
• Máxima verosimilitud
• Métodos bayesianos

o continuos (i.e. tamaño corporal)

• Máxima verosimilitud
• Constrastes filogenéticos independientes

Arbol de tribu Arini (Psitácidos) + datos morfologicos/habitat

library(ape)
library(phytools)

#cargar datos 
arini.datos.comp <-
readRDS(gzcon(url("http://marceloarayasalas.weebly.com/uploads/2/5/5/2/25524573/arini.datos.comp.rds")))

#definir arbol y datos
arb.arini <- arini.datos.comp[[1]]
datos.arini <- arini.datos.comp[[2]]
row.names(datos.arini) <- datos.arini$Species

#empatar datos y arbol
#install.packages("picante")
library(picante)
arini <- match.phylo.data(arb.arini, data = datos.arini[,-1])

Fuente: Medina-García, A., M. Araya-Salas, and T. F. Wright. "Does vocal learning accelerate acoustic diversification? Evolution of contact calls in Neotropical parrots."Journal of Evolutionary Biology 28.10 (2015): 1782-1792.

## Loading required package: maps

## 
##  # maps v3.1: updated 'world': all lakes moved to separate new #
##  # 'lakes' database. Type '?world' or 'news(package="maps")'.  #

## Loading required package: vegan

## Loading required package: permute

## Loading required package: lattice

## This is vegan 2.3-5

## Loading required package: nlme

Arbol de tribu Arini (Psitácidos) + datos morfologicos/habitat

head(arini$data, 3)

##                            Duartion.s Kurtosis Spectral.Entropy
## Anodorhynchus_hyacinthinus  0.5121600 28.53519        0.9002164
## Anodorhynchus_leari         0.3843178 31.78945        0.8734356
## Aratinga_auricapillus       0.3504207 24.49549        0.9136763
##                            Fund.freq.KHz Peak.freq.KHz Freq.range
## Anodorhynchus_hyacinthinus     0.9781165      1.801403   2.021715
## Anodorhynchus_leari            0.6286564      1.389516   2.077279
## Aratinga_auricapillus          1.5487699      3.978212   3.059284
##                            Peak.time.s Vegetation.index.EVI Body.Mass.g
## Anodorhynchus_hyacinthinus   0.3415090            0.5740557        1331
## Anodorhynchus_leari          0.4170942            0.3746777         940
## Aratinga_auricapillus        0.3882459            0.5242632         130
##                            Bill.length.mm
## Anodorhynchus_hyacinthinus          89.10
## Anodorhynchus_leari                 68.50
## Aratinga_auricapillus               23.75

Arbol de tribu Arini

plot.phylo(arini$phy, no.margin = TRUE, cex = 0.6)

Arbol de tribu Arini + tamaño corporal

#unir nombre de especies y datos
tc <- arini$data$Body.Mass.g
names(tc) <- row.names(arini$data)

#arbol con datos mapeados
dotTree(arini$phy, log(tc), standardize = T, 
        fsize = 0.6)

Arbol de tribu Arini + tamaño corporal

Arbol de tribu Arini + tamaño corporal y frecuencia (tono) del canto

#extraer rasgo
ff <- datos.arini$Fund.freq.KHz

#arbol con datos mapeados
dotTree(arb.arini, cbind(log(tc), log(ff)), standardize = T, 
        fsize = 0.6)

Arbol de tribu Arini + tamaño corporal y frecuencia (tono) del canto

library(stats4)

set.seed(1001)

x <- rnorm(100, mean = 3, sd = 2)

mean(x)

## [1] 2.998305

sd(x)

## [1] 2.288979

LL <- function(mu, sigma) {
R = suppressWarnings(dnorm(x, mu, sigma))

-sum(log(R))
}

a <- mle(LL, start = list(mu = 1, sigma=1))

b <- mle(LL, start = list(mu = 6, sigma=1))

#resultados
summary(a)

## Maximum likelihood estimation
## 
## Call:
## mle(minuslogl = LL, start = list(mu = 1, sigma = 1))
## 
## Coefficients:
##       Estimate Std. Error
## mu    2.998305  0.2277506
## sigma 2.277506  0.1610438
## 
## -2 log L: 448.4039

#resultados
summary(b)

## Warning in sqrt(diag(object@vcov)): NaNs produced

## Maximum likelihood estimation
## 
## Call:
## mle(minuslogl = LL, start = list(mu = 6, sigma = 1))
## 
## Coefficients:
##        Estimate Std. Error
## mu    -45.57025   20.64938
## sigma 227.30254        NaN
## 
## -2 log L: 1273.62

AIC: criterio de información de Akaike

#AIC
AIC(a)

## [1] 452.4039

AIC(b)

## [1] 1277.62

Máxima verosimilitud

recons.tc <- ace(x = tc, phy = arini$phy, type="continuous",
                        model="BM", method="ML")

## Warning in sqrt(diag(solve(h))): NaNs produced

#opcion en phytools 
recons.tc2 <- anc.ML(x = tc, tree = arini$phy, model = "BM")

names(recons.tc)

## [1] "loglik" "ace"    "sigma2" "CI95"   "call"

recons.tc$ace

##         52         53         54         55         56         57 
##  302.69177  305.46116  308.56286  356.15472 1103.81236  317.85829 
##         58         59         60         61         62         63 
##  254.04993  146.20062  128.44232  313.60348  320.22699  429.15766 
##         64         65         66         67         68         69 
##  523.08365  996.05995 1034.28206 1179.45780  654.80821  257.25683 
##         70         71         72         73         74         75 
##  305.94344  298.24280  212.18114  197.76131  185.02759  165.77928 
##         76         77         78         79         80         81 
##  296.36116  239.17990  204.69842  190.07049  161.21915  145.03591 
##         82         83         84         85         86         87 
##  129.68731  207.30883  154.87416  122.13610  122.95280   66.98481 
##         88         89         90         91         92         93 
##   85.32799   59.33394  116.41653  203.25818   70.99488  299.00488 
##         94         95         96         97         98         99 
##  304.55897  384.51017  178.55032  277.37492  240.28174  203.43312 
##        100        101 
##  275.77550  138.70191

Máxima verosimilitud

phenogram(arini$phy,log(c(tc,recons.tc$ace)), fsize = 0.5, 
          ylab = "log(tamaño corporal)")

contMap(tree = arini$phy, x = tc, plot=T, fsize = 0.5)

Máxima verosimilitud

## Optimizing the positions of the tip labels...

Máxima verosimilitud

Parsimonia

Ejemplo:

Algoritmo de optimización con 2 reglas aplicadas de forma descendente y ascendente:

• REGLA 1: si los clados hermanos (o especies) comparten el estado de un rasgo asignar ese estado al ancestro.
• REGLA 2: si los clados hermanos no comparten el mismo estado del rasgo asignar la union de los diferentes estados.

Parsimonia

Descendente:

Fuente: Cunningham, Clifford W., Kevin E. Omland, and Todd H. Oakley. "Reconstructing ancestral character states: a critical reappraisal."Trends in Ecology & Evolution 13.9 (1998): 361-366.

Parsimonia

Ascendente:

Parsimonia

• Optimización final

Crear rasgo discreto (Tamaño corporal)

#Converter en discreto
tc.dis <- cut(log(tc), breaks = 3)

#renombrar estados
levels(tc.dis) <- c("peq", "med", "gran")

#añadir nombres de especies
names(tc.dis) <- row.names(arini$data)

#definir colores
colors<-setNames(c("blue","red", "green"),c("peq","med", "gran"))

#mapear sobre arbol
dotTree(arini$phy, tc.dis, colors = colors, fsize = 0.5)

Crear rasgo discreto (Tamaño corporal)

• La evolución de rasgos discretos implica transiciones entre los estados del rasgo
• Las tasas de las transiciones se representa en un matriz (de transiciones)

• Cada matriz puede representar un modelo distinto
• El ajuste de los modelos puede ser comparado

• Matrices comunes
  • Tasas iguales (ER)

• Matrices comunes
  • Todas diferentes (ARD)

• Matrices comunes
  • Meristico

Comparación de modelos

# install.packages("geiger")
library(geiger)

#Crear una tabla para guardar resultados
disc.res <- data.frame(model=c("ER","ARD","simetrico", "meristico", "usuario"),
                      AICc=numeric(4),params=numeric(4))
# modelo ER
ER_mod <- fitDiscrete(phy = arini$phy,dat =  tc.dis, model="ER")

# guardar resultados
disc.res[1,-1]<- c(AICc=ER_mod$opt$aicc,ER_mod$opt$k)

# modelo ARD
ARD_mod <- fitDiscrete(phy = arini$phy,dat =  tc.dis, model="ARD")

# guardar resultados
disc.res[2,-1]<- c(AICc=ARD_mod$opt$aicc,ARD_mod$opt$k)

Comparación de modelos

# modelo simetrico
SYM_mod <- fitDiscrete(phy = arini$phy, tc.dis,
                         model= "SYM")

disc.res[3,-1] <- c(AICc=SYM_mod$opt$aicc, SYM_mod$opt$k)

#modelo meristico
tc.ord <- as.numeric(tc.dis) 

names(tc.ord) <- names(tc.dis)

meris_mod <- fitDiscrete(phy = arini$phy, tc.ord,
                         model= "meristic")

disc.res[4,-1]<- c(AICc=meris_mod$opt$aicc, meris_mod$opt$k)

Comparación de modelos

#Resultados ordenados por AICc
disc.res[order(disc.res$AICc),]

##       model     AICc params
## 4 meristico 64.82801      2
## 3 simetrico 67.08865      3
## 1        ER 68.03036      1
## 2       ARD 71.41346      6

Graficar reconstrucción

#simular 100 veces
arini.disc100 <- make.simmap(arini$phy, tc.dis, nsim=100)

# graficar arboles
par(mfrow=c(10,10))
null <- sapply(arini.disc100,plot,colors=colors,lwd=1,ftype="off")

## make.simmap is sampling character histories conditioned on the transition matrix
## 
## Q =
##              peq          med         gran
## peq  -0.01750374  0.017503736  0.000000000
## med   0.01750374 -0.024063670  0.006559934
## gran  0.00000000  0.006559934 -0.006559934
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       peq       med      gran 
## 0.3333333 0.3333333 0.3333333

## Done.

Graficar reconstrucción

Graficar reconstrucción

par(mfrow=c(1,1))

arini100 <- summary(arini.disc100,plot=T)

plot(arini100,fsize=0.6,ftype="i", colors = colors)

Graficar reconstrucción

Cambios de estados durante evolucion

arini100

## 100 trees with a mapped discrete character with states:
##  gran, med, peq 
## 
## trees have 12.07 changes between states on average
## 
## changes are of the following types:
##      gran,med gran,peq med,gran med,peq peq,gran peq,med
## x->y     0.39        0     2.19    6.85        0    2.64
## 
## mean total time spent in each state is:
##              peq         med       gran    total
## raw  215.9181289 319.9970559 67.7563332 603.6715
## prop   0.3576749   0.5300847  0.1122404   1.0000

1. Estime los valores ancestrales del largo del pico ("Bill.length.mm") de los llamados de Arini y grafíquelos.

2. Estime los valores ancestrales de la densidad de vegetación ("Vegetation.index.EVI") como un rasgo discreto (primero debe convertirlo a un rasgo discreto). Usando 2 tipos (modelos) de matrices de transición.

lp <- arini$data$Bill.length.mm
names(lp) <- row.names(arini$data)

#reconstruccion
recons.lp <- ace(x = lp, phy = arini$phy, type="continuous",
                        model="BM", method="ML")

#grafico
recons.tc2 <- anc.ML(x = tc, tree = arini$phy, model = "BM")

contMap(tree = arini$phy, x = lp, plot=T, fsize = 0.5)

## Warning in sqrt(diag(solve(h))): NaNs produced

#Converter en discreto
veg.dis <- cut(log(arini$data$Vegetation.index.EVI), breaks = 3)

#renombrar estados
levels(veg.dis) <- c("abierto", "medio", "cerrado")

#añadir nombres de especies
names(veg.dis) <- row.names(arini$data)

#definir colores
colors<-setNames(c("blue","red", "green"),c("abierto", "medio", "cerrado"))

#mapear sobre arbol
dotTree(arini$phy, veg.dis, colors = colors, fsize = 0.5)

#Crear una tabla para guardar resultados
# modelo ER
ER_mod <- fitDiscrete(phy = arini$phy,dat =  veg.dis, model="ER")

# modelo ARD
ARD_mod <- fitDiscrete(phy = arini$phy,dat =  veg.dis, model="ARD")

#AICc
cbind(ARD_mod$opt,ER_mod$opt)

##        [,1]         [,2]      
## q12    2.051529e-06 0.01508599
## q13    7.298082     0.01508599
## q21    0.03780435   0.01508599
## q23    1.713935e-10 0.01508599
## q31    1.722051     0.01508599
## q32    0.1222678    0.01508599
## lnL    -35.62071    -44.79298 
## method "subplex"    "Brent"   
## k      6            1         
## aic    83.24142     91.58597  
## aicc   85.15051     91.6676