script.Rmd

---
title: "ÉVOLUTION MOLÉCULAIRE COMPARÉE DES EXONS NUCLÉAIRES ORTHOLOGUES DE PRIMATES ET DE RONGEURS"
author: "Allyson MOUREAUX, Elian STROZYK et Mattéo TRAISSAC-MONTAHUT"
lang: fr
output:
  pdf_document: 
    toc: true
  html_document: default
date: "2022-12-15"
geometry: margin=2.5cm
fontsize: 12pt
classoption: a4
header-includes:
    - \usepackage{setspace}\onehalfspacing
    - \usepackage{lscape}
    - \newcommand{\blandscape}{\begin{landscape}}
    - \newcommand{\elandscape}{\end{landscape}}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, message=FALSE, warning=FALSE)
```

```{=tex}
\newpage
\onehalfspacing
```

# Introduction

Des arbres de probabilité ont été extraits de la plateforme OrthoMaM. Ces arbres de probabilité ont été construits à partir d'alignements de séquences exoniques similaires de plusieurs mammifères. Ce projet a pour but d'étudier la vitesse d'évolution des séquences exoniques ayant un pourcentage de G+C~3~ similaire. Comme la souris et l'Homme ont divergé depuis 90 Ma et que le pourcentage de G+C~3~ présente des variations importantes en général nous nous attendons à ce que la vitesse d'évolution soit différente. Nous souhaitons également voir si le pourcentage de G+C~3~ influence le taux d'évolution.

# Estimation du taux d'évolution pour chacune des six branches

Nous importons nos données dans une liste et créons des vecteurs de noms d'espèces pour récupérer les identifiants des n\oe{}uds à partir de `getMRCA()`.

```{r}
library(ape)
library(tidyr)

#Importation des arbres qui se trouvent dans ./trees
#création d'une liste pour les stocker
temp = list.files(path="./trees/",
                  pattern = ".rootree$"
                  )
setwd("./trees")
trees = lapply(temp, read.tree)

# Creation de vecteurs correspondant à chaque espèce et son taxon
R = "Rattus_norvegicus"; S = "Mus_musculus";
H = "Homo_sapiens"; M = "Macaca_fascicularis"

species = c(R, S, H, M); rodents = c(R, S)
primates = c(H, M); mouse = c(S, S)
human = c(H, H); macaque = c(M, M)
rat = c(R, R)
# getMRCA() a besoin d'un vecteur contenant au moins 2 éléments
```

Chaque espèce a un identifiant de n\oe{}ud différent selon l'arbre. Nous récupèrons donc ces identifiants que nous stockons dans des vecteurs. Egalement, nous créons une liste qui va comporter la matrice des distances entre chaque MRCA pour chaque arbre.

```{r}
# Récupération des MRCA (ID des noeuds) pour chaque
# arbre et pour chaque espèce/ancêtre
mrca_species = c(); mrca_rodents = c()
mrca_primates = c(); mrca_human = c()
mrca_mouse = c(); mrca_rat = c()
mrca_macaque = c()
tab_dist = c()

# Récupération de la matrice des distances pour chaque arbre
for (i in 1:length(trees)) {
  mrca_species[i] = getMRCA(trees[[i]], species)
  mrca_rodents[i] = getMRCA(trees[[i]], rodents)
  mrca_primates[i] = getMRCA(trees[[i]], primates)
  mrca_human[i] = getMRCA(trees[[i]], human)
  mrca_mouse[i] = getMRCA(trees[[i]], mouse)
  mrca_macaque[i] = getMRCA(trees[[i]], macaque)
  mrca_rat[i] = getMRCA(trees[[i]], rat)
  tab_dist[[i]] = dist.nodes(trees[[i]]) # matrice de distance
}
```

A partir des distances de matrice, nous récupérons les données de distance avec les identifiants des n\oe{}uds qui nous intéressent. Pour récupérer le taux d'évolution, il suffit de diviser la distance entre deux MRCA par le temps depuis l'événement de spéciation, ces temps ont été récupérés à partir de l'arbre schématique du sujet.

```{r}
#Création de vecteurs vides qui vont contenir les distances entre deux MRCA
# qui vont ensuite être convertis en taux d'évolution
r_rodents.mouse = c(); r_rodents.rat = c()
r_primates.human = c(); r_primates.macaque = c()
r_primates.species = c(); r_rodents.species = c()

