################################################################################
### TALLER FUNDAMENTOS DE R                                                  ###
### EJERCICIO 7.2: Sequencias aleatoreas (y in poco sobre gráficos)          ###
###                                                                          ###
### Center for Conservation and Sustainable Development                      ###
### Missouri Botanical Garden                                                ###
### Website: rbasicsworkshop.weebly.com                                      ###
################################################################################

### INTRODUCCIÓN ###############################################################
# La evolución de caracteres continuos se describe comúnmente mediante el modelo
# Browniano de evolución (Felsenstein, J. 1985. Phylogenies and the comparative 
# method. Am.Nat. 125:1-15; Felsenstein, J. 1988. Phylogenies and quantitative 
# characters. Annu. Rev. Ecol. Syst. 19:445-471). En este modelo evolutivo, 
# cambios en un caracter phenotipico durante un intervalo de tiempo dado es 
# determinado por una distribución normal con media cero y varianza proporcional 
# a la tasa de evolución.


## TAREA 1 ##
# Basado en el modelo Browniano de evolución, genere valores de 1000 cambios 
# evolutivos con una tasa de evolución = 1. Utilice la función "rnorm" y 
# almacene los valores de cambio evolutivo en un vector llamado 
# "change.character.1".

change.character.1 <- rnorm(1000, mean=0, sd=1)


## TAREA 2 ##
# Examine los valores en "change.character.1" usando la función "hist".

hist(change.character.1)


## TAREA 3 ##
# Arriba se generó valores de cambio evolutivo en un caracter que fueron 
# guardados en "change.character.1". Ahora va a generar los valores del
# caracter a través del tiempo. Para ello, asuma que el carácter ancestral
# tiene un valor de cero, y utilice la función "cumSum" para acumular el cambio 
# evolutivo a través del tiempo. Si no ha utilizado la función "cumSum" 
# anteriormente, vea la respectiva página de ayuda. Utilice la función "plot" 
# para examinar visualmente los valores del carácter a través del tiempo. Si no 
# recuerda el uso de la función "plot", vea la página de ayuda de "plot" 
# escribiendo en la consola "?plot" o "help(plot)."

phenotype.1 <- c(0, cumsum(change.character.1))
plot(phenotype.1, ylab="Phenotype", type="l", xlab="Evolutionary time", cex.lab=1.5, cex.axis=1.5)


## TAREA 4 ##
# Usando el modelo de la evolución Browniano, genere 1.000 valores de cambio 
# evolutivo con una tasa de evolución = 2. Utilice la función "rnorm" y almacene 
# los valores en un vector llamado "change.character.2".


## TAREA 5 ##
# Examine los valores en "change.character.2" usando la función "hist".


## TAREA 6 ##
# Compare "change.character.1" y "change.character.2" usando la función "summary".


## TAREA 7 ##
# La gigura 1 de Revell et al. (2008. Phylogenetic Signal, Evolutionary Process, 
# and Rate. Systematic Biology 57: 591-601; disponible en el sitio web de taller)
# muestra una filogenia de tres especies (A, B y C), y simulaciones de 
# evolución de un carácter a lo largo de una filogenia de acuerdo con un modelo 
# Browniano y con una tasa de evolución = 1. Considere la evolución del caracter
# para las especies A.
# Primero, defina el tiempo evolutivo que, de acuerdo con la Fig. 1 en Revell 
# et al., es igual a 2 unidades de tiempo (con las unidades no especificadas, 
# así que esto podría ser de 2 millones de años). Luego, construya un vector 
# usando la función "seq", y pongalo en un objecto llamado "evolution.time":

evolution.time <- seq(0.001, 2, 0.001)

# Ahora, defina la tasa de evolución en un objeto denominado "evolution.rate".
# Según Revell et al., la tasa de evolución por unidad de tiempo es 1. En el 
# paso anterior dividimos cada unidad de tiempo en 1000 partes, por lo que la 
# tasa de evolución que vamos a usar es:

evolution.rate <- 1/1000

# Genere valores de cambio fenotípico para la especie A en la Fig. 1 de Revell 
# et al. La especie A evoluciona durante 2 unidades de tiempo. Dividimos cada 
# unidad de tiempo en 1000 partes, por lo que necesitamos para generar los 
# valores de 2000 de cambios evolutivos de la especie A:

change.phenotype.A <- rnorm(2000, mean = 0, sd = evolution.rate^0.5)

# Genere valores de fenotipo para la especie A (suponiendo un fenotipo ancestral 
# = 0), y use la función "plot" para inspeccionar visualmente los valores:

phenotype.A <- cumsum(change.phenotype.A)
plot(evolution.time, phenotype.A, type="l")

# Tenga en cuenta el uso de argumento "type" de la función "plolt" Lea 
# información sobre este argumento en la página de ayuda.


## TAREA 8 ##
# Genere valores de cambio fenotípico para la especie B en la Figura 1 de Revell 
# et al., y almacene los valores en un vector llamado "change.phenotype.B".


## TAREA 9 ##
# Genere valores fenotípicos para las especies de B en la Figura 1 de Revell et 
# al., y almacene los resultados en un vector llamado "phenotype.B". Supongamos 
# que el fenotipo ancestral = 0.


