c-figure-03.Rmd

---
title: "Panels of Figure 3: Comparative genomics of the 21 NR-genomes of OsHV-1 isolated from three farming areas"
author: "Delmotte jean"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output:
  html_document:
    toc: TRUE
    code_folding: "hide"
    theme: united
    highlight: tango
    number_sections: true
editor_options: 
  chunk_output_type: console
---

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

# Fonction to install / load package if it's not here
ipak <- function(pkg){
  new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
  if (length(new.pkg)) 
    install.packages(new.pkg, dependencies = TRUE)
  sapply(pkg, require, character.only = TRUE)
}

ipak(unique(
  c("data.table", "tidyverse", "ape", "RColorBrewer", "pheatmap", "limma", "glue", "kableExtra", "plotly", "tidytree", "broom",
    "hrbrthemes", "viridis", "ggsci", # palette
    "cowplot", "scales", "maptools", "gggenes") # plot
  ) )
```

```{r Base_Path}
base_path <- "~/Documents/OshV-1-molepidemio" # Base location for the folder
```

# Preparation and import of data

Abbreviations : 
- GC: comparative genomics
- VC: variant calling 

```{r Color_analysis, echo=FALSE, warning = FALSE, message = FALSE}
## Convention
# Color Location Brest:blue | LT:Red | Thau:Green
# show_col(pal_aaas("default")(10))
Brest_color = pal_aaas("default")(10)[1]
LT_color = pal_aaas("default")(10)[2]
Thau_color = pal_aaas("default")(10)[3]

# Color Polymorphisms 
# show_col(pal_aaas("default")(10))
SNP = pal_aaas("default")(10)[2]
INDEL = pal_jco("default")(10)[3]

# Order
plot_order_samples <- tibble(ID_EXPERIMENT = c("Brest_2018_NSI_broyage_ind10_noPCR", 
                                               "Brest_2018_NSI_broyage_ind2_noPCR", 
                                               "Brest_2018_NSI_broyage_ind4_noPCR", 
                                               "Brest_2018_NSI_broyage_ind6_noPCR", 
                                               "Brest_2018_NSI_broyage_ind9_noPCR", 
                                               "LT_2018_NSI_broyage_ind10_noPCR", 
                                               "LT_2018_NSI_broyage_ind1_noPCR", 
                                               "LT_2018_NSI_broyage_ind3_noPCR", 
                                               "LT_2018_NSI_broyage_ind4_noPCR", 
                                               "LT_2018_NSI_broyage_ind7_noPCR", 
                                               "LT_2018_NSI_broyage_ind8_noPCR", 
                                               "LT_2018_NSI_broyage_ind9_noPCR", 
                                               "Thau_2018_NSI_broyage_ind10_noPCR", 
                                               "Thau_2018_NSI_broyage_ind1_noPCR", 
                                               "Thau_2018_NSI_broyage_ind3_noPCR", 
                                               "Thau_2018_NSI_broyage_ind4_noPCR",
                                               "Thau_2018_NSI_broyage_ind5_noPCR", 
                                               "Thau_2018_NSI_broyage_ind6_noPCR", 
                                               "Thau_2018_NSI_broyage_ind7_noPCR", 
                                               "Thau_2018_NSI_broyage_ind8_noPCR", 
                                               "Thau_2018_NSI_broyage_ind9_noPCR") ) %>%
  tidyr::separate(ID_EXPERIMENT, c("LOCATION","YEARS","FAM","TYPE","CONDITION","WGS_PCR"), remove = FALSE) %>%
  tidyr::unite(LOCATION, CONDITION, col="Samples",sep="_",remove = FALSE) %>%
  dplyr::left_join(., tibble(LOCATION = c("Brest","LT","Thau"),
       ORDER= c(1,2,3))) %>%
  dplyr::mutate(ORDER_samples = stringr::str_replace_all(CONDITION,"ind",""),
                ORDER_samples = as.numeric(ORDER_samples)) %>%
  tidyr::unite(LOCATION, YEARS, col="VAR",sep="_",remove = FALSE) %>%
  tidyr::unite(LOCATION, CONDITION, col="Samples",sep="_",remove = FALSE) %>%
  arrange(ORDER, ORDER_samples) %>%
  dplyr::mutate(Samples = stringr::str_replace_all(Samples,"_"," "))
```

```{r TABLE_03}
# importation Comparative genomics
comparative_genomic <- data.table::fread(glue::glue("{base_path}/results/Tables/nucmer_numref.snps"),
                  skip = 3,
         sep = "\t",
         fill=TRUE,
         quote=FALSE)
