# PACKAGES
library(tidyverse)

# LOAD DATA FROM SQL DATA BASE
load("~/marmot-quantitative-genetics/1_DATA_ANALYSIS/3_Data/statut_dominance.rds")  # table statut social
load("~/marmot-quantitative-genetics/1_DATA_ANALYSIS/3_Data/pedigree_data.rds")  # pédigrée (pour portée)
load("~/marmot-quantitative-genetics/1_DATA_ANALYSIS/3_Data/pheno.rds")  # morpho
load("~/marmot-quantitative-genetics/1_DATA_ANALYSIS/3_Data/individu_bdd.rds")  # informations individuelles
load("~/marmot-quantitative-genetics/1_DATA_ANALYSIS/3_Data/territoire_ref.rds")  # infos territoire

# AJOUT SEXE (à ajouter avec BDD)
index_temp <- match(statut_dominance$individual_id, individu_bdd$individu_id)
statut_dominance$sexe <- individu_bdd$sexe[index_temp]
head(statut_dominance)
##   individual_id territoire_id statut_social_id type_attribution annee territoire sexe
## 1            15           749                2          deduite  1989          c    f
## 2            16           749                2          deduite  1989          c    m
## 3             1             2                2          capture  1990          b    f
## 4            11             1                2          capture  1990          a    m
## 5            12             1                2          capture  1990          a    m
## 6            13             1                2          capture  1990          a    m


COMPLETER LES INDIVIDUS DOMINANTS MANQUANTS

Les données de présence de la table SQL statut_dominance (contenant les présences par territoire de chaque indiviu, via capture ou observation) sont probablement incomplètes et quelques données peuvent probablement être récupérées. Je vois 2 cas de figures où je pourrais potentiellement récupérer de l’information :


Reconstruction avec les informations du pédigrée

# FILTRER LES DOMINANCE
statut_dominance_filtered <- 
  statut_dominance |> 
  filter(!duplicated(paste0(individual_id, "_", territoire_id, "_", statut_social_id))) |> 
  arrange(individual_id, territoire_id) |> 
  filter(statut_social_id == 1)

# CREER UN ID UNIQUE TERRITOIRE_ID (i.e. annee x territoire) x INDIVIDU ID
statut_key <- paste0(statut_dominance_filtered$individual_id, "_", statut_dominance_filtered$territoire_id)

# MERES MANQUANTE PAR TERRITOIRE ID
meres_manquantes <-
  pedigree_data |>
  filter(!is.na(mother)) |>
  mutate(mere_territoire_key = paste0(mother, "_", territoire_id)) |>
  filter(!mere_territoire_key %in% statut_key) |>
  select(individual_id = mother, territoire_id) |>
  distinct()

# PERES MANQUANTE PAR TERRITOIRE ID
peres_manquants <- pedigree_data |>
  filter(!is.na(father)) |>
  mutate(pere_territoire_key = paste0(father, "_", territoire_id)) |>
  filter(!pere_territoire_key %in% statut_key) |>
  select(individual_id = father, territoire_id) |>
  distinct()

# COMBINER LES PARENTS MANQUANTS
parents_manquants <- rbind(meres_manquantes, peres_manquants) |>
  distinct()

# GET LES INFORMATIONS DU TERRITOIRE
parents_manquants <- parents_manquants |>
  left_join(territoire_ref, by = "territoire_id")

# GET SEXE
sexe_ref <- statut_dominance |>
  filter(!duplicated(individual_id)) |>
  select(individual_id, sexe) |> 
  mutate(individual_id = as.character(individual_id))

parents_manquants <- parents_manquants |>
  left_join(sexe_ref, by = "individual_id")

# NOUVELLES LIGNES A AJOUTER
nouvelles_lignes_pedigree <-
  parents_manquants |>
  filter(!is.na(territoire_id), !is.na(annee_territoire)) |>
  mutate(
    statut_social_id = 1,
    type_attribution = "deduite_pedigree",
    annee = annee_territoire
  ) |>
  select(individual_id, territoire_id, statut_social_id, type_attribution, annee, territoire, sexe)


# AJOUTER AUX DATAFRAME
statut_dominance <- rbind(
  statut_dominance,
  nouvelles_lignes_pedigree
) |>
  filter(!duplicated(paste0(individual_id, "_", territoire_id, "_", statut_social_id))) |>
  arrange(individual_id, territoire_id)

