--- word_document: fig_caption: yes author: "Laboratorio de Evolución" date: "21/08/2024" output: pdf_document: fig_caption: yes word_document: default title: "PCM 2" html_document: default editor_options: markdown: wrap: sentence --- **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** () el cual es un entorno y lenguaje de programación enfocado al análisis estadístico. Se recomienda el uso **RStudio** (), 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: ***tidyverse***: una colección de paquetes que facilitan el análisis de datos. ***broom***: para formatear salida de modelos de ajuste de estos datos ***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") library("tinytex") library("tidyverse") library("formattable") ``` **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} # para usar directorio donde se encuentra script library(rstudioapi) current_path = rstudioapi::getActiveDocumentContext()$path setwd(dirname(current_path )) getwd() # activamos paquetes library("phytools") library("tidyverse") library("formattable") # Leemos los datos, tenemos pares de spp. en distinta elevacion (categoria) # más la elevación ocupada + P50 de cada spp Datos <- read_tsv("Datos.tsv") # notar cuantas columnas tengo y que tipo de variable es (categorica y numérica) Datos # veo rapidamente los datos Datos %>% count(Elev.cat) # un conteo de ocurrencias en la variable categoría de elevación ``` **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} # Análisis exploratorios class(Datos) # que tipo de objeto es? Puede ser vector, list, matrix, data frame, etc... # Datos %>% formattable() # imprimo en pantalla toda la tabla formateada summary(Datos) # Resumen de número y tipo de variables y nro de observaciones # cual es la media de P50 y Elevación en este dataset summarise(Datos, mean(P50)) # veo un estadístico de una variable en particular summarise(Datos, mean(Elevacion)) # veo un estadístico de una variable en particular Datos %>% count(Familia) # cuantas observaciones de cada Familia Datos %>% arrange(desc(Elevacion)) # ordenamos de mayor a menor, por Elevación Datos %>% arrange(P50) # ordenamos x P50 # Datos %>% slice_max(Elevacion) # siempre puedo buscar ciertos valores # Datos %>% filter(Elevacion >= 3000) # siempre puedo filtrando con una condición # ahora podemos encadenar alguno de estos comandos # ¿cual es la media en las dos categorñías de elevacion? Datos %>% group_by(Elev.cat) %>% summarise(mean(Elevacion)) Datos %>% select(Especie, P50) # asi puedo ver facilmente dos columnas hist(Datos$Elevacion) hist(Datos$P50) # dentro de una familia de interes Datos %>% filter(Familia == "Trochilidae") %>% summarise(mean(P50)) Datos %>% filter(Familia == "Trochilidae") %>% group_by(Elev.cat) %>% summarise(mean(Elevacion)) Datos %>% filter(Familia == "Trochilidae") %>% group_by(Elev.cat) %>% summarise(mean(P50)) # en dos familias Datos %>% filter(Familia == "Anatidae" | Familia == "Trochilidae" ) %>% ggplot( aes(Elevacion, P50, group=Familia)) + geom_point(aes(color= Familia)) + theme_minimal() # Exploro relación P50 vs. Elevación (categórica) ggplot(Datos, aes(Elev.cat, P50)) + # boxplt geom_boxplot(outlier.colour = "red", outlier.shape = 1) + # vemos los outliers en rojo theme_classic() + ggtitle("P50 versus Elevación") + xlab("Elevación") ggplot(Datos, aes(Elev.cat, P50)) + geom_boxplot(outlier.colour = "red", outlier.shape = 2) + geom_jitter(width = 0.5) + # hacemos esto solo para ver mejor todos los valores theme_classic() + ggtitle("P50 versus Elevación") + xlab("Elevación") # Ahora vemos 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. ggplot(Datos, aes(Elevacion, P50)) + geom_point() + geom_smooth(method='lm') + theme_classic() + ggtitle("Relación entre P50 de HbA y elevación") Regresion <- lm(Datos$P50~Datos$Elevacion) summary.lm(Regresion) # como alternativa podemos obtener una tabla ya formateada o lista # broom::tidy(Regresion) %>% formattable() # tabla formateada # broom::glance(summary(Regresion)) %>% %>% formattable() # tabla formateada # Graficamos los valores de P50 y elevación, agregando la línea de ajuste # ...obtenida de la regresión lineal. # ``` Datos %>% filter(Familia == "Anatidae" | Familia == "Trochilidae" ) %>% ggplot( aes(Elevacion, P50, group=Familia)) + geom_point(aes(color= Familia)) + geom_smooth(method='lm') + theme_classic() + ggtitle("Relación entre P50 de HbA y elevación ") **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, # y leemos otro archivo con los valores de P50 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 #Paso necesario para que los datos y los taxa terminales se asocien correctamente names(VP50) <- row.names(P50_2) names(Velevacion) <- row.names(elevacion_2) enframe(VP50) # acá corroboramos como asociamos nombre spp. ~ valor enframe(Velevacion) # para la otra variable #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 enframe(ContrasteVP50) # aca lo veo como 'tabla' #ContrasteVelevacion enframe(ContrasteVelevacion) # aca lo veo como 'tabla' # vemos los contrastes filogenéticos calculados y también su varianza asociada # ContrasteVP50.var # ContrasteVelevacion.var enframe(ContrasteVP50.var) # la varianza P50 enframe(ContrasteVelevacion.var) # la varianza Elevacion # 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) broom::tidy(summary(RegresionP50_elevacion)) # %>% formattable() # tabla formateada broom::glance(summary(RegresionP50_elevacion)) # %>% formattable() # tabla formateada broom::glance(summary(Regresion)) # %>% formattable() # tabla formateada # broom::glance(summary(Regresion)) %>% formattable() # tabla anterior #visualizamos la relación entre los contrastes y agregamos una línea de tendenica a la regresión # Contraste <- bind_cols(enframe(ContrasteVP50), enframe(ContrasteVelevacion)) %>% rename(CVP50 = `value...2`, CElevacion = `value...4`) ggplot(Contraste, aes(CElevacion, CVP50)) + geom_point() + geom_smooth(method='lm') + theme_classic() + ggtitle("Relación entre contrastes de P50 de HbA y de elevación") ``` ```{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) #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") #modificamos la leyenda del gráfico plot(Dist_P50,fsize=c(0.5,1),outline=FALSE,lwd=c(3,7),leg.txt="P50") ``` **Pregunta 4** ¿Según los resultados obtenidos, la historia filogenética compartida de las especies parece haber influido en la evolucion de P50? ¿La inercia filogenetica (o tendencia en la historia del rasgo P50) parece haber impedido la adapatacion en algunos pares de spp?