---
title: "Práctico Filogenias 2"
author: "Laboratorio de Evolución"
date: "09/11/2020"
output:
  html_document: default
  word_document: default
  pdf_document: default
---
## Métodos Comparativos Filogenéticos


**Resumen del problema**

Es bien sabido que la presión parcial de oxígeno disminuye con la altura.  Las personas no habituadas "se apunan" o sienten malestar al trasladarse a zonas de montaña, especialmente al momento de realizar actividades físicas exigentes.  Un individuo que se traslada desde una zona baja hasta una de altura, experimenta a lo largo del tiempo una serie de ajustes fisiológicos para paliar al menos parte de esos desajustes.

Por lo tanto, cabe preguntarse si hay cambios genéticos que han sido favorecidos por la selección natural y que permiten a especies y poblaciones que viven a altas elevaciones adaptarse mejor a dichas condiciones.  Tanto en humanos como en especies animales (y muchas otras), hay estudios orientados a identificar estos cambios.

Puesto que la afinidad de la sangre por el oxígeno es un factor clave para la vida en altura, y dicha afinidad depende fuertemente de las características de la hemoglobina (recordemos que la estructura cuaternaria de la hemoglobina combina dos cadenas de tipo alfa y dos de tipo beta en un tetrámero, en torno a un núcleo de hierro), esta proteína ha sido el blanco de muchos estudios.  

En un trabajo reciente, Natarajan et al. (2016) se plantean identificar algunas de las adaptaciones de la hemoglobina para la vida en la altura.  La hipótesis de trabajo es que la selección natural pudo haber favorecido cambios en las características de la hemoglobina de las especies asociados a la elevación en la que vive cada una. 

Algunas de las ideas del artículo son:
- Realizar un estudio comparando múltiples especies de aves, procurando elegir pares de especies cercanamente relacionadas, de modo que una de las especies de cada par viva en tierras altas y otra en las tierras bajas cercanas.
- Para cada una de estas especies, aislar la hemoglobina y estudiar su afinidad con el oxígeno en el laboratorio.

El estudio incluye un análisis de los cambios en las secuencias de las hemoglobinas, que usaremos más adelante en el curso.  Por el momento, extraemos del artículo los siguientes datos para cada una de las 56 especies estudiadas:

1. La elevación en la que viven (específicamente la altura de la localidad en la que fueron estudiadas).
2. La afinidad de su alfa globina por el oxígeno, tomada en condiciones controladas de laboratorio (en un medio que aproxima las condiciones en sangre), resumida por el valor conocido como P50.
3. El árbol filogenético que se utiliza en el artículo.  Este árbol se toma como la mejor hipótesis disponible de las relaciones entre las especies.

**Preguntas 1**: Obtener de internet una curva de afinidad de la hemoglobina con el oxígeno.  ¿Qué es el valor P50?  ¿Cómo debería cambiar el P50 para que la hemoglobina tenga mayor afinidad con el oxígeno, por ejemplo en una especie de altura?

**Activando los paquetes de R**