Vérification et statistiques

  • Nombre de lignes (individual_id × territoire_id) ajoutées : 117

  • Nombre de territoire_id uniques concernés: 84



Reconstruction avec les séquences \(t\) à \(t+x\)

# OBTENIR PAR TERRITOIRE, LES SERIES TEMPORELLES DE CHAQUE INDIVIDU
statut_dominance_filtered <- 
  statut_dominance_filtered |> 
  mutate(key = paste0(individual_id, territoire, annee))

# GET GUESSABLE SEQUENCES
statut_dominance_filtered_completed <- 
  statut_dominance_filtered |> 
  group_by(individual_id, territoire) |> 
  reframe(annee = min(annee):max(annee),
          key = paste0(individual_id, territoire, annee)) |> 
  ungroup() |> 
  filter(!key %in% statut_dominance_filtered$key) |> 
  select(-key)

# MERGE LES HISTOIRES RECONSTRUITES AU DATAFRAME
df_temp <- data.frame(
  individual_id = statut_dominance_filtered_completed$individual_id,
  territoire_id = NA,
  statut_social_id = 1,
  type_attribution = "deduite_sequence",
  annee = statut_dominance_filtered_completed$annee,
  territoire = statut_dominance_filtered_completed$territoire,
  sexe = NA
) 

statut_dominance_completed <- rbind(statut_dominance, df_temp)

# COMPLETER LA COLONNE TERRITOIRE ID
index_temp <- match(
  x = paste0(statut_dominance_completed$territoire, statut_dominance_completed$annee), 
  table = paste0(territoire_ref$territoire, territoire_ref$annee)
)
statut_dominance_completed$territoire_id <- territoire_ref$territoire_id[index_temp]

# REFILTRER LES DUPLICATIONS
statut_dominance <- 
  statut_dominance_completed |> 
  filter(!duplicated(paste0(individual_id, "_", territoire_id, "_", statut_social_id))) |> 
  arrange(individual_id, territoire_id) 

# AJOUT SEXE
index_temp <- match(statut_dominance$individual_id, individu_bdd$individu_id)
statut_dominance$sexe <- individu_bdd$sexe[index_temp]


Vérification et statistiques

  • Nombre de lignes (individual_id × territoire_id) ajoutées : 60

  • Nombre de territoire_id uniques concernés: 47



STATISTIQUES DESCRIPTIVES

On commence par extraire le couple de dominant par territoire par année depuis 1990.

# FILTER DATAFRAME TO KEEP ONLY DOMINANT
individu_dominant <- 
  statut_dominance |>
  filter(statut_social_id %in% c(1)) |>
  arrange(territoire, annee) |> 
  group_by(territoire, annee) |>
  mutate(individual_id = list(unique(individual_id)),
         sexe = list(sexe))

# GET COUPLE
individu_dominant$couple <- sapply(individu_dominant$individual_id, paste, collapse = "-")

# QUICK VIEW
head(individu_dominant)
## # A tibble: 6 × 8
## # Groups:   territoire, annee [3]
##   individual_id territoire_id statut_social_id type_attribution annee territoire sexe      couple
##   <list>                <int>            <dbl> <chr>            <chr> <chr>      <list>    <chr> 
## 1 <chr [1]>                 1                1 capture          1990  a          <chr [1]> 7     
## 2 <chr [3]>                 6                1 capture          1991  a          <chr [3]> 32-4-7
## 3 <chr [3]>                 6                1 deduite          1991  a          <chr [3]> 32-4-7
## 4 <chr [3]>                 6                1 capture          1991  a          <chr [3]> 32-4-7
## 5 <chr [2]>                15                1 deduite          1992  a          <chr [2]> 32-4  
## 6 <chr [2]>                15                1 capture          1992  a          <chr [2]> 32-4


Nombre de couple unique

length(unique(individu_dominant$couple))  
## [1] 463


Nombre de territoires x années suivi


Nombre de dominant connu par territoire

Total

individu_dominant |>
  filter(!duplicated(territoire_id)) |>
  reframe(table(n_dom))
##   table(n_dom)
## 1          115
## 2          651
## 3           95
## 4           12
## 5            1


Par années



IDENTIFIER LES CAS POTENTIELS D’INFANTICIDE

L’infanticide chez la marmotte Alpine

Généralités sur l’infanticide chez la marmotte

  • Changement de dominance mâle comme femelles entraine l’infanticide
  • Tous les marmottons et les yearlings meurent en cas de changement de dominance
  • Plus la saison avance, moins les changements de dominances sont fréquents (Rebecca)


