---
title: "Étude de la démographie récente chez une espèce endémique en déclin : la centaurée de la Clape"
author : "Océane EYCHENNE"
date : "Stage M2 EE 2024 - ISEM"
output: 
  flexdashboard::flex_dashboard:
    theme: flatly
    orientation : rows 
    social : menu
    source_code : embed
    vertical_layout : scroll
---

```{r setup, include=FALSE}
library(flexdashboard)
library(htmltools)
library(reticulate)
library(ggplot2)
library(plotly)
library(tidyr)
library(scales)
library(plotly)
library(viridis)
```

# ECHANTILLONNAGE {data-icon="fa-seedling"}

## Row {.tabset}

### Fiche biologique de l'espèce {data-height="1600"}

::: {style="text-align:center;"}
![](/Users/oceane/Desktop/Stage/APPENDIX/Fiche_Biologique_2.png){width="1050px" height="650px"}
:::

------------------------------------------------------------------------

### Lieux d'échantillonnage

Cette étude a été réalisée sur 5 des sous-populations de la centaurée de la Clape. Des graines ont été récoltées en populations naturelles en juin 2013 par Eric Imbert selon un maximum de trois graines par individu et uniquement sur des plantes accessibles. Parmi l’ensemble des graines échantillonnées, **96** individus ont été utilisés pour obtenir des données génétiques avec un objectif de **30** individus pour les trois grandes populations (**Au**, **E1** et **E2**) et **3** individus pour les deux petites populations (**Pe** et **Po**).

::: {style="text-align:center;"}
![](/Users/oceane/Desktop/Stage/FIGURES/Figure_Sampling.png){width="750px" height="500px"}
:::

------------------------------------------------------------------------

# BIO-INFORMATIQUE {data-icon="fa-desktop"}

## Row {.tabset}

### Pipeline {data-height="5000"}

::: {style="text-align:center;"}
![](/Users/oceane/Desktop/Stage/FIGURES/Pipeline.png){width="800px" height="850px"}
:::

------------------------------------------------------------------------

### Avant-propos

Afin de réussir à obtenir des génomes phasés, une succession d'étapes de traitement bio-informatique doit être réalisée grâce à la pipeline [Harpy *v.0.4.0*](https://pdimens.github.io/harpy/) qui permet de traiter les données de type haplotagging.

::: {style="text-align:center;"}
![](/Users/oceane/Desktop/Stage/APPENDIX/FIGURES_BIO_INFO/Pipeline_HARPY.png){width="550px" height="150px"}
:::

Avant de lancer la pipeline Harpy, deux étapes préliminaires sont à réalisées :

(1) Conversion des données =\> passer du format BCL (données Illumina) au format *.fastq* : BCL Convert

```{bash, eval=FALSE, echo = TRUE}
# téléchargement du logiciel BCL Convert 

# conversion des fichiers
bcl-convert --bcl-input-directory /mnt/DATA3/DEVOCGEN/RAW/Ccory/Desmarais --output-directory /mnt/DATA3/DEVOCGEN/RAW/Ccory/Desmarais_E_02_bclconvert_out
```

Ainsi, on obtient des fichiers au format *.fastq.qz* : ex : Undetermined_S0_L001_I1_002

(2) Création d'un environnement mamba (conda) =\> meilleure rapidité d'exécution

```{bash, eval=FALSE, echo = TRUE}
# création de l'environnement harpy
mamba create -n harpy -c bioconda -c conda-forge harpy

# activation de l'environnement harpy 
mamba avtivate harpy
```

**NB** : Un clic sur les noms de fichier *.html* donne accès à un lien contenant les *.html* de sortie visibles par téléchargement.

------------------------------------------------------------------------

### Harpy : Démultiplexage

L'étape de demultiplexage permet, dans un premier temps, d'assigner chaque molécule (read) à une combinaison de barcodes (A, B, C et D). Puis, dans un second temps, d'assigner la combinaison de barcodes à un individu grâce au barcode D qui représente le code de l'individu.

**NB** : il faut préalablement télécharger la dépendance falco qui permet de créer les sorties multiQC

```{bash, eval=FALSE, echo = TRUE}
conda install -c bioconda falco
```

**NB** : nos individus ont été indexés sur le barcode D, or dans Harpy l'indexage est réalisé sur le barcode C, ceci pose alors problème pour réussir le démultiplexage des individus. Il faut donc modifier le script du snakemake DemuxGen1 (/miniforge3/envs/harpy/bin) afin de faire passer l'information de l'individu non pas sur le barcode C mais sur le barcode D :

```{bash, eval=FALSE, echo = TRUE}
# modification des ligne l.86 et l.100 de la commande zgrep dans le fichier snakemake "DemuxGen1"
zgrep -A3 "A..C..B..{params}" {input} | grep -v "^--$" | gzip -q > {output} || touch {output}
```

**LANCEMENT DE LA COMMANDE DE DEMULTIPLEXAGE :**

```{bash, eval=FALSE, echo = TRUE}
harpy demultiplex --threads 40 --file Undetermined_L001_R1_001.fastq.gz --samplesheet Sample_Sheet_Ccory.txt
```

Compte rendu de sortie de l'étape de démultiplexage :

Un rapport multiQC est généré à la fin de l'étape de demultiplexage ([*demultiplex.QC*](https://cloud.isem-evolution.fr/nextcloud/index.php/s/MpDxmfcY4SDXFfr?path=%2FCcory%2FDEMULTIPLEX)) :

::: {style="text-align:center;"}
![](/Users/oceane/Desktop/Stage/APPENDIX/FIGURES_BIO_INFO/Demultiplex_Status_Checks.png){width="800px" height="500px"}
:::

Un deuxième rapport généré lors de l'étape de contrôle de format des fichiers ([*filecheck.fastq.html*](https://cloud.isem-evolution.fr/nextcloud/index.php/s/MpDxmfcY4SDXFfr?path=%2FCcory%2FDEMULTIPLEX)) est également accesible.

------------------------------------------------------------------------

### Harpy : Contrôle de qualité

Ensuite, une étape de filtres de contrôle de qualité des séquences permet d'éliminer les adaptateurs de séquençage, les séquences d’index ou encore les régions de mauvaise qualité.

**LANCEMENT DE LA COMMANDE DE CONTRÔLE DE QUALITÉ :**

```{bash, eval=FALSE, echo = TRUE}
harpy qc --directory ./Demultiplex/Undetermined_S0_L002/ --threads 40 --extra-params '-c'
```

Le paramètre -c permet d'activer la correction des bases dans les régions chevauchées (pour les données Paired-End).

**NB** : traitement de l'individu D35 qui présente un excès de barcode invalide (A00B00C00D35) qui sont dû à un artéfact de bio-informatique (les 6 premières bases de l'index i5 sont identiques à 1 nucléotide près au barcode D35) =\> certaines molécules sont donc attribuées à cet individu alors qu'elles ne lui appartiennent pas, il faut donc éliminer toutes les molécules présentant cette combinaison de barcode :

```{bash, eval=FALSE, echo = TRUE}
# créer un nouveau fichier bam sans les molécules avec la combinaison A00C00B00D35
samtools view -h Ccory_E2_13_195_A.bam | grep -v 'A00C00B00D35' | samtools view -bS -o Ccory_E2_13_195_A_filter.bam

# créer le fichier d'index associé au nouveau fichier bam 
samtools view -h Ccory_E2_13_195_A_filter.bam
```

Compte rendu de sortie de l'étape de contrôle de qualité :