## TAREA 10 ##
# Utilice el código de abajo para inspeccionar visualmente la evolución de 
# fenotipos para las especies A y B como en la Fig. 1 de Revell et al. Tenga en 
# cuenta el uso de argumentos "type", "col" y "ylab". En la página de ayuda para 
# la función "plot" leer la información sobre argumentos "type" y "col". Se dará 
# cuenta de que no hay información sobre el argumento
# "col". Sin embargo, en la sección de "argumentos" (en la página de ayuda para 
# la función "plot") hay tres puntos ("...") que indica que hay más argumentos 
# en la página de ayuda para la función de "par". Ir a la página de ayuda para
# función "par" y leer la información sobre el argumento "col". Tenga en cuenta 
# también el argumento "cex" en la página de ayuda para la función de "points".

plot(evolution.time, phenotype.A, type="l", col="gray70", ylab="Phenotype") 
points(0,0, col="green", pch=19, cex=1.5) #ancestral phenotype
points(evolution.time, phenotype.B, type="l", col="red")


## TAREA 11 ##
# Puede haber un problema con la gráfica. Puede verlo? Al menos en algunas 
# iteraciones de este código, los límites de los ejes verticales no incluyen 
# todos los valores de la serie temporal de fenotipos de la especies B. Use el
# código abajo para resolver este problema.

# Note el uso de la función "cbind". Lea la página de ayuda para la función "cbind"
phenotypes.A.B <- cbind(phenotype.A, phenotype.B) 
phenotypes.A.B[1:5,] #examine the first five rows of "phenotypes.A.B".
dim(phenotypes.A.B)
max(phenotypes.A.B)
min(phenotypes.A.B)

plot(evolution.time, phenotype.A, type="l", col="gray70", ylab="Phenotypes", 
    ylim=c(min(phenotypes.A.B), max(phenotypes.A.B)))
points(0,0, col="green", pch=19, cex=1.5)
points(evolution.time, phenotype.B, type="l", col="red")
text(evolution.time[length(phenotype.A)], phenotype.A[length(phenotype.A)], labels="A", cex=3) #adds species name at the end of phenotypic series
text(evolution.time[length(phenotype.B)], phenotype.B[length(phenotype.B)], labels="B", cex=3) #adds species name at the end of phenotypic series


## TAREA 12 ##
# Genere valores para el cambio evolutivo de la especie C en la Fig. 1 de Revell 
# et al., y almacene los valores en un vector llamado "change.phenotype.C". Es
# importante darse cuenta que las especies B y C sólo han evolucionado de forma 
# independiente solamente durante una unidad de tiempo. Examine cuidadosamente 
# lafilogenia en la Figura 1 de Revell et al.

change.phenotype.C <- c(change.phenotype.B[1:1000], rnorm(1000, mean = 0, sd = evolution.rate^0.5))


## TAREA 13 ##
# Genere valores fenotípicos para la especie C en la Fig. 1 de Revell et al., y 
# la pongalos en un vector llamado "phenotype.C". Suponga que el fenotipo 
# ancestral = 0.


## TAREA 14 ##
# Use y conside cuidadosamente el siguiente código para inspeccionar visualmente 
# los valores fenotípicos para la especie A, B y C en la figura. 1 de Revell et 
# al.:

phenotypes.A.B.C <- cbind(phenotype.A, phenotype.B, phenotype.C)
phenotypes.A.B.C[1:5,]
dim(phenotypes.A.B.C)
max(phenotypes.A.B.C)
min(phenotypes.A.B.C)

plot(evolution.time, phenotype.A, type="l", col="gray70", ylab="phenotype", ylim=c(min(phenotypes.A.B.C), max(phenotypes.A.B.C)))
points(0,0, col="green", pch=19, cex=1.5)
points(evolution.time, phenotype.B, type="l", col="red")
points(evolution.time, phenotype.C, type="l", col="blue")
text(evolution.time[length(phenotype.A)], phenotype.A[length(phenotype.A)], labels="A", cex=3) 
text(evolution.time[length(phenotype.B)], phenotype.B[length(phenotype.B)], labels="B", cex=3) 
text(evolution.time[length(phenotype.C)], phenotype.C[length(phenotype.C)], labels="C", cex=3)
#ADD una línea vertical en el momento del evento de especiación separar especies B y C, como en la Fig. 1b de Revell et al ..
abline(v=1, lty=3) 




















################################################################################
### SOLUCIONES PARA TAREAS #####################################################
################################################################################

## TAREA 1 ##
# Ya tiene las respuestas


## TAREA 2 ##
# Ya tiene las respuestas


## TAREA 3 ##
# Ya tiene las respuestas


## TAREA 4 ##
change.character.2 <- rnorm(1000, mean=0, sd=2)


## TAREA 5 ##
hist(change.character.2)


## TAREA 6 ##
summary(change.character.1)
summary(change.character.2)


## TAREA 7 ##
# Ya tiene las respuestas


## TAREA 8 ##
change.phenotype.B <- rnorm(2000, mean = 0, sd = evolution.rate^0.5)


## TAREA 9 ##
phenotype.B <- cumsum(change.phenotype.B)


## TAREA 10 ##
# Ya tiene las respuestas


## TAREA 11 ##
# Ya tiene las respuestas


## TAREA 12 ##
# Ya tiene las respuestas


## TAREA 13 ##
phenotype.C <- cumsum(change.phenotype.C)


## TAREA 14 ##
# Ya tiene las respuestas