Occurrence de l’infanticide

  • Avant l’arrivée sur le terrain (de l’émergence au début de terrain)
  • Pendant le terrain
  • Après le terrain (de l’après terrain à l’entrée en hibernation)



Quantifier les changements de dominance

Dans un premier temps, nous devons caractériser les changements de dominance dont dépendent totalement les occurrences d’infanticide. D’un point de vue opérationnel, je considérerai un changement de dominance lorsque :

  • Pour une année \(t\), le nombre de dominant est supérieur à 3
  • Le couple de dominant à \(t\) est différent de celui à \(t-1\)
  • Le couple de dominant passe de 1 à 2 (i.e. un des deux n’étaient pas connus)


NOTE DE CODE
Il faut gérer les cas de figures comme à la ligne 8-9 car si “coexistence” de e.g. >3 dominants dans les données (i.e. dominance a été détectée pendant le terrain), alors automatiquement l’année \(t+1\) sera > considérée comme “changement de dominance” alors que le couple était déjà formée à \(t\).

Je rajoute donc la condition “si à l’année \(t\) il y a >3 dominants (n_dom > 3) ET que le couple à \(t+1\) est un sous-ensemble de celui de \(t\), alors ce n’est pas un changement”. Pour cela, j’utilise le package stringr.



Changement de dominance

# COMPUTE LE CHANGEMENT DE DOMINANCE PAR ANNEE PAR TERRITOIRE
individu_dominant <- 
  individu_dominant |>
  filter(!duplicated(territoire_id)) |>
  complete(territoire, annee, fill = list(NA)) |>   # on complète les années manquantes par NAs sinon génère artificellement des changements de dominances 
  arrange(territoire, annee) |> 
  group_by(territoire) |>
  mutate(
    cur = strsplit(couple, "-"),
    before = strsplit(lag(couple), "-"),
    is_na = is.na(couple),
    
    # Tous les dominants de t+1 étaient déjà présents à t
    inclus = !is.na(lag(couple)) &
             mapply(function(x, y) all(y %in% x), before, cur),
    
    changement_dom = if_else(is_na, NA, if_else(couple != lag(couple) &
                                                  !inclus, T, F, missing = F))
  ) |>
  select(-cur, -before, -inclus, -is_na) |> 
  as.data.frame()


La somme totale de changement de dominance :

## [1] 283


Par années



Changement de dominance par sexe

# IDENTIFIER LE TYPE DE CHANGEMENT DE DOMINANCE

# INFORMATIONS PAR SEXE PAR CHANGEMENT DE DOMINANCE
individu_dominant_sexe <- 
  individu_dominant |>
  arrange(territoire, annee) |>
  group_by(territoire) |>
  mutate(
    # Calculer AVANT de passer en rowwise
    sexes_couple_avant = lag(sexe),
    individual_id_avant = lag(individual_id)
  ) |>
  ungroup()

# IDENTIFIER TYPE DE CHANGEMENT (couple entier, male, femelle)
individu_dominant_sexe <- 
  individu_dominant_sexe |>
  rowwise() |>  # permet d'appliquer le mutate par ligne
  mutate(
    type_changement = if_else(
      changement_dom == TRUE,
      {
        # Extraire les individus et sexes avant et après
        ids_avant = unlist(individual_id_avant)
        ids_apres = unlist(individual_id)
        sexes_av = unlist(sexes_couple_avant)
        sexes_ap = unlist(sexe)
        
        # Gérer le cas où ids_avant est NULL (première année)
        if (is.null(ids_avant) || length(ids_avant) == 0) {
          NA_character_
        } else {
          # Séparer par sexe
          males_avant = ids_avant[sexes_av == "m"]
          males_apres = ids_apres[sexes_ap == "m"]
          femelles_avant = ids_avant[sexes_av == "f"]
          femelles_apres = ids_apres[sexes_ap == "f"]
          
          # Vérifier si chaque sexe a changé
          changement_male = !all(males_apres %in% males_avant)
          changement_femelle = !all(femelles_apres %in% femelles_avant)
          
          # Classifier le type
          if (changement_male & changement_femelle) {
            "couple_entier"
          } else if (changement_male) {
            "male"
          } else if (changement_femelle) {
            "femelle"
          } else {
            NA_character_
          }
        }
      },
      NA_character_
    )
  ) |>
  ungroup()