colnames(comparative_genomic) <-c("POS", "REF", "ALT","POS_ALT", "BUFF","DIST","R", "Q", "FRM","TAGS", "gREF","ID_EXPERIMENT")

comparative_genomic <- comparative_genomic %>%
dplyr::mutate(UniqueID= glue::glue("{POS}_{REF}>{ALT}")) %>%
    as_tibble() %>%
  group_by(UniqueID, ALT, ID_EXPERIMENT) %>%
  mutate(ALT = paste(ALT, collapse = "")) %>%
  data.frame() %>%
  select(-POS_ALT, -UniqueID, -DIST) %>%
  distinct() 
```


# Analysis

```{r}
# Get the genomic position of comparative genomic modification
CG_POS <- comparative_genomic %>%
  select(POS) %>%
  distinct() %>%
  pull()


# In the select, we don't used a regex like `contain()` because we want to be sure of the order
Table_02_binary_matrix <- comparative_genomic %>%
  dplyr::mutate(UniqueID= glue::glue("{POS}_{REF}>{ALT}"),
                AF = 1) %>%
  data.frame() %>%
  dplyr::select(UniqueID,
                ID_EXPERIMENT,
                AF) %>%
  distinct() %>%
  tidyr::pivot_wider(names_from = ID_EXPERIMENT,
                     values_from = AF,
                     values_fill = list(AF = 0)) %>%
  as.data.frame() %>% 
  select(UniqueID,
         NR_genome_Brest_2018_NSI_broyage_ind10_noPCR,
         NR_genome_Brest_2018_NSI_broyage_ind2_noPCR,
         NR_genome_Brest_2018_NSI_broyage_ind4_noPCR,
         NR_genome_Brest_2018_NSI_broyage_ind6_noPCR,
         NR_genome_Brest_2018_NSI_broyage_ind9_noPCR,
         NR_genome_LT_2018_NSI_broyage_ind10_noPCR,
         NR_genome_LT_2018_NSI_broyage_ind1_noPCR,
         NR_genome_LT_2018_NSI_broyage_ind3_noPCR,
         NR_genome_LT_2018_NSI_broyage_ind4_noPCR,
         NR_genome_LT_2018_NSI_broyage_ind7_noPCR,
         NR_genome_LT_2018_NSI_broyage_ind8_noPCR,
         NR_genome_LT_2018_NSI_broyage_ind9_noPCR,
         NR_genome_Thau_2018_NSI_broyage_ind10_noPCR,
         NR_genome_Thau_2018_NSI_broyage_ind1_noPCR,
         NR_genome_Thau_2018_NSI_broyage_ind3_noPCR,
         NR_genome_Thau_2018_NSI_broyage_ind4_noPCR,
         NR_genome_Thau_2018_NSI_broyage_ind5_noPCR,
         NR_genome_Thau_2018_NSI_broyage_ind6_noPCR,
         NR_genome_Thau_2018_NSI_broyage_ind7_noPCR,
         NR_genome_Thau_2018_NSI_broyage_ind8_noPCR,
         NR_genome_Thau_2018_NSI_broyage_ind9_noPCR) %>%
  tibble::column_to_rownames("UniqueID") %>%
  as.matrix()
```

## Figure 3 B Analysis of all genomic variations at the conchyliculture area level

```{r}
Matrix_SET <- comparative_genomic %>%
  dplyr::mutate(ID_EXPERIMENT = stringr::str_replace_all(ID_EXPERIMENT,"NR_genome_","" )) %>%
  tidyr::separate(ID_EXPERIMENT, c("LOCATION","YEARS","FAM","STEP","CONDITION","WGS_PCR"), remove = FALSE) %>%
  tidyr::unite(LOCATION, YEARS, CONDITION, col="Samples",sep="_",remove = FALSE) %>%
  dplyr::mutate(SET = glue::glue("{LOCATION}_{YEARS}"),
                SET = as.character(SET)) %>%
  dplyr::select(ID_EXPERIMENT, SET) %>%
  dplyr::distinct() %>% table()

Table_XXXXXXX<-Table_02_binary_matrix %*% Matrix_SET %>%
  as.data.frame() %>%
  rownames_to_column("UniqueID")

# impact on the genome
Table_XXXXXXX %>%
  tidyr::pivot_longer(-UniqueID) %>%
  dplyr::filter(value >0) %>%
  separate(UniqueID, c("POS","CHANGE"), sep = "_", remove = FALSE) %>%
  mutate(Samples = stringr::str_replace_all(name,"_","")) %>%
  ggplot(aes(x = POS, y=value)) +
    geom_bar(stat = "identity")

