# Sauf mention contraire, ce contenu est mis à disposition selon les termes de la licence
# Creative Commons Attribution - Partage dans les mêmes conditions 4.0 International (CC BY-SA 4.0)
# ---------------------------------------------------------------------------
# Gaëlle Lelandais <gaelle.lelandais@universite-paris-saclay.fr>
# D'après une idée de Judith Legrand :).
# ---------------------------------------------------------------------------

#---------------------
# Objectifs de la séance
#---------------------
# Présenter un exemple de simulation d'un processus biologique 
# (ici la dérive génétique), en utilisant le logiciel R.
# Objectif de distinguer, dans les lignes de code informatique
# les parties qui sont MODIFIABLES de celles qui ne le sont pas. 
# Ainsi, les simulations pourront être réalisées, en changeant les paramètres
# tels que la fréquence du caractère d'intérêt dans la population, la taille
# de la population ou le nombre d'indididus prélevés à chaque génération.

#---------------------
# La dérive génétique
#---------------------
# Rappel : La dérive génétique est la modification aléatoire de la fréquence
# des allèles au sein d’une population au cours de générations
# successives.


#---------------------
# Rappels de la séance précédente, nous avons vu comment :
# 1) Créer une population de taille n, dans laquelle les individus
#    portent un caractère d'intérêt "rouge" avec une fréquence f.
# 2) Compter le nombre d'indididus de couleur "rouge" et réaliser une
#    représentation graphique de type "barplot".
# 3) Choisir au hasard N individus dans la population.
# 4) Recréer une population de taille n à partir de N individus.
# Nous allons maintenant réaliser une simulation de dérive génétique...
#---------------------

# Partie 1---------------------
# Calculer la fréquence des individus "rouge" après une
# première génération.
#---------------------

# On commence par créer la population initiale, celle ci est composée de
# n individus avec une fréquence f0 de "rouge".
n = 100
f0 = 0.7
pop0 = sample(c("rouge", "autre"), size = n, 
               replace = T, prob = c(f0, 1-f0))

# Un sous ensemble N d'individus est sélectionné aléatoirement, 
# ceux ci reforment une population de taille n.
N = 10
ech = sample(pop0, N, replace = F) 
pop = rep(ech, n/N)
# Attention : Pour que cela fonctionne, il faut n >>> N.

# La fréquence du caractère "rouge" est calculée 
# (dans la nouvelle population) :
nbR = sum(pop == "rouge")
f = nbR/n
print(f)

# Questions : 
# La fréquence du caractère "rouge" dans la nouvelle population est-elle
# la même que la fréquence de la population initiale (f0) ? Pourquoi ?
# La fréquence du caractère "rouge" dans la nouvelle population est-elle
# la même que celle de la personne assise à côté de vous ? Pourquoi ?

# Partie 2---------------------
# Suivre l'évolution de la fréquence du caractère rouge
# au cours de G générations successives. Représenter graphique le résultat.
#---------------------

# On a déjà calculé la fréquence pour la première génération.
freq = f

# On choisit un nombre de générations
G = 10

# On répète la partie 1 pour un nombre G de générations (utilisation d'une
# boucle "for").
for(i in 2:G){
  # sélection des N individus pour la génération suivante
  ech = sample(pop, N, replace = F) 
  # création d'une nouvelle population
  pop = rep(ech, n/N)
  # calcul de la fréquence des individus "rouges"
  nbR = sum(pop == "rouge")
  f = nbR/n
  # cette fréquence est conservée pour la suite
  freq = c(freq, f)
}

# Les valeurs des fréquences ont été gardées, au fur et à mesure 
# des générations 
print(freq)

# On peut réaliser une représentation graphique de type "plot".
plot(1:G, freq, 
     main = "Evolution de la fréquence du caractère rouge",
     xlab = "Génération",
     ylab = "Frequence du caractère",
     ylim = c(0, 1),
     type = "b")

# /!\/!\/!\/!\ 
# Travail à réaliser : 
# /!\/!\/!\/!\

# -- Réaliser cette simulation plusieurs fois. 
# Changer 1) la valeur de la fréquence in initiale f0, et 2) le nombre N
# d'individus utilisés pour passer d'une génération à l'autre.
# -- Adapter les titres des graphiques et réaliser des sauvegardes
# sous la forme de fichiers PDF.

# Partie 3---------------------
# Observation de l'extinction/installation 
# du caractère "rouge" dans la population.
#---------------------

# Si le nombre de générations est suffisamment grand, le caractère "rouge"
# va s'éteindre (f = 0) ou bien s'installer (f = 1). 
# Nous disposons de toutes les briques de code pour réaliser
# une simulation complète.

# Liste des paramètres
n  = 100
f0 = 0.7
N  = 50
G  = 100

