Purpose: Analyse des profils génétiques 2016 à 2022
Author: Pierre-Alexandre QUITTET
Contact: pierre-alexandre.quittet@cefe.cnrs.fr
Code created: 15 septembre 2025
Last updated: 07 octobre 2025

1 SESSION INFORMATION

sessionInfo()

## R version 4.5.0 (2025-04-11 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26100)
## 
## Matrix products: default
##   LAPACK version 3.12.1
## 
## locale:
## [1] LC_COLLATE=French_France.utf8  LC_CTYPE=French_France.utf8   
## [3] LC_MONETARY=French_France.utf8 LC_NUMERIC=C                  
## [5] LC_TIME=French_France.utf8    
## 
## time zone: Europe/Paris
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## loaded via a namespace (and not attached):
##  [1] digest_0.6.37     R6_2.6.1          fastmap_1.2.0     xfun_0.52        
##  [5] cachem_1.1.0      knitr_1.50        htmltools_0.5.8.1 rmarkdown_2.29   
##  [9] lifecycle_1.0.4   cli_3.6.5         sass_0.4.10       jquerylib_0.1.4  
## [13] compiler_4.5.0    rstudioapi_0.17.1 tools_4.5.0       evaluate_1.0.4   
## [17] bslib_0.9.0       yaml_2.3.10       rlang_1.1.6       jsonlite_2.0.0

1.1 RANDOM INFORMATION

Bibl31 a foiré pour l’analyse de 2023_PK5, j’ai tout mis en NA
Bibl1 a foiré pour l’analyse de 2020, j’ai tout mis en NA

2 PACKAGE, DATA AND FUNCTIONS

# PACKAGES
library(RODBC)  # connection to DBMS
library(tidyverse)
library(xlsx)

# IMPORT DATA
dir = "C:/Users/quittet/Documents/marmot-quantitative-genetics"
load(paste0(dir, "/1_DATA_ANALYSIS/3_Data/pheno.rds"))
load(paste0(dir, "/1_DATA_ANALYSIS/3_Data/mat_info_id.rds"))
load(paste0(dir, "/1_DATA_ANALYSIS/3_Data/pedigree_data.rds"))

# LOAD FUNCTIONS
dir <- "C:/Users/quittet/Documents/marmot-quantitative-genetics/PEDIGREE_GENETIQUE/CONSTRUCTION_PROFIL_GENETIQUE"
source(paste0(dir, "/CERVUS_functions.R"))

3 DATA IMPORTATION