Table_XXXXXXX %>%
  tidyr::pivot_longer(-UniqueID) %>%
  dplyr::mutate(value = dplyr::case_when(value > 0 ~ 1,
                                      TRUE ~ 0)) %>%
  dplyr::mutate(name = stringr::str_replace_all(name,"_2018",""))  %>%
  tidyr::pivot_wider(names_from = name,
                     values_from = value,
                     values_fill = list(value = 0)) %>%
  column_to_rownames(var = "UniqueID") %>%
  # filter(Brest == 1 & LT == 1 & Thau == 1)
  as.matrix() %>%
  vennDiagram()
```

mafft \
    --thread 4 \
    --auto \
    oshv_wRef.fna \
    > oshv_wRef.faa

Note: Thau lagoon (Thau) and the LT (Marennes Olérons) were exchanged on the article using adobe illustrator.

## Figure 3 C Analysis of all genomic variations at the individual level

```{r}
#### Upsetplot Fig3B
set_order <- plot_order_samples %>%
  dplyr::mutate(Samples = stringr::str_replace_all(Samples,"Brest","Br"),
               Samples = stringr::str_replace_all(Samples,"LT","LT"),
               Samples = stringr::str_replace_all(Samples,"Thau","Th"),
               Samples = stringr::str_replace_all(Samples,"_"," ") ) %>%
  dplyr::select(Samples) %>% pull()

Figure_03_B_comparative_genomic_nt <-  comparative_genomic %>%
  dplyr::mutate(UniqueID= glue::glue("{POS}_{REF}>{ALT}")) %>%
    as_tibble() %>%
  group_by(UniqueID, ALT, ID_EXPERIMENT) %>%
  mutate(ALT = paste(ALT, collapse = "")) %>%
  data.frame() %>%
  select(-UniqueID) %>%
  distinct() %>%
 dplyr::mutate(ID_EXPERIMENT = stringr::str_replace_all(ID_EXPERIMENT,"NR_genome_",""),
               ID_EXPERIMENT = stringr::str_replace_all(ID_EXPERIMENT,"_NSI_broyage_"," "),
               ID_EXPERIMENT = stringr::str_replace_all(ID_EXPERIMENT,"_noPCR",""),
               ID_EXPERIMENT = stringr::str_replace_all(ID_EXPERIMENT,"_2018",""),
               ID_EXPERIMENT = stringr::str_replace_all(ID_EXPERIMENT,"Brest","Br"),
               ID_EXPERIMENT = stringr::str_replace_all(ID_EXPERIMENT,"LT","LT"),
               ID_EXPERIMENT = stringr::str_replace_all(ID_EXPERIMENT,"Thau","Th"),
   UniqueID= glue::glue("{POS}_{REF}>{ALT}"),
               AF = 1) %>%
 dplyr::select(UniqueID,
               ID_EXPERIMENT,
               AF) %>%
 tidyr::pivot_wider(names_from = ID_EXPERIMENT,
                    values_from = AF,
                    values_fill = list(AF = 0)) %>%
 as.data.frame() %>%
  tibble::column_to_rownames("UniqueID")  %>%
  as.matrix() %>%
  as.data.frame() %>%
  UpSetR::upset(nsets = 30,
               nintersects = NA,
               query.legend = "top",
               set_size.show = TRUE,
               sets = rev(set_order),
               keep.order = TRUE)
```

```{r Figure_03_B_comparative_genomic_nt_export, message=FALSE}
ggsave("Figure_03_B_comparative_genomic_nt.jpeg",
       plot = Figure_03_B_comparative_genomic_nt,
       device = "jpeg",
       path = glue::glue("{base_path}/results/Figures/jpeg/Figure_3"),
       width = 180, units = "mm", dpi = 300)

ggsave("Figure_03_B_comparative_genomic_nt.tiff",
       plot = Figure_03_B_comparative_genomic_nt,
       device = "tiff",
       path = glue::glue("{base_path}/results/Figures/tiff/Figure_3"),
       width = 180, units = "mm", dpi = 300)