# Récupération des distances entre 2 noeuds
# Les distances entre noeuds suivent la relation de Chasles
# puis division par le temps depuis l'événement de spéciation
for (i in 1:length(mrca_human)) {
  r_rodents.mouse[i] = (tab_dist[[i]][mrca_rodents[i], mrca_mouse[i]])/15
  r_rodents.rat[i] = (tab_dist[[i]][mrca_rodents[i], mrca_rat[i]])/15
  r_primates.human[i] = (tab_dist[[i]][mrca_primates[i], mrca_human[i]])/30
  r_primates.macaque[i] = (tab_dist[[i]][mrca_primates[i], mrca_macaque[i]])/30
  r_primates.species[i] = (tab_dist[[i]][mrca_primates[i], mrca_species[i]])/60
  r_rodents.species[i] = (tab_dist[[i]][mrca_rodents[i], mrca_species[i]])/75
} 
```

Voir **annexe 1** pour la distribution des taux d'évolution des six espèces sous forme de boîtes à moustache.

# Distribution des taux d'évolution

## Médiane, quantiles à 5% et 95% de chacune des six branches

```{r}
#création des variables qui vont contenir la médiane et les quantiles
taux_devolution = list("RS" = r_rodents.mouse, "RR" = r_rodents.rat, 
                       "RH" = r_primates.human, "RM" = r_primates.macaque, 
                       "RPRI" = r_primates.species, "RROD" = r_rodents.species
                       )
#création de vecteurs vides
mediane = c(); quantile5 = c(); quantile95 = c()

#application de la fonction médiane quantiles5 et quantiles95
#à la liste taux_devolution qui contient tous les taux d'évolutions
mediane = unname(unlist(lapply(taux_devolution, median)))
quantile5 = unname(unlist(lapply(taux_devolution, quantile, p = c(0.05))))
quantile95 = unname(unlist(lapply(taux_devolution, quantile, p = c(0.95))))

#création d'un dataframe résumant les informations
tab_r = data.frame(names(taux_devolution), quantile5, mediane, quantile95)
knitr::kable(tab_r)
```

## Histogramme des taux d'évolution de l'humain et de la souris

```{r, out.width="75%", fig.cap="Histogramme du taux d'évolution pour \\textit{Homo sapiens} et \\textit{Mus musculus}.", fig.align="center"}
library(ggplot2); library(reshape2); library(dplyr)

# Création d'un tableau large indicé qui est transformé en tableau long
tab_hist = data.frame(r_rodents.mouse, r_primates.human, 
                      id = 1:length(r_rodents.mouse))
tab_hist = melt(tab_hist, id.vars="id")

# Histogramme représentant le taux d'évolution de GC3 chez l'Homme et la souris
ggplot(tab_hist, aes(x=value, color=variable)) +
  geom_histogram(aes(y=..count..), fill="white", alpha=0.7,
                 position="dodge", bins=30)  +
  theme_bw() + 
  xlab(expression(Taux~"d'évolution"~("Nombre de substitutions par site" %.%Ma^{-1}))) +   
  scale_color_manual(labels = c(expression(italic(Mus~musculus)), 
                                expression(italic(Homo~sapiens))),
                     values = c("blue", "red"))