## # A tibble: 3 × 4
##   type_changement     n proportion pourcentage
##   <chr>           <int>      <dbl>       <dbl>
## 1 couple_entier      59      0.211        21.1
## 2 femelle            95      0.339        33.9
## 3 male              126      0.45         45



Quantifier l’infanticide

Détecter l’infanticide

  • Cas 1 : Présence de juvéniles + observation directe de changement de dominance
  • Cas 2 : Présence de marmottons à \(t\) + changement de dominance à \(t+1\)
  • Cas 3 : Présence de yearling à \(t\) + changement de dominance à \(t+1\)


Potentiels effets confondants

  • Absence de juvéniles (marmottons + yearlings) meurent pendant l’hibernation + changement de dominance pré-terrain
  • Portée manquée + changement de dominance pré-terrain
  • Bruce effect / pregnancy block - infanticide pendant la gestation par la propre mère en cas de changement de mâle dominant


On considérera donc comme potentiel infanticide les changements de dominance associés à la présence de marmotton ou de yearlings.


Présence de marmottons

Je vais d’abord construire la variable presence_pup par année*territoire.

# CONSTRUCT VARIABLE <PRESENCE OF LITTER>
# ... Si le territoire ID est dans le pédigrée, alors une portée cette année...
individu_dominant$presence_pup <-
  ifelse(individu_dominant$territoire_id %in% pedigree_data$territoire_id,
         T,
         F)


Présence de yearlings

Maintenant, on doit ajouter la variable presence_yearling, aussi victime de l’infanticide en cas de changement de dominance.

# IDENTIFIER LES YEARLINGS PAR ANNEE PAR TERRITOIRE
# GET EMERGENCE DATE
statut_dominance <- 
  statut_dominance |>
  full_join(pheno[c("individual_id", "age_obs", "annee_emergence")],
            by = "individual_id",
            multiple = "first") |> 
  rename(annee_emergence_pheno = annee_emergence)

statut_dominance <- 
  statut_dominance |> 
  full_join(pedigree_data[c("individual_id", "cohort")],
            by = "individual_id",
            multiple = "first") |> 
  rename(annee_emergence_pedigree = cohort)

statut_dominance <- statut_dominance |> 
  mutate(
    annee_emergence = coalesce(  # find the first non-missing value
      annee_emergence_pedigree,
      annee_emergence_pheno
    )
  ) |> 
  select(-c(annee_emergence_pheno, annee_emergence_pedigree))

# COMPUTE AGE
statut_dominance$age_annee <- 
  as.numeric(as.character(statut_dominance$annee)) -   as.numeric(as.character(statut_dominance$annee_emergence)) +
  statut_dominance$age_obs

# GET ALL YEARLINGS
yearling_df <- statut_dominance |> 
  filter(age_annee == 1) |> 
  distinct(territoire_id)

# ADD THE VARIABLE presence_yearling
yearling_df <- statut_dominance |>
  filter(age_annee == 1) |>
  distinct(individual_id, .keep_all = TRUE)

individu_dominant$presence_yearling <- 
  individu_dominant$territoire_id %in% yearling_df$territoire_id

# ON NE GARDE QU'UNE MESURE PAR TERRITOIRE ID
dominance_per_territory <- 
  individu_dominant |> 
  filter(!duplicated(territoire_id))


Présence de juvéniles

Je construis la variable presence_juvenile qui regroupe la présence de marmottons OU de yearlings.

dominance_per_territory$presence_juvenile <- dominance_per_territory$presence_pup | dominance_per_territory$presence_yearling


Potentiel infanticide

Maintenant, on construit la variable potential_infanticide où les potentiels infanticides sont tous les cas de figure où un changement de dominance coincide avec la présence de juvéniles.

dominance_per_territory$potential_infanticide <-
  dominance_per_territory$changement_dom & (dominance_per_territory$presence_juvenile)


Statistiques descriptives

Nombre total de territoire-année avec au moins un marmotton

##   n_presence_pup prop_presence_pup
## 1            631         0.7211429


Nombre total de territoire-année avec au moins un yearling

##   n_presence_pup prop_presence_pup
## 1            338         0.3862857


Nombre total de territoire-année avec au moins un juvénile

##   n_presence_pup prop_presence_pup
## 1            687         0.7851429


Nombre total de territoire-année avec potentiel infanticide

##   n_presence_pup prop_presence_pup
## 1            198         0.2262857