Figure_03_B_comparative_genomic_nt
dev.copy2pdf(file=glue::glue("{base_path}/results/Figures/eps_svg/Figure_3/Figure_03_B_comparative_genomic_nt.pdf"),out.type="cairo", width=18)
writeLines(glue::glue("pdftops -eps  {base_path}/results/Figures/eps_svg/Figure_3/Figure_03_B_comparative_genomic_nt.pdf"))
rm(Figure_03_B_comparative_genomic_nt)
```

Note: The figure 3B has been reworked on adobe illustrator

To understand how to dichotomize the sets we have based ourselves on the following table:

| Name         | Explanation of genomic changes or polymorphisms                                     | Sets                  |
|--------------|-------------------------------------------------------------------------------------|-----------------------|
| singletons   | Variation find in only one samples                                                  | n=1                   |
| Shared intra | Variation shared by at least 2 samples and not between all samples of the same area | 1 < n < N_area        |
| Common intra | Variation find in all samples of the same area                                      | n=N_area              |
| Shared inter | Variation shared by at least 2 samples and not between all samples between area     | 1 < n < N_experiment  |
| Common inter | Variation find in samples of all area                                               | n=N_experiment        |
| Global       | Variation which are find in all samples in all area                                 | n=N_area+N_experiment |


```{r}
N_Brest=5
N_LT=7
N_Thau=9

GC_table <- Table_XXXXXXX %>%
dplyr::mutate(Sets =
                  case_when(
                    Brest_2018  == 1 & Brest_2018 < N_Brest & Thau_2018 == 0 & LT_2018 == 0 ~ "Br-sing", 
                    Brest_2018  == 0 & Thau_2018 == 0 & LT_2018 == 1  & LT_2018 < N_LT~ "Lt-sing",
                    Brest_2018  == 0 & Thau_2018 == 1 & Thau_2018 < N_Thau & LT_2018 == 0 ~ "Th-sing",
                    Brest_2018  > 1 & Brest_2018 < N_Brest & Thau_2018 == 0 & LT_2018 == 0 ~ "Br-shared_intra", 
                    Brest_2018  == 0 & Thau_2018 > 1 & Thau_2018 < N_Thau & LT_2018 == 0 ~ "Th-shared_intra",
                    Brest_2018  >= 1  & Thau_2018 >= 1 & LT_2018 == 0 ~ "Br Th-shared_inter",
                    Brest_2018  >= 1  & Thau_2018 == 0 & LT_2018 >= 1 ~ "Br Lt-shared_inter",
                    Brest_2018  == 0  & Thau_2018 >= 1 & LT_2018 >= 1 ~ "Lt Th-shared_inter",
                    Brest_2018  == 0 & Thau_2018 == N_Thau & LT_2018 == 0 ~ "Th-common_intra",
                    Brest_2018  == N_Brest & Thau_2018 == 0 & LT_2018 == 0 ~ "Br-common_intra",
                    Brest_2018  == 0 & Thau_2018 == 0 & LT_2018 == N_LT ~ "Lt-common_intra",
                    Brest_2018  == 0 & Thau_2018 == 0 & LT_2018 > 1  & LT_2018 < N_LT~ "Lt-shared_intra",
                    Brest_2018  >= 1  & Thau_2018 >= 1 &  LT_2018 >=1 ~ "Common_inter",
                    Brest_2018  == N_Brest  & Thau_2018 == N_Thau &  LT_2018 == N_LT ~ "Global"
                    )
                ) %>%
  dplyr::mutate(DATASET = "GC")

Global_GC_table <- GC_table %>%
  left_join(., Table_02_binary_matrix %>%
              as.data.frame() %>%
              rownames_to_column(var = "UniqueID")) %>%
  as_tibble()


GC_table %>%
  left_join(.,Table_02_binary_matrix %>%
  as.data.frame() %>%
  rownames_to_column("UniqueID")) #%>%
  # data.table::fwrite(glue::glue("{base_path}/results/Tables/Table02_comparative_genomic_matrix.csv"),
  # sep = "\t",
  # quote=FALSE,
  # row.names = FALSE,
  # col.names = TRUE)

# Global overview
comparative_genomic %>% 
  group_by(ID_EXPERIMENT) %>%
  summarise(N= n()) %>%
  dplyr::mutate(ID_EXPERIMENT = stringr::str_replace_all(ID_EXPERIMENT,"NR_genome_","" )) %>%
  tidyr::separate(ID_EXPERIMENT, c("LOCATION","YEARS","FAM","STEP","CONDITION","WGS_PCR"), remove = FALSE) %>%
  tidyr::unite(LOCATION, YEARS, CONDITION, col="Samples",sep="_",remove = FALSE) %>%
  dplyr::mutate(Samples = stringr::str_replace_all(Samples,"_2018_"," " )) %>%
  group_by(LOCATION) %>%
  summarise(MEAN_N = mean(N),
            SD_N = sd(N))