On importe les données des différentes tables de la base de données marmotte (http://umr5558-marmottes.univ-lyon1.fr/phpPgAdmin-master/) et les profils lus sur Genemapper.

channel <- odbcConnect(dsn = "PostgreSQL30")  # DSN name given during configuration for BDD access

3.1 Table “parentee_genetique”

parentee_genetique <- sqlQuery(channel = channel,
                               query =  paste0("
SELECT
    parentee_genetique.individu_id AS individual_id,
    parentee_genetique.mere_id AS mother,
    parentee_genetique.pere_id AS father,
    parentee_genetique.pere_lod AS father_lod,
    parentee_genetique.pere_delta AS father_delta,
    date_part('year', individu.date_emergence) AS annee_emergence,
    individu.sexe,
    capture.territoire_id
FROM
    parentee_genetique
    LEFT JOIN individu ON parentee_genetique.individu_id = individu.id
    LEFT JOIN capture ON parentee_genetique.individu_id = capture.individu_id 
    LEFT JOIN territoire ON territoire.id = capture.territoire_id
    LEFT JOIN zone_etude ON territoire.zone_etude_id = zone_etude.id

WHERE
    zone_etude.nom::text = 'sassiere'::text
                  "))

# Formatage
parentee_genetique <- parentee_genetique[!duplicated(parentee_genetique$individual_id), ]
parentee_genetique <- parentee_genetique[colnames(parentee_genetique) != "territoire_id"]
parentee_genetique <- parentee_genetique |>  
mutate(father = if_else(father == "unknow", NA_character_, father))  # replace wrong levels of "father

3.2 Table “analyse_microsat”

microsatellite_data <- sqlQuery(channel = channel,
                             query =  paste0("
SELECT
    analyse_microsat.id AS id_file,
    analyse_microsat.individu_id AS individual_id,
    analyse_microsat.file_genet AS fsa,
    analyse_microsat.marqueur AS marker,
    analyse_microsat.allele_1,
    analyse_microsat.allele_2,
    date_part('year', individu.date_emergence) AS annee_emergence,
    individu.sexe,
    capture.territoire_id
FROM
    analyse_microsat
    LEFT JOIN individu ON analyse_microsat.individu_id = individu.id
    LEFT JOIN capture ON analyse_microsat.individu_id = capture.individu_id 
    LEFT JOIN territoire ON territoire.id = capture.territoire_id
    LEFT JOIN zone_etude ON territoire.zone_etude_id = zone_etude.id

WHERE
    zone_etude.nom::text = 'sassiere'::text
                  "))     

# FORMATAGE
microsatellite_data <- microsatellite_data[!duplicated(microsatellite_data$id_file), ]  # remove useless doublons
microsatellite_data <- microsatellite_data[colnames(microsatellite_data) != "territoire_id"]

# standardized name of marker (e.g. ms46 = MS46)
microsatellite_data$marker <- toupper(microsatellite_data$marker)

Vérification : On vérifie si on a bien 17 marqueurs

length(unique(microsatellite_data$marker))

## [1] 17

3.3 Profils génétiques lus sur Genemapper

# AJOUT DES PROFILS LUS SUR GENEMAPPER (2016 to 2022)
# read xlsx
pg_2017 <- read.table(file = paste0(dir, "/LECTURE_PROFIL/txt/LECTURE_PROFIL_2017.txt"), header = T, sep = ",", strip.white = T)[1:11]  # pg = profil genetique
pg_2020 <- read.table(file = paste0(dir, "/LECTURE_PROFIL/txt/LECTURE_PROFIL_2020.txt"), header = T, sep = ",", strip.white = T)[1:11]  
pg_2022 <- read.table(file = paste0(dir, "/LECTURE_PROFIL/txt/LECTURE_PROFIL_2022.txt"), header = T, sep = ",", strip.white = T)[1:11] 
pg_2023_PK1_PK2 <- read.table(file = paste0(dir, "/LECTURE_PROFIL/txt/LECTURE_PROFIL_2023_PK1_PK2.txt"), header = T, sep = ",", strip.white = T)[1:11]  
pg_2023_PK3_PK4 <- read.table(file = paste0(dir, "/LECTURE_PROFIL/txt/LECTURE_PROFIL_2023_PK3_PK4.txt"), header = T, sep = ",", strip.white = T)[1:11]  
pg_2023_PK5 <- read.table(file = paste0(dir, "/LECTURE_PROFIL/txt/LECTURE_PROFIL_2023_PK5.txt"), header = T, sep = ",", strip.white = T)[1:11] 

# add year of analysis for future stats
pg_2017$year_analysis <- "2017"
pg_2020$year_analysis <- "2020"
pg_2022$year_analysis <- "2022"
pg_2023_PK1_PK2$year_analysis <- "2023_1"
pg_2023_PK3_PK4$year_analysis <- "2023_2"
pg_2023_PK5$year_analysis <- "2023_3"

# merge datasets
pg <- rbind(pg_2017, pg_2020, pg_2022, pg_2023_PK1_PK2, pg_2023_PK3_PK4, pg_2023_PK5)
rm(list = c("pg_2017", "pg_2020", "pg_2022", "pg_2023_PK1_PK2", "pg_2023_PK3_PK4", "pg_2023_PK5"))
str(pg)

## 'data.frame':    8772 obs. of  12 variables:
##  $ Sample.File  : chr  "1_1360_E07_028.fsa" "1_1360_E07_028.fsa" "1_1360_E07_028.fsa" "1_1360_E07_028.fsa" ...
##  $ Sample.Name  : chr  "1_1360" "1_1360" "1_1360" "1_1360" ...
##  $ Marker       : chr  "Ma091" "Bibl1" "MS6" "Bibl20" ...
##  $ Allele.1     : chr  "173" "97" "158" "208" ...
##  $ Allele.2     : chr  "175" "107" "" "216" ...
##  $ Allele.3     : logi  NA NA NA NA NA NA ...
##  $ Size.1       : num  173.6 98.9 158.3 206.8 107.4 ...
##  $ Size.2       : num  176 108 NA 215 109 ...
##  $ Size.3       : logi  NA NA NA NA NA NA ...
##  $ Height.1     : int  10140 8785 24372 7420 14562 7050 19148 NA 6816 7807 ...
##  $ Height.2     : int  6435 5704 NA 4567 11281 4379 12818 NA 4576 5098 ...
##  $ year_analysis: chr  "2017" "2017" "2017" "2017" ...

# reformat SampleName columns
pg$Sample.Name <- gsub(pattern = "1_|2_|_2", x = pg$Sample.Name, replacement = "")

# remove white space and replace by NA
pg[] <- lapply(pg, function(x) {
  if (is.character(x)) gsub(" ", "", x) else x
})

pg[pg == ""] <- NA

# rename colnames
colnames(pg) <- c(
  "fsa",
  "individual_id",
  "marker",
  "allele_1",
  "allele_2",
  "allele_3",
  "size_1",
  "size_2",
  "size_3",
  "height_1",
  "height_2",
  "year_analysis"
)

# standardized name of marker (e.g. ms46 = MS46)
pg$marker <- toupper(pg$marker)

# for later
pg_raw <- pg

Vérification : On vérifie que l’on a bien 17 marqueurs :

length(unique(pg$marker))

## [1] 17

Vérification : On calcule le nombre d’individu génotypé entre 2016 et 2022 :

length(unique(pg$individual_id))

## [1] 483

D. Allainé avait importé dans la BDD les lectures automatiques GENOSCREEN. Elles sont pleines de fautes car elles n’ont pas été lues directement et beaucoup de marqueurs ne sont pas correctement calibrés. Je supprime ces données.

microsatellite_data <- microsatellite_data[!microsatellite_data$individual_id %in% pg$individual_id,]

3.4 Table “individu”

Les informations de cette table permettent d’avoir le sexe, l’année d’émergence et la plus ancienne capture de tous les individus.

individu <- sqlQuery(
  channel = channel,
  query = paste0(
    "
WITH captures_uniques AS (
    SELECT
        c.*,
        ROW_NUMBER() OVER (
            PARTITION BY c.individu_id
            ORDER BY c.date_capture ASC
        ) AS rn
    FROM capture c
)
SELECT
    i.id AS individual_id,
    date_part('year', i.date_emergence) AS annee_emergence,
    i.sexe,
    cu.territoire_id,
    date_part('year', cu.date_capture) AS annee_capture
FROM individu i
LEFT JOIN captures_uniques cu
       ON i.id = cu.individu_id AND cu.rn = 1
LEFT JOIN territoire t ON t.id = cu.territoire_id
LEFT JOIN zone_etude z ON z.id = t.zone_etude_id
WHERE z.nom::text = 'sassiere'::text
    "
  ))

individu <- individu[!duplicated(individu$individual_id),]  # remove duplicated row
individu <- individu[colnames(individu) != "territoire_id"]

# get sexe for GM data
index <- match(x = pg$individual_id, individu$individual_id)
pg$sexe <- individu$sexe[index]

# get cohort for GM data
index <- match(x = pg$individual_id, pheno$individual_id)
pg$annee_emergence <- pheno$annee_emergence[index]

# get first capture year for social pedigree
index <- match(x = pedigree_data$individual_id, individu$individual_id)
pedigree_data$first_capture <- individu$annee_capture[index]

Vérification : pas de NA dans les sexes des individus

Notre jeu de données contient 1 individu(s) avec sexe == NA…

head(pg[index_temp,])  # due à individual_id : Tnegatif... je supprime ces lignes

##                       fsa individual_id marker allele_1 allele_2 allele_3
## 170  1_1718_2_I23_088.fsa          1718  BIBL1       97     <NA>       NA
## 653    1_1776_O17_065.fsa          1776   MS53      140      142       NA
## 1136   2_1737_N07_020.fsa          1737  MA002      279     <NA>       NA
## 1619   2_1788_J21_087.fsa          1788   ST10      118      120       NA
## 2102   1_F06_1891_038.fsa          1891  BIBL4     190a     <NA>       NA
## 2585   2_C07_1895_059.fsa          1895   ST10      118      132       NA
##      size_1 size_2 size_3 height_1 height_2 year_analysis sexe annee_emergence
## 170   98.74     NA     NA    12836       NA          2017    f            2016
## 653  141.22 143.44     NA     7950     5449          2017    m            2016
## 1136 279.79     NA     NA     2292       NA          2017    m            2016
## 1619 118.98 120.88     NA     5752     4529          2017    f            2016
## 2102 191.37     NA     NA    20059       NA          2020    m            2018
## 2585 118.55 132.19     NA    33980     7878          2020    m            2018

pg <- pg[pg$individual_id != "Tnegatif",]

Ces observations sont issus de l’individu Tnegatif qui doit être un teste de méthode, on supprime cette observation.

Vérification : le nombre d’individus avec des années d’émergence NA

Parmi les individus nouvellement génotypés, 29 n’ont pas de date d’émergence, probablement des individus capturés adultes. Les ID de ces individus sont :

##  [1] 1417 1606 1625 1710 1730 1790 1791 1843 1844 1845 1847 1850 1882 1913 1943
## [16] 1958 1964 1965 1969 1989 2078 2079 2080 2166 2167 2168 2169 2171 2205

4 DATA CLEANING AND FORMATTING

Nombre d’individu génotypé dans la BDD : 1442

4.1 Formattage

On va homogénéiser le codage des homozygotes, identifier et supprimer les allèles et marqueurs non présents dans la référence de Aurélie Cohas (2008), en cas de re-test de certains individus, déterminer quel profil sera conservé, et enfin calculer le nombre d’individus n’ayant pas été génotypés.

4.1.1 RECODE DES HOMOZYGOTES

Cette étape est nécessaire car certains haplotypes sont notés 155/0 ou 155/155 ou 155/NA. On homogénise le tout.

# Fonction d'homogénéisation des homozygotes
recoding_homozygotis <- function(x){
  if(!is.na(x[1]) && x[1] != 0 && (is.na(x[2]) || (!is.na(x[2]) && x[2] == 0))){
    x <- c(x[1], x[1])
  } else if (!is.na(x[2]) && x[2] != 0 && (is.na(x[1]) || (!is.na(x[1]) && x[1] == 0))){
    x <- c(x[2], x[2])
  } else if (is.na(x[1]) && (!is.na(x[2]) && x[2] == 0)){
    x <- c(NA, NA)
  } else {
    x <- c(x[1], x[2])
  }
  return(x)
}

# application de la fonction sur les données de la BDD
microsatellite_data[c("allele_1", "allele_2")] <- t(apply(microsatellite_data[c("allele_1", "allele_2")],
                                                          1,
                                                          recoding_homozygotis, simplify = T))

# ... sur les données lues sur Genemapper
pg[c("allele_1", "allele_2")] <- t(apply(pg[c("allele_1", "allele_2")],
                                                          1,
                                                          recoding_homozygotis, simplify = T))

4.1.2 FORMATAGE DU NOM

On homogénise le noms des allèles en enlevant les lettres e.g. 180a 180b deviennent 180 et les “?” des allèles incertains.

pg$allele_1 <- gsub(pattern = "a|b|c|\\?", x = pg$allele_1, replacement = "")
pg$allele_2 <- gsub(pattern = "a|b|c|\\?", x = pg$allele_2, replacement = "")

4.1.3 ALLELES DE REFERENCE

On ne conserve que les allèles présents dans l’appendix de Cohas et al. (2008). On précise NA pour les autres.

# ALLELES DE REFERENCES
cohas_alleles <- list(
  "BIBL1" = c("95", "97", "101", "103", "107", "109", "0"),
  "BIBL18" = c("132", "137", "143", "145", "147", "149", "0"),
  "BIBL20" = c("206", "208", "216", "218", "220", "222", "0"),
  "BIBL31" = c("157", "159", "161", "163", "0"),
  "BIBL4" = c("175", "178", "188", "190", "192", "0"),
  "MA002" = c("265", "271", "279", "281", "283", "0"),
  "MA018" = c("296", "298", "0"),
  "MA066" = c("231", "233", "241", "0"),
  "MA091" = c("159", "167", "169", "171", "173", "175", "177", "179", "188", "0"),
  "MS41" = c("184", "186", "0"),
  "MS45" = c("107", "109", "111", "0"),
  "MS47" = c("176", "180", "182", "184", "186", "188", "190", "0"),
  "MS53" = c("132", "140", "142", "144", "148", "0"),
  "MS56" = c("104", "106", "108", "110", "0"),
  "MS6" = c("142", "158", "160", "0"),
  "ST10" = c("116", "118", "120", "130", "132", "134", "136", "0")
)

On fait le bilan des allèles problématiques par marqueur :

## 
##  ____________________
## MARKER : BIBL1
##  ____________________  
##  * Number problematic loci : 0 
##  * Problematic alleles :   
## 
##  ____________________
## MARKER : BIBL18
##  ____________________  
##  * Number problematic loci : 1 
##  * Problematic alleles :  133 143 
## 
##  ____________________
## MARKER : BIBL20
##  ____________________  
##  * Number problematic loci : 0 
##  * Problematic alleles :   
## 
##  ____________________
## MARKER : BIBL31
##  ____________________  
##  * Number problematic loci : 0 
##  * Problematic alleles :   
## 
##  ____________________
## MARKER : BIBL4
##  ____________________  
##  * Number problematic loci : 3 
##  * Problematic alleles :  179 186 192 
## 
##  ____________________
## MARKER : MA002
##  ____________________  
##  * Number problematic loci : 7 
##  * Problematic alleles :  277 267 275 281 279 271 
## 
##  ____________________
## MARKER : MA018
##  ____________________  
##  * Number problematic loci : 0 
##  * Problematic alleles :   
## 
##  ____________________
## MARKER : MA066
##  ____________________  
##  * Number problematic loci : 3 
##  * Problematic alleles :  233 231 235 239 
## 
##  ____________________
## MARKER : MA091
##  ____________________  
##  * Number problematic loci : 0 
##  * Problematic alleles :   
## 
##  ____________________
## MARKER : MS41
##  ____________________  
##  * Number problematic loci : 0 
##  * Problematic alleles :   
## 
##  ____________________
## MARKER : MS45
##  ____________________  
##  * Number problematic loci : 0 
##  * Problematic alleles :   
## 
##  ____________________
## MARKER : MS47
##  ____________________  
##  * Number problematic loci : 5 
##  * Problematic alleles :  178 174 184 180 176 
## 
##  ____________________
## MARKER : MS53
##  ____________________  
##  * Number problematic loci : 3 
##  * Problematic alleles :  138 142 140 
## 
##  ____________________
## MARKER : MS56
##  ____________________  
##  * Number problematic loci : 0 
##  * Problematic alleles :   
## 
##  ____________________
## MARKER : MS6
##  ____________________  
##  * Number problematic loci : 8 
##  * Problematic alleles :  158 156 168 
## 
##  ____________________
## MARKER : ST10
##  ____________________  
##  * Number problematic loci : 2 
##  * Problematic alleles :  120 128 122 132

Vérification : Nombre d’observations (i.e. donnée pour un marqueur donné) supprimées

Nombre d’observation(s) supprimée(s) : 32

On regarde les informations relatives à ces observations supprimées :

##    id_file individual_id                fsa marker annee_emergence sexe
## 1     3920          1009               <NA> BIBL18              NA    f
## 2    32943          1327   1327_A15_064.fsa   ST10            2011    f
## 3    33129          1339   1339_E17_075.fsa  MA002              NA    m
## 4    34071          1393 1-1393_O06_017.fsa  MA002              NA    m
## 5    35545          1481  S1481_M19_068.fsa   MS53            2013    m
## 6    34228          1403  S1403_C01_013.fsa  BIBL4              NA    m
## 7    34579          1425  S1425_M05_019.fsa   ST10              NA    m
## 8    34775          1436  S1436_C09_045.fsa  MA002              NA    f
## 9    34639          1428  S1428_A24_096.fsa  MA002              NA    m
## 10   35773          1496  S1496_K23_086.fsa  MA002            2013    f
## 11   34292          1408  S1408_K01_005.fsa    MS6              NA    f
## 12   37477          1596 1-1596_M07_020.fsa  MA002            2014    f
## 13   37116          1575   1575_J19_072.fsa   MS53            2014    f
## 14   36111          1516   1516_D05_029.fsa   MS47            2014    m
## 15   36519          1540   1540_D11_046.fsa   MS47            2014    f
## 16   36366          1531   1531_B09_047.fsa   MS47            2014    m
## 17   36483          1538   1538_P09_033.fsa    MS6            2014    m
## 18   41369          1623 1_1623_K01_005.fsa   MS47              NA    m
## 19   41002          1624 2_1624_N01_003.fsa  MA066              NA    m
## 20   41462          1624 1_1624_M01_003.fsa   MS53              NA    m
## 21   41658          1637 1_1637_G05_025.fsa    MS6            2015    m
## 22   41380          1635 1_1635_A05_031.fsa   MS47            2015    f
## 23   41024          1647 2_1647_L07_022.fsa  MA066            2015    m
## 24   41672          1651 1_1651_C09_045.fsa    MS6            2015    m
## 25   40661          1652 1_1652_E09_043.fsa  BIBL4            2015    m
## 26   41684          1663 1_1663_K11_038.fsa    MS6            2015    m
## 27   40675          1666 1_1666_A13_063.fsa  BIBL4            2015    f
## 28   40872          1679 2_1679_N15_052.fsa  MA002            2015    m
## 29   41704          1683 1_1683_E17_075.fsa    MS6            2015    f
## 30   41066          1689 2_1689_B19_080.fsa  MA066            2015    m
## 31   41711          1690 1_1690_C19_078.fsa    MS6            2015    f
## 32   41712          1691 1_1691_E19_076.fsa    MS6            2015    f

On dirait qu’une grosse partie est de 2014 ou 2015… Cela est sûrement lié à l’analyse de 2016, qui parait bizarre. J’ai la sensation que cette analyse est une lecture automatique de Genoscreen aussi…

4.1.4 SUPPRESSION MARQUEUR MA001

Le marqueur n’est pas présent dans les analyses de parenté des papiers de Aurélie. On le supprime de la liste.

# REMOVE MA0001 NON USED IN PARENTAGE ANALYSIS
microsatellite_data <- microsatellite_data[!microsatellite_data$marker == "MA001",]
pg <- pg[!pg$marker == "MA001",]

4.1.5 VERIFICATIONS FINALES

Vérification : Nom unique des allèles pour vérifier qu’il n’y a pas de bizareries

##  [1]  95 145 216 157 190 279 296 231 169 186 109 182 132 158 106 116 101 149 159
## [20] 281 175 184 140  97 143 208 188 298 173 107 108 118 271 147 241 167 180 142
## [39] 104 218 120 176 130 161 134  NA 137 171 265 178 192 111 222 206 233 163 177
## [58] 136 160 144 103 179 220 283 110 148

##  [1] "173" "97"  "158" "208" "106" "184" "140" NA    "95"  "216" "180" "142"
## [13] "188" "101" "108" "176" "206" "104" "190" "167" "218" "175" "171" "182"
## [25] "159" "186" "132" "178" "169" "107" "192" "220" "179" "231" "116" "298"
## [37] "157" "279" "145" "296" "147" "109" "130" "271" "143" "118" "281" "120"
## [49] "111" "265" "161" "163" "241" "233" "134" "149" "103" "144" "160" "148"
## [61] "136" "283" "110" "177"

Vérification : Tous les allèles présents dans la BDD sont bien présents dans les allèles de références d’Aurélie Cohas (2008)

Nombre d’allèle(s) non-présent(s) : 0

Vérification : Tous les allèles des données lues sur Genemapper sont bien présents dans les allèles de références d’Aurélie Cohas (2008)

Nombre d’allèle(s) non-présent(s) : 0

Vérification : Tous les allèles présents de la référence de A. Cohas sont dans mes allèles extraits de mes analyses GM

## BIBL182 BIBL206 
##   "137"   "222"

On a les allèles 137 et 222 qui ne sont pas présents, mais qui correspondent à des allèles qui ne sont plus utilisés car plus présent dans bin Genemapper (i.e. allèles de références pour GM).

Vérification : Après toutes ces modifications, on recompte le nombre de marqueur

length(unique(microsatellite_data$marker))  # BDD

## [1] 16

length(unique(pg$marker))  # GM

## [1] 16

4.2 Genotype re-test

Certains individus ont été génotypé plusieurs fois, et les résultats des deux analyses sont parfois différents. Pour choisir quelle analyse prendre, je vais d’abord voir, NA non-inclus, si les résultats par marqueurs test vs retest sont congruents. Si oui, je conserve tous allèles non-NAs.

4.2.1 MERGING DATAFRAME BDD ET GENEMAPPER

On merge les deux dataframes car certains individus ont pu être re-testé dans le batch que j’ai analysé.

# artificially add unique id_file to pg
pg$id_file <- paste0("000", 1:nrow(pg))

# ne garder que les colonnes en commun mais en conservant year_analysis
microsatellite_data$year_analysis <- NA
pg <- pg[colnames(pg) %in% colnames(microsatellite_data)]

# reorder columns
pg <- pg[colnames(microsatellite_data)]

# merging
pg_all <- rbind(pg, microsatellite_data)

4.2.2 EXTRACT RE-TESTED INDIVIDUALS

index_temp_retest <- which(duplicated(paste0(pg_all$individual_id, pg_all$marker)))  # row with duplicated marker x individual id
id_temp_retest <- unique(pg_all[index_temp_retest, "individual_id"])

Vérification : Nombre d’individu retesté

Nombre d’individu re-testé : 71.

Leur ID est :

##  [1]  137  832  849  850  864  888  904 1005 1021 1061 1278 1279 1296 1297 1300
## [16] 1303 1311 1312 1316 1319 1327 1336 1345 1352 1355 1359 1361 1375 1376 1381
## [31] 1382 1390 1399 1405 1455 1470 1481 1520 1527 1551 1575 1598 1714 1715 1716
## [46] 1717 1718 1733 1763 1770 1778 1830 1845 1853 1859 1894 1897 1899 1945 1946
## [61] 1947 1952 1957 1960 1961 1990 1994 1999 2046 2124 2138

Vérification : Répartition des re-test par analyse Genemapper

## 
##   2017   2020   2022 2023_1 2023_2 2023_3 
##     10      4     12     11      6     10

Vérification : Outlier(s) entre taille de l’allèle et allèle (e.g. quelle taille en Pb de l’allèle 180 ?)

Des outliers semblent exister en 2017 sur le marqueur BIBL4. On vérifie ces points :

##                  fsa individual_id marker allele_1 allele_2
## 1 1_1732_A07_032.fsa          1732  BIBL4      178      188
## 2 1_1780_G19_074.fsa          1780  BIBL4      178      190

Ils correspondent à l’allèle 178, qui est très rare…! Tout est bon.

4.2.3 MISMATCH TEST RE-RETEST

4.2.3.1 Etape 1**

On extrait les individus pour lesquels on a un mismatch test-retest (i.e. où le test-retest ne donne pas le même résultat)

mismatch_retest_id <- c()  # vector to stock individuals with different test/re-test results
list_error <- list()
mismatch_marker <- c()  # store mismatched marker
mismatch_analysis <- c()  # store analysis that Id did that failed

for(i in id_temp_retest){
  data_temp <- pg_all[pg_all$individual_id == i, ]
  data_temp$haplotype <- paste0(data_temp$marker, "_", data_temp$allele_1, "_", data_temp$allele_2)
  
  # check if at least one mismatch test/re-test is present
  test <- 
    data_temp |> 
    filter(!is.na(allele_1) | !is.na(allele_2)) |> 
    group_by(marker) |> 
    mutate(error = length(unique(haplotype)) != 1) |> 
    arrange(marker) |> 
    as.data.frame() |> 
    filter(error)
  
  # count longest sequence
  nb_na_seq <- 
    data_temp |> 
    group_by(fsa) |> 
    reframe(seq_non_na = n() - sum(is.na(allele_1)))  # longest profil sequence without NAs (i.e. most reliable)

  test <- 
    test |> 
    right_join(y = nb_na_seq, by = "fsa") |> 
    filter(error)
  
  # count number of duplicated test
  haplotype_count <- 
    test |> 
    filter(!is.na(allele_1)) |> 
    count(haplotype)
 
  test <- 
    test |> 
    right_join(y = haplotype_count, by = "haplotype") |> 
    filter(error)
  
  colnames(test)[ncol(test)] <- "n_test"
  
  
  # output
  if(any(test$error)){
    
    # store id
    pb_id <- unique(test$individual_id)
    pb_marker <- unlist(unique(test$marker))
    pb_analysis <-  unlist(unique(test$year_analysis))
    
    mismatch_retest_id <- c(mismatch_retest_id, pb_id)
    mismatch_marker <- c(mismatch_marker, pb_marker)
    mismatch_analysis <- c(mismatch_analysis, pb_analysis)
    
    # store data_temp to manually check results
    list_error[[pb_id]] <- test
    
  }
  }


# export df with all error
df_error <- do.call(rbind, list_error)
# write.xlsx(df_error, file = "test_retest_error_to_check.xlsx")

Vérification : Nombre de mismatch test-retest

Nombre d’individu avec au moins un mismatch test-retest : 16 (23% des individus re-testés)

Vérification : Répartition des mismatchs test-retest par marqueur

## mismatch_marker
##  BIBL1 BIBL18 BIBL20 BIBL31  BIBL4  MA002  MA018  MA066  MA091   MS41   MS45 
##      3      3      2      3      3      2      3      2      2      1      2 
##   MS47   MS53   MS56    MS6   ST10 
##      3      4      1      2     10

Les mismatchs test-retests se trouvent surtout sur le marqueur ST10.

Vérification : Répartition des mismatchs test-retest par analyse Genemapper

## mismatch_analysis
##   2017   2020   2022 2023_1 2023_2 2023_3 
##      1      2      7      6      2      8

4.2.3.2 Etape 2**

On sélectionne la meilleure analyse par marqueur selon les critères suivants (par ordre d’importance)

Qualité de la séquence - Allèles correspondant au profil avec la plus longue séquence sans NAs
Haplotype le plus fréquent - Si séquence sans NA égale, prendre l’haplotype le plus fréquent
Analyse la plus récente - Si toutes les fréquences sont égales, prendre la séquence la plus récente

df_error$year_analysis[df_error$year_analysis %in% c("2023_1", "2023_2", "2023_3")] <- 2023  # pour tirer au plus récent
to_remove <- c()

for(ind_temp in unique(df_error$individual_id)){
  marker_pb <- unique(df_error[df_error$individual_id == ind_temp, "marker"])
  
  for(marker_temp in marker_pb){
    df_filtered <- df_error[df_error$individual_id == ind_temp & df_error$marker == marker_temp, ]
    
      # a. Garder séquence(s) la/les plus longue(s)
      df_temp <- df_filtered |> 
        filter(seq_non_na == max(seq_non_na, na.rm = TRUE))
      
      # b. Si plusieurs séquences équivalentes, garder haplotype le plus fréquent
      if(nrow(df_temp) > 1){
        freq_tab <- table(df_temp$haplotype)
        max_freq <- max(freq_tab)
        best_haplos <- names(freq_tab[freq_tab == max_freq])
        
        df_temp <- df_temp |> 
          filter(haplotype %in% best_haplos)
      }
      
      # c. Si encore ex-aequo, garder la plus récente analyse
      if(nrow(df_temp) > 1 & !all(is.na(df_temp$year_analysis))){
        df_temp <- df_temp |> 
          filter(year_analysis == max(year_analysis, na.rm = TRUE))
      }
      
      # Sélection finale : garder le premier (au cas d’égalité stricte)
      df_keep <- df_temp[1, ]
      
      # Marquer tous les autres comme "à supprimer"
      ids_remove <- setdiff(df_filtered$id_file, df_keep$id_file)
      to_remove <- c(to_remove, ids_remove)
  }
}

# On supprime les lignes avec les données problématiques ou redondantes
pg_all <- pg_all[!pg_all$id_file %in% to_remove, ]

Le problème de cette méthode est que on va crée des chimères de profils génétiques issus de plusieurs génotypage (i.e. fichiers .fsa).

4.2.4 MATCH TEST-RETEST

Pour les individus sans mismatch test/retest, on supprime simplement les doublons. En cas de NA pour un marqueur donnée, on conserve la donnée non-NA.

for(i in id_temp_retest){  # [!id_temp_retest %in% mismatch_retest_id]
  
  data_temp <- pg_all[pg_all$individual_id %in% i, ]
  id_file_initial <- data_temp$id_file
  
  # on optimise les données gardées i.e. si NAs dans test et pas NAs dans retest, on garde retest
  id_file_to_keep <- 
    data_temp |> 
    group_by(marker) %>%
    # Si le marker a au moins une ligne complète, supprimer les NAs
    filter(
      !(is.na(allele_1) & is.na(allele_2) & any(!is.na(allele_1) | !is.na(allele_2)))
    ) %>%
    filter(!duplicated(marker)) |>   # on enlève les doublons car on sait qu'ils sont identiques (test précédent)
    ungroup() |> 
    dplyr::select(id_file) |> 
    as.data.frame() |> 
    unlist()
  
  # remove the non wanted id file
  id_file_to_remove <- id_file_initial[!id_file_initial %in% id_file_to_keep]
  
  # filter dataframe
  pg_all <- pg_all[!pg_all$id_file %in% id_file_to_remove, ]
}

Même problème que plus haut, cette méthode est que on va crée des chimères de profils génétiques issus de plusieurs génotypage (i.e. fichier .fsa).

Vérification : Après cette méthode, on s’assure qu’on a plus d’individus avec plusieurs données pour un marqueur donné

Nombre d’individu(s) dupliqué(s) : 0

4.3 Non-genotyped individuals

On calcule le nombre d’individu n’ayant pas été génotypés.

# we filtered by updated years 
individus_filtered <- individu[individu$annee_capture %in% c(1990:2022), ]  # annee_capture == annee de la première capture
pg_filtered <- pg_all[!duplicated(pg_all$individual_id),]

# nb d'individus non-génotypées
non_genotyped <- (!individus_filtered$individual_id %in% pg_filtered$individual_id)

Nombre d’individu(s) non-génotypé(s) : 118 / 2042 (5.78%)

Vérification : Répartition des non-génotypés dans les années

## 
## 1990 1991 1992 1996 1997 2000 2003 2011 2012 2013 2014 2016 2017 2018 2019 
##   16   11    6    8    2    1   13    1    7    6    1    1   41    1    3

La plupart des individus non-génotypés sont capturés pour la première fois en 2017, suivi de 1990 et 2003. D’après Corinne, les individus de 2017 non-génotypés correspondent à des oublis.

Vérification : Informations sur ces individus de 2017

Les ID de ces individus capturés en 2017 sont :

##  [1] 1793 1794 1795 1797 1798 1799 1800 1801 1803 1805 1806 1807 1808 1809 1810
## [16] 1811 1812 1813 1815 1817 1818 1819 1820 1821 1822 1823 1824 1825 1826 1829
## [31] 1831 1832 1833 1834 1836 1837 1838 1839 1840 1841 1842

La plupart (37) sont des marmottons (i.e. date d’émergence en 2017 et date de première capture en 2017) et 4 sont des individus capturés adultes.

Vérification : Nombre d’individus non-génotypés qui sont parents dans le pédigrée social

Ces individus correspondraient aux individus prioritaires pour les analyses car ils apporteraient de l’information de parentée.

Nombre d’individu(s) non-génotypé(s) parent(s) : 7.

Leur ID sont :

## [1]    6   45   55 1605 1818 1834 1842

Vérification : Répartition des parents non-génotypés selon l’année de capture

## 
## 1990 1991 1992 2014 2017 
##    1    1    1    1    3

La plupart sont des parents de 2017… et 3 d’entre eux sont de très vieux individus (1990, 1991 et 1992).

5 FORMATTING DATA FOR CERVUS

On formatte les données comment dans CERVUS pour avoir un tableau avec en ligne les individus, en colonne leur profil génétique au loci donné (i.e. 2 par allèles e.g. BIBL1a / BIBL1b) et des 0 à la place des NA.

5.1 Formatage

# FORMAT DF
pg_cervus <- pg_all |>
  arrange(individual_id, desc(id_file)) |>
  filter(!duplicated(paste0(marker, individual_id))) |>
  
  #formatting for cervus : ID / SEXE / marker_a / marker_b
  pivot_wider(id_cols = c(individual_id, sexe, annee_emergence),
              names_from = marker,
              values_from = c(allele_1, allele_2),
              names_sep = "",
              names_vary = "slowest"
  )

pg_cervus <- as.data.frame(pg_cervus)

# RENAME COL FOR CERVUS
former_colnames <- colnames(pg_cervus)[4:ncol(pg_cervus)] 
suffixe <- rep(x = c("a", "b"), length(former_colnames)/2)
new_names <- paste0(gsub(pattern = "allele_1|allele_2", replacement = "", x = former_colnames), suffixe)
colnames(pg_cervus)[4:ncol(pg_cervus)] <- new_names


# REPLACE NA PER 0 (cf Cervus tutorial)
pg_cervus[4:ncol(pg_cervus)] <- apply(pg_cervus[4:ncol(pg_cervus)], 2, function(x) replace_na(as.numeric(as.character(x)), 0))

# REMOVE BLANK SPACE
pg_cervus[] <- lapply(pg_cervus, function(x) {
  if (is.character(x)) gsub(" ", "", x) else x
})  # pg[] permet de garder la structure en data.frame

pg_cervus$individual_id <- as.character(pg_cervus$individual_id)  # IMPORTANT

head(pg_cervus)

##   individual_id sexe annee_emergence MA066a MA066b MA018a MA018b MA002a MA002b
## 1           100    m            <NA>    231    241    296    298    271    279
## 2          1000    f            2006    231    231    298    298    281    281
## 3          1001    f            2006    231    241    298    298    279    281
## 4          1002    f            2006    231    231    298    298    279    281
## 5          1003    m            2006    231    231    298    298    281    281
## 6          1004    f            2006    231    231    298    298    281    281
##   BIBL4a BIBL4b BIBL31a BIBL31b BIBL20a BIBL20b BIBL18a BIBL18b BIBL1a BIBL1b
## 1    175    188     157     159     208     208     137     149     97    101
## 2    190    192     157     157     208     216     143     145     95    101
## 3    188    190     157     159     208     216     143     145    101    109
## 4    188    190     157     159     208     216     145     147    101    109
## 5    188    190     157     159     208     218     143     145    101    107
## 6    188    190     157     159     208     218     145     149    101    107
##   ST10a ST10b MS56a MS56b MS6a MS6b MS53a MS53b MS47a MS47b MS45a MS45b MS41a
## 1   130   134   106   108  158  158   140   140   182   184   109   109   186
## 2   116   118   108   108  158  158   140   142   186   186   107   109   186
## 3   116   132   106   108  158  158   140   140   176   182   109   109   186
## 4   116   118   106   108  158  158   132   142   176   182   109   109   186
## 5   118   130   106   108  158  158   140   140   186   186   107   109   184
## 6   118   130   108   108  158  158   132   140   182   186   109   109   184
##   MS41b MA091a MA091b
## 1   186    167    175
## 2   186    175    175
## 3   186    159    175
## 4   186    159    175
## 5   186    175    175
## 6   184    173    175

Vérification : Présence de NA au lieu de 0 dans les allèles ?

Nombre de NA dans le tableau : 0

Vérification : Respect de l’ordre des allèles ?

Nombre de cas où allèle 1 < allèle 2 : 0

Vérification : Nombre de profils ne contenant que des 0s

Nombre d’individus sans données : 0

5.2 Vérification des fréquences alléliques

On regarde par marqueur et par années d’analyse si les fréquences alléliques à chaque marqueur ne dévie pas trop de ceux attendus spécifiés dans Cohas et al. (2008) et ainsi repérer d’éventuels outliers. Pour cela, on va comparer les déviations à l’attendues des fréquences alléliques pour chaque marqueur pour chque analyse Genemapper.

5.2.1 Importation des références

# IMPORTATION DES REFERENCES
references_freq_wide <- read.csv(paste0(dir, "/COHAS_ALLELES_FREQUENCY.csv"), sep=";", header=TRUE, stringsAsFactors = FALSE)

# Transformation en format long
references_freq <- pivot_longer(
  references_freq_wide,
  cols = everything(),
  names_to = c("Marqueur", ".value"),
  names_sep = "_"
) |> as.data.frame()
colnames(references_freq) <- c("marker", "allele", "freq")
references_freq$marker <- toupper(references_freq$marker)
references_freq <- na.omit(references_freq)

# Vérifier le résultat
head(references_freq, n = 5)

##   marker allele freq
## 1  BIBL1     95 0.15
## 2 BIBL18    132 0.01
## 3 BIBL20    206 0.01
## 4 BIBL31    157 0.50
## 5  BIBL4    175 0.13

5.2.2 Proportion allélique par jeu de données

On calcule les proportions alléliques pour les données de la BDD et de chaque analyse Genemapper.

###########################
#' # Overall proportion
###########################
freq_df <- data.frame(matrix(character(), nrow = nrow(pg_cervus)*2, ncol = 16))
colnames(freq_df) <- unique(gsub(pattern = "a|b", replacement = "", x = colnames(pg_cervus[5:ncol(pg_cervus)])))

j= 1
for(i in seq(4, (ncol(pg_cervus)), by = 2)){
  freq_df[,j] <- c(pg_cervus[,i], pg_cervus[,(i+1)])
  j = j + 1
}

freq_df[freq_df == 0] <- NA  # on remplace les 0 par NAs pour ne pas bruité les proportions
freq_df <- freq_df[names(cohas_alleles)]  # on réordonner dans le même ordre que le tableau de référence

prop_alleles_bdd <- apply(freq_df, 2, function(x) round(prop.table(table(x)), 2))  # fréquence des allèles par marqueur


###########################
#' # 2016 - 2022 proportion
###########################
# Test des proportions par années d'analyse
marker_names <- unique(references_freq$marker)
prop_alleles_gm <- vector("list", 6)
pg_prop <- pivot_longer(data = pg, cols = c("allele_1", "allele_2"), names_to = "position", values_to = "allele")
for(analyse_temp in unique(pg$year_analysis)){
  pg_prop_filtered <- pg_prop[pg_prop$year_analysis == analyse_temp, ]
  pg_prop_filtered[pg_prop_filtered == 0] <- NA  # on remplace les 0 par NAs pour ne pas bruité les proportions

  for(marker_temp in marker_names){
    prop_alleles_gm[[analyse_temp]][[marker_temp]] <- 
      round(prop.table(table(pg_prop_filtered[pg_prop_filtered$marker == marker_temp, "allele"])), 2)
  }
}

5.2.3 Déviations des proportions attendues par analyse

# BDD
marker_names <- unique(references_freq$marker)
analysis_tested <- "bdd"
references_freq[analysis_tested] <- NA

for(marker_temp in marker_names){
  data_temp <- prop_alleles_bdd[[marker_temp]]
  alleles_temp <- names(data_temp)
  freq_temp <- data_temp
  
  for(i in 1:length(alleles_temp)){
    bool_temp <- references_freq$marker == marker_temp & 
                      references_freq$allele == alleles_temp[i] 
    
    if(sum(bool_temp, na.rm = T) > 0){
          references_freq[bool_temp, analysis_tested] <- freq_temp[i] - references_freq[bool_temp, "freq"]
    }
  }
}

# GENEMAPPER ANALYSIS
for(analysis_tested in c("2017", "2020", "2022", "2023_1", "2023_2", "2023_3")){
  
  references_freq[analysis_tested] <- NA

for(marker_temp in marker_names){
  data_temp <- prop_alleles_gm[[analysis_tested]][[marker_temp]]
  alleles_temp <- names(data_temp)
  freq_temp <- data_temp
  
  for(i in 1:length(alleles_temp)){
    bool_temp <- references_freq$marker == marker_temp & 
                      references_freq$allele == alleles_temp[i] 
    
    if(sum(bool_temp, na.rm = T) > 0){
          references_freq[bool_temp, analysis_tested] <- freq_temp[i] - references_freq[bool_temp, "freq"]
    }
  }
}
}

Plus de bruit dans les données que je viens de lire mais aussi probablement car celle de la BDD sont moyennées sur 20 ans…

5.3 Save

# # SAVE 
# dir <- "C:/Users/quittet/Documents/marmot-quantitative-genetics/PEDIGREE_GENETIQUE/CERVUS/"
# write.table(pg_cervus, file = paste0(dir, "ANALYSIS_0_BDD_CERVUS/genotype.txt"), row.names = F, sep = "\t")  # analyse 0
# write.table(pg_cervus, file = paste0(dir, "ANALYSIS_1_BDD_CERVUS/genotype.txt"), row.names = F, sep = "\t")  # analyse 1
# write.table(pg_cervus, file = paste0(dir, "ANALYSIS_2_BDD_CERVUS/genotype.txt"), row.names = F, sep = "\t")  # analyse 2

6 IDENTIFICATION DES INDIVIDUS A TESTER AVEC CERVUS

Les individus candidats sont les individus répondant aux critères suivants :

Génotypé entre 2016 et 2022
Capturés marmotton (i.e. date d’émergence connu)
Sans mère génétique précisé dans la BDD

#' extract id of all new tested individuals
id_tested_2016_2022 <- unique(pg$individual_id)

# ajout année de capture la plus ancienne
index_temp <- match(pg_cervus$individual_id, individu$individual_id)
pg_cervus$oldest_capture <- individu$annee_capture[index_temp]
pg_cervus <- pg_cervus[c(1,2,3,ncol(pg_cervus), 4:(ncol(pg_cervus) - 1))]  # sort columns

# EXTRAIRE LES ID DES INDIVIDUS A TESTER AVEC CERVUS
#' ... (*) qui ont été génotypé entre 2016 et 2022
#' ... (*) mère sociale connue
#' ... (*) sans parents génétiques dans la BDD
id_pups_to_test <- pg_cervus[pg_cervus$individual_id %in% id_tested_2016_2022 &  # génotypé en 2016 - 2022
                               !is.na(pg_cervus$annee_emergence) &  # date d'émergence connue
                               !pg_cervus$individual_id %in% parentee_genetique$individual_id[!is.na(parentee_genetique$mother)],  # sans mère génétiques dans la BDD
                             "individual_id"]

Nombre d’individu(s) à tester avec CERVUS : 437

Vérification : Nombre de retest parmi les individus à analyser

Nombre d’individu ayant au moins une donnée issu d’un retest : 32.

Leur ID est :

##  [1] 1319 1520 1598 1714 1715 1716 1717 1718 1733 1763 1770 1778 1830 1845 1853
## [16] 1859 1894 1897 1899 1945 1946 1947 1952 1957 1960 1961 1990 1994 1999 2046
## [31] 2124 2138

Vérification : Informations sur ces individus re-testés

##      individual_id annee_emergence sexe annee_capture genemapper_analysis
## 357           1319            2011    m          2011              2023_3
## 577           1520            2014    f          2014              2023_1
## 662           1598            2014    m          2014                2017
## 790           1714            2016    m          2016                2017
## 791           1715            2016    m          2016                2017
## 792           1716            2016    f          2016                2017
## 793           1717            2016    m          2016                2017
## 794           1718            2016    f          2016                2017
## 811           1733            2016    f          2016                2017
## 844           1763            2016    f          2016                2017
## 852           1770            2016    f          2016                2017
## 860           1778            2016    f          2016                2017
## 918           1830            2017    f          2017              2023_1
## 934           1845              NA    m          2018                2022
## 943           1853            2018    m          2018                2020
## 949           1859            2018    m          2018                2020
## 988           1894            2018    f          2018                2020
## 991           1897            2018    f          2018                2022
## 993           1899            2018    f          2018                2020
## 1045          1945            2019    m          2019                2022
## 1046          1946            2019    f          2019                2022
## 1047          1947            2019    m          2019                2022
## 1053          1952            2019    m          2019                2022
## 1058          1957            2019    m          2019                2022
## 1062          1960            2019    f          2019                2022
## 1063          1961            2019    m          2019                2022
## 1095          1990            2019    m          2019                2022
## 1099          1994            2019    m          2019                2022
## 1104          1999            2019    m          2019                2022
## 1158          2046            2020    m          2020              2023_1
## 1245          2124            2021    m          2021              2023_1
## 1260          2138            2021    f          2021              2023_1

Les retests ont surtout eu lieu sur des marmottons en 2017 et 2022.

7 RAW PARENTAGE COMPARISON

L’idée ici est de tester si la mère a toujours au moins un de ses allèles en commun avec le marmotton. On utilise pas CERVUS ici.

7.1 Mother test

On teste ici uniquement avec la mère. On regarde si au moins une allèle de l’offspring par locus est présent chez la mère. L’idée étant qu’on connait l’identité de la mère, qui est la femelle dominante. Donc sauf cas rarissime, on ne devrait pas avoir d’erreur.

Vérification : Tous les individus à tester sont présents dans le pédigrée

Nombre d’individu à tester non-présent dans le pédigrée : 1.

Les ID sont :

## [1] "2077"

Seul l’individu 2077 n’est pas présent. Cela est du au fait que c’est un individu qui a perdu sa bague et qui est 1967 en réalité, après comparaison de son génotype (voir fichier famille). On corrige :

pg_cervus <- pg_cervus[!pg_cervus$individual_id == 1967, ]  # on conserve les data de 1967 qui ont plus de marqueur
pg_cervus$individual_id[pg_cervus$individual_id == 2077] <- "1967"  # on renomme
id_pups_to_test <- id_pups_to_test[!id_pups_to_test == "2077"]  # on supprime 2077 des pups id à tester

7.1.1 Mère sociale des individus

On filtre le pédigrée pour extraire les parents sociaux des individus à tester. On trie le pédigrée selon les critères suivants :

Individus génotypés entre 2016 et 2022
Mère sociale connue
Mère sociale a été génotypé

Vérification : Nombre d’individu avec mère sociale NA

Individus à tester avec mère == NA : 25.

Leur ID est :

ped_temp[is.na(ped_temp$mother), "individual_id"]

##  [1] "1626" "2177" "2170" "2165" "2147" "2146" "2140" "2139" "2030" "2029"
## [11] "2028" "2027" "1909" "1890" "1885" "1884" "1883" "1846" "1758" "1757"
## [21] "1756" "1738" "1722" "1721" "1720"

Vérification : Cohorte des individus avec mère sociale NA

table(ped_temp[is.na(ped_temp$mother), "cohort"], useNA = "always")

## 
## 2013 2016 2018 2020 2021 <NA> 
##    1    7    5    4    5    3

Ces bizareries ont été envoyées à Rebecca pour correction.

Vérification : Nombre de mère non-génotypée

Nombre de mère(s) non-génotypée(s) : 2

Leur ID est :

## [1] "1818" "1605"

Nombre d’individu(s) avec mère non-génotypée : 10

7.1.2 Comparaison des génotypes

On trie le pédigrée selon les critères cités plus haut :

ped_tested <- pedigree_data[pedigree_data$individual_id %in% id_pups_to_test &  # id génotypé entre 2016-2022
                              !is.na(pedigree_data$mother) &  # mère sociale connue
                              !pedigree_data$mother %in% id_non_genotyped, ]  # mère génotypée

Vérification : Nombre d’individu testé effectivement

test <- nrow(ped_tested)

Nombre d’individu(s) répondant à tous les critères de filtre : 401

On teste le match entre le profil de la mère et des offsprings à l’aide de la fonction match_genotype() disponible dans le code CERVUS_fonctions.r :

res_match <- apply(X = ped_tested, MARGIN = 1, FUN = function(x) match_genotype(x)$parent_match)

Associations correctes : 311 / 401 (78%)

7.1.3 Exploration des mismatchs

A noter que l’ensemble des parentées sociales ont été re-vérifiées à partir du fichier famille quand il y avait au moins un mismatch…

# GET MISMATCH INDIVIDUAL INFORMATION
index_temp <- as.numeric(which(!res_match))
index_temp <- match(ped_tested$individual_id, individu$individual_id)  # match cohort
ped_tested["annee_emergence"] <- individu[index_temp, "annee_emergence"]
ped_tested[["match_mother"]] <- res_match
ped_tested["n_loci_compared"] <- apply(X = ped_tested, MARGIN = 1, FUN = function(x) match_genotype(x)$n) 
ped_tested["n_error"] <- apply(X = ped_tested, MARGIN = 1, FUN = function(x) match_genotype(x)$nb_mismtach_parent_offspring)

7.1.3.1 Nombre d’erreur et nombre de loci comparé**

On calcule ici le nombre de cas où un individu a un mismatch avec sa mère. On précise également le nombre de loci comparable entre les deux (en excluant les NA).

## Nombre de marqueurs avec mismatch
##   0   1   2   4 
## 311  84   5   1

## Nombre de loci comparé
##   4   6   7   8   9  10  11  12  13  14  15  16 
##   2   4   6   8   4  13  19  36  55  83 103  68

L’immense majorité des mismatchs est situé seulement sur un marqueur. 5 individus ont 2 mismatchs avec leur mère et un individu a 4 mismatchs.

7.1.3.2 Nombre de mère unique avec mismatch**

Nombre de mère unique avec au moins un mismatch avec leur offspring : 20.

Leur ID est :

##  [1] 1292 1298 1315 1408 1432 1443 1475 1520 1560 1564 1569 1586 1626 1635 1657
## [16] 1868 1897 1913 1925 1941

7.1.3.3 Répartition des erreurs dans les analyses**

## 
##   2017   2020   2022 2023_1 2023_2 2023_3 
##      3     27     18     13     16     13

7.1.3.4 Répartition des erreurs parmi les marqueurs**

## 
##  BIBL1 BIBL18 BIBL20 BIBL31  BIBL4  MA002  MA066   MS41   MS47   ST10 
##      2      1      2     19      1      2      6     10      2     53

La grande majorité des erreurs provient du marqueur ST10 et BIBL31. J’ai vérifié les profils génétiques sur Genemapper un à un quand je le pouvais, et parfois ces mismatchs sont simplement incompréhensibles (e.g. les deux profils sont clairs, mais différents) ou bien il s’agit d’un décalage de \(\pm\) 1 avec la mère dans la BDD, probablement lié à une mauvaise lecture dans le passé. Il s’agit très certainement de problème de manipulation. Pour ne pas perdre trop de temps, je précise NA pour tous les cas où il y a mismatch.

# on match avec l'année d'analyse
index_temp <- match(df_error$individual_id, pg$individual_id)
df_error$year_analysis <- pg[index_temp, "year_analysis"]  # get analysis corresponding to mismatch
index_temp <- match(df_error$mother_id, pg$individual_id)
df_error$year_analysis_mother <- pg[index_temp, "year_analysis"]  # get analysis corresponding to mismatch

# on match le nombre d'erreur total avec la mère de l'offspring
index_temp <- match(df_error$individual_id, ped_tested$individual_id)
df_error$n_error_total <- ped_tested$n_error[index_temp]


# COMPREHENSION DES ERREURS...
#' ... Certaines erreurs viennent d'un décalage avec la mère. On regarde l'haplotype
#' ... du père pour voir de qui vient le décalage (de moi ou de la bdd)...
index_temp <- match(df_error$individual_id, pedigree_data$individual_id)
df_error$father <- pedigree_data[index_temp, "father"]

# fonction pour ajouter l'haplotype du père au DF
#' ... l'idée est de détecter de qui vient le décalage à l'allèle donné.
get_father_haplo <- function(i){
  father_temp <- df_error[i, c("father")]
  marker_temp <- df_error[i, c("marker")]
  df_temp <- df_error[i, c("marker", "off_haplo", "mother_haplo")]

  haplo_father <- unlist(pg_cervus[pg_cervus$individual_id == father_temp, grepl(paste0("\\b", marker_temp, "\\b"), gsub("a|b", "", colnames(pg_cervus)))])
  haplo_father <- paste0(haplo_father[1], "/", haplo_father[2])

  df_temp$haplo_father <- haplo_father

  if(length(haplo_father) > 0){
    return(df_temp)
  } else{
    return("Father non-génotypé")
  }
}

df_error$father_haplo <- NA
for(i in 1:nrow(df_error)){
  df_error[i, "father_haplo"] <- get_father_haplo(i)[, "haplo_father"]
}

# ON SUPPRIME LES MARQUEURS DES OFFSPRINGS AVEC UN MISMATCH AVEC LA MERE
#'... même si parfois l'erreur vient clairement du génotype de la mère... 
#'... ça me parait la solution la plus conservatrice.
for(i in 1:nrow(df_error)){
  if(df_error[i, "n_error_total"] == 1){
    marker_temp <- df_error[i, "marker"]
    id_temp <- df_error[i, "individual_id"]
    pg_cervus[pg_cervus$individual_id == id_temp, grepl(paste0("\\b", marker_temp, "\\b"), x =  gsub("a|b", "", colnames(pg_cervus)))] <- 0
  }
}

7.1.3.5 Reste des erreurs…**

On représente les mismatchs avec plus d’une erreur avec leur mère :

##    individual_id marker mother_id off_haplo mother_haplo year_analysis
## 1           1835 BIBL20      1432   218/218      208/208          2022
## 2           1835   MS47      1432   186/186      176/184          2022
## 3           1933  MA066      1475   233/233      231/231          2020
## 4           1933   ST10      1475   116/116      118/132          2020
## 5           1948  MA002      1564   279/279      281/281          2022
## 6           1948   MS41      1564   184/184      186/186          2022
## 7           2075  BIBL4      1569   175/188      190/190        2023_1
## 8           2075 BIBL18      1569   147/147      143/149        2023_1
## 9           2075   ST10      1569   118/132      120/134        2023_1
## 10          2075   MS47      1569   180/180      176/182        2023_1
## 11          1827  BIBL1      1635     97/97      101/107        2023_3
## 12          1827   MS41      1635   186/186      184/184        2023_3
## 13          1894  BIBL1      1635    97/109      101/107          2020
## 14          1894   MS41      1635   186/186      184/184          2020
##    year_analysis_mother n_error_total father father_haplo
## 1                  <NA>             2   1730      208/218
## 2                  <NA>             2   1730      186/186
## 3                2023_3             2   1526      231/231
## 4                2023_3             2   1526      116/116
## 5                  <NA>             2   1440      271/279
## 6                  <NA>             2   1440      184/186
## 7                  <NA>             4   1686      175/188
## 8                  <NA>             4   1686      143/143
## 9                  <NA>             4   1686      120/134
## 10                 <NA>             4   1686      180/180
## 11                 <NA>             2   1309       97/109
## 12                 <NA>             2   1309      186/186
## 13                 <NA>             2   1309       97/109
## 14                 <NA>             2   1309      186/186

Ces mismatchs sont répartis entre 5 mères.

Leur ID est :

## [1] 1432 1475 1564 1569 1635

Après avoir regardé à la main, les mères 1432, 1475, 1564 et 1635 n’ont que très peu de loci comparé avec leur offsprings… (respectivement 8, 9, 9 et 8/7). Cela signifie que soit l’offspring soit la mère a eu un génotypage qui a raté.

En revanche, la mère 1569 a 12 allèles comparés… Il est possible que ce ne soit pas la mère de cette portée.

Analyse des profils génétiques 2016 à 2022

Pierre-Alexandre QUITTET

2025-10-03

1 SESSION INFORMATION

1.1 RANDOM INFORMATION

2 PACKAGE, DATA AND FUNCTIONS

3 DATA IMPORTATION

3.1 Table “parentee_genetique”

3.2 Table “analyse_microsat”

3.3 Profils génétiques lus sur Genemapper

3.4 Table “individu”

4 DATA CLEANING AND FORMATTING

4.1 Formattage

4.1.1 RECODE DES HOMOZYGOTES

4.1.2 FORMATAGE DU NOM

4.1.3 ALLELES DE REFERENCE

4.1.4 SUPPRESSION MARQUEUR MA001

4.1.5 VERIFICATIONS FINALES

4.2 Genotype re-test

4.2.1 MERGING DATAFRAME BDD ET GENEMAPPER

4.2.2 EXTRACT RE-TESTED INDIVIDUALS

4.2.3 MISMATCH TEST RE-RETEST

4.2.3.1 Etape 1**

4.2.3.2 Etape 2**

4.2.4 MATCH TEST-RETEST

4.3 Non-genotyped individuals

5 FORMATTING DATA FOR CERVUS

5.1 Formatage

5.2 Vérification des fréquences alléliques

5.2.1 Importation des références

5.2.2 Proportion allélique par jeu de données

5.2.3 Déviations des proportions attendues par analyse

5.3 Save

6 IDENTIFICATION DES INDIVIDUS A TESTER AVEC CERVUS

7 RAW PARENTAGE COMPARISON

7.1 Mother test

7.1.1 Mère sociale des individus

7.1.2 Comparaison des génotypes

7.1.3 Exploration des mismatchs

7.1.3.1 Nombre d’erreur et nombre de loci comparé**

7.1.3.2 Nombre de mère unique avec mismatch**

7.1.3.3 Répartition des erreurs dans les analyses**

7.1.3.4 Répartition des erreurs parmi les marqueurs**

7.1.3.5 Reste des erreurs…**