```

Nous observons un taux d'évolution moyen plus important chez la souris que chez l'Homme **(figure 1)**. Nous pouvons émettre l'hypothèse que le pourcentage de G+C~3~ chez l'Homme évolue moins que celui de la souris puisque ces mutations seraient délétères. Les mutations à cet endroit de l'exon sont certainement contre seléctionnées chez l'Homme tandis que chez la souris elles seraient bénéfiques et se fixeraient donc avec plus de probabilité. De plus, l'Homme ayant un temps de génération plus long que la souris le taux d'évolution serait donc plus lent. Nous savons que l'horloge moléculaire peut varier entre lignée surtout si les temps de générations sont différents. Nous observons un taux d’évolution intra-espèces hétérogènes. Cela peut potentiellement être dû à des mutations qui sont liées aux conditions environnementales de l’espèce pour assurer sa survie.

# Synthèse et perspectives

Les quantiles et les médianes pour la souris et l'Homme entre les différents groupes semblent tous concorder puisqu'ils ont le même ordre de grandeur (**annexe 3**). Pour l'Homme, la dispersion semble bien plus faible que chez la souris pour laquelle les quantiles sont plus éloignés de la médiane que chez l'Homme.

D'après les tendances (**annexe 3**), nous observons que lorsque le pourcentage moyen de G+C~3~ augmente chez l'Homme le taux d'évolution augmente. Le contraire est observé chez la souris. Mais d'après le R², il n'y a pas de corrélation entre le taux d'évolution et le pourcentage de G+C~3~ étant donné qu'il est trop proche de 0 (**annexe 3**).

Pour vérifier si cela est généralisable il faudrait plus de réplicats. Ici nous n'en avons que 13. Il faudrait également comparer entre d'autre espèces avec des temps de divergence différents pour confirmer nos conclusions.

\newpage

# Annexes

## Annexe 1 : taux d'évolution pour différentes espèces

```{r, echo=FALSE, fig.cap="Boîtes à moustache des taux d'évolution de différentes espèces depuis leur ancêtre d'intérêt."}

# ANNEXE3 : facilitation de visualisation 

#Taux d'évolution pour chaque espèce
library(ggplot2)
library(reshape2)
tab_visu = data.frame(r_rodents.mouse, r_rodents.rat, r_rodents.species,
                      r_primates.human, r_primates.macaque, r_primates.species)
tab_visu = melt(tab_visu) # transformation d'un tableau large en tableau long (individu-type)

ggplot(tab_visu, aes(x=variable,y=value,fill=variable))+
  geom_boxplot(alpha = 0.7, outlier.shape = NA) +
  theme_bw() + 
  ylab(expression(Taux~"d'évolution"~("Nombre de substitutions par site" %.%Ma^{-1}))) +   
  geom_point(position=position_jitterdodge(), aes(colour = variable), size = 0.9) +
  scale_x_discrete(labels=expression("Rodentia-Mus musculus", "Rodentia-Rattus norvegicus", "Rodentia-Euarchontoglires", "Primates-Homo sapiens", "Primates-Macaca fascicularis", "Primates-Euarchontoglires")) +
  theme(legend.position="none") +
  theme(axis.title.y=element_blank()) +
  coord_flip()
```

\newpage

## Annexe 2 : tableau récapitulatif des résultats pour différents pourcentages de G+C~3~

![Tableau récapitulatif](tableau_de_ses_morts.png)

\newpage

## Annexe 3 : nuage de points du tableau récapitulatif des résultats

```{r, echo=FALSE, fig.cap="Nuage de points représentant le pourcentage G+C\\textsubscript{3} en fonction du taux d'évolution chez \\textit{Mus musculus} et \\textit{Homo sapiens}."}
# ANNEXE3 : facilitation de visualisation 
#graphe commun

library(ggplot2)
library(ggpubr)
library(ggpmisc)

dt = read.csv(file='csv.Rfriendly_complet.csv', sep=",", dec=".")
dt$espece=factor(dt$espece)
dt$quantile=factor(dt$quantile)

ggplot(data = dt, aes(x = centre_classe, y = valeur)) +
  geom_point(mapping = aes(x = centre_classe, y = valeur, color = quantile,
                           shape = espece), size = 2.5, alpha=0.8) +
  theme_bw() + 
  ylab(expression(Taux~"d'évolution"~("Nombre de substitutions par site" %.%Ma^{-1}))) +   
  xlab(expression(G+C[3]~moyen~'(%)')) +
  scale_shape_manual(labels = c(expression(italic(Homo~sapiens)), expression(italic(Mus~musculus))), values = c('circle', 'triangle')) +
  scale_color_discrete(labels=c("5%", "Médiane", "95%")) +
  labs(color="Quantile", shape="Espèce") +
  stat_poly_line(color="black") +
  stat_poly_eq(aes(label.y=0.5, label = after_stat(rr.label))) +
  facet_grid(~ espece)
  #   geom_smooth(method = "lm", se = FALSE, aes(color=quantile))
```

La droite de régression linéaire a été générée en poolant les trois quantiles mais pas les espèces.