GC_table %>%
  dplyr::filter(stringr::str_detect(Sets, "Common")) %>%
  dplyr::mutate(sum = Brest_2018 + LT_2018 +Thau_2018 ) %>%
  group_by(sum) %>%
  summarise(N=n())

GC_table %>%
  dplyr::filter(stringr::str_detect(Sets, "sing"))

# Br - LT
GC_table %>%
  dplyr::filter(stringr::str_detect(Sets, "Br Lt")) %>%
  left_join(., Table_02_binary_matrix %>%
              as.data.frame() %>%
              rownames_to_column(var = "UniqueID")) %>%
  tidyr::pivot_longer(cols = c(-UniqueID, -Brest_2018, -LT_2018, -Thau_2018, -Sets, -DATASET),
                      names_to = "ID_EXPERIMENT",
                      values_to = "Pres_abs") %>%
  dplyr::filter(Pres_abs > 0)  %>%
  group_by(ID_EXPERIMENT) %>%
  summarise(N=n())

GC_table %>%
  dplyr::filter(stringr::str_detect(Sets, "Br Lt")) %>%
  left_join(., Table_02_binary_matrix %>%
              as.data.frame() %>%
              rownames_to_column(var = "UniqueID")) %>%
  tidyr::pivot_longer(cols = c(-UniqueID, -Brest_2018, -LT_2018, -Thau_2018, -Sets, -DATASET),
                      names_to = "ID_EXPERIMENT",
                      values_to = "Pres_abs") %>%
  dplyr::filter(Pres_abs > 0)  %>%
  group_by(ID_EXPERIMENT) %>%
  summarise(N=n()) %>%
  tidyr::separate(ID_EXPERIMENT, c("NR","genome","LOCATION","YEARS","NSI","BROYAGE", "IND","PCR"), sep = "_", remove = TRUE) %>%
  dplyr::select(-NR, -genome, -YEARS, -NSI, -BROYAGE, -PCR) %>%
  unite(LOCATION, IND, col="Samples", sep = " ", remove=FALSE) %>%
  ggplot(aes(x = Samples, y = N, fill = LOCATION)) +
  geom_bar(stat = "identity") +
  geom_text(aes(label=N), hjust=-1) +
  theme_minimal()+
    scale_fill_manual(values=list(Brest = Brest_color, 
                                LT = LT_color)) +
  coord_flip()

# Lt - Th
GC_table %>%
  dplyr::filter(stringr::str_detect(Sets, "Lt Th")) %>%
  left_join(., Table_02_binary_matrix %>%
              as.data.frame() %>%
              rownames_to_column(var = "UniqueID")) %>%
  tidyr::pivot_longer(cols = c(-UniqueID, -Brest_2018, -LT_2018, -Thau_2018, -Sets, -DATASET),
                      names_to = "ID_EXPERIMENT",
                      values_to = "Pres_abs") %>%
  dplyr::filter(Pres_abs > 0)  %>%
  group_by(ID_EXPERIMENT) %>%
  summarise(N=n()) #%>%
  # data.table::fwrite(glue::glue("{base_path}/results/Tables/Fig03A_analysis_Lt_Th.csv"),
  # sep = "\t",
  # quote=FALSE,
  # row.names = FALSE,
  # col.names = TRUE)

GC_table %>%
  dplyr::filter(stringr::str_detect(Sets, "Lt Th")) %>%
  left_join(., Table_02_binary_matrix %>%
              as.data.frame() %>%
              rownames_to_column(var = "UniqueID")) %>%
  tidyr::pivot_longer(cols = c(-UniqueID, -Brest_2018, -LT_2018, -Thau_2018, -Sets, -DATASET),
                      names_to = "ID_EXPERIMENT",
                      values_to = "Pres_abs") %>%
  dplyr::filter(Pres_abs > 0)  %>%
  group_by(ID_EXPERIMENT) %>%
  summarise(N=n()) %>%
  tidyr::separate(ID_EXPERIMENT, c("NR","genome","LOCATION","YEARS","NSI","BROYAGE", "IND","PCR"), sep = "_", remove = TRUE) %>%
  dplyr::select(-NR, -genome, -YEARS, -NSI, -BROYAGE, -PCR) %>%
  unite(LOCATION, IND, col="Samples", sep = " ", remove=FALSE) %>%
  ggplot(aes(x = Samples, y = N, fill = LOCATION)) +
  geom_bar(stat = "identity") +
  geom_text(aes(label=N), hjust=-1) +
  theme_minimal()+
    scale_fill_manual(values=list(Thau = Thau_color, 
                                LT = LT_color)) +
  coord_flip()