Un premier rapport est généré à la fin de l'étape de contrôle de qualité sur l'analyse des reads ([*qc.report.html*](https://cloud.isem-evolution.fr/nextcloud/index.php/s/MpDxmfcY4SDXFfr?path=%2FCcory%2FQC)) :

::: {style="text-align:center;"}
![](/Users/oceane/Desktop/Stage/APPENDIX/FIGURES_BIO_INFO/QC_Filtered_Reads.png){width="800px" height="450px"}
:::

Un deuxième rapport sur l'analyse des barcodes ([*barcode.summary.html*](https://cloud.isem-evolution.fr/nextcloud/index.php/s/MpDxmfcY4SDXFfr?path=%2FCcory%2FQC)) est également accesible :

::: {style="text-align:center;"}
```{r}
plotdata <- read.csv("/Users/oceane/Desktop/Stage/APPENDIX/DATA_BIO_INFO/QC.csv") %>%
  select(Sample, ValidBarcodes, InvalidBarcodes) %>% 
  rename(Valid = ValidBarcodes, Invalid = InvalidBarcodes)
p <-  plotdata %>%
  pivot_longer(-Sample, names_to = "Barcode", values_to = "Count") %>%
  ggplot(aes(y = Sample, x = Count, fill = Barcode)) +
  geom_col() +
  theme_minimal() +
  theme(panel.grid.major.y = element_blank(), axis.title.y = element_blank()) +
  scale_fill_manual(values = c(Valid = "#afd5e7", Invalid = "#c02554")) +
  scale_x_continuous(labels = scales::comma) +
  labs(title = "Barcode Count and Validity")
ggplotly(p, height = 750)
```
:::

De plus, un troisème type de rapport qui contient de nombreuses informations pour chaque individu est accesible pour l'individu Ccory_AU_13_131_A ([*Ccory_AU_13_131_A.html*](https://cloud.isem-evolution.fr/nextcloud/index.php/s/MpDxmfcY4SDXFfr?path=%2FCcory%2FQC)).

------------------------------------------------------------------------

### Harpy : Alignement

Par la suite, les séquences doivent être alignées sur le génome de référence de la centaurée (Ccory_REF.fasta) avec l’algorithme BWA-BEM.

**NB** : il faut préalablement télécharger la dépendance Rmarkdown (de R) qui permet de créer des fichiers de sorties

```{bash, eval=FALSE, echo = TRUE}
conda install R
```

**NB** : l'alignement ne fonctionnent pas correctement sur tous les individus. Plus précisément, le script ne prend pas systématiquement les brins R et F pour chaque individu. En effet, de manière aléatoire (à cause d'une erreur dans le script), il prend deux fois le même brin. Il faut donc modifier le snakemake align-bwa (/miniforge3/envs/harpy/bin/) :

```{bash, eval=FALSE, echo = TRUE}
# modification de la ligne l.20 dans le script snakemake "align-bwa.smk"
lst = sorted(glob.glob(seq_dir + "/" + wildcards.sample + "*"))
```

**LANCEMENT DE LA COMMANDE D'ALIGNEMENT :**

```{bash, eval=FALSE, echo = TRUE}
harpy align bwa --genome Ccory_REF.fasta --directory ./QC --threads 40
```

Là alors, il est possible de visualiser les sorties bwa sur [IGV](https://github.com/igvteam/igv) :

(1) charger le génome de référence

(2) importer le fichier *.bam* de chaque individu

(3) visualiser les reads sous la forme liée avec l'outil *BX-tag*

Compte rendu de sortie de l'étape d'alignement :

Un premier rapport contenant des statistiques sur les barcodes est généré pour chaque individu, celui de l'individu Ccory_AU_13_131_A ([*Ccory_AU_13_131_A.bxstats*](https://cloud.isem-evolution.fr/nextcloud/index.php/s/MpDxmfcY4SDXFfr?path=%2FCcory%2FALIGN)) est accesible :

::: {style="display:flex; justify-content:center;"}
```{r}

infile = "/Users/oceane/Desktop/Stage/APPENDIX/DATA_BIO_INFO/Alignement_Individu_Ccory_AU_13_131_A.bxstats"
tb <- read.table(infile, header = T, sep = "\t") %>% select(-start, -end)
tb$valid <- tb$bx
tb[!(tb$valid %in% c("noBX", "invalidBX")), "valid"] <- "validBX"
tb$valid <- gsub("BX", " BX", tb$valid)
p <- tb %>% filter(valid == "valid BX", length_inferred > 10) %>%
    ggplot(aes(x = length_inferred)) +
        geom_histogram(binwidth = 2500, fill = "#70a7d9", boundary = 0.5) + #, color = "white") +
        labs(title = "Inferred Haplotagged Molecule Lengths", subtitle = "Ccory_AU_13_131_A.bxstats") +
        xlab("Inferred Length of Molecule (2,500 bp bins)") +
        ylab("Number of occurrences") +
        theme_light() +
        scale_y_continuous(labels = scales::comma) +
        scale_x_continuous(labels = scales::comma) +
        theme(
            panel.grid.minor.y = element_blank(),
        )
ggplotly(p, height = 500)
```
:::

Un deuxième rapport contenant des informations sur la couverture des séquences est généré pour chaque individu, celui de l'individu Ccory_AU_13_131_A ([*Ccory_AU_13_131_A.cov*](https://cloud.isem-evolution.fr/nextcloud/index.php/s/MpDxmfcY4SDXFfr?path=%2FCcory%2FALIGN)) est également accesible :

::: {style="display:flex; justify-content:center;"}
```{r}

covfile = "/Users/oceane/Desktop/Stage/APPENDIX/DATA_BIO_INFO/Alignement_Individu_Ccory_AU_13_131_A.cov"
gencov <- read.table(covfile, header = F)
gencov <- gencov[, c(1, 2, 3, 5)]
colnames(gencov) <- c("contig", "start", "end", "depth")
sname <- gsub("\\.\\w?\\w?\\w?.?cov.gz", "", basename(covfile))
q99 <- quantile(gencov$depth, 0.99)
distplot <- ggplot(gencov, aes(x = depth)) +
  geom_histogram(binwidth = 10, color = "#8484bd" , fill = "#9393d2", linewidth = 0.2) +
  theme_minimal() +
  scale_y_continuous(labels = scales::comma) +
  scale_x_continuous(labels = scales::comma) +
  xlim(NA, q99) +
  labs(title = "Distribution of Alignment Depths") +
  xlab("Depth") +
  ylab("Count")
ggplotly(distplot,  height = 500) 
```
:::

Un troisème rapport global ([*bwa.stats*](https://cloud.isem-evolution.fr/nextcloud/index.php/s/MpDxmfcY4SDXFfr?path=%2FCcory%2FALIGN)) est également accesible.

------------------------------------------------------------------------

### Harpy : Appel de variants

Une fois les séquences alignées, les fichiers d’alignement (au format *.bam*) sont utilisés pour identifier les variants, c’est-à-dire appeler les SNPs en utilisant le ‘variant caller’ mpileup multiple (c'est-à-dire, qui prend en compte la structure sous forme de différentes populations) développés par bcftools

**LANCEMENT DE LA COMMANDE D'APPEL DES VARIANTS (MPILEUP) :**

```{bash, eval=FALSE, echo = TRUE}
harpy snp mpileup --threads 80 --genome /mnt/Data3/DEVOCGEN/REFERENCE_GENOME/Ccory_REF.fasta --directory Align/bwa --windowsize 250000
```

La valeur du paramètre --windowsize (par défaut fixé à 50000) doit être augmenté car le génome de la centaurée est grand (765Mb).

Compte rendu de sortie de l'étape d'appel des variants :

Un premier rapport contenant des informations sur les données brutes ([*variants.raw.html*](https://cloud.isem-evolution.fr/nextcloud/index.php/s/MpDxmfcY4SDXFfr?path=%2FCcory%2FVARIANTS)) est accesible :

::: {style="display:flex; justify-content:center;"}
```{r}
dataL <- read.table("/Users/oceane/Desktop/Stage/APPENDIX/DATA_BIO_INFO/Alignement.csv", sep = ",", header = TRUE)
psc = dataL
tidy_psc <- pivot_longer(psc[,c(1,3,6,7,8,9,10,11)], -Sample , names_to = "Metric", values_to = "Count")
pscplot <- ggplot(data = tidy_psc, mapping = aes(x = Count, y = Sample, color = Metric)) +
  geom_point(size = 2) +
  labs(title = "Individual Statistics") +
  theme_bw() +
  scale_x_continuous(n.breaks=9) +
  xlab("Count") +
  scale_color_brewer(palette = "Dark2") +
  theme(panel.grid.minor = element_blank(), panel.grid.major.y = element_blank(), axis.title.y = element_blank())

ggplotly(pscplot, height = 500) %>% partial_bundle()

```
:::

Un deuxième rapport sur les données normalisées ([*variants.normalized.html*](https://cloud.isem-evolution.fr/nextcloud/index.php/s/MpDxmfcY4SDXFfr?path=%2FCcory%2FVARIANTS)), est également accesible.

------------------------------------------------------------------------

### Harpy : Imputation

Une fois l’étape d’appel de variants terminée, l’imputation des génotypes peut être réalisée. Cette étape, réalisée par STITCH va permettre d’imputer les informations manquantes dans les données génétiques.

Avant tout, une filtration des SNP est réalisée afin d'éliminer les SNP avec Minimum Allele Count (MAC)\<2, ces SNPs sont considérés comme non informatifs pour l'imputation des génomes :

```{bash, eval=FALSE, echo = TRUE}
bcftools filter -e "MAC<2" variants.raw.bcf -o /mnt/Data3/DEVOCGEN/RAW/Ccory/Desmarais_E_02_bclconvert_out/Variants_filtered/variants_filtered_2.bcf -O b
```

Tout d'abord, il faut créer un fichier de paramètres STITCH avec les valeurs choisies pour chacun des pramètres :

```{bash, eval=FALSE, echo = TRUE}
# création d'un fichier de paramètres stitch
harpy stitchparams -o Ccory_STITCH_params_50000_10_5_80 
```

Le choix des différents paramètres a été optimisé en testant plusieurs valeurs de paramètres puis en évaluant leurs performances en visualisant la distribution de l'info score (ci-dessous). Les paramètres optimisant la performance de l’imputation sont method = diploid, bxlimit = 50000, K = 30, S = 5, nGen = 30.

**LANCEMENT DE LA COMMANDE D'IMPUTATION DES GÉNOMES :**

```{bash, eval=FALSE, echo = TRUE}
harpy impute -v /mnt/Data3/DEVOCGEN/RAW/Ccory/Desmarais_E_02_bclconvert_out/Variants/mpileup/variants.raw.bcf -p /mnt/Data3/DEVOCGEN/RAW/Ccory/Desmarais_E_02_bclconvert_out/Ccory_STITCH_params_50000_10_5_80 -t 80 -d /mnt/Data3/DEVOCGEN/RAW/Ccory/Desmarais_E_02_bclconvert_out/Align/bwa
```

Compte rendu de sortie de l'étape d'imputation des génotypes :

Un premier rapport contenant des informations sur les données brutes ([*variants.imputed.html*](https://cloud.isem-evolution.fr/nextcloud/index.php/s/MpDxmfcY4SDXFfr?path=%2FCcory%2FIMPUTE)) est accesible à l'échelle populationnelle :

::: {style="text-align:center;"}
![](/Users/oceane/Desktop/Stage/APPENDIX/FIGURES_BIO_INFO/Impute_Info_Score.png){width="700px" height="350px"}
:::

Ce rapport contient également des informations sur les données brutes à l'échelle individuelle :

::: {style="text-align:center;"}
```{r, fig.height=8}
comparefile = "/Users/oceane/Desktop/Stage/APPENDIX/DATA_BIO_INFO/impute.compare.stats"
gcts <- read.table(comparefile, header = F)[, c(2, 23, 25, 24, 27)]
names(gcts) <- c("sample","Homozygote Ref", "Heterozygote",	"Homozygote Alt", "Missing")
gcts_long <- pivot_longer(gcts, -1, names_to = "Conversion", values_to = "Count")

gct_plot = ggplot(gcts_long, aes(y =  sample, x = Count, fill = Conversion)) +
              geom_col(position = "stack") +
              scale_x_continuous(expand = c(0,0)) +
              ylab("") +
              theme_minimal() +
              theme(panel.grid.major.y = element_blank(), axis.ticks.y = element_line(size = 0.5)) +
              scale_fill_manual(
                labels = c("→ Homozygous Alt", "→ Heterozygote", "→ Missing", "→ Homozygous Ref"),
                values = c("#93d2fa", "#83afa8", "#424043", "#ffd75f")
              ) +
              labs(title = "Missing Genotype Imputations")

ggplotly(gct_plot, heigth = 800)
```
:::

------------------------------------------------------------------------

### Harpy : Phasage

Il est maintenant possible de se servir des données imputées afin de pahser les génomes. Cette étape, réalisée par HAPCUT2, permet de mettre en phase les génotypes en haplotypes.

Avant de phaser les données, une étape de filtration des données est nécessaire afin de ne garder que :

(1) les SNPs avec un INFO SCORE \> 0.5

```{bash, eval=FALSE, echo = TRUE}
bcftools view -i 'INFO/INFO_SCORE > 0.7' /mnt/Data3/DEVOCGEN/RAW/Ccory/Desmarais_E_02_bclconvert_out/Impute/modeldiploid_usebxTrue_bxlimit50000_k30_s1_ngen30/variants.imputed.bcf -o variants.imputed_filtered_info_score.bcf -O b 
```

(2) les SNPs avec un taux de missingness \< 0.5

```{bash, eval=FALSE, echo = TRUE}
vcftools --bcf /mnt/Data3/DEVOCGEN/RAW/Ccory/Desmarais_E_02_bclconvert_out/Impute/modeldiploid_usebxTrue_bxlimit50000_k30_s1_ngen30/variants.imputed_filtered_info_score.bcf --max-missing 0.8 --recode-bcf --out variants.imputed_filtered
```

(3) les SNPs avec un Minimum Allele Count \> 2

```{bash, eval=FALSE, echo = TRUE}
bcftools filter -e "MAC<2" /mnt/Data3/DEVOCGEN/RAW/Ccory/Desmarais_E_02_bclconvert_out/Impute/modeldiploid_usebxTrue_bxlimit50000_k30_s1_ngen30/variants.imputed.bcf -o variants.imputed_filtered_info_score.bcf -O b 
```

**LANCEMENT DE LA COMMANDE D'IMPUTATION DES GÉNOMES :**

```{bash, eval=FALSE, echo = TRUE}
harpy phase --threads 50 --vcf /mnt/Data3/DEVOCGEN/RAW/Ccory/Desmarais_E_02_bclconvert_out/Impute/modeldiploid_usebxTrue_bxlimit50000_k30_s1_ngen30/variants.imputed_filtered.bcf --directory /mnt/Data3/DEVOCGEN/RAW/Ccory/Desmarais_E_02_bclconvert_out/Align/bwa
```

Compte rendu de sortie de l'étape d'appel des variants :

Un rapport contenant des informations sur les génotypes phasés ([*phase.html*](https://cloud.isem-evolution.fr/nextcloud/index.php/s/MpDxmfcY4SDXFfr?path=%2FCcory%2FPHASE)) est accesible :

::: {style="display:flex; justify-content:center;"}
```{r, fig.height=2.7}
file = "/Users/oceane/Desktop/Stage/APPENDIX/DATA_BIO_INFO/blocks.summary"
df <- read.table(
  file, header = T,
  colClasses = c("factor","factor","integer","integer", "integer")
 )
df_plot = ggplot(df, aes(block_length, after_stat(count))) + 
  geom_density(color = "#859f9a", fill = "#bfe4dc") +
  theme_minimal() +
  theme(
    panel.grid.minor.y = element_blank(),
    panel.grid.minor.x = element_blank()
    ) +
   scale_y_continuous()  +
   scale_x_log10(
    breaks = trans_breaks("log10", function(x) 10^x),
    labels = trans_format("log10", math_format(10^.x))
  ) +
  xlab("Haplotype Length (log scale)") +
  ylab("Count") +
  labs(title = "Distribution of Haplotype Lengths")

ggplotly(df_plot, height = 500)
```
:::

L'ensemble des étapes du pipeline Harpy sont maintenant terminées !

Avant de pouvoir utiliser le fichier *.bcf* pour obtenir les segements IBD, plusieurs étapes de traitement des données vont être réalisés :

**(1)** Obtention de statistique de phasage avec [*PhaseME*](https://github.com/sinamajidian/phaseme)

Cette étape va permettre d'obtenir des statistiques de phasage pour chaque chromosome dont la statistique N50 qui est intéressante pour avoir une idée de la qualité du phasage. PhaseMe ne travaille que sur 22 chromosomes numérotés de 1 à nchr, il faut donc modifier le *.bcf* pour l'adapter à ce critère :

```{bash, eval=FALSE, echo = TRUE}
# téléchargement des scripts phaseme.py et stats.py sur le Github 

# split du .bcf en 3 sous-jeux de données : commande suivante pour les chromosomes allant de 258 à 10026
bcftools filter -i 'CHROM="Super-Scaffold_258" || CHROM="Super-Scaffold_294" || CHROM="Super-Scaffold_584" || CHROM="Super-Scaffold_607" || CHROM="Super-Scaffold_623" || CHROM="Super-Scaffold_700" || CHROM="Super-Scaffold_1212" || CHROM="Super-Scaffold_1225" || CHROM="Super-Scaffold_2367" || CHROM="Super-Scaffold_100002" || CHROM="Super-Scaffold_100004" || CHROM="Super-Scaffold_100006" || CHROM="Super-Scaffold_100007" || CHROM="Super-Scaffold_100009" || CHROM="Super-Scaffold_100010" || CHROM="Super-Scaffold_100016" || CHROM="Super-Scaffold_100017" || CHROM="Super-Scaffold_100018" || CHROM="Super-Scaffold_100018" || CHROM="Super-Scaffold_100024" || CHROM="Super-Scaffold_1000025" || CHROM="Super-Scaffold_100026"' variants.phased.vcf -o variants.phased_258_100026.vcf -O z

# création de fichiers .txt contenant les correspondances de nom entre le format "Super-Scaffold_258" et le nouveau nom au format "1" pour chacun des sous-jeux de données : exemple du format pour les chromosomes allant de 258 à 10026 : corresp_258_100026.txt
Super-Scaffold_258    1
Super-Scaffold_294    2
Super-Scaffold_584    3
...
Super-Scaffold_100026   22

# renommer les noms des Super-Scaffolds dans le .bcf
bcftools annotate --rename-chrs corresp_258_100026.txt variants.phased_258_100026.vcf -O z -o variants.phased_rename_258_100026.vcf
```

**LANCEMENT DE LA COMMANDE PHASE ME :**

```{bash, eval=FALSE, echo = TRUE}
# lancement de la commande phaseME
python3 phaseme.py improver variants.phased_rename_258_100026.bcf out_258_100026

# merger dans Calc les 3 fichiers de sortie quality.csv en 1 seul
```

**(2)** Filtration des chromosomes 

L'étape précédente a permis de mettre en évidence que certains Super-Scaffolds avaient des valeurs de N50 très faible (entre 0 et 5 Kb). En regardant leurs longueurs, il s'agit de Super-Scaffolds généralement de petite longueur (< 5 Mb). Les Super-Scaffolds avec une taille inférieur à 5Mb vont donc être éliminée du jeu de données étant donné que des segments IBD de minimum 1 Mb sont visés. 

```{bash, eval=FALSE, echo = TRUE}
# élimination des chromosomes inférieurs à 5 Mb
bcftools filter -e 'CHROM="Super-Scaffold_100052" ||CHROM="Super-Scaffold_585" || CHROM="Super-Scaffold_100025" || CHROM="Super-Scaffold_137" || CHROM="Super-Scaffold_100033" || CHROM="Super-Scaffold_100087" || CHROM="Super-Scaffold_28" || CHROM="Super-Scaffold_100046" ||CHROM="Super-Scaffold_19" || CHROM="Super-Scaffold_35" || CHROM="Super-Scaffold_96" || CHROM="Super-Scaffold_100095" || CHROM="Super-Scaffold_258" || CHROM="Super-Scaffold_100060" || CHROM="Super-Scaffold_100183" || CHROM="Super-Scaffold_100161" || CHROM="Super-Scaffold_233" || CHROM="Super-Scaffold_100143" || CHROM="Super-Scaffold_190" || CHROM="Super-Scaffold_100215" || CHROM="Super-Scaffold_100166" || CHROM="Super-Scaffold_294"' variants.phased.vcf -o variants.phased_without_scaffold.vcf -O z
```

**(3)** Phasage à l'échelle chromosomique avec Shapeit4

En sortie de HapCut2, les haplotypes ne sont pas phasés à l'échelle chromosomique, c'est-à-dire que les haplotypes phasés ne sont pas phasés entre eux, cette étape réalisée par SHAPEIT4, permet donc un phasage "chromosome-scale" :

```{bash, eval=FALSE, echo = TRUE}
# installation de la version 4.2.1 de Shapeit 
mamba install -c bioconda shapeit4=4.2.1 
```

Le logiciel travaille sur les chromosomes 1 par 1, des boucles pour traiter les chromosomes à la suite, les indexer et les merger ont donc été générées : 

**LANCEMENT DE LA COMMANDE SHAPE IT 4 :**

```{bash, eval=FALSE, echo = TRUE}
# création d'une boucle pour automatiser la commande shapeit4 sur chaque chromosome 1 par 1 : shapeit.sh 
awk '{print}' "$1" /mnt/Data3/DEVOCGEN/Ccory/Desmarais_E_02_bclconvert_out/List_Chr.csv | while read chromosome; do     shapeit4 --input /mnt/Data3/DEVOCGEN/Ccory/Desmarais_E_02_bclconvert_out/Phase/variants.phased_without_scaffold.vcf --sequencing --region $chromosome --use-PS 0.0001 --thread 50 --log phased_shapeit.log --output /mnt/Data3/DEVOCGEN/Ccory/Desmarais_E_02_bclconvert_out/Phase/SHAPEIT/variants.phased_shapeit_$chromosome.bcf; 
  done

# lancement de la boucle shapeit.sh 
bash shapeit.sh 

# création d'une boucle pour indexer tous les fichiers shapeit : shapeit_index.sh 
for fichier in /mnt/Data3/DEVOCGEN/Ccory/Desmarais_E_02_bclconvert_out/Phase/SHAPEIT/variants.phased_shapeit/variants.phased_shapeit_Super-Scaffold_*.bcf; do 
  echo "Indexation du fichier {fichier}..."
  bcftools index ${fichier}
  done 
  
# lancement de la boucle shapeit_index.sh
bash shapeit_index.sh

# création d'une boucle pour concaténer tous les fichiers shapeit en 1 seul : shapeit_concat.sh 
fichiers_a_fusionner = ""
for fichier in /mnt/Data3/DEVOCGEN/Ccory/Desmarais_E_02_bclconvert_out/Phase/SHAPEIT/variants.phased_shapeit/variants.phased_shapeit_Super-Scaffold_*.bcf; do 
  fichiers_a_fusionner = "${fichiers_a_fusionner} ${fichier}"
  done 
echo "Fusion des fichiers BCF"
bcftools concat -O z -o "/mnt/Data3/DEVOCGEN/Ccory/Desmarais_E_02_bclconvert_out/Phase/SHAPEIT/variants.phased_shapeit/variants.phased_shapeit.bcf" ${fichiers_a_fusionner}

# lancement de la boucle shapeit_concat.sh
bash shapeit_concat.sh
```

------------------------------------------------------------------------

### Segments IBD : IBD-Caller

Les données sont maintenant phasés à l'échelle chromosomique, les segments IBD vont pouvoir être obtenus. Pour cela, 2 méthodes différentes vont être utilisées : 

(1) méthode IBD-Caller avec FastSMC / HapNe

(2) méthode basé sur l'ARG avec ts_infer/ts_date

Pour la première méthode FastSMC, un fichier de l'évolution de la taille de population dans les dernières décennies et des fichiers *.hap.gz* pour chaque Super-Scaffolds sont nécessaires. Il est possible d'obtenir ces fichiers via le logiciel HapNe en utilisant la méthode "ld". 

```{bash, eval=FALSE, echo = TRUE}
# création d'un environnement mamba : HapNe
mamba create -n HapNe 

# téléchargement de plink2 nécessaire à HapNe
mamba install -c bioconda plink2

# téléchargement de HapNe
pip install hapne
```

Ce logiciel a spécialement été conçu pour traiter des données de génétique humaine, il y a donc un certains nombres de modifications à effectuer afin de pouvoir faire tourner le logiciel correctement. Tout d'abord, un des fichiers à confectionner est un fichier de configuration : 

```{bash, eval=FALSE, echo = TRUE}
# création du fichier Config
[CONFIG]
vcf_file = /mnt/Data3/DEVOCGEN/IBD/HapNe/VCF/variants.phased_rename_without_scaffold_E2
map = /mnt/Data3/DEVOCGEN/IBD/HapNe/MAPS/genetic_map_@.txt
pseudo_diploid = FALSE
output_folder = /mnt/Data3/DEVOCGEN/IBD/HapNe
genome_build = grch37
method = ld
```

Un des fichiers à avoir également est une carte génétique pour chaque chromosome avec 3 colonnes (Position génétique | Taux de recombinaison | Position physique) : 
```{bash, eval=FALSE, echo = TRUE}
# création des cartes génétiques de chaque chromosome
bcftools view variants.phased_shapeit.vcf | grep 'Super-Scaffold_100002' | bcftools query -f '%POS\t%INFO/CM\n' variants.phased_shapeit.vcf | awk -F '\t' '{print $1, 1, $2}' > genetic_map_Super-Scaffold_100002.txt
```

Un fichier *.vcf* avec des noms de chromosomes numérotés de 1 à nchr : 
```{bash, eval=FALSE, echo = TRUE}
# création d'un fichier de correspondance .txt au format :
Super-Scaffold_100002   1
Super-Scaffold_100004   2
Super-Scaffold_100006   3
... 
Super-Scaffold_75   34

# renommer les chromosomes dans le .vcf
bcftools annotate --rename-chrs corresp_without_scaffold.txt variants.phased_without_scaffold.vcf -O z -o variants.phased_rename_without_scaffold.vcf.gz
```

En parallèle, il faut également modifier des parties des différents scripts et fichiers dans "/miniforge/envs/FastSMC/lib/python3.12/site-packages/hapne" : 

Tout d'abord, il faut modifier le fichier regions_grch37.txt :

-   CHR : nom des Super-Scaffolds préalablement renommé de 1 à nchr

-   NAME : un identifiant pour chaque Super-Scaffold : ex: SS100002

-   FROM_BP : la position du premier SNP du Super-Scaffold

-   TO_BP : la position du dernier SNP du Super-Scaffold

Les 2 dernières valeurs sont obtenues grâce à ces commandes : 
```{bash, eval=FALSE, echo = TRUE}
bcftools view variants.phased_shapeit.vcf | grep '^Super-Scaffold_100002' | less  # => position du premier SNP du Super-Scaffold_100002
bcftools view variants.phased_shapeit.vcf | grep '^Super-Scaffold_100002'         # => position du dernier SNP du Super-Scaffold_100002
```

Le script Python utils.py : 
```{bash, eval=FALSE, echo = TRUE}
# modification de la ligne 59 
years_per_gen = 5     # à adapter en fonction du temps de génération de l'espèce

# modification de la ligne 161
n = n.to_numpy()
```

Le script Python tools.py : 
```{bash, eval=FALSE, echo = TRUE}
# modification des valeurs représentant le nombre de chromosomes 
for ii in range(34)     # ligne 110
if not (0 <= index < 34)    # ligne 116
for ii in range(34)     # ligne 162
if not (0 <= index < 34)    # ligne 168
for ii in range(34)     # ligne 212
if not (0 <= index < 34)    # ligne 218

# modification de la commande plink pour spliter le vcf => ajouter après la ligne 153
--allow-extra-chr \ 
--autosome-num 34 \ 
--export haps bgz"   # permet de générer les fichiers .hap.gz
```

Le script HapNe.py : 
```{bash, eval=FALSE, echo = TRUE}
# modification des lignes 174 et 176
self.sig2 = ((1 + yearly_rate_growth) ** 5 -1)**2     # à modifier en fonction du temps de génération de l'espèce
self.butter_time = self.t_max - 4   # à adapter en conséquence du terme précédent 
```

**LANCEMENT DE LA COMMANDE HAPNE :**

```{bash, eval=FALSE, echo = TRUE}
python HapNe.py --config_file ./Config
```

------------------------------------------------------------------------

### Segments IBD : IBD-caller

Pour faire tourner FastSMC, il faut comme fichiers d'entrées les haplotypes phasés *.hap.gz*, les cartes génétiques *.map* générées précédemment, ainsi qu'un fichier 
*.decodingQuantities.gz* qui est obtenu grâce au script Decoding_quantities disponibles dans le package ASMC : 

```{bash, eval=FALSE, echo = TRUE}
# création d'un environnement mamba : ASMC 
conda create -n ASMC

# téléchargement du package ASMC qui contient à fois ASMC, FASTSMC et DECODING_QUANTITIES
pip install asmc 
```

Le fichier *.decodingQuantities.gz* peut être créer à partir d'un fichier de l'évolution de la démographie dans les dernières générations*.demo*, d'un fichier de discrétisation du temps *.disc* et d'un fichier de fréquences des allèles *.frq*.

Les fichiers *.demo* et *.disc* doivent être créés manuellement : 
```{bash, eval=FALSE, echo = TRUE}
# création du fichier .demo au format : stable.demo
0   2000
50000   2000 

# création du fichier .disc au format crée à partir du fichier 30-100-2000_CEU.disc disponible sur le github: 0_50846.disc
0.0
10.0
20.0
30.0
40.0
... 
50846.1
```

Le fichier *.frq* doit être au format : 

-   CHROM : nom des Super-Scaffolds préalablement renommé de 1 à nchr

-   SNP : un identifiant pour chaque SNP : ex: 1_SNP131

-   REF : l'allèle de référence

-   ALT : l'allèle alternatif 

-   MAF : la fréquence de l'allèle minoritaire 

-   NCHR : le nombre d'haplotypes portant l'allèle

Toutes ces informations peuvent être obtenus grâce aux logiciels BCFtools et VCFtools : 

```{bash, eval=FALSE, echo = TRUE}
# obtention des colonnes REF et ALT 
bcftools query -f '%REF\t%ALT\n' variants.phased_rename_without_scaffold.vcf -o REF_ALT.txt

# obtention des colonnes CHROM, MAF et NCHR 
vcftools --vcf variants.phased_rename_without_scaffold.vcf --freq --out MAF.freq
awk '{print $6}' MAF.frq | sed 's/[CGTA]://' > MAF.txt

# obtention de la colonne SNP
awk 'NR > 1 {print $1"_"$2}' MAF.frq > ID_SNP.txt

# merger toutes les colonnes pour créer le fichier 
paste ID_SNP.txt REF_ALT.txt MAF.frq MAF.txt| awk 'BEGIN {OFS="\t"} {print $4, $1, $2, $3, $7, $10}' > frequences.frq
```

**LANCEMENT DE LA COMMANDE DECONDING_QUANTITIES :**

```{bash, eval=FALSE, echo = TRUE}
import pathlib

from asmc.preparedecoding import *

files_dir = (pathlib.Path('/mnt/DATA3/DEVOCGEN/IBD/FastSMC/DECODING_QUANTITIES') / 'test' / 'regression').resolve()

demo_file = str(files_dir / 'stable.demo')
disc_file = str(files_dir / '0_50846.disc')
freq_file = str(files_dir / 'frequences.frq')

# Calcul avec un fichier .disc
dq = prepare_decoding(
    demography=demo_file,
    discretization=disc_file,
    frequencies=freq_file,
    samples=60, 
)
```

Il est maintenant possible de faire tourner FastSMC avec comme fichiers d'entrées les haplotypes phasés *.hap.gz*, les cartes génétiques *.map*, les deconding quantities *.decodingQuantities.gz* générés précédemment. 

**LANCEMENT DE LA COMMANDE FAST SMC :**

------------------------------------------------------------------------

### Segments IBD : ARG

Pour la deuxième méthode basé sur l'utilisation de l'information contenu dans l'ARG (Graphe de Recombinaison Ancestral), le couple de logiciel ts-infer/ts-date sont utilisées. 

```{bash, eval=FALSE, echo = TRUE}

# téléchargement des 2 logiciels 

```

Le fichier d'entrée pour ts-infer est un *.bcf* phasé et orienté, il faut donc orienté notre *.bcf* grâce à des génomes de référence de groupe externe. Dans notre cas, les génomes de *Centaurea solstitialis*, *Cynara cardunculus* et *Silybum marianum*. 

**LANCEMENT DE LA COMMANDE TS-INFER :**

------------------------------------------------------------------------

# SIMULATIONS {data-icon="fa-sitemap"}

## Row {.tabset}

### Pipeline {data-height="15000"}

::: {style="text-align:center;"}
![](/Users/oceane/Desktop/Stage/FIGURES/Pipeline_Simulations.png){width="650px" height="800px"}
:::

------------------------------------------------------------------------

### Obtention des segments IBD

(1) Simulation en backward (coalescence) avec [*msprime*](https://tskit.dev/msprime/docs/stable/intro.html)

(2) Simulation en forward avec [*SLIM*](https://github.com/MesserLab/SLiM)

(3) Obtention des segments IBD between et within avec [*tskit*](https://tskit.dev/tskit/docs/stable/introduction.html)

```{python, eval=FALSE, echo = TRUE}
import msprime, tskit, pyslim, random, subprocess, os
import numpy as np
import pandas as pd
import matplotlib
from matplotlib import pyplot as plt


############################################################
################ SIMULATION PAR COALESCENCE ################
############################################################

demography = msprime.Demography()

# initialisation des 5 sous-populations
demography.add_population(name="p0", initial_size=1000)
demography.add_population(name="p1", initial_size=200)
demography.add_population(name="p2", initial_size=200)
demography.add_population(name="p3", initial_size=200)
demography.add_population(name="p4", initial_size=200)
demography.add_population(name="p5", initial_size=200)
all_pop = ["p1", "p2", "p3", "p4", "p5"]

# évènement de coalescence en une seule population 
demography.add_population_split(time=10000, derived=all_pop, ancestral="p0")

# définition du taux de migration
# tirer 10 valeurs de m équitablement réparties dans la plage de données [0.00037, 0.012] 
'''
borne_inf = 0.00037			# définition de la borne inférieure 
borne_sup = 0.012			# définition de la borne supérieure 
nb_valeurs_m = 10			# définition du nombre de valeurs souhaitées 
increment_m = (borne_sup - borne_inf) / (nb_valeurs_m-1)  # calcul de la valeur de l'incrément 
vecteur_m = [borne_inf + i * increment_m for i in range(nb_valeurs_m)]  # création du vecteur avec toutes les valeurs de m calculé en ajoutant la valeur de l'incrément à chaque fois 
m = round(vecteur_m[9], 5)		# récupérer une valeur de m et l'arrondir à 3 chiffres après la virgule : la valeur [] est à changer pour chaque simulation

# tirer les valeurs de m équitablement réparties dans chaque intervalle log10 
intervalle_10_4 = [0.0004, 0.0005, 0.0006, 0.0007, 0.0008, 0.0009]     		   # créer l'intervalle de 10-4 
intervalle_10_3 = [0.001, 0.002, 0.003, 0.004, 0.005, 0.006, 0.007, 0.008, 0.009]  # créer l'intervalle de 10-3 
intervalle_10_2 = [0.01]							   # créer l'intervalle de 10-2
m = round(intervalle_10_3[8], 5)	# récupérer une valeur de m et l'arrondir à 5 chiffres après la virgule : la valeur [] est à changer pour chaque simulation
'''
# tirer les valeurs de m équitablement réparties dans les intervalles 'vides' en terme de temps  
intervalle_1 = [0.000416, 0.000432, 0.000448, 0.000464, 0.000480]     		# créer l'intervalle 1 
intervalle_2 = [0.00094, 0.0011, 0.00124, 0.00138, 0.00152]  			# créer l'intervalle 2 							   
m = round(intervalle_2[4], 6)	# récupérer une valeur de m et l'arrondir à 6 chiffres après la virgule : la valeur [] est à changer pour chaque simulation

demography.set_symmetric_migration_rate(all_pop, m) 	# définir le taux de migration comme symétrique entre toutes les paires de populations

demography.sort_events()

# simulation par coalescence avec msprime
ts = msprime.sim_ancestry(
    samples={"p1":200, "p2":200, "p3":200, "p4":200, "p5":200, },
    demography=demography,
    sequence_length = 85000000, 
    recombination_rate = 1e-8, 
)

# acquisition des metadata
ts = pyslim.annotate(ts, model_type="WF", tick=1)

# ajout des mutations (nécessaire car calcul des FST dans SLiM)
ts = msprime.sim_mutations(ts, rate=1e-8, model=msprime.SLiMMutationModel(type=1))

# création du tree sequence obtenu par coalescence 
ts.dump("coal.trees")


############################################################
################### SIMULATION DANS SLIM ###################
############################################################

# run du script SLiM
subprocess.run(["slim", "Script.slim"])

# renommer le tree sequence de sortie de SLiM en incluant la valeur de l'ID, du t et du m de la simulation
dossier = '/DATA/sdb1/OCEANE'			  	# définir le dossier dans lequel il y a le fichier de tree sequence de sortie SLiM
valeur_ID = []
valeur_t = []
for nom_fichier in os.listdir(dossier):		
	if 'ID_t' in nom_fichier:		  	# si le nom de fichier continent l'expression "ID_t"
		split = nom_fichier.split('ID_t') 	# spliter le nom de fichier au niveau de l'expression "ID_t" (qui est propre au fichier de tree sequence de sortie SLiM)
		valeur_1 = split[0]		 	# récupérer la première valeur du split qui contient "forw_valeurID"
		valeur_1_split = valeur_1.split('_')	# spliter cette première valeur au niveau du "_" 
		valeur_ID = int(valeur_1_split[1])	# récupérer l'entier de la valeur de l'ID 
		valeur_2 = split[1]			# récupérer la deuxième valeur du split qui contient "valeurt.trees"
		valeur_2_split = valeur_2.split('.')	# spliter cette première valeur au niveau du "."
		valeur_t = int(valeur_2_split[0])	# récupérer l'entier de la valeur de t
		nouveau_nom = f'forw_ID{valeur_ID}_m{m}_t{valeur_t}.trees'  # créer le nouveau nom en rajoutant la valeur de m de la simulation
		os.rename(nom_fichier, nouveau_nom)	# renommer le fichier avec son nouveau nom


############################################################
################# SAMPLE DU TREE SEQUENCE ##################
############################################################

# chargement du tree sequence obtenu après simulation dans SLiM
ts_forw = tskit.load(f'forw_ID{valeur_ID}_m{m}_t{valeur_t}.trees')

# sous-échantillonnage de 20 individus dans chacune des 5 sous-populations
sample_ind = [] 				# création d'un vecteur pour stocker les individus échantillonnés
list_ind = ts_forw.individuals()  		# récupération dans un vecteur de tous les individus de la simulation
arrays_dict = {}  				# création d'un dictionnaire pour stocker les tableaux numpy
rng = np.random.default_rng()			# création d'un générateur de nombres aléatoires 
for i in range(1,6):   				# boucle sur les 5 populations
	array_name = f"list_id_{i}" 		# définition du nom du tableau contenant les individus de la population i 
	arrays_dict[array_name] = np.array([])  # création d'un nouveau tableau numpy pour la population i 
	arrays_dict[array_name] = arrays_dict[array_name].astype(int)	# conversion du tableau numpy en entiers
	for y in range(0, 1000):		# boucle sur les 1000 individus de la simulation
		ind = list_ind[y]		# création d'un vecteur pour stocker les informations de chaque individu y 
		if ind.metadata['subpopulation'] == i:	# si la population de l'individu y est égal à la population i, alors : 
			arrays_dict[array_name] = np.append(arrays_dict[array_name], ind.id)	# remplir le tabelau de la pop i avec l'ID de l'individu 
	keep_inds = rng.choice(arrays_dict[array_name], 20, replace = False)	# tirer aléatoirement 20 individus dans le tableau de la pop i 
	sample_ind.extend(keep_inds)		# stocker ce tirage dans le vecteur de sous_ech 

# extraire les noeuds des individus échantillonnés 
keep_nodes = []             # création d'un vecteur pour stocker les noeuds conservés
for i in sample_ind:        # boucle sur les individus échantillonnés  
	keep_nodes.extend(ts_forw.individual(i).nodes)

# simplifier le tree sequence original pour ne garder que l'information des individus échantillonnés  
ts_forw_simpl = ts_forw.simplify(keep_nodes) 

# création du tree sequence obtenu après simplification (avec info de m et de t) 
ts_forw_simpl.dump(f'forw_ID{valeur_ID}_m{m}_t{valeur_t}_simpl.trees')


############################################################
###################### MÉTAPOPULATION ######################
############################################################

'''
### Obtention des segments IBD ###

# chargement du tree sequence obtenu après simplification dans SLiM
ts_forw_simpl = tskit.load(f'forw_ID{valeur_ID}_m{m}_t{valeur_t}_simpl.trees')

# afficher les statistiques du tree sequence
print(f"Loaded a simulated tree sequence of {ts_forw_simpl.num_samples} haploid genomes from {ts_forw_simpl.num_individuals} diploid individuals:")
print(f"containing {ts_forw_simpl.num_sites} variable sites over {ts_forw_simpl.sequence_length} base pairs")

# obtenir les segements IBD depuis le tree sequence simplifié 
segments = ts_forw_simpl.ibd_segments(store_pairs=True, store_segments=True)

Pair = []
Left = np.array([])
Right = np.array([])
Node = np.array([])
Time = np.array([])

for k, v in segments.items():
    kLeft = v.left 
    kRight = v.right
    kNode = v.node
    order = kLeft.argsort()
    l = kLeft[order]
    n = kNode[order]
    r = kRight[order]
    i = np.where(np.convolve(n + (n[0] == 0), [-1, 1]) != 0)[0][:-1]
    kPair = [k, ]*len(i)
    kLeft = l[i]
    kRight = np.append(r[i[1:]-1], r[len(r)-1])
    kNode = n[i]   
    Pair = Pair + kPair
    Left = np.append(Left, kLeft)
    Right = np.append(Right, kRight)
    Node = np.append(Node, kNode)
    Time = np.append(Time, ts.tables.nodes.time[kNode])
Length = np.subtract(Right,Left)
dfIBD = pd.DataFrame(columns = ['Pair', 'Left', 'Right', 'Length', 'Node', 'Time'])
dfIBD['Pair'] = Pair
dfIBD['Left'] = Left
dfIBD['Right'] = Right
dfIBD['Length'] = Length
dfIBD['Node'] = Node
dfIBD['Time'] = Time

Num_IBD = len(dfIBD.Pair)
print(f"Recorded {Num_IBD} IBD segments")

dfIBD.to_csv('out_IBD_{:.3f}_smplf.csv'.format(m))


### Visualisation des segments IBD ###

# création d'un histogramme du nombre de segements IBD en fonction de leur longueur
res = plt.hist(dfIBD["Length"], range = (0, ts_forw_simpl.sequence_length), bins=40, color = 'thistle', edgecolor = 'white')
plt.xlabel('Longueur')1202421
plt.ylabel('Nombre')
plt.title('IBD')
plt.yscale('log')
plt.close()
print(f"Generated IBD segments length distribution plot")

# ajout de la courbe analytique calculé à partir de Palamara 2012
number_segments = len(dfIBD)
limits = list(map(int,res[1]))
m = np.add(limits[:-1],limits[1:])/2
M = m.tolist()

C = []
for i in range(len(limits)-1):
    u = limits[i]*10**(-6)
    v = (limits[i+1]-1)*10**(-6)
    proba = 50**2*(((50+ts_forw_simpl.num_individuals*u)**(-2))-((50+ts_forw_simpl.num_individuals*v)**(-2)))
    C = C + [proba * number_segments]

plt.hist(dfIBD["Length"], range = (0, ts_forw_simpl.sequence_length), bins=40, color = 'thistle', edgecolor = 'white')
plt.xlabel('Length')
plt.ylabel('Number')
plt.title('IBD')
plt.scatter(M, C, marker = '+', color = 'purple')
plt.yscale('log')
plt.savefig("/DATA/sdb1/OCEANE/IBD_length_smplf.png")
plt.show()
'''

############################################################
#################### BETWEEN VS. WITHIN ####################
############################################################

# créer des sous-tableaux contenant les haplotypes appartenant à la même population 
pop_arrays_dic = {}			# création d'un dictionnaire pour stocker les tableaux numpy des populations
for i in range(0,5):			# boucle sur les 5 populations
	a = i+1				# ajout de 1 à la valeur de i pour correspondre à la bonne valeur de population (car i=1 correspond à la pop 1) 
	array_name = f"pop_{a}" 	# définition du nom du tableau contenant les haplotypes de la population i 
	pop_arrays_dic[array_name] = np.array([])	# création d'un nouveau tableau numpy pour la population i
	pop_arrays_dic[array_name] = pop_arrays_dic[array_name].astype(int)		# conversion du tableau numpy en entiers
	pop_arrays_dic[array_name] = np.append(pop_arrays_dic[array_name], ts_forw_simpl.samples(population = a)) # remplir le tabelau de la pop i avec les IDs des haplotypes appartenant à la population i 


### Obtention des segments IBD ###

## BETWEEN ## 
# obtenir les segements IBD depuis le tree sequence simplifié en mode BETWEEN 
segments_between = ts_forw_simpl.ibd_segments(between=[pop_arrays_dic['pop_1'], pop_arrays_dic['pop_2'], pop_arrays_dic['pop_3'], pop_arrays_dic['pop_4'], pop_arrays_dic['pop_5']], store_segments = True, store_pairs=True)

# obtenir les valeurs des coordonnées (gauche et droite), le noeud ainsi que l'âge de l'ancêtre commun de chacun des segments IBD between
Pair = []
Left = np.array([])
Right = np.array([])
Node = np.array([])
Time = np.array([])

for k, v in segments_between.items():
    kLeft = v.left 
    kRight = v.right
    kNode = v.node
    order = kLeft.argsort()
    l = kLeft[order]
    n = kNode[order]
    r = kRight[order]
    i = np.where(np.convolve(n + (n[0] == 0), [-1, 1]) != 0)[0][:-1]
    kPair = [k, ]*len(i)
    kLeft = l[i]
    kRight = np.append(r[i[1:]-1], r[len(r)-1])
    kNode = n[i]   
    Pair = Pair + kPair
    Left = np.append(Left, kLeft)
    Right = np.append(Right, kRight)
    Node = np.append(Node, kNode)
    Time = np.append(Time, ts_forw_simpl.tables.nodes.time[kNode]) 
Length = np.subtract(Right,Left)
dfIBD_between = pd.DataFrame(columns = ['Pair', 'Left', 'Right', 'Length', 'Node', 'Time'])
dfIBD_between['Pair'] = Pair
dfIBD_between['Left'] = Left
dfIBD_between['Right'] = Right
dfIBD_between['Length'] = Length
dfIBD_between['Node'] = Node
dfIBD_between['Time'] = Time

Num_IBD_between = len(dfIBD_between.Pair)
print(f"Recorded {Num_IBD_between} IBD segments")

dfIBD_between.to_csv(f'out_IBD_between_ID{valeur_ID}_m{m}_t{valeur_t}_simpl.csv')

## WITHIN ## 
# obtenir les segements IBD depuis le tree sequence simplifié en mode WITHIN pour la population 1
segments_within_1 = ts_forw_simpl.ibd_segments(within=pop_arrays_dic['pop_1'], store_segments = True, store_pairs=True)

# obtenir les valeurs des coordonnées (gauche et droite), le noeud ainsi que l'âge de l'ancêtre commun de chacun des segments IBD within
Pair = []
Left = np.array([])
Right = np.array([])
Node = np.array([])
Time = np.array([])

for k, v in segments_within_1.items():
    kLeft = v.left 
    kRight = v.right
    kNode = v.node
    order = kLeft.argsort()
    l = kLeft[order]
    n = kNode[order]
    r = kRight[order]
    i = np.where(np.convolve(n + (n[0] == 0), [-1, 1]) != 0)[0][:-1]
    kPair = [k, ]*len(i)
    kLeft = l[i]
    kRight = np.append(r[i[1:]-1], r[len(r)-1])
    kNode = n[i]   
    Pair = Pair + kPair
    Left = np.append(Left, kLeft)
    Right = np.append(Right, kRight)
    Node = np.append(Node, kNode)
    Time = np.append(Time, ts_forw_simpl.tables.nodes.time[kNode])
Length = np.subtract(Right,Left)
dfIBD_within_1 = pd.DataFrame(columns = ['Pair', 'Left', 'Right', 'Length', 'Node', 'Time'])
dfIBD_within_1['Pair'] = Pair
dfIBD_within_1['Left'] = Left
dfIBD_within_1['Right'] = Right
dfIBD_within_1['Length'] = Length
dfIBD_within_1['Node'] = Node
dfIBD_within_1['Time'] = Time

Num_IBD_within_1 = len(dfIBD_within_1.Pair)
print(f"Recorded {Num_IBD_within_1} IBD segments")

dfIBD_within_1.to_csv(f'out_IBD_within_1_ID{valeur_ID}_m{m}_t{valeur_t}_simpl.csv')

# obtenir les segements IBD depuis le tree sequence simplifié en mode WITHIN pour la population 2
segments_within_2 = ts_forw_simpl.ibd_segments(within=pop_arrays_dic['pop_2'], store_segments = True, store_pairs=True) 

# obtenir les valeurs des coordonnées (gauche et droite), le noeud ainsi que l'âge de l'ancêtre commun de chacun des segments IBD within
Pair = []
Left = np.array([])
Right = np.array([])
Node = np.array([])
Time = np.array([])

for k, v in segments_within_2.items():
    kLeft = v.left 
    kRight = v.right
    kNode = v.node
    order = kLeft.argsort()
    l = kLeft[order]
    n = kNode[order]
    r = kRight[order]
    i = np.where(np.convolve(n + (n[0] == 0), [-1, 1]) != 0)[0][:-1]
    kPair = [k, ]*len(i)
    kLeft = l[i]
    kRight = np.append(r[i[1:]-1], r[len(r)-1])
    kNode = n[i]   
    Pair = Pair + kPair
    Left = np.append(Left, kLeft)
    Right = np.append(Right, kRight)
    Node = np.append(Node, kNode)
    Time = np.append(Time, ts_forw_simpl.tables.nodes.time[kNode])
Length = np.subtract(Right,Left)
dfIBD_within_2 = pd.DataFrame(columns = ['Pair', 'Left', 'Right', 'Length', 'Node', 'Time'])
dfIBD_within_2['Pair'] = Pair
dfIBD_within_2['Left'] = Left
dfIBD_within_2['Right'] = Right
dfIBD_within_2['Length'] = Length
dfIBD_within_2['Node'] = Node
dfIBD_within_2['Time'] = Time

Num_IBD_within_2 = len(dfIBD_within_2.Pair)
print(f"Recorded {Num_IBD_within_2} IBD segments")

dfIBD_within_2.to_csv(f'out_IBD_within_2_ID{valeur_ID}_m{m}_t{valeur_t}_simpl.csv')

# obtenir les segements IBD depuis le tree sequence simplifié en mode WITHIN pour la population 3
segments_within_3 = ts_forw_simpl.ibd_segments(within=pop_arrays_dic['pop_3'], store_segments = True, store_pairs=True) # population 3

# obtenir les valeurs des coordonnées (gauche et droite), le noeud ainsi que l'âge de l'ancêtre commun de chacun des segments IBD within
Pair = []
Left = np.array([])
Right = np.array([])
Node = np.array([])
Time = np.array([])

for k, v in segments_within_3.items():
    kLeft = v.left 
    kRight = v.right
    kNode = v.node
    order = kLeft.argsort()
    l = kLeft[order]
    n = kNode[order]
    r = kRight[order]
    i = np.where(np.convolve(n + (n[0] == 0), [-1, 1]) != 0)[0][:-1]
    kPair = [k, ]*len(i)
    kLeft = l[i]
    kRight = np.append(r[i[1:]-1], r[len(r)-1])
    kNode = n[i]   
    Pair = Pair + kPair
    Left = np.append(Left, kLeft)
    Right = np.append(Right, kRight)
    Node = np.append(Node, kNode)
    Time = np.append(Time, ts_forw_simpl.tables.nodes.time[kNode])
Length = np.subtract(Right,Left)
dfIBD_within_3 = pd.DataFrame(columns = ['Pair', 'Left', 'Right', 'Length', 'Node', 'Time'])
dfIBD_within_3['Pair'] = Pair
dfIBD_within_3['Left'] = Left
dfIBD_within_3['Right'] = Right
dfIBD_within_3['Length'] = Length
dfIBD_within_3['Node'] = Node
dfIBD_within_3['Time'] = Time

Num_IBD_within_3 = len(dfIBD_within_3.Pair)
print(f"Recorded {Num_IBD_within_3} IBD segments")

dfIBD_within_3.to_csv(f'out_IBD_within_3_ID{valeur_ID}_m{m}_t{valeur_t}_simpl.csv')

# obtenir les segements IBD depuis le tree sequence simplifié en mode WITHIN pour la population 4
segments_within_4 = ts_forw_simpl.ibd_segments(within=pop_arrays_dic['pop_4'], store_segments = True, store_pairs=True) # population 4

# obtenir les valeurs des coordonnées (gauche et droite), le noeud ainsi que l'âge de l'ancêtre commun de chacun des segments IBD within
Pair = []
Left = np.array([])
Right = np.array([])
Node = np.array([])
Time = np.array([])

for k, v in segments_within_4.items():
    kLeft = v.left 
    kRight = v.right
    kNode = v.node
    order = kLeft.argsort()
    l = kLeft[order]
    n = kNode[order]
    r = kRight[order]
    i = np.where(np.convolve(n + (n[0] == 0), [-1, 1]) != 0)[0][:-1]
    kPair = [k, ]*len(i)
    kLeft = l[i]
    kRight = np.append(r[i[1:]-1], r[len(r)-1])
    kNode = n[i]   
    Pair = Pair + kPair
    Left = np.append(Left, kLeft)
    Right = np.append(Right, kRight)
    Node = np.append(Node, kNode)
    Time = np.append(Time, ts_forw_simpl.tables.nodes.time[kNode])
Length = np.subtract(Right,Left)
dfIBD_within_4 = pd.DataFrame(columns = ['Pair', 'Left', 'Right', 'Length', 'Node', 'Time'])
dfIBD_within_4['Pair'] = Pair
dfIBD_within_4['Left'] = Left
dfIBD_within_4['Right'] = Right
dfIBD_within_4['Length'] = Length
dfIBD_within_4['Node'] = Node
dfIBD_within_4['Time'] = Time

Num_IBD_within_4 = len(dfIBD_within_4.Pair)
print(f"Recorded {Num_IBD_within_4} IBD segments")

dfIBD_within_4.to_csv(f'out_IBD_within_4_ID{valeur_ID}_m{m}_t{valeur_t}_simpl.csv')

# obtenir les segements IBD depuis le tree sequence simplifié en mode WITHIN pour la population 5
segments_within_5 = ts_forw_simpl.ibd_segments(within=pop_arrays_dic['pop_5'], store_segments = True, store_pairs=True) # population 5

# obtenir les valeurs des coordonnées (gauche et droite), le noeud ainsi que l'âge de l'ancêtre commun de chacun des segments IBD within
Pair = []
Left = np.array([])
Right = np.array([])
Node = np.array([])
Time = np.array([])

for k, v in segments_within_5.items():
    kLeft = v.left 
    kRight = v.right
    kNode = v.node
    order = kLeft.argsort()
    l = kLeft[order]
    n = kNode[order]
    r = kRight[order]
    i = np.where(np.convolve(n + (n[0] == 0), [-1, 1]) != 0)[0][:-1]
    kPair = [k, ]*len(i)
    kLeft = l[i]
    kRight = np.append(r[i[1:]-1], r[len(r)-1])
    kNode = n[i]   
    Pair = Pair + kPair
    Left = np.append(Left, kLeft)
    Right = np.append(Right, kRight)
    Node = np.append(Node, kNode)
    Time = np.append(Time, ts_forw_simpl.tables.nodes.time[kNode])
Length = np.subtract(Right,Left)
dfIBD_within_5 = pd.DataFrame(columns = ['Pair', 'Left', 'Right', 'Length', 'Node', 'Time'])
dfIBD_within_5['Pair'] = Pair
dfIBD_within_5['Left'] = Left
dfIBD_within_5['Right'] = Right
dfIBD_within_5['Length'] = Length
dfIBD_within_5['Node'] = Node
dfIBD_within_5['Time'] = Time

Num_IBD_within_5 = len(dfIBD_within_5.Pair)
print(f"Recorded {Num_IBD_within_5} IBD segments")

dfIBD_within_5.to_csv(f'out_IBD_within_5_ID{valeur_ID}_m{m}_t{valeur_t}_simpl.csv')

# fusionner en un seul tableau tous les sous-tableaux des différentes populations within
out_IBD_within = pd.concat([dfIBD_within_1, dfIBD_within_2, dfIBD_within_3, dfIBD_within_4, dfIBD_within_5], ignore_index=True)
out_IBD_within.to_csv(f'out_IBD_within_ID{valeur_ID}_m{m}_t{valeur_t}_simpl.csv', index = False)


### Visualisation des segments IBD ###

# création d'un histogramme du nombre de segements IBD en fonction de leur longueur
plt.figure()
	plt.hist(out_IBD_within["Length"], range = (0, 85000000), bins=60, color = '#66CC99', edgecolor = 'white', label = "Within", alpha = 0.7)
	plt.hist(out_IBD_between["Length"], range = (0, 85000000), bins=60, color = '#009999', edgecolor = 'white', label = "Between")
	plt.xlabel('Longueur')
	plt.ylabel('Nombre')
	plt.title(f'Distribution en longueur des segments IBD simulés', fontsize = 16, loc = 'center')
	plt.suptitle(f'ID{ID} m{m} t{t}', fontsize = 12)
	plt.yscale('log')
	plt.legend(loc = 'upper right', frameon = False)
	plt.savefig(f'IBD_length_vs__ID{ID}_m{m}_t{t}_simpl.png')
	print(f"Generated IBD segments length distribution ID{ID}_m{m}_t{t} plot")
```

La simulation dans SLiM intégrée dans le script précédent lance le script suivant dans SLiM :

```{python, eval=FALSE, echo = TRUE}
// initialisation des paramètres 
initialize(){
	initializeTreeSeq();
	// définition des constantes
	defineConstant("m", 0.0);         // taux de migration nul
	defineConstant("Ne", 200);        // taille des dèmes
	defineConstant("R", 1e-8);        // taux de recombinaison moyen 
	
	// définition des mutations 
	initializeMutationRate(0);
	initializeMutationType("m1", 0.5, "f", 0.0);   // mutation neutre et co-dominante
	
	//  définition du génome
	initializeGenomicElementType("g1", m1, 1.0);
	initializeGenomicElement(g1, 0, 84999999);
	initializeRecombinationRate(R);
}

// importation du tree-sequence obtenu par coalescence 
1 early(){
	sim.readFromPopulationFile("/DATA/sdb1/OCEANE/coal.trees");
	defineConstant("simID", getSeed());
}

// définition des taux de migration 
1: late(){
	subpopCount = 5;
	for (i in 1:subpopCount)
		for (j in 1:subpopCount)
			if (i !=j)
				sim.subpopulations[i-1].setMigrationRates(j, m);
}

// calcul à chaque génération du FST entre toutes les paires de pop 
1:10000 late(){
	subpopCount = 5;
	for (i in 1:subpopCount)
		for (j in 1:subpopCount)
			if (i != j)
				FST = calcFST(sim.subpopulations[i-1].genomes, sim.subpopulations[j-1].genomes);
	mean_FST = mean(FST);
	// lorsque la valeur de FST dépasse 0.25 (FST actuel dans la métapop) => arrêter la simu 
	if (mean_FST >= 0.25){
		defineConstant("t", sim.cycle);
		cat("sim_ID == " + simID + " time == " + t);
		sim.simulationFinished();
		// enregistrer le tree sequence avec comme nom : "forw + simID + t + .trees"
		sim.treeSeqOutput("/DATA/sdb1/OCEANE/" + "forw_" + simID + "ID_t" + t + ".trees");
	}
}
```

------------------------------------------------------------------------

### Calcul de la statistique de différenciation

Afin d’avoir une statistique représentative de la différence entre la distribution des segments IBD en ‘within’ et en ‘between’, une statistique FST a été calculée pour chaque simulation sur l’intervalle de temps récent allant de 0 à 300 générations (subdivisé en pas de 5).

```{python, eval=FALSE, echo = TRUE}
import pandas as pd
import matplotlib
from matplotlib import pyplot as plt

# récupération du tableau contenant les valeurs de l'ID, de t et de m de la simulation 
Params = pd.read_csv("/DATA/sdb1/OCEANE/ID_Migration_Time_2.csv")

# création d'une boucle pour réaliser le calcul de FST sur chaque simulation 
i = 0
for i in range(0,26):
	# récupération de la valeur de l'ID de la simulation
	Valeurs_ID = Params['ID']
	ID = Valeurs_ID[i]
	# récupération de la valeur de m de la simulation
	Valeurs_m = Params['Migration']
	m = Valeurs_m[i]
	m = round(m, 6)
	# récupération de la valeur de t de la simulation
	Valeurs_t = Params['Time']
	t = Valeurs_t[i]

	# définition des intervalles de 5 compris entre 0 et 300 
	intervals = range(0, 301, 5)


	############################################################
	################### CALCUL PROBA BETWEEN ###################
	############################################################

	# chargement du tableau between 
	IBD_between = pd.read_csv(f'/DATA/sdb1/OCEANE/OUT_BETWEEN_2/out_IBD_between_ID{ID}_m{m}_t{t}_simpl.csv')

	# création d'une nouvelle colonne 'Interval' indiquant l'intervalle correspondant à la valeur de Time 
	IBD_between['Interval'] = pd.cut(IBD_between['Time'], bins = intervals)

	# calcul de la somme des longueurs pour chacunes des paires et chacune des intervalles 
	sum_length_between = IBD_between.groupby(['Pair', 'Interval'])['Length'].sum().reset_index()

	# remplacer les valeurs NaN par des 0.0
	sum_length_between['Length'] = sum_length_between['Length'].fillna(0.0)

	# calcul de la probabilité between à chacune des intervalles et pour chacune des paires 
	sum_length_between['Proba_between'] = sum_length_between['Length'] / 85000000

	# faire la moyenne de l'ensemble des paires sur chaque intervalle 
	mean_proba_between = sum_length_between.groupby('Interval')['Proba_between'].mean()

	mean_proba_between = mean_proba_between.fillna(0.0)


	############################################################
	################### CALCUL PROBA WITHIN ####################
	############################################################

	# chargement du tableau within 
	IBD_within = pd.read_csv(f'/DATA/sdb1/OCEANE/OUT_WITHIN_2/out_IBD_within_ID{ID}_m{m}_t{t}_simpl.csv')

	# création d'une nouvelle colonne 'Interval' indiquant l'intervalle correspondant à la valeur de Time 
	IBD_within['Interval'] = pd.cut(IBD_within['Time'], bins = intervals)

	# calcul de la somme des longueurs pour chacunes des paires et chacune des intervalles 
	sum_length_within = IBD_within.groupby(['Pair', 'Interval'])['Length'].sum().reset_index()

	# remplacer les valeurs NaN par des 0.0
	sum_length_within['Length'] = sum_length_within['Length'].fillna(0.0)

	# calcul de la probabilité within à chacune des intervalles et pour chacune des paires 
	sum_length_within['Proba_within'] = sum_length_within['Length'] / 85000000

	# faire la moyenne de l'ensemble des paires sur chaque intervalle 
	mean_proba_within = sum_length_within.groupby('Interval')['Proba_within'].mean()

	mean_proba_within = mean_proba_within.fillna(0.0)


	############################################################
	######################## CALCUL FST ########################
	############################################################

	# calculer la valeur de FST = (pw - pb) / (1-pb)
	mean_proba = pd.concat([mean_proba_between, mean_proba_within], axis=1)
	mean_proba['FST'] = (mean_proba['Proba_within'] - mean_proba['Proba_between']) / (1 - mean_proba['Proba_between'])

	# exporter le tableau au format .csv pour créer le plot sur Rstudio
	Time = intervals[1:61]		# prendre les valeurs des intervalles 
	mean_proba['Time'] = Time 	# rajouter une colonne au tableau contenant la valeur supérieure de chaque intervalle 
	mean_proba['Migration'] = m	# rajouter une colonne au tableau contenant la valeur de migration de la simulation 
	mean_proba.to_csv(f'Proba_FST_ID{ID}_m{m}_t{t}.csv')

	print(f"Generated Proba_FST_ID{ID}_m{m}_t{t}.csv")
```

------------------------------------------------------------------------

### Distribution en longueur des segments IBD (between vs. within)

Pour chaque simulation, un histogramme de la distribution en longueur des segments IBD comparant les segments between et les segments within a été obtenu. Les 2 extrêmes de valeur de migration (m) sont présentées ici :

::: {style="display:flex; justify-content:center;"}
::: {style="flex:1; padding-left:10px;"}
![](/Users/oceane/Desktop/Stage/SIMULATIONS/IBD_length_vs_ID1795848085227_m0.00037_t7_simpl.png){width="600px" height="400px"}
:::

::: {style="flex:1; padding-right:8px;"}
![](/Users/oceane/Desktop/Stage/SIMULATIONS/IBD_length_vs_ID1879781714091_m0.012_t213_simpl.png){width="600px" height="400px"}
:::
:::

L'ensemble des graphiques pour chacune des valeurs de m sont accesibles en téléchargement via un clic [ici](https://cloud.isem-evolution.fr/nextcloud/index.php/s/MpDxmfcY4SDXFfr?path=%2FCcory%2FSIMULATIONS%2FPLOT%20LENGTH%20IBD).

------------------------------------------------------------------------

### Espace des possibles du couple migration-temps

Les valeurs du couple m-t pour toutes les simulations ont été représentées sous la forme d'un graphique illustrant l'espace des possibles du couple m-t :

::: {style="display:flex; justify-content:center;"}
```{r}
data <- read.csv("/Users/oceane/Desktop/Stage/APPENDIX/SIMULATIONS/ID_Migration_Time.csv")

plot_m_t = ggplot(data, aes(x = Migration, y = Time)) + theme_classic() + 
  geom_point(colour = "red") +
  labs(title = "", y = "Nombre de générations (t)", x = "Taux de migration (m)") + 
  scale_x_log10(limits = c(0.0001, 0.1), breaks = c(0.0001, 0.001, 0.01, 0.1)) + 
  theme(plot.title = element_text(size = 32, face = "bold", hjust = 0.5),
        axis.title = element_text(size = 25, hjust = 0.5),
        axis.text = element_text(size = 15))
ggplotly(plot_m_t, height = 500, width = 1050)
```
:::

------------------------------------------------------------------------

### Distribution du FST en fonction du temps

Pour chaque simulation, un graphique de la distribution des valeurs de FST en fonction du temps a été obtenu. Les 2 extrêmes de valeur de migration (m) sont présentées ici :

::: {style="display:flex; justify-content:center;"}
::: {style="flex:1; padding-left:10px;"}
```{r}

Proba_FST_ID1795848085227_m0.00037_t7 = read.csv("/Users/oceane/Desktop/Stage/APPENDIX/SIMULATIONS/FST_Time_0.00037.csv")

plot_m_0.00037 = ggplot(Proba_FST_ID1795848085227_m0.00037_t7, aes(y = FST, x = Time)) + theme_classic() + 
  geom_point(colour = "#2A788EFF") +
  labs(title = "", subtitle = "ID1795848085227  m0.00037  t7",  y = "FST", x = "Temps") + 
  theme(plot.title = element_text(size = 32, face = "bold", hjust = 0.5),
        plot.subtitle = element_text(size = 20, hjust = 0.5),
        axis.title = element_text(size = 25, hjust = 0.5),
        axis.text = element_text(size = 15))

ggplotly(plot_m_0.00037, height = 500)
```
:::

::: {style="flex:1; padding-right:10px;"}
```{r}

Proba_FST_ID1879781714091_m0.012_t213 = read.csv("/Users/oceane/Desktop/Stage/APPENDIX/SIMULATIONS/FST_Time_0.012.csv")

plot_m_0.012 = ggplot(Proba_FST_ID1879781714091_m0.012_t213, aes(y = FST, x = Time)) + theme_classic() + 
  geom_point(colour = "#2A788EFF") +
  labs(title = "", subtitle = "ID1879781714091  m0.012  t213",  y = "", x = "Temps") + 
  theme(plot.title = element_text(size = 32, face = "bold", hjust = 0.5),
        plot.subtitle = element_text(size = 20, hjust = 0.5),
        axis.title = element_text(size = 25, hjust = 0.5),
        axis.text = element_text(size = 15))

ggplotly(plot_m_0.012, height = 500)
```
:::
:::

------------------------------------------------------------------------

### Points de cassure de la distribution du FST

Les 2 représentations précédentes ont ensuite été traités pour fitter sur les distributions 3 modèles linéaires afin de trouver les 2 points de cassures entre ces modèles. Le premier point de cassure représentant le moment où la valeur de FST subit un décrochage, ce qui s'apparente alors au moment où la migration a cessé dans le passé. Les 2 extrêmes sont une nouvelle fois présentées ici :

::: {style="display:flex; justify-content:center;"}
::: {style="flex:1; padding-left:10px;"}
![](/Users/oceane/Desktop/Stage/APPENDIX/SIMULATIONS/seg_0.00037.png){width="600px" height="400px"}
:::

::: {style="flex:1; padding-right:8px;"}
![](/Users/oceane/Desktop/Stage/APPENDIX/SIMULATIONS/seg_0.012.png){width="600px" height="400px"}
:::
:::

------------------------------------------------------------------------

### Distribution du temps estimé en fonction du temps observé

Ainsi, le temps estimé à partir des valeurs de FST (premier point de cassure) a pu être représenté en fonction du temps observé dans les simulations :

::: {style="justify-content:center;"}
```{r}
data_time_all = read.csv("/Users/oceane/Desktop/Stage/APPENDIX/SIMULATIONS/data_time_all.csv")

plot_time_all = ggplot(data_time_all, aes(y = Est., x = Time)) + theme_classic() + 
  geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed", size= 0.8, alpha = 0.8) + 
  geom_smooth(method = "lm", colour = "#39568CFF", fill = "#39568CFF", alpha = 0.2) + 
  geom_point(colour = "#39568CFF", size = 2.5) +
  geom_errorbar(aes(ymin = Est. - St.Err, ymax = Est. + St.Err), width = 0.1, colour = "#39568CFF") + 
  annotate("text", x = 190, y = -20, label = "r = 0.801 (p < 0.001)", col = "#39568CFF", size = 7) +
  labs(title = "",  y = "Temps estimé (t.est)\n ", x = "Temps (t)") + 
  theme(axis.title = element_text(size = 25, hjust = 0.5),
        axis.text = element_text(size = 20))
ggplotly(plot_time_all, height = 500, width = 1050)
```
:::

------------------------------------------------------------------------

### Distribution en "3D"

Et pour finir, un graphique en "3D" a permis de représenter les valeurs de FST en fonction des 2 paramètres d'intérêts, à savoir le temps et la migration.

::: {style="text-align:center;"}
```{r}
Proba_FST_all = read.csv("/Users/oceane/Desktop/Stage/APPENDIX/SIMULATIONS/Proba_FST_all.csv")
Time_est = read.csv("/Users/oceane/Desktop/Stage/APPENDIX/SIMULATIONS/Time_est.csv")
x_breaks = c(0.0001, 0.001, 0.01, 0.1)

plot_3D = ggplot(Proba_FST_all, aes(y = Time, x = Migration, color = FST)) + theme_classic() + 
    geom_point(shape = 15, size = 3) + 
    geom_point(data = Time_est, aes(x = Migration, y = Est.), shape = 16, size = 3, col = "red") + 
    labs(title = "", y = "Temps (t)\n", x = "Migration (m)", color = "FST") + 
    scale_color_viridis(option = "D", direction = -1) +
    scale_x_log10(breaks = x_breaks) + 
    theme(axis.title = element_text(size = 25, hjust = 0.5),
          legend.title = element_text(size = 20, hjust = 0.5), 
          legend.text = element_text(size = 15), 
          axis.text = element_text(size = 20))
ggplotly(plot_3D, height = 500, width = 1050)
```
:::

------------------------------------------------------------------------