# Population initiale
pop0 = sample(c("rouge", "autre"), size = n, 
              replace = T, prob = c(f0, 1-f0))

# Première génération
ech = sample(pop0, N, replace = F) 
pop = rep(ech, n/N)
nbR = sum(pop == "rouge")
f = nbR/n

# Pour suivre de la fréquence au cours des générations suivantes
freq = f

# Pour garder la génération d'extinction/stabilisation
date_ext = NULL

for(i in 2:G){
  ech = sample(pop, N, replace = F) 
  pop = rep(ech, n/N)
  nbR = sum(pop == "rouge")
  f = nbR/n
  freq = c(freq, f)
  
  # Si la fréquence est égale à 0 ou 1, il y a 
  # extinction/stabilisation
  if(((f == 0) | (f == 1)) & is.null(date_ext)){
    date_ext = i
  }
}

plot(1:G, freq, 
     main = "Evolution de la fréquence du caractère rouge",
     xlab = "Génération",
     ylab = "Frequence du caractère",
     ylim = c(0, 1),
     type = "b")

abline(v = date_ext, col = "red")  

# Partie 4---------------------
# Mise en fonction de la simulation 
#---------------------

# Création de la fonction
simulation_derive <- function(n, f0, N, G, graph = T){
  
  # Population initiale
  pop0 = sample(c("rouge", "autre"), size = n, 
                replace = T, prob = c(f0, 1-f0))
  
  # Première génération
  ech = sample(pop0, N, replace = F) 
  pop = rep(ech, n/N)
  nbR = sum(pop == "rouge")
  f = nbR/n
  
  # Suivi de la fréquence au cours des générations suivantes
  freq = f

  # Pour garder la date d'extinction/stabilisation
  date_ext = NA
  
  for(i in 2:G){
    ech = sample(pop, N, replace = F) 
    pop = rep(ech, n/N)
    nbR = sum(pop == "rouge")
    f = nbR/n
    freq = c(freq, f)
    
    # Si la fréquence est égale à 0 ou 1, il y a 
    # extinction/stabilisation
    if(((f == 0) | (f == 1)) & is.na(date_ext)){
      date_ext = i
    }
  }
  
  if(graph == T){
    plot(1:G, freq, 
         main = "Evolution de la fréquence du caractère rouge",
         xlab = "Génération",
         ylab = "Frequence du caractère",
         ylim = c(0, 1),
         type = "b")
    
    # Si l'extinction/stabilisation a eu lieu, elle est notée sur le graphique
    abline(v = date_ext, col = "red")
  }
  
  # la date d'extinction/satabilisation du caractère est renvoyée
  return(date_ext)
}

# Utilisation de la fonction
simulation_derive(n = 100, f0 = 0.5, N = 20, G = 100)

# Partie 5 ---------------------
# Etude de la date d'extinction, en fonction de f0 et N. 
# La taille de la population est augmentée à 10000 !
# Les simulations sont répétées 100 fois.
#---------------------

# Nombre de répétitions
Nsimul = 100

# Valeurs du paramètre testé
val = c(10, 50, 100, 500, 1000)

# Table pour garder les résultats
table_res = NULL

# Boucle for sur les différentes valeurs de paramètres
for(i in val){
  dates = NULL
  for(s in 1:Nsimul){
    dates = c(dates, 
              simulation_derive(n = 10000, f0 = 0.7, N = i, G = 100,
                                graph = F))
  }
  table_res = cbind(table_res, dates)
}

# Ajout des noms de colonnes
colnames(table_res) = val

# Graphique des résultats
boxplot(table_res, main = "Effet de la valeur de N, f0 = 0.7",
        ylab = "Date extinction")

# /!\/!\/!\/!\ 
# Travail à réaliser : 
# /!\/!\/!\/!\

# 1) Tester l'effet de N sur la date d'extinction, pour f0 = 0.4
# et f0 = 0.8 (nombre de génération G = 100)
# 2) Tester l'effet du nombre de générations G sur la date d'extinction, 
# pour f0 = 0.4 et f0 = 0.8
# 3) Réaliser des représentations graphiques, les exporter sous la forme d'un
#    fichier PDF.

# Bilan de la séance :
# 
# Au cours de cette séance, les logiciels R et Rstudio ont été utilisés.
# Leur utilisation vous a permis de réaliser une simulation d'un processus
# biologique. Vous avez pu tester les conséquences des changements des 
# valeurs de paramètres de la simulation, réaliser des représentation graphiques
# et sauvegarder vos résultats.
#
# L'étape suivant (si utile), consiste à étudier la programmation R. Des 
# ressources sont disponibles ici :
# https://gaellelelandais.github.io/teaching/initiation-R-base/ 
# ###################