Para esta actividad práctica usaremos, **R** (https://cran.r-project.org/) el cual es un entorno y lenguaje de programación enfocado al análisis estadístico. 
Se recomienda el uso **RStudio** (https://rstudio.com/), que es un entorno de desarrollo integrado (IDE), en conjunción con **R**, lo cual facilita el manejo de datos y la realización de los análisis. 

Utilizaremos los siguientes paquetes de **R**, los cuales deben ser previamente instalados:

_**R Markdown**_: para formatear el archivo de salida de los análisis.  
_**knitr**_: para convertir los archivos de R markdown en varios formatos.  
_**phytools**_: varias funciones para análisis filogenéticos, principalmente orientado a la biología comparada.  

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library("markdown")
library("knitr")
library("phytools")
```

**Trabajando con los datos originales**

Mediante el comando "setwd(dir)" podemos indicarle a R en donde están ubicados los archivos con los cuales vamos a estar trabajando, siendo "dir" la ruta hacia nuestro directorio (e.g. C:\\User\\Desktop\\R).También podemos ubicar el directorio de trabajo de R mediante "getwd()" y colocar ahí nuestros archivos de interés.

El bloque de código siguiente lee solamente los datos de elevación y los de P50 de la HbA. Se hace una exploración de esos datos.

```{r}
# Leemos las dos tablas de datos y combinarlas en un solo "data frame".
# Notar que en esta parte son archivos nuevos con los nombres de las variables cambiados
  P50<-read.table("HbA_KCLHIP_P50.txt")
  elevacion<-read.table("Elevacion_2.txt", header = TRUE)
  Datos = cbind(elevacion, P50) # acá se combinan los datos de las dos tablas
  Datos$Elev.cat <- c("Alta", "Baja") # agrego variable categórica Elevación a cada par de spp. 
  table(Datos$Elev.cat) # chequeo que quedan 28 pares   
  print(Datos)
  rm(P50, elevacion)
```

**Preguntas 2**:  La tabla de datos incluye varios pares de especies de un mismo género. Elegir algunos de esos pares para discutir:

1. ¿De qué especies se trata?  Averiguar algo de los nombres comunes, familias a las que pertenecen, distribución geográfica.
2.¿Qué tendencias se observan al examinar en varios de esos pares la relación entre altura y P50?    

**Examinando la relación entre elevación y P50**

En el siguiente bloque se obtiene el coeficiente de correlación entre las dos variables de interés, se aplica un modelo de regresión lineal entre dichas variables, y se grafican los valores junto con la línea de tendencia obtenida en la regresión.
```{r}
attach(Datos) # Fijamos las Variables de Datos (ya no es necesario llamarlas a través del comando Datos$x)

# Análisis exploratorios
str(Datos) # Resumen de número y tipo de variables y nro de observaciones 
Datos$Elevacion <- as.numeric(as.character(Datos$Elevacion))

summary(Familia) # cuantas observaciones de cada Familia
summary(P50)
hist(P50)
summary(Elevacion)
hist(Elevacion)

# Exploro relación P50 vs. Elevación (categórica)  
boxplot(P50~Elev.cat)

# ...y la relación P50 vs. Elevación (numérica)
# Calculamos el coeficiente de correlación entre las dos variables
  R = cor(Datos$Elevacion,Datos$P50)
  
# Aplicar un modelo lineal (lm) de regresión. Se define a P50 como la variable dependiente y a Elevacion como la variable independiente. 
  Regresion <- lm(Datos$P50~Datos$Elevacion)
  summary.lm(Regresion)

# Graficamos los valores de P50 y elevación, agregando la línea de ajuste obtenida de la regresión lineal. 
  plot(Datos$P50~Datos$Elevacion, main = "Relación entre P50 de HbA y elevación")
  lines(Datos$Elevacion, fitted(Regresion), col="blue")
    
```

**Preguntas 3** 
1. ¿Cuál es el valor del coeficiente de correlación?  Qué sugiere respecto a la relación entre las variables (signo, intensidad)?
2. Examinar la tabla de regresión.  ¿Qué fracción de la varianza está explicada por el modelo? 
3. Identificar y explicar los parámetros relevantes del modelo y el p-valor asociado a cada uno.  ¿Qué sugieren estos resultados sobre la hipótesis de trabajo.
4. Examinar la gráfica, incluyendo la predicción obtenida por regresión lineal. ¿Qué sugiere la gráfica y cómo se relaciona con los puntos discutidos más arriba?

Las especies, y por lo tanto sus rasgos, no pueden ser tratadas como variables independientes desde el punto de vista estadístico dado que comparten una historia evolutiva en común.Este problema se hace aún más evidente en taxa cercanamente emparentados. La respuesta a tal problema fue propuesta por Felsenstein (1985) mediante el cálculo de Contrastes Filogenéticos Independientes (CFI), mediante el cual se transforma a los rasgos analizados en variables independientes, empleando la filogenia de las especies como marco de análisis. En la siguiente sección, vamos a abordar este problema y analizaremos nuevamente la correlación entre la P50 y la Altura, pero esta vez corrigiendo la falta de independencia de los datos, utilizando los CFI como nuevas variables, independientes de la historia evolutiva compartida entre las diferentes especies de aves.


```{r}

#leemos el árbol filogenético
tree<-read.tree("tree.tre")

#enraizamos el árbol en su punto medio
treeR<-midpoint.root(tree)

#visualizamos el árbol enraizado
plotTree(treeR, edge.width=1, ftype="i", fsize=0.7)

#para simplificar los comandos, leemos otro archivo que contiene solo las especies y las máximas alturas registradas, así como otro archivo con los valores de P50 (se puede corroborar y si no se modificó para la etapa anterior, se puede usar el mismo en ambas partes)
elevacion_2<-read.table("Elevacion.txt")
P50_2<-read.table("HbA_KCLHIP_P50.txt") 

#Creamos un vector numérico de las variables para calcular los contrastes
VP50 <- P50_2$P50
Velevacion <- elevacion_2$Elevacion

#Para que los datos y los taxa terminales se asocien correctamente
names(VP50) <- row.names(P50_2)
names(Velevacion) <- row.names(elevacion_2)

#calculamos los contrastes filogenéticamente independientes
ContrasteVP50 <- pic(VP50, treeR)
ContrasteVelevacion <- pic(Velevacion, treeR)

#si queremos extraer además de los contrastes, su varianza asociada 
ContrasteVP50.var <- pic(VP50, treeR, var.contrasts=TRUE)
ContrasteVelevacion.var <- pic(Velevacion, treeR, var.contrasts=TRUE)

#Ahora visualizamos los contrastes filogenéticos calculados anteriormente
ContrasteVP50
ContrasteVelevacion

# vemos los contrastes filogenéticos calculados y también su varianza asociada
ContrasteVP50.var
ContrasteVelevacion.var

#vemos los contrastes asociados a cada nodo, a partir del archivo que también contiene el cálculo de las varianzas. Los contrastes asociados a P50 aparecerán en azul y los asociados a la Elevación aparecerán en verde.

plot(treeR, no.margin=TRUE,edge.width=1 ,cex=c(0.6,0.6))
nodelabels(round(ContrasteVelevacion.var [,1], 3), adj = c(1, -0.4), frame="n", col = "blue", cex=0.6)

#volvemos a plotear la filogenia solo para que nos se superpongan ambos contrastes
plot(treeR, no.margin=TRUE, edge.width=1 ,cex=c(0.6,0.6))
nodelabels(round (ContrasteVP50.var [,1], 3), adj = c(1, 1.1), frame="n", cex=0.6, col = "darkgreen") 
```

```{r}
#Analizamos la evolución correlacionada de ambos rasgos a partir de los archivos sin el cálculo de varianzas
RegresionP50_elevacion <- lm(ContrasteVP50~ContrasteVelevacion)

#Visualizamos los parámetros de la regresión
summary.lm(RegresionP50_elevacion)

#visualizamos la relación entre los contrastes y agregamos una línea de tendenica a la regresión
plot(ContrasteVP50~ContrasteVelevacion, main = "Relación entre contrastes de P50 de HbA y de elevación")
lines(ContrasteVelevacion, fitted(RegresionP50_elevacion), col="blue")
```

```{r}
#Veamos ahora cómo se distribuyen los caracteres sobre la filogenia
#primero definimos el nombre de los datos
P50_2<-setNames(P50_2[,1],rownames(P50_2))
elevacion_2<-setNames(elevacion_2[,1],rownames(elevacion_2))

#y luego obtenemos su distribución sobre la filogenia propuesta
#para el valor de P50
Dist_P50<-contMap(treeR,P50_2,fsize=c(0.5,1),outline=FALSE)

#modificamos la leyenda del gráfico
plot(Dist_P50,fsize=c(0.5,1),outline=FALSE,lwd=c(3,7),leg.txt="P50")

#para las altitudes máximas registradas para cada taxón
Dist_elevacion<-contMap(treeR,elevacion_2,fsize=c(0.5,1),outline=FALSE)

#nuevamente, modificamos la leyenda del gráfico
plot(Dist_elevacion,fsize=c(0.5,1),outline=FALSE,lwd=c(3,7),leg.txt="Elevacion")

```


**Pregunta 4 ** 

¿Según los resultados obtenidos, la historia filogenética compartida de las especies parece haber influido en la adaptación de la hemoglobina a la hipoxia?