```

## Location of the positions within the OsHV-1 genome

```{r}
Consensus_size = 186356
var_CG_across_NR_consensus <- Table_02_binary_matrix %>%
  as.data.frame() %>%
  rownames_to_column("UniqueID") %>%
  pivot_longer(-UniqueID) %>%
  tidyr::separate(UniqueID, c("POS","MODIFICATION"), sep = "_", remove = FALSE) %>%
  dplyr::mutate(POS = as.numeric(POS)) %>%
  dplyr::filter(value != 0) %>%
  tidyr::separate(name, c("NR","genome","LOCATION","YEARS","NSI","BROYAGE", "IND","PCR"), sep = "_", remove = TRUE) %>%
  dplyr::select(-NR, -genome, -YEARS, -NSI, -BROYAGE, -PCR) %>%
  unite(LOCATION, IND, col="Samples", remove=FALSE) %>%
  left_join(.,GC_table) %>%
  dplyr::mutate(Specificite = 
                  case_when(
                    stringr::str_detect(Sets,"_full") ~ "Specific",
                    stringr::str_detect(Sets," ") ~ "Shared",
                    stringr::str_detect(Sets,"Common") ~ "Common",
                    TRUE ~ "normal"
                    ),
                 Region = case_when(
                  POS < 164287 ~ "UL",
                  POS > 164288 & POS < 171692 ~ "TRL",
                  POS > 171693 & POS < 173202 ~ "X",
                  POS > 173203 & POS < 182984 ~ "IRL",
                  TRUE ~ "US")) %>%
  # dplyr::filter(Specificite != "Common") %>% 
  left_join(., plot_order_samples)


var_CG_across_NR_consensus  %>%
ggplot(aes(x = POS, y = forcats::fct_reorder2(Samples, ORDER_samples, ORDER), fill = LOCATION )) +
  geom_point(aes(color = LOCATION, shape= Specificite))+
    xlim(0, Consensus_size) +
  theme_minimal()+
  geom_vline(xintercept = Consensus_size, linetype = "dotdash", size = 0.5, color = "#5C5C5C") +
  geom_vline(xintercept = 0, size = 0.25, color = "black") +
  # geom_hline(yintercept = 0, size = 1, color = "black") +
    # facet_grid(Samples~.) +
  theme(axis.line = element_line(size = 0.5, colour = "black")) +
    scale_color_manual(values=list(Brest = Brest_color, 
                                LT = LT_color,
                                Thau = Thau_color)) +
  scale_shape_manual(values=c(1,4, 3))

# Table_02_binary_matrix %>%
#   as.data.frame() %>%
#   rownames_to_column("UniqueID") %>%
#   pivot_longer(-UniqueID) %>%
#   tidyr::separate(UniqueID, c("POS","MODIFICATION"), sep = "_", remove = FALSE) %>%
#   dplyr::mutate(POS = as.numeric(POS)) %>%
#   dplyr::filter(value != 0) %>%
#   tidyr::separate(name, c("NR","genome","LOCATION","YEARS","NSI","BROYAGE", "IND","PCR"), sep = "_", remove = TRUE) %>%
#   dplyr::select(-NR, -genome, -YEARS, -NSI, -BROYAGE, -PCR) %>%
#   unite(LOCATION, IND, col="Samples", remove=FALSE) %>%
#   ggplot(aes(x = POS, fill = LOCATION )) +
#   geom_histogram(bins = 1000) +
#     xlim(0, Consensus_size) +
#   theme_minimal()+
#   geom_vline(xintercept = Consensus_size, linetype = "dotdash", size = 0.5, color = "#5C5C5C") +
#   geom_vline(xintercept = 0, size = 0.25, color = "black") +
#   # geom_hline(yintercept = 0, size = 1, color = "black") +
#     # facet_grid(LOCATION~.) +
#   theme(axis.line = element_line(size = 0.5, colour = "black")) +
#     scale_fill_manual(values=list(Brest = Brest_color, 
#                                 LT = LT_color,
#                                 Thau = Thau_color))
```

In future research I encourage the author to use the notions of [jaccard index](https://en.wikipedia.org/wiki/Jaccard_index) for the analysis of inter and intra sensembles variations.

```{r}
sessionInfo()
```