---
title: "Parse pelican results"
author: "Carine"
date: "2025"
output: 
  html_document:
    toc: true # table of content true
    number_sections: true  ## if you want number sections at each table header
    theme: united  # many options for theme, this one is my favorite.
    highlight: tango  # specifies the syntax highlighting style
    code_folding: hide
    toc_float: true
    toc_depth: 5
editor_options: 
  markdown: 
    wrap: sentence
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, cache = T, warning = FALSE, message = FALSE) 
library(tidyverse)
library(ggplot2)
library(cowplot)


output_figures <- "output_parse_pelican_results/"
dir.create(output_figures, recursive = T,showWarnings =  F)

total_output_figures <- str_c(output_figures,"/total/")
dir.create(total_output_figures, recursive = T,showWarnings =  F)

murinae_output_figures <- str_c(output_figures,"/murinae/")
dir.create(murinae_output_figures, recursive = T,showWarnings =  F)

```

# Load data



## Species tree and attributes

```{r}
library(ggtree)
library(treeio)
library(ape)


total_nhx_file <- "../data/trees_and_metadata/local_online_species_tree_rep_8_20200605.raxml.annotated_habitat_filtered_rotated.nhx"
murinae_nhx_file <- "../data/trees_and_metadata/local_online_species_tree_rep_8_20200605.raxml.annotated_habitat_filtered_murinae_rotated.nhx"

total_nhx_tree <- read.nhx(total_nhx_file)
murinae_nhx_tree <- read.nhx(murinae_nhx_file)

total_tree_plot <- ggtree(total_nhx_tree,ladderize=F)
total_tips_order <- get_taxa_name(total_tree_plot)

murinae_tree_plot <- ggtree(murinae_nhx_tree,ladderize=F)
murinae_tips_order <- get_taxa_name(murinae_tree_plot)

```

### Visualize metadata

```{r}
species_attributes <- read_tsv("../data/trees_and_metadata/attributes_Seq_dataonline_20220630.tsv")
species_attributes <- species_attributes %>% filter(! species_attributes$leaf_name %in% c("Arvicanthis_niloticus_lmo","Nannomys_cot","Peromyscus_leucopus"))

species_attributes_l <- species_attributes %>% dplyr::select(leaf_name,total,ancient_trans,recent_trans,murinae,Aride) %>%
  pivot_longer(-c(leaf_name,Aride), names_to = "dataset",values_to = "presence") %>%
  mutate(presence = ifelse(presence=="no","no","yes"),
         Aride2 = ifelse(presence=="no","abs",Aride),
         dataset = gsub("_trans","",dataset),
         dataset=factor(dataset,level = c("total","ancient","recent","murinae"))
         )

total_species_attributes <- species_attributes %>% filter(total != "no") %>%
  mutate(leaf_name_f = factor(leaf_name, level=total_tips_order),
         label_r = leaf_name) %>%
  relocate(label_r)

murinae_species_attributes <- species_attributes %>% filter(murinae != "no") %>% 
  mutate(leaf_name_f = factor(leaf_name, level=murinae_tips_order),
         label_r = leaf_name) %>%
  relocate(label_r)
```


```{r, fig.height=10}
library(ggnewscale)
library(ggExtra)
library(ggtreeExtra)
#BiocManager::install("ggtreeExtra")

species_attributes <- species_attributes %>% 
  mutate(leaf_name_f = factor(leaf_name, level=total_tips_order),
         label_r = leaf_name) %>%
  relocate(label_r)

p_tree <- total_tree_plot %<+%
  species_attributes  +
  geom_treescale(x = 0, y = 1) +
  geom_tiplab(aes(fill = Aride),
              geom = "label",
              alpha=0.5,
              align=TRUE,
              size=3,
              linesize=.3)  + geom_rootedge(rootedge = 0.01) + 
  geom_tippoint(aes(shape = total, color = Aride), size = 4) +
  scale_color_manual(values = c("yes"="orange","no"="darkgreen")) +
  scale_fill_manual(values = c("yes"="orange","no"="darkgreen","abs"="grey90")) +
  scale_shape_manual(name = "Expression", values = c("yes"=15,"only_s"=16)) +
  geom_fruit(data =species_attributes_l,
             geom=geom_tile,
             mapping=aes(y=leaf_name, x=dataset,fill=Aride2),
             color = "black",
             offset = 0.20,
             size = 0.02,
             axis.params=list(
                         axis="x",
                         title = "Dataset",
                         text.angle = 45,
                         text.size = 2,
                         line.size = 0,
                         vjust = 0.5
                     ))

p_tree

```


### Add transitions


```{r, fig.height=10}
species_attributes_with_transition <- species_attributes %>%
  dplyr::select(label_r, Aride) %>%
  mutate(transition = case_when(
    label_r %in% c("Mus_macedonicus") ~ "T1",
    label_r %in% c("Nannomys_indutus") ~ "T2",
    label_r %in% c("Mastomys_coucha","Mastomys_erythroleucus","Mastomys_kollmanspergeri") ~ "T3",
    label_r %in% c("Apodemus_mystacinus") ~ "T4",
    label_r %in% c("Micaelamys_namaquensis") ~ "T5",
    label_r %in% c("Arvicanthis_niloticus_sen") ~ "T6",
    label_r %in% c("Acomys_dimidiatus") ~ "T7",
    label_r %in% c("Meriones_unguiculatus") ~ "T8",
    label_r %in% c("Peromyscus_eremicus") ~ "T9",
    label_r %in% c("Cricetulus_griseus_crigri","Mesocricetus_auratus") ~ "T10",
    label_r %in% c("Saccostomus_campestris") ~ "T11",
    label_r %in% c("Eospalax_fontanierii","Nannospalax_galili") ~ "T12",
    label_r %in% c("Jaculus_jaculus") ~ "T13",
    label_r %in% c("Chaetodipus_baileyi") ~ "T14",
    label_r %in% c("Dipodomys_spectabilis2","Dipodomys_ordii") ~ "T15",
    label_r %in% c("Fukomys_micklemi","Fukomys_damarensis","Fukomys_mechowii") ~ "T16",
    label_r %in% c("Heterocephalus_glaber_male") ~ "T17",
    label_r %in% c("Chinchilla_lanigera") ~ "T18",
    label_r %in% c("Octodon_degus") ~ "T19",
    label_r %in% c("Urocitellus_parryii") ~ "T20",
    label_r %in% c("Spermophilus_dauricus") ~ "T21",
    T ~ "M"
  )
  )
```


```{r, fig.height=6, fig.width=3}
murinae_transition <- murinae_species_attributes %>%
  mutate(leaf_name_f = factor(leaf_name, levels = rev(murinae_tips_order)),
         label_r = leaf_name) %>%
  arrange(leaf_name_f) %>%
  mutate (y = 1:n()) %>%
  dplyr::select(label_r, Aride, leaf_name_f,y) %>%
  left_join(species_attributes_with_transition) %>%
  arrange(-y)


# repérer les positions où la transition change
murinae_limites_y <- murinae_transition %>% 
  filter(transition != "M") %>%
  group_by(transition,Aride) %>%
  summarise(min_y = min(y)-0.4,
            max_y = max(y)+0.4,
            sp_l = paste0(label_r, collapse = ","))

murinae_df_nodes <- murinae_tree_plot$data %>% left_join(murinae_transition)


# Filtrer les nœuds selon une métadonnée, par exemple "group"
murinae_nodes_of_interest <- murinae_df_nodes %>%
  filter(Aride == "yes") %>%
  group_by(transition) %>%
  summarise(node_MRCA =  ifelse(length(label) > 1,
                                getMRCA(as.phylo(murinae_nhx_tree), tip =  label),
                                node)
  )

murinae_p <- murinae_tree_plot %<+%
  murinae_species_attributes  +
  geom_treescale(x=0, y=1) +
  geom_tiplab(aes(fill=Aride),
              geom = "label",
              alpha=0.5,
              align=TRUE,
              size=5,
              linesize=.1,
              hjust = 1,
              offset = 0.1)  +
  geom_rootedge(rootedge = 0.001) + 
  geom_tippoint(aes(shape = total, color = Aride), size = 5) +
  scale_color_manual(values = c("yes"="orange","no"="darkgreen")) +
  scale_fill_manual(values = c("yes"="orange","no"="darkgreen","abs"="grey90")) +
  scale_shape_manual(name = "Expression", values = c("yes"=15,"only_s"=16)) + 
  theme(legend.position="top") +
  guides(fill=guide_legend(nrow=1,byrow=TRUE),
         color=guide_legend(override.aes=list(fill=NA))) +
  facet_wrap(~"Murinae") +
  geom_hilight(node = murinae_nodes_of_interest$node_MRCA, fill = c("grey85"),to.bottom = T,  xmax = Inf, type="auto", size = 1.8, color = "white"
              )

murinae_p
```



```{r, fig.height=12, fig.width=5}
total_transition <- total_species_attributes %>%
  mutate(leaf_name_f = factor(leaf_name, levels = rev(total_tips_order)),
         label_r = leaf_name) %>%
  arrange(leaf_name_f) %>%
  mutate (y = 1:n()) %>%
  dplyr::select(label_r, Aride, leaf_name_f,y) %>%
  left_join(species_attributes_with_transition) %>%
  arrange(-y)


# repérer les positions où la transition change
total_limites_y <- total_transition %>% 
  filter(transition != "M") %>%
  group_by(transition,Aride) %>%
  summarise(min_y = min(y)-0.4,
            max_y = max(y)+0.4,
            sp_l = paste0(label_r, collapse = ","))

total_df_nodes <- total_tree_plot$data %>% left_join(total_transition)


# Filtrer les nœuds selon une métadonnée, par exemple "group"
total_nodes_of_interest <- total_df_nodes %>%
  filter(Aride == "yes") %>%
  group_by(transition) %>%
  summarise(node_MRCA =  ifelse(length(label) > 1,
                                getMRCA(as.phylo(total_nhx_tree), tip =  label),
                                node)
  )

total_p <- total_tree_plot %<+%
  total_species_attributes  +
  geom_treescale(x=0, y=1) +
  geom_tiplab(aes(fill=Aride),
              geom = "label",
              alpha=0.5,
              align=TRUE,
              size=5,
              linesize=.1,
              hjust = 1,
              offset = 0.2)  +
  geom_rootedge(rootedge = 0.001) + 
  geom_tippoint(aes(shape = total, color = Aride), size = 3) +
  scale_color_manual(values = c("yes"="orange","no"="darkgreen")) +
  scale_fill_manual(values = c("yes"="orange","no"="darkgreen","abs"="grey90")) +
  scale_shape_manual(name = "Expression", values = c("yes"=15,"only_s"=16)) + 
  theme(legend.position="top") +
  guides(fill=guide_legend(nrow=1,byrow=TRUE),
         color=guide_legend(override.aes=list(fill=NA))) +
  facet_wrap(~"Total_seq") +
  geom_hilight(node = total_nodes_of_interest$node_MRCA, fill = c("grey85"),to.bottom = T,  xmax = Inf, type="auto", size = 1.8, color = "white"
              )

total_p
```

```{r, fig.height=10, fig.width=6}
trees_list <- list(Murinae = murinae_p,
                   Total = total_p)

```

## Get NCBI ID (ENTREZ ID) to run enrichment on Reactome database

```{r, fig.height=10}
library(biomaRt)
#mus_mart <-  biomaRt::useEnsembl(
#  biomart = "genes",
#  dataset = "mmusculus_gene_ensembl",
#  host = "http://jan2020.archive.ensembl.org"
#) #version 99

mus_mart <-  biomaRt::useEnsembl(
  biomart = "genes",
  dataset = "mmusculus_gene_ensembl",
  host = "https://apr2020.archive.ensembl.org"
) #version 100


MGI2ENTREZ <- biomaRt::getBM(
  attributes=c("mgi_symbol","entrezgene_id"),
  mart = mus_mart) %>%
  dplyr::filter(!is.na(entrezgene_id))

MGI2ENTREZ_ok <- MGI2ENTREZ %>%
   dplyr::filter(!duplicated(mgi_symbol))

MGI2ENTREZ_ok[MGI2ENTREZ_ok$mgi_symbol == "Gm17657","entrezgene_id"] = 102640495 #https://www.ncbi.nlm.nih.gov/gene/102640495
MGI2ENTREZ_ok[MGI2ENTREZ_ok$mgi_symbol == "Pabpn1","entrezgene_id"] = 54196 #https://www.ncbi.nlm.nih.gov/gene/54196
```

## load expr data

```{r, eval = T}
DE_output_dir="../data/expression/"

total_count_df<- read.table("../data/expression/dataset_tot_final/Table_resBasemean.tsv", row.names = 1)

total_DE_score_df <- read_tsv(paste0(DE_output_dir,"dataset_tot_final/Table_DE_results_cond_season.tsv"),
                              col_names = c("MGI","baseMean","log2FoldChange","lfcSE","stat","pvalue","padj","is_DE"), skip = 1) %>%
  mutate(dataset = "Total",
         log2FoldChange =-log2FoldChange)
  


murinae_DE_score_df <- read_tsv(paste0(DE_output_dir,"dataset_muridae_final/Table_DE_results_cond_season.tsv"),
                             col_names = c("MGI","baseMean","log2FoldChange","lfcSE","stat","pvalue","padj","is_DE"), skip = 1) %>%
  mutate(dataset = "Murinae",
         log2FoldChange =-log2FoldChange)
      

ancient_DE_score_df <- read_tsv(paste0(DE_output_dir,"dataset_ancient_final/Table_DE_results_cond_season.tsv"),
                              col_names = c("MGI","baseMean","log2FoldChange","lfcSE","stat","pvalue","padj","is_DE"), skip = 1) %>%
  mutate(dataset = "Ancient",
         log2FoldChange =-log2FoldChange)

table(ancient_DE_score_df$is_DE)
table(ancient_DE_score_df$is_DE, ancient_DE_score_df$log2FoldChange>0)

recent_DE_score_df <- read_tsv(paste0(DE_output_dir,"dataset_recent_final/Table_DE_results_cond_season.tsv"),
                              col_names = c("MGI","baseMean","log2FoldChange","lfcSE","stat","pvalue","padj","is_DE"), skip = 1) %>%
  mutate(dataset = "Recent",
         log2FoldChange =-log2FoldChange)

table(recent_DE_score_df$is_DE)
table(recent_DE_score_df$is_DE, recent_DE_score_df$log2FoldChange>0)

gene_DE_score_df <- bind_rows(total_DE_score_df,murinae_DE_score_df) %>%
  left_join(MGI2ENTREZ_ok, by = c ("MGI"="mgi_symbol")) %>%
  mutate(SYMBOL=MGI) %>%
  dplyr::rename(ENTREZID="entrezgene_id")
dim(gene_DE_score_df)

table(is.na(gene_DE_score_df$ENTREZID))
table(is.na(gene_DE_score_df$SYMBOL))

gene_DE_score_df %>% filter(is_DE == "FDR<0.1") %>%
   ggplot(aes(x=dataset)) +
   geom_bar(fill = "grey") +
   geom_text(stat = "count", aes(label = ..count.., y = ..count..), vjust = -0.5) +
   theme_bw () +
   geom_hline(yintercept = 800, size = 0)
```




## get alignement attributes per dataset

```{r}
total_true_trim_aa_fasta_dir = "../data/aa_alignment_hmmcleaner_3min_sp_renamed_onlygap_removed_seuil_0.1_filtered_sequences/"
murinae_true_trim_aa_fasta_dir = "../data/aa_alignment_hmmcleaner_3min_sp_renamed_onlygap_removed_seuil_0.1_only_in_dataset_murinae/"


total_genes_attributes <- tibble(
  dir_aa_path = c(total_true_trim_aa_fasta_dir) ) %>%
  group_modify(~ { tibble(fasta_fn=list.files(.x$dir_aa_path, full.names = T))}) %>%
  bind_rows() %>%
  mutate(alignment = basename(fasta_fn) %>% str_remove(".fasta"),
         dataset = "Total") %>%
  separate(alignment,into ="MGI", sep = "_", remove = F)

murinae_genes_attributes <- tibble(
  dir_aa_path = c(murinae_true_trim_aa_fasta_dir) ) %>%
  group_modify(~ { tibble(fasta_fn=list.files(.x$dir_aa_path, full.names = T))}) %>%
  bind_rows() %>%
  mutate(alignment = basename(fasta_fn) %>% str_remove(".fasta"),
         dataset = "Murinae") %>%
  separate(alignment,into ="MGI", sep = "_", remove = F)

## identify duplicated MGI ...
total_duplicated_MGI <- total_genes_attributes %>% group_by(dataset,MGI) %>%
 summarise(n=n()) %>%
  filter(n > 1)


total_genes_attributes <- total_genes_attributes %>%
  filter(! MGI  %in% total_duplicated_MGI$MGI )
  


murinae_duplicated_MGI <- murinae_genes_attributes %>% group_by(dataset,MGI) %>%
 summarise(n=n()) %>%
  filter(n > 1)


murinae_genes_attributes <- murinae_genes_attributes %>%
  filter(! MGI  %in% murinae_duplicated_MGI$MGI )
  


duplicated_MGI <- bind_rows(total_duplicated_MGI, murinae_duplicated_MGI)


## add nb species mesic/xeric per gene

mesic_species <- species_attributes %>% filter(Aride == "no") %>% pull(leaf_name)
xeric_species <- species_attributes %>% filter(Aride == "yes") %>% pull(leaf_name)


get_nb_x_M_species <- function(fasta_tmp_fn) {
  #fasta_tmp_fn = "../data/input_ali/aa_alignment_hmmcleaner_3min_sp_renamed_onlygap_removed_seuil_0.1/Prkar1b_4PTSB.renamed.cleaned.aa.fasta"
  #get_nb_x_M_species("../data/input_ali/aa_alignment_hmmcleaner_3min_sp_renamed_onlygap_removed_seuil_0.1/Prkar1b_4PTSB.renamed.cleaned.aa.fasta")

  fasta_tmp <- readLines(fasta_tmp_fn)
  
  sequence_names <- grep("^>",fasta_tmp, value = T)
  sequence_names <- gsub("^>","",sequence_names)
  
  nb_mesic_species <- sum(sequence_names %in% mesic_species)
  nb_xeric_species <- sum(sequence_names %in% xeric_species)
  nb_seqs <- length(sequence_names)
  
  tibble(nb_mesic_species = nb_mesic_species,
         nb_xeric_species = nb_xeric_species,
         nb_seqs = nb_seqs)
}


total_genes_attributes_tmp2 <- total_genes_attributes %>%
  rowwise() %>%
  mutate(nb_seq_df = get_nb_x_M_species(fasta_fn)) %>%
  unnest(nb_seq_df) %>%
  ungroup()


murinae_genes_attributes_tmp2 <- murinae_genes_attributes %>%
  rowwise() %>%
  mutate(nb_seq_df = get_nb_x_M_species(fasta_fn)) %>%
  unnest(nb_seq_df) %>%
  ungroup()


table(file.exists(murinae_genes_attributes_tmp2$fasta_fn))
table(file.exists(total_genes_attributes_tmp2$fasta_fn))

ggplot(murinae_genes_attributes_tmp2, aes(x = nb_mesic_species+0.5, y = nb_xeric_species+0.5)) +
  geom_bin2d(binwidth = 1) +
  scale_fill_gradientn(colours= rev(rainbow(5))) +
  geom_abline(slope = 1, linetype = "dashed", color = "grey") +
  theme_bw() + facet_wrap(~dataset, scale = "free")

ggplot(total_genes_attributes_tmp2, aes(x = nb_mesic_species+0.5, y = nb_xeric_species+0.5)) +
  geom_bin2d(binwidth = 1) +
  scale_fill_gradientn(colours= rev(rainbow(5))) +
  geom_abline(slope = 1, linetype = "dashed", color = "grey") +
  theme_bw() + facet_wrap(~dataset, scale = "free")
```



```{r}
# add expression data
total_genes_attributes_ok <- total_genes_attributes_tmp2 %>%
  left_join(gene_DE_score_df %>% dplyr::select(MGI,baseMean,log2FoldChange,padj,dataset))  %>%
  dplyr::rename(padj_DEseq2 = "padj",
                log2FoldChange_DEseq2 = "log2FoldChange",
                baseMean_DEseq2 = "baseMean")

murinae_genes_attributes_ok <- murinae_genes_attributes_tmp2 %>%
  left_join(gene_DE_score_df %>% dplyr::select(MGI,baseMean,log2FoldChange,padj,dataset))  %>%
  dplyr::rename(padj_DEseq2 = "padj",
                log2FoldChange_DEseq2 = "log2FoldChange",
                baseMean_DEseq2 = "baseMean")

```


```{r fig.height=4, fig.width=8}
gene_use_in_exp_summary <- bind_rows(total_genes_attributes_ok, murinae_genes_attributes_ok)  %>% 
  filter( !is.na(padj_DEseq2)) %>%
  group_by(dataset) %>%
  summarise(min_X = min(nb_xeric_species, na.rm = T),
            mean_X = mean(nb_xeric_species, na.rm = T),
            med_X = median(nb_xeric_species, na.rm = T),
            max_X = max(nb_xeric_species, na.rm = T),
            min_M = min(nb_mesic_species, na.rm = T),
            mean_M = mean(nb_mesic_species, na.rm = T),
            med_M = median(nb_mesic_species, na.rm = T),
            max_M = max(nb_mesic_species, na.rm = T))

## apply filter to use only gene with good expr to be sure to have good alignment.

threshold_df <- tibble( dataset = c("Ancient","Murinae","Recent","Total"),
        min_X = c(4,6,4,25),
        min_M = c(20,9,20,20))

threshold_df <- tibble( dataset = c("Murinae","Total"),
        min_X = c(7,25),
        min_M = c(9,20))
```



```{r fig.height=4, fig.width=8}
ggplot(total_genes_attributes_ok , aes(x = nb_mesic_species+0.5, y = nb_xeric_species+0.5)) +
  geom_bin2d(binwidth = 1) +
  scale_fill_gradientn(colours= rev(rainbow(5))) +
  geom_abline(slope = 1, linetype = "dashed", color = "grey") +
  geom_hline(data=threshold_df %>% filter(dataset == "Total"), aes(yintercept =min_X)) +
  geom_vline(data=threshold_df %>% filter(dataset == "Total"), aes(xintercept =min_M)) +
  theme_bw() + facet_wrap(~dataset+is.na(padj_DEseq2), scale = "free")



total_genes_attributes_ok_filtered <- total_genes_attributes_ok %>%
  left_join(threshold_df) %>%
  mutate(is_nearly_complete = nb_xeric_species >= min_X & nb_mesic_species >= min_M,
         has_duplicated_MGI =  paste0(MGI, dataset)  %in% paste0(duplicated_MGI$MGI, duplicated_MGI$dataset))

dim(total_genes_attributes_ok_filtered)

 table(total_genes_attributes_ok_filtered$is_nearly_complete)
 table(is.na(total_genes_attributes_ok_filtered$padj_DEseq2))
 table(total_genes_attributes_ok_filtered$is_nearly_complete)


table(EXP=!is.na(total_genes_attributes_ok_filtered$padj_DEseq2), SEQ=total_genes_attributes_ok_filtered$is_nearly_complete)


total_genes_attributes_ok_filtered <- total_genes_attributes_ok_filtered %>% filter(!is.na(padj_DEseq2), is_nearly_complete)

```


```{r fig.height=4, fig.width=8}
ggplot(murinae_genes_attributes_ok , aes(x = nb_mesic_species+0.5, y = nb_xeric_species+0.5)) +
  geom_bin2d(binwidth = 1) +
  scale_fill_gradientn(colours= rev(rainbow(5))) +
  geom_abline(slope = 1, linetype = "dashed", color = "grey") +
  geom_hline(data=threshold_df %>% filter(dataset == "Murinae"), aes(yintercept =min_X)) +
  geom_vline(data=threshold_df %>% filter(dataset == "Murinae"), aes(xintercept =min_M)) +
  theme_bw() + facet_wrap(~dataset+is.na(padj_DEseq2), scale = "free")


murinae_genes_attributes_ok_filtered <- murinae_genes_attributes_ok %>%
  left_join(threshold_df) %>%
  mutate(is_nearly_complete = nb_xeric_species >= min_X & nb_mesic_species >= min_M,
         has_duplicated_MGI =  paste0(MGI, dataset)  %in% paste0(duplicated_MGI$MGI, duplicated_MGI$dataset))

dim(murinae_genes_attributes_ok_filtered)

table(murinae_genes_attributes_ok_filtered$is_nearly_complete)
table(is.na(murinae_genes_attributes_ok_filtered$padj_DEseq2))
table(murinae_genes_attributes_ok_filtered$is_nearly_complete)


table(EXP=!is.na(murinae_genes_attributes_ok_filtered$padj_DEseq2), SEQ=murinae_genes_attributes_ok_filtered$is_nearly_complete)


murinae_genes_attributes_ok_filtered <- murinae_genes_attributes_ok_filtered %>% filter(!is.na(padj_DEseq2), is_nearly_complete)
```

## Pelican results

```{r}
total_pelican_tsv = "../output_pelican_2025_Carine/total/aridity_2cats_aa_alignment_hmmcleaner_3min_sp_renamed_onlygap_removed_seuil_0.1_filtered_sequences_M1_with_simu/"
murinae_pelican_tsv = "../output_pelican_2025_Carine/murinae/aridity_2cats_aa_alignment_hmmcleaner_3min_sp_renamed_onlygap_removed_seuil_0.1_filtered_sequences_M1_with_simu/"
```

```{r}
total_true_all_sites_fn <- str_c(total_pelican_tsv,"all_sites.tsv")
total_true_all_sites <- read_tsv(total_true_all_sites_fn) %>% filter(alignment %in% (total_genes_attributes_ok_filtered %>% filter(dataset == "Total", is_nearly_complete) %>% pull(alignment))) %>%
  mutate(dataset = "Total",
         dp = "T",
         aagtr_pval = ifelse(is.na(aagtr_pval), 1, aagtr_pval),
         log10_aagtr_pval = log10(aagtr_pval)
         )


murinae_true_all_sites_fn <- str_c(murinae_pelican_tsv,"all_sites.tsv")
murinae_true_all_sites <- read_tsv(murinae_true_all_sites_fn) %>% filter(alignment %in% (murinae_genes_attributes_ok_filtered %>% filter(dataset == "Murinae", is_nearly_complete) %>% pull(alignment))) %>%
  mutate(dataset = "Murinae",
         dp = "M",
         aagtr_pval = ifelse(is.na(aagtr_pval), 1, aagtr_pval),
         log10_aagtr_pval = log10(aagtr_pval)
         )

score_all_sites_all_dataset_df <- bind_rows(total_true_all_sites,murinae_true_all_sites)

score_all_sites_all_dataset_df %>%
  ggplot(aes(x = dataset, y =log10_aagtr_pval )) +
  geom_violin() +
  theme_bw()  + facet_wrap(~dataset, scale = "free_x")
```

# Create spike alignement based on observed AA frequencies


## get score per pos per dataset

```{r}
#get best pval per MGI

total_score_all_tested_sites_df <- total_true_all_sites %>%
  group_by(alignment,dataset) %>%
  mutate(ali_length = max(site, na.rm=T),
  nb_non_constant_site = sum(tag != "Constant site") ) %>%
  ungroup() %>%
  group_by(dataset) %>%
  filter(tag == "Test") %>%
  mutate(rank_aagtr_all_sites = rank(aagtr_pval),
         log10_aagtr_pval = log10(aagtr_pval)
         ) %>%
  arrange(rank_aagtr_all_sites)


total_best_pval_per_MGI_ok <- 
  total_score_all_tested_sites_df %>%
  group_by(alignment,dataset) %>%
  slice_min(log10_aagtr_pval, n = 1, with_ties=F) %>%
  rename(
         best_site = 'site',
         nseq_best_site = "nseq",
         naa_best_site = "naa",
         best_aagtr_pval = "aagtr_pval")  %>%
  ungroup() %>%
  group_by(dataset) %>%
  mutate(rank_best_sites = rank(best_aagtr_pval),
         log10_best_aagtr_pval = log10(best_aagtr_pval)
         ) %>%
  arrange(rank_best_sites) %>%
  left_join( bind_rows(total_genes_attributes_ok_filtered,murinae_genes_attributes_ok_filtered) %>% dplyr::select(alignment,MGI,fasta_fn, dataset)) %>%
  dplyr::select(alignment,MGI,nb_non_constant_site,nseq_best_site,naa_best_site,best_aagtr_pval,best_site,dataset,dp,rank_aagtr_all_sites,ali_length,nb_non_constant_site,log10_best_aagtr_pval,fasta_fn) %>%
    ungroup() 

murinae_score_all_tested_sites_df <- murinae_true_all_sites %>%
  group_by(alignment,dataset) %>%
  mutate(ali_length = max(site, na.rm=T),
  nb_non_constant_site = sum(tag != "Constant site") ) %>%
  ungroup() %>%
  group_by(dataset) %>%
  filter(tag == "Test") %>%
  mutate(rank_aagtr_all_sites = rank(aagtr_pval),
         log10_aagtr_pval = log10(aagtr_pval)
         ) %>%
  arrange(rank_aagtr_all_sites)


murinae_best_pval_per_MGI_ok <- 
  murinae_score_all_tested_sites_df %>%
  group_by(alignment,dataset) %>%
  slice_min(log10_aagtr_pval, n = 1, with_ties=F) %>%
  rename(
         best_site = 'site',
         nseq_best_site = "nseq",
         naa_best_site = "naa",
         best_aagtr_pval = "aagtr_pval")  %>%
  ungroup() %>%
  group_by(dataset) %>%
  mutate(rank_best_sites = rank(best_aagtr_pval),
         log10_best_aagtr_pval = log10(best_aagtr_pval)
         ) %>%
  arrange(rank_best_sites) %>%
  left_join( bind_rows(total_genes_attributes_ok_filtered,murinae_genes_attributes_ok_filtered) %>% dplyr::select(alignment,MGI,fasta_fn, dataset)) %>%
  dplyr::select(alignment,MGI,nb_non_constant_site,nseq_best_site,naa_best_site,best_aagtr_pval,best_site,dataset,dp,rank_aagtr_all_sites,ali_length,nb_non_constant_site,log10_best_aagtr_pval,fasta_fn) %>%
    ungroup() 
  

```


## observe n_sites/min_pval relation

```{r fig.height=6, fig.width=8}
relation_min_aagtr_pval_nsites <-  bind_rows(total_best_pval_per_MGI_ok,murinae_best_pval_per_MGI_ok) %>%
  ggplot(aes(x=nb_non_constant_site, y = log10_best_aagtr_pval)) +
  theme_bw() +
  geom_point(color="grey") +
  geom_density2d() +                      # points
  geom_smooth(method = "lm", se = FALSE, color = "black")  +
  facet_wrap(~dataset)

relation_min_aagtr_pval_nsites

```



```{r}

get_seq_from_fasta_fn_pos <- function(fasta_fn, site) {
  
  #fasta_fn = score_all_sites_all_dataset_df_top_X_add_genes_att$fasta_fn[1]
  #site = score_all_sites_all_dataset_df_top_X_add_genes_att$site[1]
  fasta_df <- phylotools::read.fasta(fasta_fn)
  fasta_df$AA <- str_sub(fasta_df$seq.text,site,site)
  
  tibble(seqnames = fasta_df$seq.name,
         AA = fasta_df$AA)
  
}


total_best_pval_per_MGI_ok_add_genes_att <- 
  total_best_pval_per_MGI_ok %>%
  dplyr::select(alignment,MGI,best_site,best_aagtr_pval,dataset,fasta_fn,dp) %>%
  arrange(best_aagtr_pval) %>% mutate(MGI_site = paste0(MGI,"_",best_site,"_",dp),
                                 labels_simplified = paste0(MGI,"_",best_site))


murinae_best_pval_per_MGI_ok_add_genes_att <- 
  murinae_best_pval_per_MGI_ok %>%
  dplyr::select(alignment,MGI,best_site,best_aagtr_pval,dataset,fasta_fn,dp) %>%
  arrange(best_aagtr_pval) %>% mutate(MGI_site = paste0(MGI,"_",best_site,"_",dp),
                                 labels_simplified = paste0(MGI,"_",best_site))



total_best_pval_per_MGI_ok_with_seq<- 
  total_best_pval_per_MGI_ok_add_genes_att %>%
  ungroup() %>%
  rowwise() %>%
  mutate(seq_per_sp_df = list(get_seq_from_fasta_fn_pos(fasta_fn, best_site))) %>%
  unnest(seq_per_sp_df)  %>%
  left_join(species_attributes %>% dplyr::select(leaf_name, Aride), by = c("seqnames" = "leaf_name"))  %>%
  mutate(MGI_site = factor(MGI_site, levels = total_best_pval_per_MGI_ok_add_genes_att$MGI_site),
         seqnames = factor(seqnames, levels = rev(total_tips_order))) 


murinae_best_pval_per_MGI_ok_with_seq<- 
  murinae_best_pval_per_MGI_ok_add_genes_att %>%
  ungroup() %>%
  rowwise() %>%
  mutate(seq_per_sp_df = list(get_seq_from_fasta_fn_pos(fasta_fn, best_site))) %>%
  unnest(seq_per_sp_df)  %>%
  left_join(species_attributes %>% dplyr::select(leaf_name, Aride), by = c("seqnames" = "leaf_name"))  %>%
  mutate(MGI_site = factor(MGI_site, levels = murinae_best_pval_per_MGI_ok_add_genes_att$MGI_site),
         seqnames = factor(seqnames, levels = rev(murinae_tips_order))) 



```


## calcul aa freq

```{r}
compte_lettres <- function(x) {
  # Sépare la chaîne en caractères individuels
  chars <- sort(unlist(strsplit(x, "")))
  # Utilise rle pour obtenir les lettres et leur longueur consécutive
  r <- rle(chars)
  # Concatène les résultats façon "3A 13E"
  res <- paste0(r$lengths, r$values, collapse = "")
  return(res)
}

# Exemple d'utilisation
chaine <- "AAAEEEEEEEEEEEEE"
compte_lettres(chaine)  # Retourne: "3A 13E"
```


```{r}
total_best_pval_per_MGI_ok_with_seq_compte_AA <- 
  total_best_pval_per_MGI_ok %>%
  rowwise() %>%
  filter(naa_best_site==2) %>%
  mutate(seq_per_sp_df = list(get_seq_from_fasta_fn_pos(fasta_fn, best_site))) %>%
  unnest(seq_per_sp_df)  %>%
  left_join(species_attributes %>% dplyr::select(label_r,Aride), by = c("seqnames"="label_r")) %>%
  group_by(alignment,Aride) %>%
  filter(AA != "-") %>%
  summarise(aa_u = paste0(unique(sort(AA)), collapse = ""),
  ) %>%
  pivot_wider(names_from = Aride, values_from = aa_u) %>%
  rowwise() %>%
  mutate(aa_AM = paste0(no, ">" ,yes),
         new_aa = ifelse(str_length(no) == 1, str_remove(yes,no),yes),
         new_aa_l =str_length(new_aa),
         old_aa =  ifelse(str_length(no) == 2, str_remove(no,yes),no),
         old_aa_l = str_length(old_aa)
  )


total_best_pval_per_MGI_ok_with_seq_compte_AA_bilan <-
  total_best_pval_per_MGI_ok_with_seq_compte_AA %>%
  group_by(aa_AM,new_aa,old_aa,new_aa_l,old_aa_l) %>%
  summarise(n=n()) %>%
  arrange(-n)

total_df_AA_bilan_ok <- total_best_pval_per_MGI_ok_with_seq_compte_AA_bilan %>%
  filter(new_aa_l == 1, old_aa_l==1) %>%
  dplyr::select(old_aa,new_aa,n) %>%
  group_by(old_aa,new_aa) %>%
  summarise(n = sum(n)) %>%
  arrange(-n)

dim(total_df_AA_bilan_ok)

murinae_best_pval_per_MGI_ok_with_seq_compte_AA <- 
  murinae_best_pval_per_MGI_ok %>%
  rowwise() %>%
  filter(naa_best_site==2) %>%
  mutate(seq_per_sp_df = list(get_seq_from_fasta_fn_pos(fasta_fn, best_site))) %>%
  unnest(seq_per_sp_df)  %>%
  left_join(species_attributes %>% dplyr::select(label_r,Aride), by = c("seqnames"="label_r")) %>%
  group_by(alignment,Aride) %>%
  filter(AA != "-") %>%
  summarise(aa_u = paste0(unique(sort(AA)), collapse = ""),
  ) %>%
  pivot_wider(names_from = Aride, values_from = aa_u) %>%
  rowwise() %>%
  mutate(aa_AM = paste0(no, ">" ,yes),
         new_aa = ifelse(str_length(no) == 1, str_remove(yes,no),yes),
         new_aa_l =str_length(new_aa),
         old_aa =  ifelse(str_length(no) == 2, str_remove(no,yes),no),
         old_aa_l = str_length(old_aa)
  )





murinae_best_pval_per_MGI_ok_with_seq_compte_AA_bilan <-
  murinae_best_pval_per_MGI_ok_with_seq_compte_AA %>%
  group_by(aa_AM,new_aa,old_aa,new_aa_l,old_aa_l) %>%
  summarise(n=n()) %>%
  arrange(-n)

murinae_best_pval_per_MGI_ok_with_seq_compte_AA_bilan

murinae_df_AA_bilan_ok <- murinae_best_pval_per_MGI_ok_with_seq_compte_AA_bilan %>%
  filter(new_aa_l == 1, old_aa_l==1) %>%
  dplyr::select(old_aa,new_aa,n) %>%
  group_by(old_aa,new_aa) %>%
  summarise(n = sum(n)) %>%
  arrange(-n)

dim(murinae_df_AA_bilan_ok)

```



```{r}
total_species_attributes_with_transition <- species_attributes_with_transition %>% filter(label_r %in% total_species_attributes$label_r)
murinae_species_attributes_with_transition <- species_attributes_with_transition %>% filter(label_r %in% murinae_species_attributes$label_r)

total_nb_transition <- total_species_attributes_with_transition %>% filter(Aride == "yes") %>% pull(transition) %>% unique() %>% length()
murinae_nb_transition <- murinae_species_attributes_with_transition %>% filter(Aride == "yes") %>% pull(transition) %>% unique() %>% length()


total_combi_list <- unlist(lapply(1:total_nb_transition, function(k) {
  x <-  1:total_nb_transition
  comb_nk <- combn(x, k, simplify=FALSE)
  if (length(comb_nk) > 200) {
    comb_nk <- sample(comb_nk, size=200)
  }
  comb_nk  
  }
  )
  , recursive=FALSE)

murinae_combi_list <- unlist(lapply(1:murinae_nb_transition, function(k) {
  x <-  1:murinae_nb_transition
  comb_nk <- combn(x, k, simplify=FALSE)
  if (length(comb_nk) > 200) {
    comb_nk <- sample(comb_nk, size=200)
  }
  comb_nk  
  }
  )
  , recursive=FALSE)

length(total_combi_list)
length(murinae_combi_list)


# Fonction : transforme chaque liste en chaîne binaire "00001" etc.
bin_encode <- function(idx, n = 21) {
    code <- rep(0, n)
    code[idx] <- 1
    paste(code, collapse = "")
  }

#--------
generer_aa_fasta <- function(df_aa, df_species, phenotypes_target, file = "output.fasta") {
  lignes_fasta <- character(0)
  for (j in seq_len(nrow(df_species))) {
    seq_sp <- character(0)
    species <- df_species$label_r[j]
    phenotype <- df_species$transition[j]
    header <- paste0(">", species)
    for (i in seq_len(nrow(df_aa))) {
      aa1 <- df_aa$old_aa[i]
      aa2 <- df_aa$new_aa[i]
      
      # Attribution aa selon phénotype cible
      aa <- if (phenotype %in% phenotypes_target) aa1 else aa2
      seq_sp <- c(seq_sp, aa)
    }
    lignes_fasta <- c(lignes_fasta, header, paste0(seq_sp, collapse = ""))
  }
  writeLines(lignes_fasta, file)
}

# Utilisation : aa1 va aux phénotypes "Positif", aa2 aux autres

dir.create("spike_aa_data/total/", recursive = T)

for (combi_i in total_combi_list) {
  phenotypes_target=paste0("T", combi_i)
  n_t = length(combi_i)
  combi_str = bin_encode(combi_i, n=total_nb_transition)
  output_file = str_c("spike_aa_data/total/","simu_",n_t,"_",total_nb_transition,"_",combi_str,".fasta")
  generer_aa_fasta(phenotypes_target, df_aa=total_df_AA_bilan_ok, df_species=total_species_attributes_with_transition, file = output_file)

}

dir.create("spike_aa_data/murinae/", recursive = T)

for (combi_i in murinae_combi_list) {
  phenotypes_target=paste0("T", combi_i)
  n_t = length(combi_i)
  combi_str = bin_encode(combi_i, n=murinae_nb_transition)
  output_file = str_c("spike_aa_data/murinae/","simu_",n_t,"_",murinae_nb_transition,"_",combi_str,".fasta")
  generer_aa_fasta(phenotypes_target, df_aa=murinae_df_AA_bilan_ok, df_species=murinae_species_attributes_with_transition, file = output_file)

}



```



## get perfect spike pelican results


```{r}
total_pelican_simu_results <- read_tsv("../output_pelican_2025_Carine/total/aridity_2cats_aa_alignment_total_spike_simu/all_sites.tsv")
murinae_pelican_simu_results <- read_tsv("../output_pelican_2025_Carine/murinae//aridity_2cats_aa_alignment_murinae_spike_simu/all_sites.tsv")

total_pelican_simu_results <- total_pelican_simu_results %>%
  filter(grepl("^simu", alignment)) %>%
  mutate( simu_methode = "Perfect\nspike\nsimulation",
         n_t = paste0(str_split_i(alignment, "_", 2), "_", total_nb_transition),
         combi_str = str_split_i(alignment, "_", 4),
         log10_aagtr_pval = log10(aagtr_pval)
  ) %>%
  mutate(n_t = factor(n_t, levels = paste0(1:total_nb_transition,"_",total_nb_transition)))   %>%
  mutate(dataset = "Total",
         dp = "T")


murinae_pelican_simu_results <- murinae_pelican_simu_results %>%
  filter(grepl("^simu", alignment)) %>%
  mutate(simu_methode = "Perfect\nspike\nsimulation",
         n_t = paste0(str_split_i(alignment, "_", 2), "_", murinae_nb_transition),
         combi_str = str_split_i(alignment, "_", 4),
         log10_aagtr_pval = log10(aagtr_pval)
  ) %>%
  mutate(n_t = factor(n_t, levels = paste0(1:murinae_nb_transition,"_",murinae_nb_transition))) %>%
  mutate(dataset = "Murinae",
         dp = "M")


```



```{r, fig.height=10, fig.width=20}
total_pelican_simu_results_summary <- total_pelican_simu_results %>%
  group_by(simu_methode,n_t,dataset) %>%
  summarise(s_0.5 = median(log10(aagtr_pval)),
            s_0 = min(log10(aagtr_pval)),
            s_1 = max(log10(aagtr_pval)),
            )  %>%   pivot_longer(-c(simu_methode,n_t,dataset)) %>%
  mutate(q = as.numeric(str_split_i(n_t,"_",1)),
         n = as.numeric(str_split_i(n_t,"_",1)),
         q_str = paste0(n, "/", total_nb_transition)) %>%
  mutate(q2 = as.numeric(str_split_i(name,"_",2)),
         name = fct_reorder2(name, q2, q) ) %>%
  arrange(-n) 

total_pelican_simu_results_summary


murinae_pelican_simu_results_summary <- murinae_pelican_simu_results %>%
  group_by(simu_methode,n_t,dataset) %>%
  summarise(s_0.5 = median(log10(aagtr_pval)),
            s_0 = min(log10(aagtr_pval)),
            s_1 = max(log10(aagtr_pval)),
            )  %>%   pivot_longer(-c(simu_methode,n_t,dataset)) %>%
  mutate(q = as.numeric(str_split_i(n_t,"_",1)),
         n = as.numeric(str_split_i(n_t,"_",1)),
         q_str = paste0(n, "/", murinae_nb_transition)) %>%
  mutate(q2 = as.numeric(str_split_i(name,"_",2)),
         name = fct_reorder2(name, q2, q) ) %>%
  arrange(-n) 

murinae_pelican_simu_results_summary

```

## get neutral evolution pelican results


```{r}
format_no_trailing_zeros <- function(x) {
  s <- sprintf("%.10f", x)       # 10 décimales fixes
  sub("\\.?0+$", "", s)          # supprime zéros et point inutile en fin
}
```

simulate alignment under neutral model with pastek https://gitlab.in2p3.fr/pveber/pastek/-/tree/main?ref_type=heads

```{r}
total_pelican_neutral_results <- read_tsv("../output_pelican_2025_Carine/total/total_neutrally/all_sites.tsv") 
murinae_pelican_neutral_results <- read_tsv("../output_pelican_2025_Carine/murinae/murinae_neutrally/all_sites.tsv")


total_pelican_neutral_results <- total_pelican_neutral_results %>%
  mutate(
    simu_methode = "Neutral\nsimulation",
    combi_str = paste0(rep("0",total_nb_transition), collapse = ""),
    aagtr_pval = ifelse(is.na(aagtr_pval),1,aagtr_pval),
    log10_aagtr_pval = log10(aagtr_pval)
  ) %>%
  mutate(dataset = "Total",
         dp = "T")


murinae_pelican_neutral_results <- murinae_pelican_neutral_results %>%
  mutate(
    simu_methode = "Neutral\nsimulation",
    combi_str = paste0(rep("0",murinae_nb_transition), collapse = ""),
    aagtr_pval = ifelse(is.na(aagtr_pval),1,aagtr_pval),
    log10_aagtr_pval = log10(aagtr_pval)
  ) %>%
  mutate(dataset = "Murinae",
         dp = "M")

```


```{r}
total_pelican_neutral_results_summary <- total_pelican_neutral_results %>%
  group_by(simu_methode, dataset) %>%
  summarise(n_1 = max(log10(aagtr_pval)),
            n_0.5 = median(log10(aagtr_pval)),
            n_0.05  = quantile(log10_aagtr_pval, 0.05),
            n_0.01 = quantile(log10_aagtr_pval, 0.01),
            n_0.005 = quantile(log10_aagtr_pval, 0.005),
            n_0.001 = quantile(log10_aagtr_pval, 0.001),
            n_0.0005 = quantile(log10_aagtr_pval, 0.0005),
            n_0.0001 = quantile(log10_aagtr_pval, 0.0001),
            n_0.00005 = quantile(log10_aagtr_pval, 0.00005),
            n_0.00001 = quantile(log10_aagtr_pval, 0.00001),
            n_0 = min(log10(aagtr_pval)),
            ) %>%   pivot_longer(-c(simu_methode,dataset)) %>%
  mutate(q = as.numeric(str_split_i(name,"_",2)),
         q_str = format_no_trailing_zeros(q),
         name = fct_reorder(name, q) ) %>%
  arrange(-q) 


total_pelican_neutral_results_summary



murinae_pelican_neutral_results_summary  <- murinae_pelican_neutral_results %>%
  group_by(simu_methode, dataset) %>%
  summarise(n_1 = max(log10(aagtr_pval)),
            n_0.5 = median(log10(aagtr_pval)),
            n_0.05  = quantile(log10_aagtr_pval, 0.05),
            n_0.01 = quantile(log10_aagtr_pval, 0.01),
            n_0.005 = quantile(log10_aagtr_pval, 0.005),
            n_0.001 = quantile(log10_aagtr_pval, 0.001),
            n_0.0005 = quantile(log10_aagtr_pval, 0.0005),
            n_0.0001 = quantile(log10_aagtr_pval, 0.0001),
            n_0.00005 = quantile(log10_aagtr_pval, 0.00005),
            n_0.00001 = quantile(log10_aagtr_pval, 0.00001),
            n_0 = min(log10(aagtr_pval)),
  )  %>%   pivot_longer(-c(simu_methode,dataset)) %>%
  mutate(q = as.numeric(str_split_i(name,"_",2)),
         q_str = format_no_trailing_zeros(q),
         name = fct_reorder(name, q) ) %>%
  arrange(-q) 

murinae_pelican_neutral_results_summary



```


```{r}
total_neutral_all_tested_sites_df <- total_pelican_neutral_results %>%
  group_by(alignment,dataset) %>%
  mutate(ali_length = max(site, na.rm=T),
  nb_non_constant_site = sum(tag != "Constant site") ) %>%
  ungroup() %>%
  group_by(dataset) %>%
  filter(tag == "Test") %>%
  mutate(rank_aagtr_all_sites = rank(aagtr_pval),
         log10_aagtr_pval = log10(aagtr_pval)
         ) %>%
  arrange(rank_aagtr_all_sites)



total_pelican_neutral_all_tested_results_summary <- total_neutral_all_tested_sites_df %>%
  group_by(simu_methode, dataset) %>%
  summarise(n_1 = max(log10(aagtr_pval)),
            n_0.5 = median(log10(aagtr_pval)),
            n_0.05  = quantile(log10_aagtr_pval, 0.05),
            n_0.01 = quantile(log10_aagtr_pval, 0.01),
            n_0.005 = quantile(log10_aagtr_pval, 0.005),
            n_0.001 = quantile(log10_aagtr_pval, 0.001),
            n_0.0005 = quantile(log10_aagtr_pval, 0.0005),
            n_0.0001 = quantile(log10_aagtr_pval, 0.0001),
            n_0.00005 = quantile(log10_aagtr_pval, 0.00005),
            n_0.00001 = quantile(log10_aagtr_pval, 0.00001),
            n_0 = min(log10(aagtr_pval)),
            ) %>%   pivot_longer(-c(simu_methode,dataset)) %>%
  mutate(q = as.numeric(str_split_i(name,"_",2)),
         q_str = format_no_trailing_zeros(q),
         name = fct_reorder(name, q) ) %>%
  arrange(-q) 


total_pelican_neutral_all_tested_results_summary

```

```{r}
murinae_neutral_all_tested_sites_df <- murinae_pelican_neutral_results %>%
  group_by(alignment,dataset) %>%
  mutate(ali_length = max(site, na.rm=T),
  nb_non_constant_site = sum(tag != "Constant site") ) %>%
  ungroup() %>%
  group_by(dataset) %>%
  filter(tag == "Test") %>%
  mutate(rank_aagtr_all_sites = rank(aagtr_pval),
         log10_aagtr_pval = log10(aagtr_pval)
         ) %>%
  arrange(rank_aagtr_all_sites)



murinae_pelican_neutral_all_tested_results_summary <- murinae_neutral_all_tested_sites_df %>%
  group_by(simu_methode, dataset) %>%
  summarise(n_1 = max(log10(aagtr_pval)),
            n_0.5 = median(log10(aagtr_pval)),
            n_0.05  = quantile(log10_aagtr_pval, 0.05),
            n_0.01 = quantile(log10_aagtr_pval, 0.01),
            n_0.005 = quantile(log10_aagtr_pval, 0.005),
            n_0.001 = quantile(log10_aagtr_pval, 0.001),
            n_0.0005 = quantile(log10_aagtr_pval, 0.0005),
            n_0.0001 = quantile(log10_aagtr_pval, 0.0001),
            n_0.00005 = quantile(log10_aagtr_pval, 0.00005),
            n_0.00001 = quantile(log10_aagtr_pval, 0.00001),
            n_0 = min(log10(aagtr_pval)),
            ) %>%   pivot_longer(-c(simu_methode,dataset)) %>%
  mutate(q = as.numeric(str_split_i(name,"_",2)),
         q_str = format_no_trailing_zeros(q),
         name = fct_reorder(name, q) ) %>%
  arrange(-q) 


murinae_pelican_neutral_all_tested_results_summary

```


### murinae

```{r}
murinae_neutral_sites<- bind_rows(
  murinae_pelican_neutral_results %>% filter(tag == "Test") %>% slice_sample(n=20),
   murinae_pelican_neutral_results %>% filter(tag == "Test") %>% slice_min(log10_aagtr_pval,n=20),) %>%
  mutate(fasta_fn =  str_c("../output_pelican_2025_Carine/murinae/pastek/aln/", alignment,".fasta")) %>%
  mutate(MGI =alignment,
         MGI_site = paste0(alignment,"_",site,"_",dp),
         labels_simplified = paste0(alignment,"_",site),
         MGI_site = fct_reorder(MGI_site, log10_aagtr_pval)) %>%
  arrange(MGI_site)



murinae_neutral_sites_with_seq <- murinae_neutral_sites %>%
  ungroup() %>%
  rowwise() %>%
  mutate(seq_per_sp_df = list(get_seq_from_fasta_fn_pos(fasta_fn, site))) %>%
  unnest(seq_per_sp_df)  %>%
  left_join(species_attributes %>% dplyr::select(leaf_name, Aride), by = c("seqnames" = "leaf_name"))  %>%
  mutate(MGI_site = factor(MGI_site, levels = levels(murinae_neutral_sites$MGI_site)),
         MGI = factor(MGI, levels = unique(murinae_neutral_sites$MGI)),
         seqnames = factor(seqnames, levels = rev(murinae_tips_order))) 



```


```{r fig.height=10, fig.width=15}

murinae_neutral_sites_with_seq_plot <- ggplot(murinae_neutral_sites_with_seq, aes(x = MGI_site, y = seqnames, fill = AA, label= AA, col = Aride)) +
 # geom_hline(yintercept = c(0.5,3.5,7.5,6.5,11.5,12.5,15.5,14.5,13.5)) +
  geom_rect(
    data = murinae_limites_y, 
    aes(
      xmin = -Inf,
      xmax = Inf,
      ymin = min_y,
      ymax = max_y,
      fill = Aride
    ),
    alpha = 1, # Transparence pastel
    inherit.aes = FALSE
  ) +
  geom_tile(alpha=0.8,  size = 1, width = 0.8, height=0.8, col = ifelse(murinae_neutral_sites_with_seq$Aride =="yes", "grey85","white")) + geom_text(  size = 6, color = "black") +
  scale_color_manual(values = c("darkgreen","darkorange")) +
  theme_bw()  + ggforce::facet_row(~MGI, scale = "free_x", space='free', strip.position="top")  +
  scale_x_discrete(position = "top",
                   breaks = murinae_neutral_sites_with_seq$MGI_site,
                   labels = murinae_neutral_sites_with_seq$labels_simplified) +
  theme(legend.key=element_blank(),panel.background = element_rect(fill = "white"), panel.grid.major = element_blank(),
        axis.text.x = element_text(colour = "black", size = 10, angle = 45, vjust = 0.3, hjust = 0)) +
  scale_fill_manual(values = c(palette_AA_LETTER,"yes"="grey85"))  +
  theme(legend.position="top") +guides(fill=guide_legend(nrow=1,byrow=TRUE),
                                       color=guide_legend(override.aes=list(fill=NA))) 
murinae_neutral_sites_with_seq_plot



murinae_neutral_sites_with_pval_plot <-
  bind_rows(murinae_pelican_neutral_all_tested_results_summary) %>% 
  mutate(facet = "Murinae") %>%
  ggplot()+
  geom_point(data = murinae_neutral_sites, aes(x=MGI_site, y = log10_aagtr_pval))   +
  ggrepel::geom_text_repel(aes(label = q_str, y = value), x = 0.5, min.segment.length = 0, size = 2)+
  geom_hline(yintercept = c(0,-2.5), size = 0, col = "grey") +
  theme_bw() +
  geom_hline(aes(yintercept = value, col = simu_methode))  +
  scale_color_manual(values = c("Neutral\nsimulation" = "#145db6", "Purifying\nselection\nsimulation" ="#6c2da3", "Perfect\nspike\nsimulation" = "#7c0000", "Best\nObs\nSite"="black"))  +
    scale_x_discrete(position = "top",
                   breaks = murinae_neutral_sites_with_seq$MGI_site,
                   labels = murinae_neutral_sites_with_seq$labels_simplified) +
   theme(legend.key=element_blank(),panel.background = element_rect(fill = "white"), panel.grid.major = element_blank(),
        axis.text.x = element_text(colour = "black", size = 6, angle = 45, vjust = 0.3, hjust = 0)) 



fig_murinae_neutral<-  cowplot::plot_grid(
  cowplot::plot_grid(
    cowplot::plot_grid(NULL, murinae_p  , rel_heights = c(5,16), ncol=1),
    cowplot::plot_grid(murinae_neutral_sites_with_seq_plot + theme(axis.text.y=element_blank()) + ylab(""), NULL, rel_heights = c(100,3), ncol=1),
    rel_widths = c(3,19)
  ),
  cowplot::plot_grid(NULL, murinae_neutral_sites_with_pval_plot + theme(legend.position = "none"),
                     rel_widths = c(3,19),  ncol=2),
  ncol=1,rel_heights = c(10,4)
)

 
 
 ggsave(plot=fig_murinae_neutral, filename = str_c(output_figures,"A_pelican_murinae_neutral.pdf"), width = 35, height = 10)

  
```






### total

```{r}
total_neutral_sites<- bind_rows(
  total_pelican_neutral_results %>% filter(tag == "Test") %>% slice_sample(n=30),
   total_pelican_neutral_results %>% filter(tag == "Test") %>% slice_min(log10_aagtr_pval,n=20),) %>%
  mutate(fasta_fn =  str_c("../output_pelican_2025_Carine/total/pastek/aln/", alignment,".fasta")) %>%
  mutate(MGI =alignment,
         MGI_site = paste0(alignment,"_",site,"_",dp),
         labels_simplified = paste0(alignment,"_",site),
         MGI_site = fct_reorder(MGI_site, log10_aagtr_pval)) %>%
  arrange(MGI_site)



total_neutral_sites_with_seq <- total_neutral_sites %>%
  ungroup() %>%
  rowwise() %>%
  mutate(seq_per_sp_df = list(get_seq_from_fasta_fn_pos(fasta_fn, site))) %>%
  unnest(seq_per_sp_df)  %>%
  left_join(species_attributes %>% dplyr::select(leaf_name, Aride), by = c("seqnames" = "leaf_name"))  %>%
  mutate(MGI_site = factor(MGI_site, levels = levels(total_neutral_sites$MGI_site)),
         MGI = factor(MGI, levels = unique(total_neutral_sites$MGI)),
         seqnames = factor(seqnames, levels = rev(total_tips_order))) 



```


```{r fig.height=10, fig.width=15}

total_neutral_sites_with_seq_plot <- ggplot(total_neutral_sites_with_seq, aes(x = MGI_site, y = seqnames, fill = AA, label= AA, col = Aride)) +
 # geom_hline(yintercept = c(0.5,3.5,7.5,6.5,11.5,12.5,15.5,14.5,13.5)) +
  geom_rect(
    data = total_limites_y, 
    aes(
      xmin = -Inf,
      xmax = Inf,
      ymin = min_y,
      ymax = max_y,
      fill = Aride
    ),
    alpha = 1, # Transparence pastel
    inherit.aes = FALSE
  ) +
  geom_tile(alpha=0.8,  size = 1, width = 0.8, height=0.8, col = ifelse(total_neutral_sites_with_seq$Aride =="yes", "grey85","white")) + geom_text(  size = 6, color = "black") +
  scale_color_manual(values = c("darkgreen","darkorange")) +
  theme_bw()  + ggforce::facet_row(~MGI, scale = "free_x", space='free', strip.position="top")  +
  scale_x_discrete(position = "top",
                   breaks = total_neutral_sites_with_seq$MGI_site,
                   labels = total_neutral_sites_with_seq$labels_simplified) +
  theme(legend.key=element_blank(),panel.background = element_rect(fill = "white"), panel.grid.major = element_blank(),
        axis.text.x = element_text(colour = "black", size = 10, angle = 45, vjust = 0.3, hjust = 0)) +
  scale_fill_manual(values = c(palette_AA_LETTER,"yes"="grey85"))  +
  theme(legend.position="top") +guides(fill=guide_legend(nrow=1,byrow=TRUE),
                                       color=guide_legend(override.aes=list(fill=NA))) 
total_neutral_sites_with_seq_plot



total_neutral_sites_with_pval_plot <-
  bind_rows(total_pelican_neutral_all_tested_results_summary) %>% 
  mutate(facet = "Total") %>%
  ggplot()+
  geom_point(data = total_neutral_sites, aes(x=MGI_site, y = log10_aagtr_pval))   +
  ggrepel::geom_text_repel(aes(label = q_str, y = value), x = 0.5, min.segment.length = 0, size = 2)+
  geom_hline(yintercept = c(0,-3), size = 0, col = "grey") +
  theme_bw() +
  geom_hline(aes(yintercept = value, col = simu_methode))  +
  scale_color_manual(values = c("Neutral\nsimulation" = "#145db6", "Purifying\nselection\nsimulation" ="#6c2da3", "Perfect\nspike\nsimulation" = "#7c0000", "Best\nObs\nSite"="black"))  +
    scale_x_discrete(position = "top",
                   breaks = total_neutral_sites_with_seq$MGI_site,
                   labels = total_neutral_sites_with_seq$labels_simplified) +
   theme(legend.key=element_blank(),panel.background = element_rect(fill = "white"), panel.grid.major = element_blank(),
        axis.text.x = element_text(colour = "black", size = 6, angle = 45, vjust = 0.3, hjust = 0)) 



fig_total_neutral<-  cowplot::plot_grid(
  cowplot::plot_grid(
    cowplot::plot_grid(NULL, total_p  , rel_heights = c(2.1,16), ncol=1),
    cowplot::plot_grid(total_neutral_sites_with_seq_plot + theme(axis.text.y=element_blank()) + ylab(""), NULL, rel_heights = c(100,1), ncol=1),
    rel_widths = c(3,19)
  ),
  cowplot::plot_grid(NULL, total_neutral_sites_with_pval_plot + theme(legend.position = "none"),
                     rel_widths = c(3,19),  ncol=2),
  ncol=1,rel_heights = c(10,4)
)

 
 
 ggsave(plot=fig_total_neutral, filename = str_c(output_figures,"A_pelican_total_neutral.pdf"), width = 45, height = 20,limitsize = FALSE) 
```



## get observed data pelican results

### best site per ali


```{r}
total_best_pval_per_MGI_ok_summary <- total_best_pval_per_MGI_ok %>%
  mutate(simu_methode = "Best\nObs\nSite")%>%
  group_by(simu_methode, dataset) %>%
  summarise (b_1 = max(log10_best_aagtr_pval),
            b_0.5 = median(log10_best_aagtr_pval),
            b_0.4  = quantile(log10_best_aagtr_pval, 0.4),
            b_0.3  = quantile(log10_best_aagtr_pval, 0.3),
            b_0.2  = quantile(log10_best_aagtr_pval, 0.2),
            b_0.1  = quantile(log10_best_aagtr_pval, 0.1),
            b_0.05  = quantile(log10_best_aagtr_pval, 0.05),
            b_0.01 = quantile(log10_best_aagtr_pval, 0.01),
            b_0 = min(log10_best_aagtr_pval),
            )  %>%  pivot_longer(-c(simu_methode,dataset)) %>%
  mutate(q = as.numeric(str_split_i(name,"_",2)),
         q_str = format_no_trailing_zeros(q),
         name = fct_reorder(name, q) ) %>%
  arrange(-q) 

total_best_pval_per_MGI_ok_summary


murinae_best_pval_per_MGI_ok_summary  <- murinae_best_pval_per_MGI_ok%>%
  mutate(simu_methode = "Best\nObs\nSite")%>%
  group_by(simu_methode, dataset) %>%
  summarise (b_1 = max(log10_best_aagtr_pval),
            b_0.5 = median(log10_best_aagtr_pval),
            b_0.4  = quantile(log10_best_aagtr_pval, 0.4),
            b_0.3  = quantile(log10_best_aagtr_pval, 0.3),
            b_0.2  = quantile(log10_best_aagtr_pval, 0.2),
            b_0.1  = quantile(log10_best_aagtr_pval, 0.1),
            b_0.05  = quantile(log10_best_aagtr_pval, 0.05),
            b_0.01 = quantile(log10_best_aagtr_pval, 0.01),
            b_0 = min(log10_best_aagtr_pval),
            )  %>%  pivot_longer(-c(simu_methode,dataset)) %>%
  mutate(q = as.numeric(str_split_i(name,"_",2)),
         q_str = format_no_trailing_zeros(q),
         name = fct_reorder(name, q) ) %>%
  arrange(-q) 

murinae_best_pval_per_MGI_ok_summary



```


## Threshold from spike


```{r fig.height=8, fig.width=12}
total_spike_violin <- bind_rows( total_pelican_simu_results_summary %>% filter(q2==0, n <= 6),total_best_pval_per_MGI_ok_summary)  %>%
  mutate(facet = "Total") %>%
  ggplot()+
  geom_hline(aes(yintercept = value, col = simu_methode))  +
  scale_color_manual(values = c("Neutral\nsimulation" = "#145db6", "Purifying\nselection\nsimulation" ="#6c2da3", "Perfect\nspike\nsimulation" = "#7c0000", "Best\nObs\nSite"="black")) +
  facet_wrap(~facet + simu_methode, ncol = 4) +
  theme_bw() +
  ggrepel::geom_text_repel(aes(label = q_str, y = value), x = "", min.segment.length = 0, size = 2) +
  geom_point(data = tibble(dataset = "", log10_aagtr_pval = 0),aes(x=dataset, y = log10_aagtr_pval), size = 0, col = "grey") +
  geom_violin(data=total_true_all_sites %>% mutate(dataset = "Total\n  all\nsites"), aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "black") +
  geom_violin(data=total_score_all_tested_sites_df %>% mutate(dataset = "Total\n  only\n non\nconstant"), aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "black") +
      geom_violin(data=total_best_pval_per_MGI_ok %>% mutate(dataset = "Total\n  only\nbest\nsites"), aes(x=dataset, y = log10_best_aagtr_pval), alpha = 0.5, col = "black") +
  geom_violin(data=total_pelican_simu_results %>% mutate(dataset = "Total\nperfect\nspike\n1/21") %>% filter(n_t=="1_21")  %>% dplyr::select(dataset, log10_aagtr_pval) , aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#7c0000") +
  geom_violin(data=total_pelican_simu_results %>% mutate(dataset = "Total\nperfect\nspike\n2/21") %>% filter(n_t=="2_21")  %>% dplyr::select(dataset, log10_aagtr_pval) , aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#7c0000") +
  geom_violin(data=total_pelican_simu_results %>% mutate(dataset = "Total\nperfect\nspike\n3/21") %>% filter(n_t=="3_21")  %>% dplyr::select(dataset, log10_aagtr_pval) , aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#7c0000") 


total_spike_violin


```




```{r fig.height=8, fig.width=12}
murinae_spike_violin <- bind_rows(murinae_pelican_simu_results_summary %>% filter(q2==0, n <= 4),murinae_best_pval_per_MGI_ok_summary)  %>%
  mutate(facet = "Murinae") %>%
  ggplot()+
  geom_hline(aes(yintercept = value, col = simu_methode))  +
  scale_color_manual(values = c("Neutral\nsimulation" = "#145db6", "Purifying\nselection\nsimulation" ="#6c2da3", "Perfect\nspike\nsimulation" = "#7c0000", "Best\nObs\nSite"="black")) +
  facet_wrap(~facet + simu_methode, ncol = 4) +
  theme_bw() +
  ggrepel::geom_text_repel(aes(label = q_str, y = value), x = "", min.segment.length = 0, size = 2) +
  geom_point(data = tibble(dataset = "", log10_aagtr_pval = 0),aes(x=dataset, y = log10_aagtr_pval), size = 0, col = "grey") +
  geom_violin(data=murinae_true_all_sites %>% mutate(dataset = "Murinae\n  all\nsites"), aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "black") +
  geom_violin(data=murinae_score_all_tested_sites_df %>% mutate(dataset = "Murinae\n  only\n non\nconstant"), aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "black") +
    geom_violin(data=murinae_best_pval_per_MGI_ok %>% mutate(dataset = "Murinae\n  only\nbest\nsites"), aes(x=dataset, y = log10_best_aagtr_pval), alpha = 0.5, col = "black") +
  geom_violin(data=murinae_pelican_simu_results %>% mutate(dataset = "Murinae\nperfect\nspike\n1/5") %>% filter(n_t=="1_5")  %>% dplyr::select(dataset, log10_aagtr_pval) , aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#7c0000") +
  geom_violin(data=murinae_pelican_simu_results %>% mutate(dataset = "Murinae\nperfect\nspike\n2/5") %>% filter(n_t=="2_5")  %>% dplyr::select(dataset, log10_aagtr_pval) , aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#7c0000") +
  geom_violin(data=murinae_pelican_simu_results %>% mutate(dataset = "Murinae\nperfect\nspike\n3/5") %>% filter(n_t=="3_5")  %>% dplyr::select(dataset, log10_aagtr_pval) , aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#7c0000") 


murinae_spike_violin


```


```{r, fig.height=10, fig.width=12}

total_relation_min_log10_aagtr_pval_nsites_simu2 <- total_pelican_simu_results %>%
  filter(n_t %in% c("1_21","2_21","3_21","20_21","21_21")) %>%
  ggplot(aes(x=site, y = log10_aagtr_pval)) +
  theme_bw() +
  ggrastr::geom_point_rast(aes(color= combi_str), show.legend = T) +
  geom_line(aes(color= combi_str) )+
  theme(legend.position = "top") +
  facet_wrap(~ "Total" + n_t, ncol = 11, scale = "free_x") +
  geom_hline(data = total_pelican_simu_results_summary %>% filter(q2==0.5), aes(yintercept = value)) +
  geom_hline(data = total_pelican_simu_results_summary %>% filter(q2==0), aes(yintercept = value), linetype = "dashed") +
  geom_hline(data = total_pelican_simu_results_summary %>% filter(q2==1), aes(yintercept = value), linetype = "dashed") +
  theme(legend.position = "none")

total_relation_min_log10_aagtr_pval_nsites_simu2
```



```{r, fig.height=5, fig.width=15}

murinae_relation_min_log10_aagtr_pval_nsites_simu2 <- murinae_pelican_simu_results %>%
  ggplot(aes(x=site, y = log10_aagtr_pval)) +
  theme_bw() +
  ggrastr::geom_point_rast(aes(color= combi_str), show.legend = T) +
  geom_line(aes(color= combi_str) )+
  theme(legend.position = "top") +
  facet_wrap(~"Murinae" +n_t, ncol=5) +
  geom_hline(data = murinae_pelican_simu_results_summary %>% filter(q2==0.5), aes(yintercept = value)) +
  geom_hline(data = murinae_pelican_simu_results_summary %>% filter(q2==0), aes(yintercept = value), linetype = "dashed") +
  geom_hline(data = murinae_pelican_simu_results_summary %>% filter(q2==1), aes(yintercept = value), linetype = "dashed") +
  theme(legend.position = "none")


murinae_relation_min_log10_aagtr_pval_nsites_simu2
```

### murinae

```{r}
murinae_spike_sites<- bind_rows(
  murinae_pelican_simu_results %>% group_by(n_t) %>% slice_max(log10_aagtr_pval,n=4, with_ties = F),
  murinae_pelican_simu_results %>% group_by(n_t) %>% slice_min(log10_aagtr_pval,n=4, with_ties = F)) %>%
  mutate(fasta_fn =  str_c("spike_aa_data/murinae/", alignment,".fasta")) %>%
  mutate(MGI =alignment,
         MGI_site = paste0(alignment,"_",site,"_",dp),
         labels_simplified = paste0(alignment,"_",site),
         MGI_site = fct_reorder2(MGI_site, n_t,log10_aagtr_pval)) %>%
  arrange(MGI_site)



murinae_spike_sites_with_seq <- murinae_spike_sites %>%
  ungroup() %>%
  rowwise() %>%
  mutate(seq_per_sp_df = list(get_seq_from_fasta_fn_pos(fasta_fn, site))) %>%
  unnest(seq_per_sp_df)  %>%
  left_join(species_attributes %>% dplyr::select(leaf_name, Aride), by = c("seqnames" = "leaf_name"))  %>%
  mutate(MGI_site = factor(MGI_site, levels = levels(murinae_spike_sites$MGI_site)),
         MGI = factor(MGI, levels = unique(murinae_spike_sites$MGI)),
         seqnames = factor(seqnames, levels = rev(murinae_tips_order))) 



```


```{r fig.height=10, fig.width=15}

murinae_spike_sites_with_seq_plot <- ggplot(murinae_spike_sites_with_seq, aes(x = MGI_site, y = seqnames, fill = AA, label= AA, col = Aride)) +
 # geom_hline(yintercept = c(0.5,3.5,7.5,6.5,11.5,12.5,15.5,14.5,13.5)) +
  geom_rect(
    data = murinae_limites_y, 
    aes(
      xmin = -Inf,
      xmax = Inf,
      ymin = min_y,
      ymax = max_y,
      fill = Aride
    ),
    alpha = 1, # Transparence pastel
    inherit.aes = FALSE
  ) +
  geom_tile(alpha=0.8,  size = 1, width = 0.8, height=0.8, col = ifelse(murinae_spike_sites_with_seq$Aride =="yes", "grey85","white")) + geom_text(  size = 6, color = "black") +
  scale_color_manual(values = c("darkgreen","darkorange")) +
  theme_bw()  + ggforce::facet_row(~MGI, scale = "free_x", space='free', strip.position="top")  +
  scale_x_discrete(position = "top",
                   breaks = murinae_spike_sites_with_seq$MGI_site,
                   labels = murinae_spike_sites_with_seq$labels_simplified) +
  theme(legend.key=element_blank(),panel.background = element_rect(fill = "white"), panel.grid.major = element_blank(),
        axis.text.x = element_text(colour = "black", size = 10, angle = 45, vjust = 0.3, hjust = 0)) +
  scale_fill_manual(values = c(palette_AA_LETTER,"yes"="grey85"))  +
  theme(legend.position="top") +guides(fill=guide_legend(nrow=1,byrow=TRUE),
                                       color=guide_legend(override.aes=list(fill=NA))) 
murinae_spike_sites_with_seq_plot



murinae_spike_sites_with_pval_plot <-
  bind_rows(murinae_pelican_simu_results_summary %>% filter(q2==0, n <= 5))  %>%
  mutate(facet = "Murinae") %>%
  ggplot()+
  geom_point(data = murinae_spike_sites, aes(x=MGI_site, y = log10_aagtr_pval))   +
  ggrepel::geom_text_repel(aes(label = q_str, y = value), x = 0.5, min.segment.length = 0, size = 2)+
  geom_hline(yintercept = c(0,-4.5), size = 0, col = "grey") +
  theme_bw() +
  geom_hline(aes(yintercept = value, col = simu_methode))  +
  scale_color_manual(values = c("Neutral\nsimulation" = "#145db6", "Purifying\nselection\nsimulation" ="#6c2da3", "Perfect\nspike\nsimulation" = "#7c0000", "Best\nObs\nSite"="black"))  +
    scale_x_discrete(position = "top",
                   breaks = murinae_spike_sites_with_seq$MGI_site,
                   labels = murinae_spike_sites_with_seq$labels_simplified) +
   theme(legend.key=element_blank(),panel.background = element_rect(fill = "white"), panel.grid.major = element_blank(),
        axis.text.x = element_text(colour = "black", size = 6, angle = 45, vjust = 0.3, hjust = 0)) 


murinae_spike_violin_small <- bind_rows(murinae_pelican_simu_results_summary %>% filter(q2==0, n <= 5))  %>%
  mutate(facet = "Murinae") %>%
  ggplot()+
  geom_hline(aes(yintercept = value, col = simu_methode))  +
  scale_color_manual(values = c("Neutral\nsimulation" = "#145db6", "Purifying\nselection\nsimulation" ="#6c2da3", "Perfect\nspike\nsimulation" = "#7c0000", "Best\nObs\nSite"="black")) +
  theme_bw() +
  ggrepel::geom_text_repel(aes(label = q_str, y = value), x = 0.5, min.segment.length = 0, size = 2) +
  geom_hline(yintercept = c(0,-4.5), size = 0, col = "grey") +
  geom_violin(data=murinae_true_all_sites %>% mutate(dataset = "Murinae\n  all\nsites"), aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "black") +
  geom_violin(data=murinae_score_all_tested_sites_df %>% mutate(dataset = "Murinae\n  only\n non\nconstant"), aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "black") +
    geom_violin(data=murinae_best_pval_per_MGI_ok %>% mutate(dataset = "Murinae\n  only\nbest\nsites"), aes(x=dataset, y = log10_best_aagtr_pval), alpha = 0.5, col = "black") +
  geom_violin(data=murinae_pelican_simu_results %>% mutate(dataset = "Murinae\nperfect\nspike\n1/5") %>% filter(n_t=="1_5")  %>% dplyr::select(dataset, log10_aagtr_pval) , aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#7c0000") +
  geom_violin(data=murinae_pelican_simu_results %>% mutate(dataset = "Murinae\nperfect\nspike\n2/5") %>% filter(n_t=="2_5")  %>% dplyr::select(dataset, log10_aagtr_pval) , aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#7c0000") +
  geom_violin(data=murinae_pelican_simu_results %>% mutate(dataset = "Murinae\nperfect\nspike\n3/5") %>% filter(n_t=="3_5")  %>% dplyr::select(dataset, log10_aagtr_pval) , aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#7c0000") +
    geom_violin(data=murinae_pelican_simu_results %>% mutate(dataset = "Murinae\nperfect\nspike\n4/5") %>% filter(n_t=="4_5")  %>% dplyr::select(dataset, log10_aagtr_pval) , aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#7c0000") +
    geom_violin(data=murinae_pelican_simu_results %>% mutate(dataset = "Murinae\nperfect\nspike\n5/5") %>% filter(n_t=="5_5")  %>% dplyr::select(dataset, log10_aagtr_pval) , aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#7c0000") +
  theme(legend.position="")  +
  scale_x_discrete(position = "top")


murinae_spike_violin_small



fig_murinae_spike<-  cowplot::plot_grid(
  murinae_relation_min_log10_aagtr_pval_nsites_simu2,
  cowplot::plot_grid(
    cowplot::plot_grid(NULL, murinae_p  , rel_heights = c(3,16), ncol=1),
    cowplot::plot_grid(murinae_spike_sites_with_seq_plot + theme(axis.text.y=element_blank()) + ylab(""), NULL, rel_heights = c(100,3), ncol=1),
    rel_widths = c(3,19)
  ),
  cowplot::plot_grid(murinae_spike_violin_small, murinae_spike_sites_with_pval_plot + theme(legend.position = "none"),
                     rel_widths = c(3,19),  ncol=2),
  ncol=1,rel_heights = c(4,10,4)
)

 
 
 ggsave(plot=fig_murinae_spike, filename = str_c(output_figures,"A_pelican_murinae_spike.pdf"), width = 35, height = 15)

  
```


### total

```{r}

n_t_total_ok <-  paste0(c(1,2,3,4,5,6,20,21),"_",total_nb_transition)
total_spike_sites<- bind_rows(
  total_pelican_simu_results %>% filter(n_t %in% n_t_total_ok)  %>% group_by(n_t)  %>% slice_max(log10_aagtr_pval,n=3, with_ties = F),
  total_pelican_simu_results %>% group_by(n_t) %>% filter(n_t %in% n_t_total_ok)  %>% slice_min(log10_aagtr_pval,n=3, with_ties = F)) %>%
  mutate(fasta_fn =  str_c("spike_aa_data//total/", alignment,".fasta")) %>%
  mutate(is.file = file.exists(fasta_fn),
         MGI =alignment,
         MGI_site = paste0(alignment,"_",site,"_",dp),
         labels_simplified = paste0(alignment,"_",site),
         MGI_site = fct_reorder2(MGI_site, n_t,log10_aagtr_pval)) %>%
  arrange(MGI_site)

all(total_spike_sites$is.file)

total_spike_sites_with_seq <- total_spike_sites %>%
  ungroup() %>%
  rowwise() %>%
  mutate(seq_per_sp_df = list(get_seq_from_fasta_fn_pos(fasta_fn, site))) %>%
  unnest(seq_per_sp_df)  %>%
  left_join(species_attributes %>% dplyr::select(leaf_name, Aride), by = c("seqnames" = "leaf_name"))  %>%
  mutate(MGI_site = factor(MGI_site, levels = levels(total_spike_sites$MGI_site)),
         MGI = factor(MGI, levels = unique(total_spike_sites$MGI)),
         seqnames = factor(seqnames, levels = rev(total_tips_order))) 





```


```{r fig.height=10, fig.width=15}

total_spike_sites_with_seq_plot <- ggplot(total_spike_sites_with_seq, aes(x = MGI_site, y = seqnames, fill = AA, label= AA, col = Aride)) +
 # geom_hline(yintercept = c(0.5,3.5,7.5,6.5,11.5,12.5,15.5,14.5,13.5)) +
  geom_rect(
    data = total_limites_y, 
    aes(
      xmin = -Inf,
      xmax = Inf,
      ymin = min_y,
      ymax = max_y,
      fill = Aride
    ),
    alpha = 1, # Transparence pastel
    inherit.aes = FALSE
  ) +
  geom_tile(alpha=0.8,  size = 1, width = 0.8, height=0.8, col = ifelse(total_spike_sites_with_seq$Aride =="yes", "grey85","white")) + geom_text(  size = 6, color = "black") +
  scale_color_manual(values = c("darkgreen","darkorange")) +
  theme_bw()  + ggforce::facet_row(~MGI, scale = "free_x", space='free', strip.position="top")  +
  scale_x_discrete(position = "top",
                   breaks = total_spike_sites_with_seq$MGI_site,
                   labels = total_spike_sites_with_seq$labels_simplified) +
  theme(legend.key=element_blank(),panel.background = element_rect(fill = "white"), panel.grid.major = element_blank(),
        axis.text.x = element_text(colour = "black", size = 10, angle = 45, vjust = 0.3, hjust = 0)) +
  scale_fill_manual(values = c(palette_AA_LETTER,"yes"="grey85"))  +
  theme(legend.position="top") +guides(fill=guide_legend(nrow=1,byrow=TRUE),
                                       color=guide_legend(override.aes=list(fill=NA))) 
total_spike_sites_with_seq_plot



total_spike_sites_with_pval_plot <-
  bind_rows(total_pelican_simu_results_summary %>% filter(q2==0, n %in% c(1,2,3,4,5,6,20,21)))  %>%
  mutate(facet = "Total") %>%
  ggplot()+
  geom_point(data = total_spike_sites, aes(x=MGI_site, y = log10_aagtr_pval))   +
  ggrepel::geom_text_repel(aes(label = q_str, y = value), x = 0.5, min.segment.length = 0, size = 2)+
  geom_hline(yintercept = c(0,-4.5), size = 0, col = "grey") +
  theme_bw() +
  geom_hline(aes(yintercept = value, col = simu_methode))  +
  scale_color_manual(values = c("Neutral\nsimulation" = "#145db6", "Purifying\nselection\nsimulation" ="#6c2da3", "Perfect\nspike\nsimulation" = "#7c0000", "Best\nObs\nSite"="black"))  +
    scale_x_discrete(position = "top",
                   breaks = total_spike_sites_with_seq$MGI_site,
                   labels = total_spike_sites_with_seq$labels_simplified) +
   theme(legend.key=element_blank(),panel.background = element_rect(fill = "white"), panel.grid.major = element_blank(),
        axis.text.x = element_text(colour = "black", size = 3, angle = 45, vjust = 0.3, hjust = 0)) 


total_spike_violin_small <- bind_rows(total_pelican_simu_results_summary %>% filter(q2==0, n %in% c(1,2,3,4,20,21)))  %>%
  mutate(facet = "Total") %>%
  ggplot()+
  geom_hline(aes(yintercept = value, col = simu_methode))  +
  scale_color_manual(values = c("Neutral\nsimulation" = "#145db6", "Purifying\nselection\nsimulation" ="#6c2da3", "Perfect\nspike\nsimulation" = "#7c0000", "Best\nObs\nSite"="black")) +
  theme_bw() +
  ggrepel::geom_text_repel(aes(label = q_str, y = value), x = 0.5, min.segment.length = 0, size = 2) +
  geom_hline(yintercept = c(0,-4.5), size = 0, col = "grey") +
  geom_violin(data=total_true_all_sites %>% mutate(dataset = "Total\n  all\nsites"), aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "black") +
  geom_violin(data=total_score_all_tested_sites_df %>% mutate(dataset = "Total\n  only\n non\nconstant"), aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "black") +
    geom_violin(data=total_best_pval_per_MGI_ok %>% mutate(dataset = "Total\n  only\nbest\nsites"), aes(x=dataset, y = log10_best_aagtr_pval), alpha = 0.5, col = "black") +
  geom_violin(data=total_pelican_simu_results %>% mutate(dataset = "Total\nperfect\nspike\n 1/21") %>% filter(n_t=="1_21")  %>% dplyr::select(dataset, log10_aagtr_pval) , aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#7c0000") +
  geom_violin(data=total_pelican_simu_results %>% mutate(dataset = "Total\nperfect\nspike\n 2/21") %>% filter(n_t=="2_21")  %>% dplyr::select(dataset, log10_aagtr_pval) , aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#7c0000") +
  geom_violin(data=total_pelican_simu_results %>% mutate(dataset = "Total\nperfect\nspike\n 3/21") %>% filter(n_t=="3_21")  %>% dplyr::select(dataset, log10_aagtr_pval) , aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#7c0000") +
  geom_violin(data=total_pelican_simu_results %>% mutate(dataset = "Total\nperfect\nspike\n 4/21") %>% filter(n_t=="4_21")  %>% dplyr::select(dataset, log10_aagtr_pval) , aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#7c0000") +
    geom_violin(data=total_pelican_simu_results %>% mutate(dataset = "Total\nperfect\nspike\n20/21") %>% filter(n_t=="20_21")  %>% dplyr::select(dataset, log10_aagtr_pval) , aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#7c0000") +
    geom_violin(data=total_pelican_simu_results %>% mutate(dataset = "Total\nperfect\nspike\n21/21") %>% filter(n_t=="21_21")  %>% dplyr::select(dataset, log10_aagtr_pval) , aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#7c0000") +
  theme(legend.position="")  +
  scale_x_discrete(position = "top")


total_spike_violin_small



fig_total_spike<-  cowplot::plot_grid(
  total_relation_min_log10_aagtr_pval_nsites_simu2,
  cowplot::plot_grid(
    cowplot::plot_grid(NULL, total_p  , rel_heights = c(3,16), ncol=1),
    cowplot::plot_grid(total_spike_sites_with_seq_plot + theme(axis.text.y=element_blank()) + ylab(""), NULL, rel_heights = c(100,1), ncol=1),
    rel_widths = c(3,19)
  ),
  cowplot::plot_grid(total_spike_violin_small, total_spike_sites_with_pval_plot + theme(legend.position = "none"),
                     rel_widths = c(3,19),  ncol=2),
  ncol=1,rel_heights = c(6,10,4)
)

 
 
 ggsave(plot=fig_total_spike, filename = str_c(output_figures,"A_pelican_total_spike.pdf"), width = 45, height = 30,limitsize = FALSE)

  
```


```{r fig.height=8, fig.width=14}

figA <- cowplot::plot_grid(
  total_relation_min_log10_aagtr_pval_nsites_simu2,total_violin_spike,
  murinae_relation_min_log10_aagtr_pval_nsites_simu2, murinae_violin_spike,
  rel_heights = c(3,1), ncol=2
)

ggsave(plot=figA, filename = str_c(output_figures,"A_pelican_spike.pdf"), width = 14, height = 10)

```


```{r fig.height=8, fig.width=10}
total_score_all_tested_sites_df <- total_score_all_tested_sites_df  %>%
  left_join(total_genes_attributes_ok %>% dplyr::select(alignment,MGI)) %>%
  mutate(MGI_site = paste0(MGI,"_",site,"_",dp),
         labels_simplified = paste0(MGI,"_",site))

murinae_score_all_tested_sites_df <- murinae_score_all_tested_sites_df  %>%
  left_join(murinae_genes_attributes_ok %>% dplyr::select(alignment,MGI)) %>%
  mutate(MGI_site = paste0(MGI,"_",site,"_",dp),
         labels_simplified = paste0(MGI,"_",site))

combined_score_all_tested_sites_df <- total_score_all_tested_sites_df  %>%
  inner_join(murinae_score_all_tested_sites_df, by = c("alignment","site","tag","labels_simplified"), suffix = c(".total",".murinae"))

combined_score_all_tested_sites_df_top <- combined_score_all_tested_sites_df %>%
  filter(log10_aagtr_pval.total < total_pelican_neutral_all_tested_results_summary %>% filter(q == 0.001) %>% pull(value) &
          log10_aagtr_pval.murinae < murinae_pelican_neutral_all_tested_results_summary %>% filter(q == 0.001) %>% pull(value) )

dim(combined_score_all_tested_sites_df)
dim(combined_score_all_tested_sites_df_top)
dim(total_score_all_tested_sites_df)
dim(murinae_score_all_tested_sites_df)

combined_score_all_tested_sites_df_top %>% filter(MGI.total == "Fam20a") %>% arrange(aagtr_pval.total) %>% dplyr::select(labels_simplified,  MGI_site.total, log10_aagtr_pval.murinae, log10_aagtr_pval.total,log10_aagtr_pval_cut_th.total, log10_aagtr_pval_cut_th.murinae ) 

```


```{r fig.height=6, fig.width=12}

p <- combined_score_all_tested_sites_df %>%
  ggplot(aes(y=log10(aagtr_pval.murinae), x = log10(aagtr_pval.total))) +
  ggrastr::geom_point_rast(color = "grey", alpha=0.8, size = 0.5, raster.dpi = 300) +
  geom_density_2d( color = "black") +
  theme_bw() + scale_x_continuous(position = "bottom")  +
   geom_hline(data = murinae_pelican_neutral_all_tested_results_summary %>% filter(q == 0.001), aes(yintercept = value), color = "#145db6" ) +
  
  geom_vline(data = total_pelican_neutral_all_tested_results_summary %>% filter(q == 0.001), aes(xintercept = value), color = "#145db6" ) +
  
  
  
  geom_hline(data = murinae_pelican_simu_results_summary%>% filter(q2 == 0), aes(yintercept = value), linetype = "dashed", color = "darkred" , linetype = "dashed" ) +
  
  geom_vline(data = total_pelican_simu_results_summary %>% filter(q2 == 0, n <=6), aes(xintercept = value), linetype = "dashed", color = "darkred" , linetype = "dashed" ) +

  theme(legend.position = "none") +
  geom_text(data = murinae_pelican_simu_results_summary%>% filter(q2 == 0), aes(y = value, x = 0.1, label = str_replace(n_t,"_","/")), vjust = -1.5 , color = "grey15", size = 4 )+
  geom_text(data = total_pelican_simu_results_summary%>% filter(q2 == 0, n <=6), aes(y = 0.2, x = value, label = str_replace(n_t,"_","/")), hjust = -1.5 , color = "grey15" , size = 4) +
  xlab("log10(pvalue.total)") + ylab("log10(pvalue.murinae)")


 p_with_label <- p +
  ggrepel::geom_text_repel(data = combined_score_all_tested_sites_df_top %>% slice_min(log10_aagtr_pval.murinae + log10_aagtr_pval.total, n = 10)  ,
                                ylim = c(-4.3,-2), xlim = c(-4,-1.8),
                               aes(label = labels_simplified ),
                               force = 20,  segment.alpha = 0.4, size = 3,
                               box.padding = 0.8, 
                               min.segment.length = 0,
                               seed = 22, segment.col ="grey60" , col ="black" , 
                               max.overlaps = Inf)  +
  ggrepel::geom_text_repel(data = combined_score_all_tested_sites_df %>% slice_min(log10_aagtr_pval.total, n = 5) ,
                                ylim = c(-1.5,0.3), xlim = c(-4,-3),
                               aes(label = labels_simplified ),
                               force = 20,  segment.alpha = 0.4, size = 3,
                               box.padding = 0.8, 
                               min.segment.length = 0,
                               seed = 22, segment.col ="grey60" , col ="black" , 
                               max.overlaps = Inf) +
  ggrepel::geom_text_repel(data = combined_score_all_tested_sites_df %>% slice_min(log10_aagtr_pval.murinae, n = 5) ,
                                ylim = c(-4.3,-2.5), xlim = c(-1.8,0),
                               aes(label = labels_simplified ),
                               force = 20,  segment.alpha = 0.4, size = 3,
                               box.padding = 0.8, 
                               min.segment.length = 0,
                               seed = 22, segment.col ="grey60" , col ="black" , 
                               max.overlaps = Inf) 

p_with_label


ggMarginal(p_with_label, type = "violin", fill = "grey", colour = "black", alpha = 0.5,size = 10)


```

```{r fig.height=8, fig.width=14}


figC <-p_with_label

ggsave(plot=figC, filename = str_c(output_figures,"C_pelican_crossing_total_murinae.pdf"), width = 12, height = 10, unit = "cm")
ggsave(plot=p, filename = str_c(output_figures,"C_pelican_crossing_total_murinae_without_names.pdf"), width = 10/1.5, height = 6/1.5)

```


```{r}
total_bornes <- c(total_pelican_neutral_all_tested_results_summary$value, log10(0.0000000001))
total_label_bornes <- rev(c(paste0("inf_", total_pelican_neutral_all_tested_results_summary$q_str)))


total_neutral_bornes <- c(0,total_pelican_neutral_all_tested_results_summary %>% filter(q == 0.001) %>% pull(value), log10(0.0000000001))
total_neutral_label_bornes <- rev(c("sup_th","inf_th"))


murinae_bornes <- c(murinae_pelican_neutral_all_tested_results_summary$value, log10(0.0000000001))
murinae_label_bornes <- rev(c(paste0("inf_", murinae_pelican_neutral_all_tested_results_summary$q_str)))

murinae_neutral_bornes <- c(0,murinae_pelican_neutral_all_tested_results_summary %>% filter(q == 0.001) %>% pull(value), log10(0.0000000001))
murinae_neutral_label_bornes <- rev(c("sup_th","inf_th"))


total_true_all_sites  <- total_true_all_sites %>%
  mutate(
         log10_aagtr_pval_cut = cut(log10_aagtr_pval, breaks = total_bornes, include.lowest = TRUE, right = FALSE, label = total_label_bornes),
         log10_aagtr_pval_cut_th = cut(log10_aagtr_pval, breaks = total_neutral_bornes, include.lowest = TRUE, right = FALSE, label = total_neutral_label_bornes)
  )

total_score_all_tested_sites_df  <- total_score_all_tested_sites_df %>%
  mutate(
         log10_aagtr_pval_cut = cut(log10_aagtr_pval, breaks = total_bornes, include.lowest = TRUE, right = FALSE, label = total_label_bornes),
         log10_aagtr_pval_cut_th = cut(log10_aagtr_pval, breaks = total_neutral_bornes, include.lowest = TRUE, right = FALSE, label = total_neutral_label_bornes)
  )
total_pelican_simu_results  <- total_pelican_simu_results %>%
  mutate(
         log10_aagtr_pval_cut = cut(log10_aagtr_pval, breaks = total_bornes, include.lowest = TRUE, right = FALSE, label = total_label_bornes),
         log10_aagtr_pval_cut_th = cut(log10_aagtr_pval, breaks = total_neutral_bornes, include.lowest = TRUE, right = FALSE, label = total_neutral_label_bornes)
  )

total_best_pval_per_MGI_ok  <- total_best_pval_per_MGI_ok %>%
  mutate(
         log10_aagtr_pval_cut = cut(log10_best_aagtr_pval, breaks = total_bornes, include.lowest = TRUE, right = FALSE, label = total_label_bornes),
         log10_aagtr_pval_cut_th = cut(log10_best_aagtr_pval, breaks = total_neutral_bornes, include.lowest = TRUE, right = FALSE, label = total_neutral_label_bornes)
  )


murinae_true_all_sites  <- murinae_true_all_sites %>%
  mutate(
         log10_aagtr_pval_cut = cut(log10_aagtr_pval, breaks = murinae_bornes, include.lowest = TRUE, right = FALSE, label = murinae_label_bornes),
         log10_aagtr_pval_cut_th = cut(log10_aagtr_pval, breaks = murinae_neutral_bornes, include.lowest = TRUE, right = FALSE, label = murinae_neutral_label_bornes)
  )

murinae_score_all_tested_sites_df  <- murinae_score_all_tested_sites_df %>%
  mutate(
         log10_aagtr_pval_cut = cut(log10_aagtr_pval, breaks = murinae_bornes, include.lowest = TRUE, right = FALSE, label = murinae_label_bornes),
         log10_aagtr_pval_cut_th = cut(log10_aagtr_pval, breaks = murinae_neutral_bornes, include.lowest = TRUE, right = FALSE, label = murinae_neutral_label_bornes)
  )

murinae_pelican_simu_results  <- murinae_pelican_simu_results %>%
  mutate(
         log10_aagtr_pval_cut = cut(log10_aagtr_pval, breaks = murinae_bornes, include.lowest = TRUE, right = FALSE, label = murinae_label_bornes),
         log10_aagtr_pval_cut_th = cut(log10_aagtr_pval, breaks = murinae_neutral_bornes, include.lowest = TRUE, right = FALSE, label = murinae_neutral_label_bornes)
  )

murinae_best_pval_per_MGI_ok  <- murinae_best_pval_per_MGI_ok %>%
  mutate(
         log10_aagtr_pval_cut = cut(log10_best_aagtr_pval, breaks = murinae_bornes, include.lowest = TRUE, right = FALSE, label = murinae_label_bornes),
         log10_aagtr_pval_cut_th = cut(log10_best_aagtr_pval, breaks = murinae_neutral_bornes, include.lowest = TRUE, right = FALSE, label = murinae_neutral_label_bornes)
  )



# Comptage du nombre de valeurs dans chaque intervalle

table(total_score_all_tested_sites_df$log10_aagtr_pval_cut)


#inf_21_21 inf_20_21 inf_19_21 inf_18_21 inf_17_21 inf_16_21 inf_15_21 inf_14_21 inf_13_21 inf_12_21 inf_11_21 inf_10_21  inf_9_21  inf_8_21  inf_7_21  inf_6_21  inf_5_21  inf_4_21  inf_3_21  inf_2_21  inf_1_21     inf_0 
#        0         0         0         0         0         0         0         0         0         0         0         0         1         5         6        58        64       459      2048     10430     86851    288119 


table(total_score_all_tested_sites_df$log10_aagtr_pval_cut_th)
#inf_th sup_th 
#  1031 387010 
  
  
table(murinae_score_all_tested_sites_df$log10_aagtr_pval_cut)

 
#inf_5_5 inf_4_5 inf_3_5 inf_2_5 inf_1_5   inf_0 
#     12       7      60     777   12815   81405 

table(murinae_score_all_tested_sites_df$log10_aagtr_pval_cut_th)

#inf_th sup_th 
#  1237  93839 

table(murinae_true_all_sites$log10_aagtr_pval_cut_th)
#inf_th sup_th 
#  1237 959170 

table(total_true_all_sites$log10_aagtr_pval_cut_th)
#inf_th  sup_th 
#  1031  1246450 


```



```{r}
murinae_nb_site_per_cat_th <- murinae_true_all_sites %>%
  group_by(log10_aagtr_pval_cut_th) %>%
  summarise(nb_site = n(),
            percent_site = nb_site / nrow(murinae_true_all_sites) * 100 ) %>%
  mutate(cumul_nb_site = cumsum(nb_site),
         cumul_percent_site = cumsum(percent_site)) %>%
  mutate(dataset = "Murinae")


murinae_plot_nb_site_per_cat_th <- murinae_nb_site_per_cat_th %>%
  ggplot(aes(x= log10_aagtr_pval_cut_th, y =  cumul_percent_site)) +
  geom_bar(stat = "identity") +
  theme_bw() +
  geom_text(aes(label = paste0(cumul_nb_site, " (", round(cumul_percent_site,3), "%)")), y=50)+
  facet_wrap(~dataset) +
  geom_vline(xintercept = 1.5, color = "blue" , linetype = "dashed")  +
  coord_flip()

murinae_nb_gene_per_cat_th <- murinae_true_all_sites %>%
  group_by(alignment) %>%
  slice_min(log10_aagtr_pval, n = 1, with_ties = F) %>%
  group_by(log10_aagtr_pval_cut_th) %>%
  summarise(nb_gene = n(),
            percent_gene = nb_gene / length(unique(murinae_true_all_sites$alignment)) * 100 ) %>%
  mutate(cumul_nb_gene = cumsum(nb_gene),
         cumul_percent_gene = cumsum(percent_gene))  %>%
  mutate(dataset = "Murinae")


murinae_plot_nb_gene_per_cat_th <- murinae_nb_gene_per_cat_th %>%
  ggplot(aes(x= log10_aagtr_pval_cut_th, y =  cumul_percent_gene)) +
  geom_bar(stat = "identity") +
  theme_bw() +
  geom_text(aes(label = paste0(cumul_nb_gene, " (", round(cumul_percent_gene,3), "%)")), y=65)+
  facet_wrap(~dataset) +
  geom_vline(xintercept = 1.5, color = "blue" , linetype = "dashed") +
  ggtitle("Nb gene with at least 1 site (cumul sum)") +
  coord_flip()

murinae_plot_nb_site_per_cat_th
murinae_plot_nb_gene_per_cat_th
```

```{r}
murinae_nb_site_per_cat <- murinae_true_all_sites %>%
  group_by(log10_aagtr_pval_cut) %>%
  summarise(nb_site = n(),
            percent_site = nb_site / nrow(murinae_true_all_sites) * 100 ) %>%
  mutate(cumul_nb_site = cumsum(nb_site),
         cumul_percent_site = cumsum(percent_site)) %>%
  mutate(dataset = "Murinae")


murinae_plot_nb_site_per_cat <- murinae_nb_site_per_cat %>%
  ggplot(aes(x= log10_aagtr_pval_cut, y =  cumul_percent_site)) +
  geom_bar(stat = "identity") +
  theme_bw() +
  geom_text(aes(label = paste0(cumul_nb_site, " (", round(cumul_percent_site,3), "%)")), y=50)+
  facet_wrap(~dataset) +
  geom_vline(xintercept = 3.5, color = "darkred" , linetype = "dashed")  +
  coord_flip()

murinae_nb_gene_per_cat <- murinae_true_all_sites %>%
  group_by(alignment) %>%
  slice_min(log10_aagtr_pval, n = 1, with_ties = F) %>%
  group_by(log10_aagtr_pval_cut) %>%
  summarise(nb_gene = n(),
            percent_gene = nb_gene / length(unique(murinae_true_all_sites$alignment)) * 100 ) %>%
  mutate(cumul_nb_gene = cumsum(nb_gene),
         cumul_percent_gene = cumsum(percent_gene))  %>%
  mutate(dataset = "Murinae")


murinae_plot_nb_gene_per_cat <- murinae_nb_gene_per_cat %>%
  ggplot(aes(x= log10_aagtr_pval_cut, y =  cumul_percent_gene)) +
  geom_bar(stat = "identity") +
  theme_bw() +
  geom_text(aes(label = paste0(cumul_nb_gene, " (", round(cumul_percent_gene,3), "%)")), y=50)+
  facet_wrap(~dataset) +
  geom_vline(xintercept = 3.5, color = "darkred" , linetype = "dashed") +
  ggtitle("Nb gene with at least 1 site (cumul sum)") +
  coord_flip()

murinae_plot_nb_site_per_cat 
murinae_plot_nb_gene_per_cat
```


```{r fig.width=12, fig.height=5}
total_nb_site_per_cat <- total_true_all_sites %>%
  group_by(log10_aagtr_pval_cut) %>%
  summarise(nb_site = n(),
            percent_site = nb_site / nrow(total_true_all_sites) * 100 ) %>%
  mutate(cumul_nb_site = cumsum(nb_site),
         cumul_percent_site = cumsum(percent_site))  %>%
  mutate(dataset = "Total")

total_plot_nb_site_per_cat <- total_nb_site_per_cat %>%
  ggplot(aes(x= log10_aagtr_pval_cut, y =  cumul_percent_site)) +
  geom_bar(stat = "identity") +
  theme_bw() +
  geom_text(aes(label = paste0(cumul_nb_site, " (", round(cumul_percent_site,3), "%)")), y=50) +
  facet_wrap(~dataset) +
  geom_vline(xintercept = 6.5, color = "darkred" , linetype = "dashed")  +
  coord_flip()


total_nb_gene_per_cat <- total_true_all_sites %>%
  group_by(alignment) %>%
  slice_min(log10_aagtr_pval, n = 1, with_ties = F) %>%
  group_by(log10_aagtr_pval_cut) %>%
  summarise(nb_gene = n(),
            percent_gene = nb_gene / length(unique(total_true_all_sites$alignment)) * 100 ) %>%
  mutate(cumul_nb_gene = cumsum(nb_gene),
         cumul_percent_gene = cumsum(percent_gene))  %>%
  mutate(dataset = "Total")




total_plot_nb_gene_per_cat <- total_nb_gene_per_cat %>%
  ggplot(aes(x= log10_aagtr_pval_cut, y =  cumul_percent_gene)) +
  geom_bar(stat = "identity") +
  theme_bw() +
  geom_text(aes(label = paste0(cumul_nb_gene, " (", round(cumul_percent_gene,3), "%)")), y=50) +
  facet_wrap(~dataset) +
  geom_vline(xintercept = 6.5, color = "darkred" , linetype = "dashed") +
  ggtitle("Nb gene with at least 1 site (cumul sum)")   +
  coord_flip()

total_plot_nb_site_per_cat 
total_plot_nb_gene_per_cat
```



```{r}
total_nb_site_per_cat_th <- total_true_all_sites %>%
  group_by(log10_aagtr_pval_cut_th) %>%
  summarise(nb_site = n(),
            percent_site = nb_site / nrow(total_true_all_sites) * 100 ) %>%
  mutate(cumul_nb_site = cumsum(nb_site),
         cumul_percent_site = cumsum(percent_site)) %>%
  mutate(dataset = "Total")


total_plot_nb_site_per_cat_th <- total_nb_site_per_cat_th %>%
  ggplot(aes(x= log10_aagtr_pval_cut_th, y =  cumul_percent_site)) +
  geom_bar(stat = "identity") +
  theme_bw() +
  geom_text(aes(label = paste0(cumul_nb_site, " (", round(cumul_percent_site,3), "%)")), y=50)+
  facet_wrap(~dataset) +
  geom_vline(xintercept = 1.5, color = "blue" , linetype = "dashed")  +
  coord_flip()

total_nb_gene_per_cat_th <- total_true_all_sites %>%
  group_by(alignment) %>%
  slice_min(log10_aagtr_pval, n = 1, with_ties = F) %>%
  group_by(log10_aagtr_pval_cut_th) %>%
  summarise(nb_gene = n(),
            percent_gene = nb_gene / length(unique(total_true_all_sites$alignment)) * 100 ) %>%
  mutate(cumul_nb_gene = cumsum(nb_gene),
         cumul_percent_gene = cumsum(percent_gene))  %>%
  mutate(dataset = "Total")


total_plot_nb_gene_per_cat_th <- total_nb_gene_per_cat_th %>%
  ggplot(aes(x= log10_aagtr_pval_cut_th, y =  cumul_percent_gene)) +
  geom_bar(stat = "identity") +
  theme_bw() +
  geom_text(aes(label = paste0(cumul_nb_gene, " (", round(cumul_percent_gene,3), "%)")), y=65)+
  facet_wrap(~dataset) +
  geom_vline(xintercept = 1.5, color = "blue" , linetype = "dashed") +
  ggtitle("Nb gene with at least 1 site (cumul sum)") +
  coord_flip()

total_plot_nb_site_per_cat_th
total_plot_nb_gene_per_cat_th
```


```{r, fig.height=4}
murinae_nb_sites_per_ali <- 
  murinae_score_all_tested_sites_df %>%
  bind_rows(murinae_pelican_simu_results) %>%
  dplyr::select(alignment,aagtr_pval,log10_aagtr_pval,log10_aagtr_pval_cut) %>%
  filter(log10_aagtr_pval_cut %in% (paste0("inf_",murinae_pelican_neutral_all_tested_results_summary %>% filter(q==0.001) %>% pull(q_str)))) %>%
  left_join(murinae_genes_attributes_ok_filtered %>% dplyr::select(MGI,alignment,padj_DEseq2,fasta_fn,log2FoldChange_DEseq2)) %>%
  group_by(MGI,alignment,fasta_fn,padj_DEseq2) %>%
  summarise(n_site_inf_th=n()) %>%
  arrange(-n_site_inf_th) %>%
  ungroup()%>%
  mutate(
         is_spike = grepl("simu_",alignment),
         is_de = padj_DEseq2 < 0.1) %>%
  mutate(n_site_inf_th = ifelse(n_site_inf_th> 40, 40, n_site_inf_th)) %>%
  mutate(dataset = "Murinae")


murinae_plot_nb_sites_per_ali <- murinae_nb_sites_per_ali %>% filter(!is_spike) %>%
  group_by(dataset,n_site_inf_th) %>%
  summarise(nb_gene = n()) %>%
  ggplot(aes(x= n_site_inf_th, y = nb_gene)) +
  geom_bar(stat = "identity") +
  theme_bw() +
  geom_text(aes(label = nb_gene), vjust = -0.2) +
  facet_wrap(~dataset) +
  geom_vline(xintercept = 0, color = "darkred" , linetype = "dashed") 

murinae_plot_nb_sites_per_ali
```



```{r, fig.height=4}
murinae_nb_sites_per_ali_th <- 
  murinae_score_all_tested_sites_df %>%
  dplyr::select(alignment,aagtr_pval,log10_aagtr_pval,log10_aagtr_pval_cut,log10_aagtr_pval_cut_th) %>%
  filter(log10_aagtr_pval_cut_th %in% ("inf_th")) %>%
  left_join(murinae_genes_attributes_ok_filtered %>% dplyr::select(MGI,alignment,padj_DEseq2,fasta_fn,log2FoldChange_DEseq2)) %>%
  group_by(MGI,alignment,fasta_fn,padj_DEseq2) %>%
  summarise(n_site_inf_th=n()) %>%
  arrange(-n_site_inf_th) %>%
  ungroup()%>%
  mutate(
         is_de = padj_DEseq2 < 0.1) %>%
  mutate(n_site_inf_th = ifelse(n_site_inf_th> 40, 40, n_site_inf_th)) %>%
  mutate(dataset = "Murinae")


murinae_plot_nb_sites_per_ali_th <- murinae_nb_sites_per_ali_th %>% 
  group_by(dataset,n_site_inf_th) %>%
  summarise(nb_gene = n()) %>%
  ggplot(aes(x= n_site_inf_th, y = nb_gene)) +
  geom_bar(stat = "identity") +
  theme_bw() +
  geom_text(aes(label = nb_gene), vjust = -0.2) +
  facet_wrap(~dataset) +
  geom_vline(xintercept = 0, color = "blue" , linetype = "dashed") 

murinae_plot_nb_sites_per_ali_th
```


```{r, fig.height=4}
total_nb_sites_per_ali <- 
  total_score_all_tested_sites_df %>%
  bind_rows(total_pelican_simu_results) %>%
  dplyr::select(alignment,aagtr_pval,log10_aagtr_pval,log10_aagtr_pval_cut) %>%
  filter(log10_aagtr_pval_cut %in%(paste0("inf_",total_pelican_neutral_all_tested_results_summary %>% filter(q==0.001) %>% pull(q_str)))) %>%
  left_join(total_genes_attributes_ok_filtered %>% dplyr::select(MGI,alignment,padj_DEseq2,fasta_fn,log2FoldChange_DEseq2)) %>%
  group_by(MGI,alignment,fasta_fn,padj_DEseq2) %>%
  summarise(n_site_inf_th=n()) %>%
  arrange(-n_site_inf_th) %>%
  ungroup()%>%
  mutate(
         is_spike = grepl("simu_",alignment),
         is_de = padj_DEseq2 < 0.1) %>%
  mutate(n_site_inf_th = ifelse(n_site_inf_th> 40, 40, n_site_inf_th)) %>%
  mutate(dataset = "Total")


total_plot_nb_sites_per_ali <- total_nb_sites_per_ali %>% filter(!is_spike) %>%
  group_by(dataset,n_site_inf_th) %>%
  summarise(nb_gene = n()) %>%
  ggplot(aes(x= n_site_inf_th, y = nb_gene)) +
  geom_bar(stat = "identity") +
  theme_bw() +
  geom_text(aes(label = nb_gene), vjust = -0.2) +
  facet_wrap(~dataset) +
  geom_vline(xintercept = 0, color = "darkred" , linetype = "dashed") 
total_plot_nb_sites_per_ali
```




```{r, fig.height=4}
total_nb_sites_per_ali_th <- 
  total_score_all_tested_sites_df %>%
  dplyr::select(alignment,aagtr_pval,log10_aagtr_pval,log10_aagtr_pval_cut,log10_aagtr_pval_cut_th) %>%
  filter(log10_aagtr_pval_cut_th %in% ("inf_th")) %>%
  left_join(total_genes_attributes_ok_filtered %>% dplyr::select(MGI,alignment,padj_DEseq2,fasta_fn,log2FoldChange_DEseq2)) %>%
  group_by(MGI,alignment,fasta_fn,padj_DEseq2) %>%
  summarise(n_site_inf_th=n()) %>%
  arrange(-n_site_inf_th) %>%
  ungroup()%>%
  mutate(
         is_de = padj_DEseq2 < 0.1) %>%
  mutate(n_site_inf_th = ifelse(n_site_inf_th> 40, 40, n_site_inf_th)) %>%
  mutate(dataset = "Total")


total_plot_nb_sites_per_ali_th <- total_nb_sites_per_ali_th %>% 
  group_by(dataset,n_site_inf_th) %>%
  summarise(nb_gene = n()) %>%
  ggplot(aes(x= n_site_inf_th, y = nb_gene)) +
  geom_bar(stat = "identity") +
  theme_bw() +
  geom_text(aes(label = nb_gene), vjust = -0.2) +
  facet_wrap(~dataset) +
  geom_vline(xintercept = 0, color = "blue" , linetype = "dashed") 

total_plot_nb_sites_per_ali_th
```


```{r fig.width=12, fig.height=5}


fig_D <- cowplot::plot_grid(
  cowplot::plot_grid(
    total_th_spike_violin_obs,  total_plot_nb_site_per_cat,  total_plot_nb_gene_per_cat ,  total_plot_nb_sites_per_ali + ylim(c(0,450)), nrow = 1, rel_widths = c(2,1,1,1)
  ),

  cowplot::plot_grid(
    murinae_th_spike_violin_obs,  murinae_plot_nb_site_per_cat ,  murinae_plot_nb_gene_per_cat ,  murinae_plot_nb_sites_per_ali + ylim(c(0,450)), nrow = 1 , rel_widths = c(2,1,1,1)
  ),
  ncol=1
)

 ggsave(plot=fig_D, filename = str_c(output_figures,"D_define_th_using_spike.pdf"), width = 15, height = 10)

 
 fig_D2 <- cowplot::plot_grid(
  cowplot::plot_grid(
    total_th_spike_violin_obs_th,  total_plot_nb_site_per_cat_th,  total_plot_nb_gene_per_cat_th ,  total_plot_nb_sites_per_ali_th + ylim(c(0,550)), nrow = 1, rel_widths = c(2,1,1,1)
  ),
  
  cowplot::plot_grid(
    murinae_th_spike_violin_obs_th,  murinae_plot_nb_site_per_cat_th ,  murinae_plot_nb_gene_per_cat_th ,  murinae_plot_nb_sites_per_ali_th + ylim(c(0,550)), nrow = 1, rel_widths = c(2,1,1,1)
  ),
  ncol=1
)

 ggsave(plot=fig_D2, filename = str_c(output_figures,"D2_define_th_using_neutral_evolution.pdf"), width = 15, height = 10)

  
  

```





## Threshold comparaison and definition


```{r fig.height=8, fig.width=12}
total_neutral_purifying_spike_violin <- bind_rows(total_pelican_neutral_all_tested_results_summary, total_pelican_purifying_all_tested_results_summary, total_pelican_simu_results_summary %>% filter(q2==0, n <= 6),total_best_pval_per_MGI_ok_summary)  %>%
  mutate(facet = "Total") %>%
  ggplot()+
  geom_hline(aes(yintercept = value, col = simu_methode))  +
  scale_color_manual(values = c("Neutral\nsimulation" = "#145db6",
                                "Purifying\nselection\nsimulation" ="#6c2da3",
                                "Perfect\nspike\nsimulation" = "#7c0000",
                                "Best\nObs\nSite"="black")) +
  facet_wrap(~facet + simu_methode, ncol = 4, scale = "free_x") +
  theme_bw() +
  ggrepel::geom_text_repel(aes(label = q_str, y = value), x = 0.5, min.segment.length = 0, size = 2) +
  #geom_point(data = tibble(dataset = "", log10_aagtr_pval = 0),aes(x=dataset, y = log10_aagtr_pval), size = 0, col = "grey") +
  geom_violin(data=total_true_all_sites %>% mutate(dataset = "Total\n  all\nsites", simu_methode = "Best\nObs\nSite" ), 
              aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "black") +
  geom_violin(data=total_score_all_tested_sites_df %>% mutate(dataset = "Total\n  only\n non\nconstant", simu_methode = "Best\nObs\nSite" ), 
              aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "black") +
  geom_violin(data=total_best_pval_per_MGI_ok %>% mutate(dataset = "Total\n  only\nbest\nsites",simu_methode = "Best\nObs\nSite"), 
              aes(x=dataset, y = log10_best_aagtr_pval), alpha = 0.5, col = "black") +
  geom_violin(data=total_pelican_neutral_results %>% mutate(dataset = "Total\n neutral\nevolution", simu_methode = "Neutral\nsimulation") %>% dplyr::select(dataset, log10_aagtr_pval,simu_methode),
              aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#145db6") +
    geom_violin(data=total_neutral_all_tested_sites_df %>% mutate(dataset = "Total\n neutral\nevolution\nonly\n non\nconstant", simu_methode = "Neutral\nsimulation") %>% dplyr::select(dataset, log10_aagtr_pval,simu_methode),
              aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#145db6") +
  geom_violin(data=total_pelican_purifying_results %>% mutate(dataset = "Total\n purifying\nselection", simu_methode ="Purifying\nselection\nsimulation")  %>% dplyr::select(dataset, log10_aagtr_pval,simu_methode) ,
              aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#6c2da3") +
  geom_violin(data=total_purifying_all_tested_sites_df %>% mutate(dataset = "Total\n purifying\nselection\nonly\n non\nconstant", simu_methode ="Purifying\nselection\nsimulation")  %>% dplyr::select(dataset, log10_aagtr_pval,simu_methode) ,
                aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#6c2da3") +
  

  geom_violin(data=total_pelican_simu_results %>% filter(n_t=="1_21") %>% mutate(dataset = "Total\nperfect\nspike\n1/21", simu_methode = "Perfect\nspike\nsimulation")   %>% dplyr::select(dataset, log10_aagtr_pval,simu_methode) ,
              aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#7c0000") +
  geom_violin(data=total_pelican_simu_results %>% filter(n_t=="2_21")  %>% mutate(dataset = "Total\nperfect\nspike\n2/21",simu_methode = "Perfect\nspike\nsimulation")  %>% dplyr::select(dataset, log10_aagtr_pval,simu_methode) ,
              aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#7c0000") +
  geom_violin(data=total_pelican_simu_results  %>% filter(n_t=="3_21") %>% mutate(dataset = "Total\nperfect\nspike\n3/21",simu_methode = "Perfect\nspike\nsimulation")  %>% dplyr::select(dataset, log10_aagtr_pval,simu_methode) ,
              aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#7c0000") +
  geom_violin(data=total_pelican_simu_results  %>% filter(n_t=="4_21") %>% mutate(dataset = "Total\nperfect\nspike\n4/21",simu_methode = "Perfect\nspike\nsimulation")  %>% dplyr::select(dataset, log10_aagtr_pval,simu_methode) ,
              aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#7c0000") 


total_neutral_purifying_spike_violin


```




```{r fig.height=8, fig.width=12}
murinae_neutral_purifying_spike_violin <- bind_rows(murinae_pelican_neutral_all_tested_results_summary, murinae_pelican_purifying_all_tested_results_summary, murinae_pelican_simu_results_summary %>% filter(q2==0, n <= 5),murinae_best_pval_per_MGI_ok_summary)  %>%
  mutate(facet = "Murinae") %>%
  ggplot()+
  geom_hline(aes(yintercept = value, col = simu_methode))  +
  scale_color_manual(values = c("Neutral\nsimulation" = "#145db6", "Purifying\nselection\nsimulation" ="#6c2da3", "Perfect\nspike\nsimulation" = "#7c0000", "Best\nObs\nSite"="black")) +
  facet_wrap(~facet + simu_methode, ncol = 4, scale="free_x") +
  theme_bw() +
  ggrepel::geom_text_repel(aes(label = q_str, y = value), x = 0.5, min.segment.length = 0, size = 2) +
  #geom_point(data = tibble(dataset = "", log10_aagtr_pval = 0),aes(x=dataset, y = log10_aagtr_pval), size = 0, col = "grey") +
  geom_violin(data=murinae_true_all_sites %>% mutate(dataset = "Murinae\n  all\nsites", simu_methode = "Best\nObs\nSite"),
              aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "black") +
  geom_violin(data=murinae_score_all_tested_sites_df %>% mutate(dataset = "Murinae\n  only\n non\nconstant", simu_methode = "Best\nObs\nSite"),
              aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "black") +
  geom_violin(data=murinae_best_pval_per_MGI_ok %>% mutate(dataset = "Murinae\n  only\nbest\nsites",simu_methode = "Best\nObs\nSite"),
              aes(x=dataset, y = log10_best_aagtr_pval), alpha = 0.5, col = "black") +
  geom_violin(data=murinae_pelican_neutral_results %>% mutate(dataset = "Murinae\n neutral\nevolution", simu_methode = "Neutral\nsimulation") %>% dplyr::select(dataset, log10_aagtr_pval,simu_methode),
              aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#145db6") +
  geom_violin(data=murinae_neutral_all_tested_sites_df %>% mutate(dataset = "Murinae\n neutral\nevolution\nonly\n non\nconstant", simu_methode = "Neutral\nsimulation") %>% dplyr::select(dataset, log10_aagtr_pval,simu_methode),
              aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#145db6") +
  
  geom_violin(data=murinae_pelican_purifying_results %>% mutate(dataset = "Murinae\n purifying\nselection",simu_methode = "Purifying\nselection\nsimulation") %>% dplyr::select(dataset, log10_aagtr_pval,simu_methode) ,
              aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#6c2da3") +
  geom_violin(data=murinae_purifying_all_tested_sites_df %>% mutate(dataset = "Murinae\n purifying\nselection\nonly\n non\nconstant", simu_methode ="Purifying\nselection\nsimulation")  %>% dplyr::select(dataset, log10_aagtr_pval,simu_methode) ,
              aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#6c2da3") +
  geom_violin(data=murinae_pelican_simu_results  %>% filter(n_t=="1_5") %>% mutate(dataset = "Murinae\nperfect\nspike\n1/5", simu_methode = "Perfect\nspike\nsimulation")  %>% dplyr::select(dataset, log10_aagtr_pval,simu_methode),
              aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#7c0000") +
  geom_violin(data=murinae_pelican_simu_results  %>% filter(n_t=="2_5") %>% mutate(dataset = "Murinae\nperfect\nspike\n2/5", simu_methode = "Perfect\nspike\nsimulation")  %>% dplyr::select(dataset, log10_aagtr_pval,simu_methode),
              aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#7c0000") +
  geom_violin(data=murinae_pelican_simu_results %>% filter(n_t=="3_5")  %>% mutate(dataset = "Murinae\nperfect\nspike\n3/5", simu_methode = "Perfect\nspike\nsimulation")  %>% dplyr::select(dataset, log10_aagtr_pval,simu_methode),
              aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#7c0000")  +
  geom_violin(data=murinae_pelican_simu_results %>% filter(n_t=="4_5")  %>% mutate(dataset = "Murinae\nperfect\nspike\n4/5", simu_methode = "Perfect\nspike\nsimulation")  %>% dplyr::select(dataset, log10_aagtr_pval,simu_methode),
              aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#7c0000")  +
  geom_violin(data=murinae_pelican_simu_results %>% filter(n_t=="5_5")  %>% mutate(dataset = "Murinae\nperfect\nspike\n5/5", simu_methode = "Perfect\nspike\nsimulation")  %>% dplyr::select(dataset, log10_aagtr_pval,simu_methode),
              aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#7c0000") 


murinae_neutral_purifying_spike_violin


```


```{r}
fig_compa <- total_neutral_purifying_spike_violin/murinae_neutral_purifying_spike_violin

ggsave(plot=fig_compa, filename = str_c(output_figures,"/D_pelican_neutral_purifying_spike_violin.pdf"), width = 12, height = 12)
```



```{r fig.height=4, fig.width=12}
murinae_random_spike_violin <- bind_rows(  murinae_pelican_simu_results_summary %>% filter(q2==0, n <= 5))  %>%
  mutate(facet = "Murinae") %>%
  ggplot()+
  geom_hline(aes(yintercept = value, col = simu_methode), linetype = "dashed")  +
  geom_hline(data = murinae_pelican_neutral_all_tested_results_summary %>% filter(q %in% c(0.001) ), aes(yintercept = value, col = simu_methode))  +
  scale_color_manual(values = c("Neutral\nsimulation" = "#145db6", "Purifying\nselection\nsimulation" ="#6c2da3", "Perfect\nspike\nsimulation" = "#7c0000", "Best\nObs\nSite"="black")) +
  facet_grid(~"Murinae_seq" + simu_methode,  scales = "free_x", space = "free_x")+
  theme_bw() +
   ggrepel::geom_text_repel(data = murinae_pelican_neutral_all_tested_results_summary %>% filter(q %in% c(0.001) ), aes(label = q_str, y = value), x = 0.5, min.segment.length = 0, size = 4, nudge_y = -0.6, direction = "y") +

  ggrepel::geom_text_repel(aes(label = q_str, y = value), x = 0.5, min.segment.length = 0, size = 4, nudge_y = -0.2, direction = "y") +
  #geom_point(data = tibble(dataset = "", log10_aagtr_pval = 0),aes(x=dataset, y = log10_aagtr_pval), size = 0, col = "grey") +
  geom_violin(data=murinae_true_all_sites %>% mutate(dataset = "All\nsites", simu_methode = "Best\nObs\nSite"),
              aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "black") +
  geom_violin(data=murinae_score_all_tested_sites_df %>% mutate(dataset = "Only non\nconstant", simu_methode = "Best\nObs\nSite"),
              aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "black") +
  # geom_violin(data=murinae_pelican_neutral_results %>% mutate(dataset = " Random\nevolution", simu_methode = "Neutral\nsimulation") %>% dplyr::select(dataset, log10_aagtr_pval,simu_methode),
  #             aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#145db6") +
  geom_violin(data=murinae_neutral_all_tested_sites_df %>% mutate(dataset = "Negative\ncontroles", simu_methode = "Neutral\nsimulation") %>% dplyr::select(dataset, log10_aagtr_pval,simu_methode),
              aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#145db6") +
  geom_violin(data=murinae_pelican_simu_results  %>% filter(n_t=="1_5") %>% mutate(dataset = "1/5", simu_methode = "Perfect\nspike\nsimulation")  %>% dplyr::select(dataset, log10_aagtr_pval,simu_methode),
              aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#7c0000") +
  geom_violin(data=murinae_pelican_simu_results  %>% filter(n_t=="2_5") %>% mutate(dataset = "2/5", simu_methode = "Perfect\nspike\nsimulation")  %>% dplyr::select(dataset, log10_aagtr_pval,simu_methode),
              aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#7c0000") +
  geom_violin(data=murinae_pelican_simu_results %>% filter(n_t=="3_5")  %>% mutate(dataset = "3/5\nPositive controles", simu_methode = "Perfect\nspike\nsimulation")  %>% dplyr::select(dataset, log10_aagtr_pval,simu_methode),
              aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#7c0000")  +
  geom_violin(data=murinae_pelican_simu_results %>% filter(n_t=="4_5")  %>% mutate(dataset = "4/5", simu_methode = "Perfect\nspike\nsimulation")  %>% dplyr::select(dataset, log10_aagtr_pval,simu_methode),
              aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#7c0000")  +
  geom_violin(data=murinae_pelican_simu_results %>% filter(n_t=="5_5")  %>% mutate(dataset = "5/5", simu_methode = "Perfect\nspike\nsimulation")  %>% dplyr::select(dataset, log10_aagtr_pval,simu_methode),
              aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#7c0000") +
  theme(legend.position = "none") +
  xlab("") +
  ylab("log10(pvalue.murinae)")


murinae_random_spike_violin


```



```{r fig.height=4, fig.width=12}
total_random_spike_violin <-bind_rows(  total_pelican_simu_results_summary %>% filter(q2==0, n <= 5))  %>%
  mutate(facet = "total") %>%
  ggplot()+
  geom_hline(aes(yintercept = value, col = simu_methode), linetype = "dashed")  +
  geom_hline(data = total_pelican_neutral_all_tested_results_summary %>% filter(q %in% c(0.001) ), aes(yintercept = value, col = simu_methode))  +
  scale_color_manual(values = c("Neutral\nsimulation" = "#145db6", "Purifying\nselection\nsimulation" ="#6c2da3", "Perfect\nspike\nsimulation" = "#7c0000", "Best\nObs\nSite"="black")) +
  facet_grid(~"Total_seq" + simu_methode,  scales = "free_x", space = "free_x")+
  theme_bw() +
   ggrepel::geom_text_repel(data = total_pelican_neutral_all_tested_results_summary %>% filter(q %in% c(0.001) ), aes(label = q_str, y = value), x = 0.5, min.segment.length = 0, size = 4, nudge_y = -0.6, direction = "y") +

  ggrepel::geom_text_repel(aes(label = q_str, y = value), x = 0.5, min.segment.length = 0, size = 4, nudge_y = -0.2, direction = "y") +
  #geom_point(data = tibble(dataset = "", log10_aagtr_pval = 0),aes(x=dataset, y = log10_aagtr_pval), size = 0, col = "grey") +
  geom_violin(data=total_true_all_sites %>% mutate(dataset = "All\nsites", simu_methode = "Best\nObs\nSite"),
              aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "black") +
  geom_violin(data=total_score_all_tested_sites_df %>% mutate(dataset = "Only non\nconstant", simu_methode = "Best\nObs\nSite"),
              aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "black") +
  # geom_violin(data=total_pelican_neutral_results %>% mutate(dataset = " Random\nevolution", simu_methode = "Neutral\nsimulation") %>% dplyr::select(dataset, log10_aagtr_pval,simu_methode),
  #             aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#145db6") +
  geom_violin(data=total_neutral_all_tested_sites_df %>% mutate(dataset = "Negative\ncontroles", simu_methode = "Neutral\nsimulation") %>% dplyr::select(dataset, log10_aagtr_pval,simu_methode),
              aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#145db6") +
  geom_violin(data=total_pelican_simu_results  %>% filter(n_t=="1_21") %>% mutate(dataset = "1/21", simu_methode = "Perfect\nspike\nsimulation")  %>% dplyr::select(dataset, log10_aagtr_pval,simu_methode),
              aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#7c0000") +
  geom_violin(data=total_pelican_simu_results  %>% filter(n_t=="2_21") %>% mutate(dataset = "2/21", simu_methode = "Perfect\nspike\nsimulation")  %>% dplyr::select(dataset, log10_aagtr_pval,simu_methode),
              aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#7c0000") +
  geom_violin(data=total_pelican_simu_results %>% filter(n_t=="3_21")  %>% mutate(dataset = "3/21\nPositive controles", simu_methode = "Perfect\nspike\nsimulation")  %>% dplyr::select(dataset, log10_aagtr_pval,simu_methode),
              aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#7c0000")  +
  geom_violin(data=total_pelican_simu_results %>% filter(n_t=="4_21")  %>% mutate(dataset = "4/21", simu_methode = "Perfect\nspike\nsimulation")  %>% dplyr::select(dataset, log10_aagtr_pval,simu_methode),
              aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#7c0000")  +
  geom_violin(data=total_pelican_simu_results %>% filter(n_t=="5_21")  %>% mutate(dataset = "5/21", simu_methode = "Perfect\nspike\nsimulation")  %>% dplyr::select(dataset, log10_aagtr_pval,simu_methode),
              aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#7c0000") +
    geom_violin(data=total_pelican_simu_results %>% filter(n_t=="6_21")  %>% mutate(dataset = "6/21", simu_methode = "Perfect\nspike\nsimulation")  %>% dplyr::select(dataset, log10_aagtr_pval,simu_methode),
              aes(x=dataset, y = log10_aagtr_pval), alpha = 0.5, col = "#7c0000") +
  theme(legend.position = "none")  +
  xlab("") +
  ylab("log10(pvalue.total)")


total_random_spike_violin


```


```{r}
ggsave(plot=total_random_spike_violin, filename = str_c(output_figures,"/D_pelican_random_spike_violin_total.pdf"), width = 18, height = 10, units = "cm")

ggsave(plot=murinae_random_spike_violin, filename = str_c(output_figures,"/D_pelican_random_spike_violin_murinae.pdf"), width = 13, height = 10, units = "cm")

ggsave(plot=murinae_random_spike_violin/total_random_spike_violin, filename = str_c(output_figures,"/D_pelican_random_spike_violin_murinae_total.pdf"), width = 14, height = 20, units = "cm")

```

##effect of the nb of tested site per ali on best pval


```{r}

len_bornes <- c(0,100,200,300,400,500,600,Inf)
murinae_len_bornes <- quantile(murinae_best_pval_per_MGI_ok$ali_length, c(0,0.2,0.4,0.6,0.8,1))
total_len_bornes <- quantile(total_best_pval_per_MGI_ok$ali_length, c(0,0.2,0.4,0.6,0.8,1))

murinae_best_pval_per_MGI_ok <-  murinae_best_pval_per_MGI_ok %>%
  mutate(group = ifelse(log10_best_aagtr_pval < murinae_pelican_neutral_all_tested_results_summary %>% filter(q == 0.001) %>% pull(value), "sites inf\nthreshold","sites sup\nthreshold"),
         group_len =  cut(ali_length, breaks = murinae_len_bornes, include.lowest = TRUE, right = FALSE)) 

table(murinae_best_pval_per_MGI_ok$group_len)

total_best_pval_per_MGI_ok <-  total_best_pval_per_MGI_ok %>% 
  mutate(group = ifelse(log10_best_aagtr_pval < total_pelican_neutral_all_tested_results_summary %>% filter(q == 0.001) %>% pull(value), "sites inf\nthreshold","sites sup\nthreshold"),
          group_len =  cut(ali_length, breaks = total_len_bornes, include.lowest = TRUE, right = FALSE)) 

table(total_best_pval_per_MGI_ok$group_len)

murinae_var_all_sites <- murinae_best_pval_per_MGI_ok %>%
   mutate(group = "all sites\ntogether") %>%
  group_by(dataset,group) %>%
  summarise(R2_pval = summary(lm(best_aagtr_pval~ ali_length))$r.squared,
            R2_log10_pval = summary(lm(log10_best_aagtr_pval~ ali_length))$r.squared)

murinae_var_all_sites_split <- murinae_best_pval_per_MGI_ok %>%
  group_by(dataset,group) %>%
  summarise(R2_pval = summary(lm(best_aagtr_pval~ ali_length))$r.squared,
            R2_log10_pval = summary(lm(log10_best_aagtr_pval~ ali_length))$r.squared)


murinae_var_all_sites_split_len <- murinae_best_pval_per_MGI_ok %>%
  group_by(dataset,group_len) %>%
  summarise(R2_pval = summary(lm(best_aagtr_pval~ ali_length))$r.squared,
            R2_log10_pval = summary(lm(log10_best_aagtr_pval~ ali_length))$r.squared)



total_var_all_sites <- total_best_pval_per_MGI_ok %>%
  mutate(group = "all sites\ntogether") %>%
  group_by(dataset,group) %>%
  summarise(R2_pval = summary(lm(best_aagtr_pval~ ali_length))$r.squared,
            R2_log10_pval = summary(lm(log10_best_aagtr_pval~ ali_length))$r.squared)

total_var_all_sites_split <- total_best_pval_per_MGI_ok %>%
   group_by(dataset,group) %>%
   summarise(R2_pval = summary(lm(best_aagtr_pval~ ali_length))$r.squared,
            R2_log10_pval = summary(lm(log10_best_aagtr_pval~ ali_length))$r.squared)


total_var_all_sites_split_len <- total_best_pval_per_MGI_ok %>%
  group_by(dataset,group_len) %>%
  summarise(R2_pval = summary(lm(best_aagtr_pval~ ali_length))$r.squared,
            R2_log10_pval = summary(lm(log10_best_aagtr_pval~ ali_length))$r.squared)



fig_variance_log10 <- bind_rows(murinae_var_all_sites,murinae_var_all_sites_split,total_var_all_sites,total_var_all_sites_split) %>%
  ggplot(aes(x= group, y =R2_log10_pval)) +
  geom_col(fill = "grey", col = "black") +
  facet_wrap(~dataset) +
  theme_bw() +
  geom_text(aes(label = round(R2_log10_pval, digit = 3), vjust = c(1.5,-0.5,1.5,1.5,-0.5,1.5))) +
  labs(title = "summary(lm(log10(best_aagtr_pval) ~ ali_length))$r.squared", y = "")
fig_variance_log10

fig_variance <- bind_rows(murinae_var_all_sites,murinae_var_all_sites_split,total_var_all_sites,total_var_all_sites_split) %>%
  ggplot(aes(x= group, y =R2_pval)) +
  geom_col(fill = "grey", col = "black") +
  facet_wrap(~dataset) +
  theme_bw() +
  geom_text(aes(label = round(R2_pval, digit = 3), vjust = c(1.5,-0.5,1.5,1.5,-0.5,1.5))) +
  labs(title = "summary(lm(best_aagtr_pval ~ ali_length))$r.squared", y = "")
fig_variance



fig_variance_log10_len <- bind_rows(murinae_var_all_sites_split_len, total_var_all_sites_split_len) %>%
  ggplot(aes(x= group_len, y =R2_log10_pval)) +
  geom_col(fill = "grey", col = "black") +
  facet_wrap(~dataset, scale = "free_x") +
  theme_bw() +
  geom_text(aes(label = round(R2_log10_pval, digit = 3), vjust = c(1.5,-0.5,-0.5,-0.5,-0.5,1.5,-0.5,-0.5,-0.5,-0.5))) +
  labs(title = "summary(lm(log10(best_aagtr_pval) ~ ali_length))$r.squared", y = "")
fig_variance_log10_len
```

```{r fig.height=10, fig.width=15}
relation_min_aagtr_pval_nsites_group <-  bind_rows(total_best_pval_per_MGI_ok,murinae_best_pval_per_MGI_ok) %>% 
  ggplot(aes(x=ali_length, y = log10_best_aagtr_pval)) +
  theme_bw() +
  ggrastr::geom_point_rast(aes(color =  paste0(dataset,group)), alpha = 0.2) +
  #geom_density2d() +                      # points
  geom_smooth(method = "lm", se = T,  aes(group =  paste0(dataset,group)),color = "black")  +
   facet_wrap(~dataset, scale= "free_x") +
  geom_hline(data = murinae_pelican_simu_results_summary%>% filter(n != 2), aes(yintercept = med_pval), linetype = "dashed", color = "grey15" ) +
  geom_hline(data = total_pelican_simu_results_summary %>% filter(n <=9) %>% filter(n != 4), aes(yintercept = med_pval), linetype = "dashed", color = "grey15" ) +
  geom_hline(data = murinae_pelican_simu_results_summary %>% filter(n_t == "2_5"), aes(yintercept = med_pval), size = 1.5, color = "darkred" , linetype = "dashed" ) +
  geom_hline(data = total_pelican_simu_results_summary %>% filter(n_t == "4_21"), aes(yintercept = med_pval), size = 1.5, color = "darkred" , linetype = "dashed") +
  geom_text(data = murinae_pelican_simu_results_summary, x = 500, aes( y = med_pval, label = str_replace(n_t,"_","/")), vjust = 1.5 , color = "grey15", size = 4 ) +
  geom_text(data = total_pelican_simu_results_summary %>% filter(n <=9), x = 900, aes( y = med_pval, label = str_replace(n_t,"_","/")), vjust = 1.5 , color = "grey15", size = 4 )

relation_min_aagtr_pval_nsites_group


relation_min_aagtr_pval_nsites_group_len <-  bind_rows(total_best_pval_per_MGI_ok,murinae_best_pval_per_MGI_ok) %>% 
  ggplot(aes(x=ali_length, y = log10_best_aagtr_pval)) +
  theme_bw() +
  ggrastr::geom_point_rast(aes(color =  paste0(dataset,group_len)), alpha = 0.2) +
  #geom_density2d() +                      # points
  geom_smooth(method = "lm", se = T,  aes(group =  paste0(dataset,group_len)),color = "black")  +
  facet_wrap(~dataset, scale= "free_x") +
  geom_hline(data = murinae_pelican_simu_results_summary %>% filter(n_t == "2_5"), aes(yintercept = med_pval), size = 1.5, color = "darkred" , linetype = "dashed" ) +
  geom_hline(data = total_pelican_simu_results_summary %>% filter(n_t == "4_21"), aes(yintercept = med_pval), size = 1.5, color = "darkred" , linetype = "dashed") 
 
  

relation_min_aagtr_pval_nsites_group_len


relation_min_aagtr_pval_nsites <-  bind_rows(total_best_pval_per_MGI_ok,murinae_best_pval_per_MGI_ok) %>% 
  ggplot(aes(x=ali_length, y = log10_best_aagtr_pval)) +
  theme_bw() +
  ggrastr::geom_point_rast(color =  "grey", alpha = 0.4, size = 1) +
  #geom_density2d() +                      # points
  geom_smooth(method = "lm", se = T, color = "black")  +
  facet_wrap(~dataset, scale= "free_x") +
  geom_hline(data = murinae_pelican_simu_results_summary%>% filter(q2 == 0, n <=5), aes(yintercept = value), linetype = "dashed", color = "darkred" ) +
  geom_hline(data = total_pelican_simu_results_summary %>% filter(q2 == 0, n <=6), aes(yintercept = value), linetype = "dashed", color = "darkred" ) +
  geom_hline(data = murinae_pelican_neutral_all_tested_results_summary %>% filter(q == 0.001), aes(yintercept = value), size = 1, color = "#145db6"  ) +
  geom_hline(data = total_pelican_neutral_all_tested_results_summary %>% filter(q == 0.001), aes(yintercept = value), size = 1, color = "#145db6" ) +
  geom_text(data = murinae_pelican_simu_results_summary %>% filter(q2 == 0, n <=5), x = 1300, aes( y = value, label = str_replace(n_t,"_","/")), vjust = 1.5 , color = "grey15", size = 4 ) +
  geom_text(data = total_pelican_simu_results_summary%>% filter(q2 == 0, n <=6), x = 2000, aes( y = value, label = str_replace(n_t,"_","/")), vjust = 1.5 , color = "grey15", size = 4 )

relation_min_aagtr_pval_nsites
```


```{r fig.height=10, fig.width=15}
modele_total <- lm(log10_best_aagtr_pval ~ ali_length, data = total_best_pval_per_MGI_ok)
plot(modele_total)
summary(modele_total)
summarry_modele_total <- summary(modele_total)
modele_total_rsquared = summarry_modele_total$r.squared
modele_total_pval = summarry_modele_total$coefficients[,4] [2]  
# Call:
# lm(formula = log10_best_aagtr_pval ~ ali_length, data = total_best_pval_per_MGI_ok)
# 
# Residuals:
#      Min       1Q   Median       3Q      Max 
# -2.22908 -0.22501  0.02696  0.24770  1.21423 
# 
# Coefficients:
#               Estimate Std. Error t value Pr(>|t|)    
# (Intercept) -1.353e+00  1.706e-02  -79.28   <2e-16 ***
# ali_length  -3.911e-04  2.843e-05  -13.76   <2e-16 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 0.3857 on 2348 degrees of freedom
# Multiple R-squared:  0.07457,	Adjusted R-squared:  0.07417 
# F-statistic: 189.2 on 1 and 2348 DF,  p-value: < 2.2e-16


modele_murinae <- lm(log10_best_aagtr_pval ~ ali_length, data = murinae_best_pval_per_MGI_ok)
plot(modele_murinae)
summary(modele_murinae)
summarry_modele_murinae <- summary(modele_murinae)
modele_murinae_rsquared = summarry_modele_murinae$r.squared
modele_murinae_pval = summarry_modele_murinae$coefficients[,4] [2]  

# Call:
# lm(formula = log10_best_aagtr_pval ~ ali_length, data = murinae_best_pval_per_MGI_ok)
# 
# Residuals:
#      Min       1Q   Median       3Q      Max 
# -2.18511 -0.22582  0.07553  0.19370  1.03106 
# 
# Coefficients:
#               Estimate Std. Error t value Pr(>|t|)    
# (Intercept) -1.139e+00  1.656e-02  -68.81   <2e-16 ***
# ali_length  -3.647e-04  3.532e-05  -10.33   <2e-16 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 0.3568 on 2272 degrees of freedom
# Multiple R-squared:  0.04483,	Adjusted R-squared:  0.04441 
# F-statistic: 106.6 on 1 and 2272 DF,  p-value: < 2.2e-16




ggsave(plot=relation_min_aagtr_pval_nsites +
  theme(
    text = element_text(size = 8),          # taille par défaut pour tous les textes
    axis.title = element_text(size = 8),    # titre des axes
    axis.text = element_text(size = 8),     # étiquettes des axes
    strip.text = element_text(size = 8)     # titres des facettes
  ) +  ggtitle(
  str_c("Murinae modele: ", "R^2: ", modele_murinae_rsquared , " Pval: ", modele_murinae_pval, "\n",
        "Total modele: ", "R^2: ", modele_total_rsquared , " Pval: ", modele_total_pval)
), filename = str_c(output_figures,"Fsupp_correlation_best_pval_ali_length.pdf"), width = 17, height = 8, unit = "cm")
```


```{r fig.height=12, fig.width=8}
fig_F_correlation_best_pval_ali_length <- cowplot::plot_grid(
  relation_min_aagtr_pval_nsites + theme(legend.position = "top"),
  fig_variance_log10,
  relation_min_aagtr_pval_nsites_group_len+ theme(legend.position = "top"),
  fig_variance_log10_len,
  relation_min_aagtr_pval_nsites_group + theme(legend.position = "top"),
  
  fig_variance,
  ncol = 1, rel_heights = c(2,1,2,1,2,1)
)
fig_F_correlation_best_pval_ali_length

ggsave(plot=fig_F_correlation_best_pval_ali_length, filename = str_c(output_figures,"F2_correlation_best_pval_ali_length.pdf"), width = 8, height = 18)

  
```

## plot top sites per datasets


```{r}
murinae_best_sites<- 
  murinae_score_all_tested_sites_df %>%
  filter(log10_aagtr_pval <= murinae_pelican_neutral_all_tested_results_summary %>% filter(q == 0.001) %>% pull(value)) %>%
  left_join(murinae_genes_attributes_ok %>% dplyr::select(MGI,alignment,fasta_fn))  %>%
  arrange(log10_aagtr_pval) %>%
  mutate(MGI_site = paste0(MGI,"_",site,"_",dp),
         labels_simplified = paste0(MGI,"_",site)) %>%
  slice_min(log10_aagtr_pval, n= 20)


total_best_sites<- 
  total_score_all_tested_sites_df %>%
  filter(log10_aagtr_pval <= total_pelican_neutral_all_tested_results_summary %>% filter(q == 0.001) %>% pull(value)) %>%
  left_join(total_genes_attributes_ok %>% dplyr::select(MGI,alignment,fasta_fn)) %>%
  arrange(log10_aagtr_pval) %>%
  mutate(MGI_site = paste0(MGI,"_",site,"_",dp),
         labels_simplified = paste0(MGI,"_",site))   %>%
  slice_min(log10_aagtr_pval, n= 20)



murinae_best_sites_with_seq <- murinae_best_sites %>%
   ungroup() %>%
  rowwise() %>%
  mutate(seq_per_sp_df = list(get_seq_from_fasta_fn_pos(fasta_fn, site))) %>%
  unnest(seq_per_sp_df)  %>%
  left_join(species_attributes %>% dplyr::select(leaf_name, Aride), by = c("seqnames" = "leaf_name"))  %>%
  mutate(MGI_site = paste0(MGI,"_",site,"_",dp),
                                 labels_simplified = paste0(MGI,"_",site)) %>%
  mutate(MGI_site = factor(MGI_site, levels = murinae_best_sites$MGI_site),
         MGI = factor(MGI, levels = unique(murinae_best_sites$MGI)),
         seqnames = factor(seqnames, levels = rev(murinae_tips_order))) 



total_best_sites_with_seq <- total_best_sites %>%
   ungroup() %>%
  rowwise() %>%
  mutate(seq_per_sp_df = list(get_seq_from_fasta_fn_pos(fasta_fn, site))) %>%
  unnest(seq_per_sp_df)  %>%
  left_join(species_attributes %>% dplyr::select(leaf_name, Aride), by = c("seqnames" = "leaf_name"))  %>%
  mutate(MGI_site = paste0(MGI,"_",site,"_",dp),
                                 labels_simplified = paste0(MGI,"_",site)) %>%
  mutate(MGI_site = factor(MGI_site, levels = total_best_sites$MGI_site),
         MGI = factor(MGI, levels = unique(total_best_sites$MGI)),
         seqnames = factor(seqnames, levels = rev(total_tips_order))) 



combined_best_sites <- combined_score_all_tested_sites_df %>%
  filter( (log10_aagtr_pval.murinae <= murinae_pelican_neutral_all_tested_results_summary %>% filter(q == 0.001) %>% pull(value) &
             log10_aagtr_pval.total <= total_pelican_neutral_all_tested_results_summary %>% filter(q == 0.001) %>% pull(value)))  %>%
  left_join(total_genes_attributes_ok %>% dplyr::select(MGI,alignment,fasta_fn))  %>%
  arrange(log10_aagtr_pval.total) %>%
  mutate(MGI_site = paste0(MGI,"_",site,"_T"),
         labels_simplified = paste0(MGI,"_",site))  %>%
  slice_min(log10_aagtr_pval.total*log10_aagtr_pval.murinae, n= 20)
  
  
combined_best_sites_with_seq <- combined_best_sites %>%
   ungroup() %>%
  rowwise() %>%
  mutate(seq_per_sp_df = list(get_seq_from_fasta_fn_pos(fasta_fn, site))) %>%
  unnest(seq_per_sp_df)  %>%
  left_join(species_attributes %>% dplyr::select(leaf_name, Aride), by = c("seqnames" = "leaf_name"))  %>%
  mutate(MGI_site = paste0(MGI,"_",site,"_T"),
                                 labels_simplified = paste0(MGI,"_",site)) %>%
  mutate(MGI_site = factor(MGI_site, levels = combined_best_sites$MGI_site),
         MGI = factor(MGI, levels = unique(combined_best_sites$MGI)),
         seqnames = factor(seqnames, levels = rev(total_tips_order))) 





combined_best_sites_murinae <- combined_score_all_tested_sites_df_top %>%
  filter( (log10_aagtr_pval.murinae <= murinae_pelican_neutral_all_tested_results_summary %>% filter(q == 0.001) %>% pull(value) &
             log10_aagtr_pval.total <= total_pelican_neutral_all_tested_results_summary %>% filter(q == 0.001) %>% pull(value)))  %>%
  left_join(murinae_genes_attributes_ok %>% dplyr::select(MGI,alignment,fasta_fn))  %>%
  arrange(log10_aagtr_pval.total) %>%
  mutate(MGI_site = paste0(MGI,"_",site,"_T"),
         labels_simplified = paste0(MGI,"_",site))  %>%
  slice_min(log10_aagtr_pval.total+log10_aagtr_pval.murinae, n= 20)
  
  
combined_best_sites_murinae_with_seq <- combined_best_sites_murinae %>%
   ungroup() %>%
  rowwise() %>%
  mutate(seq_per_sp_df = list(get_seq_from_fasta_fn_pos(fasta_fn, site))) %>%
  unnest(seq_per_sp_df)  %>%
  left_join(species_attributes %>% dplyr::select(leaf_name, Aride), by = c("seqnames" = "leaf_name"))  %>%
  mutate(MGI_site = paste0(MGI,"_",site,"_T"),
                                 labels_simplified = paste0(MGI,"_",site)) %>%
  mutate(MGI_site = factor(MGI_site, levels = combined_best_sites_murinae$MGI_site),
         MGI = factor(MGI, levels = unique(combined_best_sites_murinae$MGI)),
         seqnames = factor(seqnames, levels = rev(total_tips_order))) 

```





```{r}
#sudo apt install libproj-dev
#BiocManager::install("ggmsa")
library(ggmsa)

scheme_AA <- ggmsa:::scheme_AA
scheme_AA["-","LETTER"] = "grey80"
scheme_AA["P","LETTER"] = "#f1cc37"
scheme_AA["NA","LETTER"] = "#ffffff"

palette_AA_LETTER = scheme_AA$LETTER
names(palette_AA_LETTER) = rownames(scheme_AA)


```


```{r fig.height=10, fig.width=15}

murinae_plot_bestX_pos_th <- ggplot(murinae_best_sites_with_seq, aes(x = MGI_site, y = seqnames, fill = AA, label= AA, col = Aride)) +
 # geom_hline(yintercept = c(0.5,3.5,7.5,6.5,11.5,12.5,15.5,14.5,13.5)) +
  geom_rect(
    data = murinae_limites_y, 
    aes(
      xmin = -Inf,
      xmax = Inf,
      ymin = min_y,
      ymax = max_y,
      fill = Aride
    ),
    alpha = 1, # Transparence pastel
    inherit.aes = FALSE
  ) +
  
  geom_tile(alpha=0.8,  size = 1, width = 0.8, height=0.8, col = ifelse(murinae_best_sites_with_seq$Aride =="yes", "grey85","white")) + geom_text(  size = 6, color = "black") +

  scale_color_manual(values = c("darkgreen","darkorange")) +
   theme_bw()  + ggforce::facet_row(~MGI, scale = "free_x", space='free', strip.position="top")  +
  scale_x_discrete(position = "top",
                   breaks = murinae_best_sites_with_seq$MGI_site,
                   labels = murinae_best_sites_with_seq$labels_simplified) +
  theme(legend.key=element_blank(),panel.background = element_rect(fill = "white"), panel.grid.major = element_blank(),
        axis.text.x = element_text(colour = "black", size = 12, angle = 90, vjust = 0.3, hjust = 0)) +
  scale_fill_manual(values = c(palette_AA_LETTER,"yes"="grey85"))  +
  theme(legend.position="top") +guides(fill=guide_legend(nrow=1,byrow=TRUE),
                                       color=guide_legend(override.aes=list(fill=NA))) 
murinae_plot_bestX_pos_th


 fig_B2_site <- cowplot::plot_grid(
    cowplot::plot_grid(NULL, murinae_p  , rel_heights = c(2.3,16), ncol=1),
    cowplot::plot_grid(murinae_plot_bestX_pos_th + theme(axis.text.y=element_blank()) + ylab(""), NULL, rel_heights = c(100,3), ncol=1),
                  rel_widths = c(10,19))
 
 ggsave(plot=fig_B2_site, filename = str_c(output_figures,"B2_pelican_murinae_best_sites.pdf"), width = 15, height = 10)

  
```





```{r fig.height=20, fig.width=20}

total_plot_bestX_pos_th <- ggplot(total_best_sites_with_seq, aes(x = MGI_site, y = seqnames, fill = AA, label= AA, col = Aride)) +
 # geom_hline(yintercept = c(0.5,3.5,7.5,6.5,11.5,12.5,15.5,14.5,13.5)) +
  geom_rect(
    data = total_limites_y, 
    aes(
      xmin = -Inf,
      xmax = Inf,
      ymin = min_y,
      ymax = max_y,
      fill = Aride
    ),
    alpha = 1, # Transparence pastel
    inherit.aes = FALSE
  ) +
  
  geom_tile(alpha=0.8,  size = 1, width = 0.8, height=0.8, col = ifelse(total_best_sites_with_seq$Aride =="yes", "grey85","white")) +
  geom_text(size = 6, color = "black") +
  
  scale_color_manual(values = c("darkgreen","darkorange")) +
  theme_bw()  + ggforce::facet_row(~MGI, scale = "free_x", space='free', strip.position="top")  +
  scale_x_discrete(position = "top",
                   breaks = total_best_sites_with_seq$MGI_site,
                   labels = total_best_sites_with_seq$labels_simplified) +
  theme(legend.key=element_blank(),panel.background = element_rect(fill = "white"), panel.grid.major = element_blank(),
        axis.text.x = element_text(colour = "black", size = 12, angle = 90, vjust = 0.3, hjust = 0)) +
  scale_fill_manual(values = c(palette_AA_LETTER,"yes"="grey85"))  +
  theme(legend.position="top") +guides(fill=guide_legend(nrow=1,byrow=TRUE),
                                       color=guide_legend(override.aes=list(fill=NA))) 
total_plot_bestX_pos_th


fig_B1_site <- cowplot::plot_grid(
    cowplot::plot_grid(NULL, total_p  + theme(legend.position = "none"), rel_heights = c(0.8,16), ncol=1),
    cowplot::plot_grid(total_plot_bestX_pos_th + theme(axis.text.y=element_blank()) + ylab("") + theme(legend.position = "none") +xlab("") ,
                       NULL, rel_heights = c(150,3), ncol=1),
                  rel_widths = c(4,6))
fig_B1_site
ggsave(plot=fig_B1_site  , filename = str_c(output_figures,"B1_pelican_total_best_sites.pdf"), width = 18, height = 24)

```





```{r fig.height=10, fig.width=15}

combined_plot_bestX_pos_murinae_th <- ggplot(combined_best_sites_murinae_with_seq, aes(x = MGI_site, y = seqnames, fill = AA, label= AA, col = Aride)) +
 # geom_hline(yintercept = c(0.5,3.5,7.5,6.5,11.5,12.5,15.5,14.5,13.5)) +
  geom_rect(
    data = murinae_limites_y, 
    aes(
      xmin = -Inf,
      xmax = Inf,
      ymin = min_y,
      ymax = max_y,
      fill = Aride
    ),
    alpha = 1, # Transparence pastel
    inherit.aes = FALSE
  ) +
  
  geom_tile(alpha=0.8,  size = 1, width = 0.8, height=0.8, col = ifelse(combined_best_sites_murinae_with_seq$Aride =="yes", "grey85","white")) +
  geom_text(  size = 3, color = "black") +

  scale_color_manual(values = c("darkgreen","darkorange")) +
   theme_bw()  + ggforce::facet_row(~MGI, scale = "free_x", space='free', strip.position="top")  +
  scale_x_discrete(position = "top",
                   breaks = combined_best_sites_murinae_with_seq$MGI_site,
                   labels = combined_best_sites_murinae_with_seq$labels_simplified) +
  theme(legend.key=element_blank(),panel.background = element_rect(fill = "white"), panel.grid.major = element_blank(),
        axis.text.x = element_text(colour = "black", size = 12, angle = 90, vjust = 0.3, hjust = 0)) +
  scale_fill_manual(values = c(palette_AA_LETTER,"yes"="grey85"))  +
  theme(legend.position="top") +guides(fill=guide_legend(nrow=1,byrow=TRUE),
                                       color=guide_legend(override.aes=list(fill=NA))) 
combined_plot_bestX_pos_murinae_th


 fig_combine_murinae_site <- cowplot::plot_grid(
    cowplot::plot_grid(NULL, murinae_p  , rel_heights = c(4,16), ncol=1),
    cowplot::plot_grid(murinae_plot_bestX_pos_th + theme(axis.text.y=element_blank()) + ylab(""), NULL, rel_heights = c(100,3), ncol=1),
                  rel_widths = c(10,19))
 
 ggsave(plot=fig_combine_murinae_site, filename = str_c(output_figures,"B3_pelican_combined_murinae_best_sites.pdf"), width = 14, height = 6)

  
```





```{r fig.height=20, fig.width=20}

combined_plot_bestX_pos_th <- ggplot(combined_best_sites_with_seq, aes(x = MGI_site, y = seqnames, fill = AA, label= AA, col = Aride)) +
 # geom_hline(yintercept = c(0.5,3.5,7.5,6.5,11.5,12.5,15.5,14.5,13.5)) +
  geom_rect(
    data = total_limites_y %>% mutate(Aride = ifelse(transition %in% c("T1","T2","T3","T4","T5"), "yes_m", "yes")), 
    aes(
      xmin = -Inf,
      xmax = Inf,
      ymin = min_y,
      ymax = max_y,
      fill = Aride
    ),
    alpha = 1, # Transparence pastel
    inherit.aes = FALSE
  ) +
  
  geom_tile(alpha=0.8,  size = 1, width = 0.8, height=0.8, col = ifelse(combined_best_sites_with_seq$Aride =="yes", "grey85","white")) + geom_text(  size = 6, color = "black") +

  scale_color_manual(values = c("darkgreen","darkorange")) +
   theme_bw()  + ggforce::facet_row(~MGI, scale = "free_x", space='free', strip.position="top")  +
  scale_x_discrete(position = "top",
                   breaks = combined_best_sites_with_seq$MGI_site,
                   labels = combined_best_sites_with_seq$labels_simplified) +
  theme(legend.key=element_blank(),panel.background = element_rect(fill = "white"), panel.grid.major = element_blank(),
        axis.text.x = element_text(colour = "black", size = 12,  angle = 90, vjust = 0.3, hjust = 0)) +
  scale_fill_manual(values = c(palette_AA_LETTER,"yes"="grey85", "yes_m" = "#9a9ca1"))  +
  theme(legend.position="top") +guides(fill=guide_legend(nrow=1,byrow=TRUE),
                                       color=guide_legend(override.aes=list(fill=NA))) 
combined_plot_bestX_pos_th




figC_site <- cowplot::plot_grid(
  cowplot::plot_grid(NULL, total_p  , rel_heights = c(1.1,16), ncol=1),
  cowplot::plot_grid(combined_plot_bestX_pos_th + theme(axis.text.y=element_blank()) + ylab(""),
                     NULL, rel_heights = c(150,1), ncol=1),
  rel_widths =  c(4,6))

ggsave(plot=figC_site, filename = str_c(output_figures,"C_pelican_crossing_total_murinae_best_sites.pdf"), width = 20, height = 20)


  
```




```{r fig.height=4, fig.width=20}
##TODO
score_all_sites_all_dataset_df_top_X_add_genes_att_with_seq_suumary <- 
  score_all_sites_all_dataset_df_top_X_add_genes_att_with_seq %>%
  group_by(dataset,MGI_site,Aride) %>%
  mutate(n_group = n()) %>%
  group_by(dataset,MGI_site,Aride,AA) %>%
  dplyr::summarise(n_cond = n(),
            perc = n_cond/unique(n_group) * 100) %>%
  mutate(MGI_site_Aride = paste0(MGI_site,"_",Aride),
         MGI_site_Aride = factor(MGI_site_Aride),
         MGI_site_Aride_n = as.numeric(MGI_site_Aride) - 1,
         MGI_site_Aride_n2 = ifelse(Aride == "yes", MGI_site_Aride_n -0.1, MGI_site_Aride_n +0.1))
  


df_values_breaks_MGI <- tibble(
  values_x = score_all_sites_all_dataset_df_top_X_add_genes_att_with_seq_suumary$MGI_site_Aride_n,
  labels_unique = score_all_sites_all_dataset_df_top_X_add_genes_att_with_seq_suumary$MGI_site,
  labels_simplified = str_remove(score_all_sites_all_dataset_df_top_X_add_genes_att_with_seq_suumary$MGI_site, "_[ARTBM]$")
) %>% unique() %>%
  group_by(labels_unique,labels_simplified) %>%
  summarise(breaks_x = mean(values_x)) %>%
  arrange(breaks_x)



ggplot(score_all_sites_all_dataset_df_top_X_add_genes_att_with_seq_suumary, aes(x=MGI_site_Aride_n2, y = perc, fill = AA, group = Aride, color = Aride)) +
  geom_bar(stat = "identity", width = 0.6 )  +
  theme_bw()  + facet_wrap(~dataset, scale = "free", strip.position="right", ncol = 1)    +
  theme(legend.key=element_blank(), 
   axis.text.x = element_text(colour = "black", size = 12, face = "bold", angle = 90, vjust = 0.3, hjust = 1)) +
  scale_x_continuous(breaks = df_values_breaks_MGI$breaks_x, expand = c(0.005,0.005) , labels =   df_values_breaks_MGI$labels_simplified) +
  scale_color_manual(values = c("darkgreen","darkorange3"))   + scale_fill_manual(values = palette_AA_LETTER)  +
  theme(legend.position="top") +guides(fill=guide_legend(nrow=1,byrow=TRUE),
                                       color=guide_legend(override.aes=list(fill=NA)))


ggsave(paste0(output_figures_site, "/top_",top_X,"_pos_per_dataset_AA_percentages.pdf"), limitsize = F, height = 4, width = 6+6 + top_X * 0.1)




```


```{r fig.height=5, fig.width=20}
#TODO

ggplot(score_all_sites_all_dataset_df_top_X_add_genes_att_with_seq, aes(x = MGI_site, group = dataset , y = -log10(aagtr_pval))) +
  geom_line(color ="grey") +
  geom_point() +
  theme_bw()  + facet_wrap(~dataset, scale = "free", strip.position="right", ncol = 1)  +
  scale_x_discrete(position = "top", breaks = score_all_sites_all_dataset_df_top_X_add_genes_att$MGI_site, labels = score_all_sites_all_dataset_df_top_X_add_genes_att$labels_simplified) + theme(legend.key=element_blank(), 
  axis.text.x = element_text(colour = "black", size = 12, face = "bold", angle = 90, vjust = 0.3, hjust = 0))

ggsave(paste0(output_figures_site, "/top_",top_X,"_pos_per_dataset_pelican_pval_per_site.pdf"), limitsize = F, height = 5, width = 12 + top_X * 0.2)


  
```


## Enrichment (GSEA)

```{r}
library(clusterProfiler)
library(ReactomePA)
library("org.Mm.eg.db")
#BiocManager::install("org.Mm.eg.db")
#BiocManager::install("ReactomePA")
```



```{r}
TERM2GENE_DEG_SIGN <-  gene_DE_score_df %>%
  filter(padj < 0.1)  %>%
  mutate(TERM=paste0(dataset,"_DEG_O1"))%>%
  dplyr::select(TERM,SYMBOL)

TERM2GENE_DEG_not_SIGN <-  gene_DE_score_df %>%
  filter(padj >= 0.1)  %>%
  mutate(TERM=paste0(dataset,"_NS_O1"))%>%
  dplyr::select(TERM,SYMBOL)

TERM2GENE_DEG_UP <-  gene_DE_score_df %>%
  filter(padj < 0.1, log2FoldChange > 0)  %>%
  mutate(TERM=paste0(dataset,"_DEG_O1_UP"))%>%
  dplyr::select(TERM,SYMBOL)

TERM2GENE_DEG_DOWN <-  gene_DE_score_df %>%
  filter(padj < 0.1, log2FoldChange < 0)  %>%
  mutate(TERM=paste0(dataset,"_DEG_O1_DOWN"))%>%
  dplyr::select(TERM,SYMBOL)

TERM2GENE_DEG <- bind_rows(TERM2GENE_DEG_SIGN,TERM2GENE_DEG_UP,TERM2GENE_DEG_DOWN,TERM2GENE_DEG_not_SIGN)
```






```{r}

MGI2ENTREZ <- biomaRt::getBM(
  attributes=c("mgi_symbol","entrezgene_id"),
  mart = mus_mart) %>%
  dplyr::filter(!is.na(entrezgene_id))

# Extract all Reactome pathways for mouse
mm_reactome <- ReactomePA::gson_Reactome(organism = "mouse")
mm_reactome_gsid2name <- mm_reactome@gsid2name

# Find the exact pathway ID for SLC-mediated transmembrane transport
pathway_name <- c("SLC-mediated transmembrane transport")
slc_pathway_id <- mm_reactome_gsid2name %>%
  filter(name %in% pathway_name )

slc_pathway_id_gene  <-mm_reactome@gsid2gene  %>%
  filter(gsid %in% slc_pathway_id )

# Get all mouse genes annotated to that Reactome pathway
slc_genes_entrez <-MGI2ENTREZ %>% filter(entrezgene_id %in% slc_pathway_id_gene$gene )


total_slc_genes_entrez <- slc_genes_entrez %>% filter(mgi_symbol %in% total_results_pelican_MGI_SYMBOL_ENTREZ_ok$SYMBOL )
murinae_slc_genes_entrez <- slc_genes_entrez %>% filter(mgi_symbol %in% murinae_results_pelican_MGI_SYMBOL_ENTREZ_ok$SYMBOL )


```


### prep Murinae 

```{r}
murinae_results_pelican_MGI_SYMBOL_ENTREZ_ok <- murinae_best_pval_per_MGI_ok %>% 
  left_join(MGI2ENTREZ_ok, by = c ("MGI"="mgi_symbol")) %>%
  mutate(SYMBOL=MGI,
         rank_gene = rank(best_aagtr_pval)) %>%
  dplyr::rename(ENTREZID="entrezgene_id") %>%
  relocate(SYMBOL,ENTREZID,best_aagtr_pval,rank_gene)

murinae_results_pelican_MGI_SYMBOL_ENTREZ_ok <- murinae_results_pelican_MGI_SYMBOL_ENTREZ_ok %>%
  mutate(for_oa = MGI %in% (murinae_nb_sites_per_ali %>% filter(!is_spike)%>% pull(MGI)),
          for_oa_th = log10_best_aagtr_pval <= (murinae_pelican_neutral_all_tested_results_summary %>% filter(q == 0.001) %>% pull(value))
)


score_to_use = "best_aagtr_pval"

murinae_rank_gene_score_df = murinae_results_pelican_MGI_SYMBOL_ENTREZ_ok %>%
    dplyr::rename(SCORE=!!rlang::parse_expr(score_to_use)) %>%
    filter(!is.na(SCORE))

murinae_gene_score_geneList_SYMBOL <- -log10(murinae_rank_gene_score_df$SCORE)
names(murinae_gene_score_geneList_SYMBOL) <- murinae_rank_gene_score_df$SYMBOL

murinae_gene_score_geneList_SYMBOL <- murinae_gene_score_geneList_SYMBOL[!is.na(murinae_gene_score_geneList_SYMBOL)]
murinae_gene_score_geneList_SYMBOL <- sort(murinae_gene_score_geneList_SYMBOL, decreasing = T)

length(murinae_gene_score_geneList_SYMBOL)

#for REACTOME use ENTREZID
murinae_rank_gene_score_df <- murinae_rank_gene_score_df %>% filter(!is.na(ENTREZID))
murinae_gene_score_geneList_with_ENTREZ <-   -log10(murinae_rank_gene_score_df$SCORE)
names(murinae_gene_score_geneList_with_ENTREZ) <- as.character(murinae_rank_gene_score_df$ENTREZID)

murinae_gene_score_geneList_with_ENTREZ <- murinae_gene_score_geneList_with_ENTREZ[!is.na(murinae_gene_score_geneList_with_ENTREZ)]
murinae_gene_score_geneList_with_ENTREZ <- sort(murinae_gene_score_geneList_with_ENTREZ, decreasing = T)

length(murinae_gene_score_geneList_with_ENTREZ)

```


### prep Total

```{r}

total_results_pelican_MGI_SYMBOL_ENTREZ_ok <- total_best_pval_per_MGI_ok %>% 
  left_join(MGI2ENTREZ_ok, by = c ("MGI"="mgi_symbol")) %>%
  mutate(SYMBOL=MGI,
         rank_gene = rank(best_aagtr_pval)) %>%
  dplyr::rename(ENTREZID="entrezgene_id") %>%
  relocate(SYMBOL,ENTREZID,best_aagtr_pval,rank_gene)
```


```{r}
total_results_pelican_MGI_SYMBOL_ENTREZ_ok <- total_best_pval_per_MGI_ok %>% 
  left_join(MGI2ENTREZ_ok, by = c ("MGI"="mgi_symbol")) %>%
  mutate(SYMBOL=MGI,
         rank_gene = rank(best_aagtr_pval)) %>%
  dplyr::rename(ENTREZID="entrezgene_id") %>%
  relocate(SYMBOL,ENTREZID,best_aagtr_pval,rank_gene)


total_results_pelican_MGI_SYMBOL_ENTREZ_ok <- total_results_pelican_MGI_SYMBOL_ENTREZ_ok %>%
  mutate(for_oa = MGI %in% (total_nb_sites_per_ali %>% filter(!is_spike)%>% pull(MGI)),
         for_oa_th = log10_best_aagtr_pval <= (total_pelican_neutral_all_tested_results_summary %>% filter(q == 0.001) %>% pull(value))
         )


score_to_use = "best_aagtr_pval"

total_rank_gene_score_df = total_results_pelican_MGI_SYMBOL_ENTREZ_ok %>%
    dplyr::rename(SCORE=!!rlang::parse_expr(score_to_use)) %>%
    filter(!is.na(SCORE))

total_gene_score_geneList_SYMBOL <- -log10(total_rank_gene_score_df$SCORE)
names(total_gene_score_geneList_SYMBOL) <- total_rank_gene_score_df$SYMBOL

total_gene_score_geneList_SYMBOL <- total_gene_score_geneList_SYMBOL[!is.na(total_gene_score_geneList_SYMBOL)]
total_gene_score_geneList_SYMBOL <- sort(total_gene_score_geneList_SYMBOL, decreasing = T)

length(total_gene_score_geneList_SYMBOL)

#for REACTOME use ENTREZID
total_rank_gene_score_df <- total_rank_gene_score_df %>% filter(!is.na(ENTREZID))
total_gene_score_geneList_with_ENTREZ <-   -log10(total_rank_gene_score_df$SCORE)
names(total_gene_score_geneList_with_ENTREZ) <- as.character(total_rank_gene_score_df$ENTREZID)

total_gene_score_geneList_with_ENTREZ <- total_gene_score_geneList_with_ENTREZ[!is.na(total_gene_score_geneList_with_ENTREZ)]
total_gene_score_geneList_with_ENTREZ <- sort(total_gene_score_geneList_with_ENTREZ, decreasing = T)

length(total_gene_score_geneList_with_ENTREZ)

```

### Total OA th

```{r}


total_oa_GO_BP_X <- enrichGO(gene     = total_results_pelican_MGI_SYMBOL_ENTREZ_ok %>% filter(for_oa_th) %>% pull(SYMBOL),
                       universe = total_results_pelican_MGI_SYMBOL_ENTREZ_ok %>%  pull(SYMBOL),
                       keyType = "SYMBOL",
                       OrgDb         = org.Mm.eg.db,
                       ont          = "BP",
                       minGSSize    = 10,
                       maxGSSize    = 500,
                       pvalueCutoff = 1,
                       qvalueCutoff = 1,
                       readable = T)

total_oa_GO_MF_X <-enrichGO(gene     = total_results_pelican_MGI_SYMBOL_ENTREZ_ok %>% filter(for_oa_th) %>% pull(SYMBOL),
                       universe = total_results_pelican_MGI_SYMBOL_ENTREZ_ok %>%  pull(SYMBOL),
                      keyType = "SYMBOL",
                      OrgDb         = org.Mm.eg.db,
                      ont          = "MF",
                      minGSSize    = 10,
                      maxGSSize    = 500,
                      pvalueCutoff = 1,
                      qvalueCutoff = 1,
                      readable = T)

total_oa_GO_CC_X <-enrichGO(gene     = total_results_pelican_MGI_SYMBOL_ENTREZ_ok %>% filter(for_oa_th) %>% pull(SYMBOL),
                       universe = total_results_pelican_MGI_SYMBOL_ENTREZ_ok %>%  pull(SYMBOL),
                      keyType = "SYMBOL",
                      OrgDb         = org.Mm.eg.db,
                      ont          = "CC",
                      minGSSize    = 10,
                      maxGSSize    = 500,
                      pvalueCutoff = 1,
                      qvalueCutoff = 1,
                      readable = T)
    
total_oa_GO_ALL_X <-enrichGO(gene     = total_results_pelican_MGI_SYMBOL_ENTREZ_ok %>% filter(for_oa_th) %>% pull(SYMBOL),
                       universe = total_results_pelican_MGI_SYMBOL_ENTREZ_ok %>%  pull(SYMBOL),
                      keyType = "SYMBOL",
                      OrgDb         = org.Mm.eg.db,
                      ont          = "ALL",
                      minGSSize    = 10,
                      maxGSSize    = 500,
                      pvalueCutoff = 1,
                      qvalueCutoff = 1,
                      readable = T)
  


total_oa_pathway_X <- enrichPathway(gene     = as.character(total_results_pelican_MGI_SYMBOL_ENTREZ_ok %>% filter(for_oa_th) %>% pull(ENTREZID)),
                              universe = as.character(total_results_pelican_MGI_SYMBOL_ENTREZ_ok %>% pull(ENTREZID)),
                              organism = "mouse",
                              minGSSize    = 10,
                              maxGSSize    = 500,
                              pAdjustMethod = "BH",
                              pvalueCutoff = 1,
                              qvalueCutoff = 1,
                              readable=T)


total_oa_DEG_custom_X <- enricher(gene =   total_results_pelican_MGI_SYMBOL_ENTREZ_ok %>% filter(for_oa_th) %>% pull(SYMBOL),
                          universe = total_results_pelican_MGI_SYMBOL_ENTREZ_ok %>%  pull(SYMBOL),
                          TERM2GENE = TERM2GENE_DEG %>%  filter(grepl("Total", TERM)),
                          minGSSize    = 10,
                          maxGSSize    = 500,
                          pAdjustMethod = "none",
                          pvalueCutoff = 1,
                          qvalueCutoff = 1)


```



```{r}
total_nb_GO_BP_sign_oa <- total_oa_GO_BP_X %>% as_tibble() %>% dplyr::filter(p.adjust < 0.1) %>% nrow()
total_nb_GO_MF_sign_oa <- total_oa_GO_MF_X %>% as_tibble() %>% dplyr::filter(p.adjust < 0.1) %>% nrow()
total_nb_GO_CC_sign_oa <- total_oa_GO_CC_X %>% as_tibble() %>% dplyr::filter(p.adjust < 0.1) %>% nrow()
total_nb_GO_ALL_sign_oa <- total_oa_GO_ALL_X %>% as_tibble() %>% dplyr::filter(p.adjust < 0.1) %>% nrow()
total_nb_GO_RA_sign_oa <- total_oa_pathway_X %>% as_tibble() %>% dplyr::filter(p.adjust < 0.1) %>% nrow()
total_nb_GS_DEG_oa <-  total_oa_DEG_custom_X %>% as_tibble() %>% dplyr::filter(p.adjust < 0.1) %>% nrow()


total_nb_signi_GS_oa_df <- tibble(DB = c("ALL GO","Reactome","DEG (no corection test multiple)"), #"BP","MF","CC",
       nb_GO_inf0.1 = c(#total_nb_GO_BP_sign,
                        #total_nb_GO_MF_sign,
                        #total_nb_GO_CC_sign,
                        total_nb_GO_ALL_sign_oa,
                        total_nb_GO_RA_sign_oa,
                        total_nb_GS_DEG_oa)) %>%
  mutate(dataset = "Total")


total_nb_signi_GS_oa_plot <- total_nb_signi_GS_oa_df %>%
  ggplot(aes(x=DB, y = nb_GO_inf0.1) ) +
  geom_col() +
  geom_text(aes(label =nb_GO_inf0.1,size = 4 ),  vjust = -1, colour = "black") + 
  theme_bw() +
  ggtitle("OA : Nb geneset signi padj < 0.1") +
  xlab("") +
  ylab("nb_inf0.1") +
  ylim(c(0,12)) +
  facet_wrap(~dataset)

total_nb_signi_GS_oa_plot
```


```{r}
scale_pval <- scale_fill_gradientn( limits = c(0,1), breaks = c(0,0.1,0.2,0.3,0.4,0.5,1), colors = c("deepskyblue4","deepskyblue2","grey","grey","grey","grey","grey","grey","grey","grey","grey"))
scale_size <- scale_size_continuous( range = c(3,10), breaks = c(10,30, 50,100,200))

# best  GS

if(total_oa_GO_BP_X %>% filter(p.adjust < 0.1) %>% nrow() > 0) {
  total_p_best_go_bp_oa <-
    dotplot( total_oa_GO_BP_X %>%filter(p.adjust < 0.1), orderBy = "pvalue", decreasing = T, label_format = 100, showCategory = 160) +
    scale_y_discrete(limits=rev) + facet_wrap(~ "Total_seq - (Pelican)"+ "GO BP") + ggtitle("OA")
} else {
  total_p_best_go_bp_oa <- NULL
}
total_p_best_go_bp_oa + scale_pval + scale_size

if(total_oa_GO_MF_X %>% filter(p.adjust < 0.1) %>% nrow() > 0) {
  total_p_best_go_mf_oa <-
    dotplot( total_oa_GO_MF_X %>%filter(p.adjust < 0.1), orderBy = "pvalue", decreasing = T, label_format = 100, showCategory = 160) +
    scale_y_discrete(limits=rev) + facet_wrap(~ "Total_seq - (Pelican)"+ "GO MF") + ggtitle("OA")
} else {
  total_p_best_go_mf_oa <- NULL
}
total_p_best_go_mf_oa + scale_pval + scale_size


if(total_oa_GO_CC_X %>% filter(p.adjust < 0.1) %>% nrow() > 0) {
  total_p_best_go_cc_oa <-
    dotplot( total_oa_GO_CC_X %>%filter(p.adjust < 0.1), orderBy = "pvalue", decreasing = T, label_format = 100, showCategory = 160) +
    scale_y_discrete(limits=rev) + facet_wrap(~ "Total_seq - (Pelican)"+ "GO CC") + ggtitle("OA")
} else {
  total_p_best_go_cc_oa <- NULL
}
total_p_best_go_cc_oa + scale_pval + scale_size



if(total_oa_GO_ALL_X %>% filter(p.adjust < 0.1) %>% nrow() > 0) {
  total_p_best_go_all_oa <- dotplot( total_oa_GO_ALL_X %>% filter(p.adjust < 0.1),orderBy = "pvalue", decreasing = T, label_format = 100, showCategory = 160) +
    scale_y_discrete(limits=rev) + facet_wrap(~ "Total_seq - (Pelican) "+ ONTOLOGY, scale = "free", ncol = 1, strip.position = "right" ) + ggtitle("OA")
} else {
  total_p_best_go_all_oa <- NULL
}

total_p_best_go_all_oa+ scale_pval + scale_size

if(total_oa_pathway_X %>% filter(p.adjust < 0.1) %>% nrow() > 0) {
total_p_best_go_ra_oa <- dotplot( total_oa_pathway_X %>% filter(p.adjust < 0.1),orderBy = "pvalue", decreasing = T, label_format = 100, showCategory = 160,) +
  scale_y_discrete(limits=rev) + facet_wrap(~ "Total_seq - (Pelican)"+ "REACTOME") + ggtitle("OA")
} else {
  total_p_best_go_ra_oa <- NULL
}
total_p_best_go_ra_oa + scale_pval + scale_size

if(total_oa_DEG_custom_X %>% filter(p.adjust < 0.1) %>% nrow() > 0) {
total_p_best_DEG_oa <- dotplot( total_oa_DEG_custom_X%>% filter(p.adjust < 0.1), orderBy = "pvalue", decreasing = T, label_format = 100, showCategory = 160,) +
  scale_y_discrete(limits=rev) + facet_wrap(~ "Total_seq - (Pelican)"+ "DEG No correction test multiple")  + ggtitle("OA")
} else {
  total_p_best_DEG_oa <- NULL
}
total_p_best_DEG_oa + scale_pval + scale_size


```




```{r,fig.width=12, fig.height=5}
total_oa_gene_DE_score_df_simplify <-total_results_pelican_MGI_SYMBOL_ENTREZ_ok %>%
  left_join(gene_DE_score_df %>% dplyr::select(SYMBOL,dataset,is_DE,log2FoldChange) %>% filter(dataset == "Total")) %>%
  dplyr::select(SYMBOL,ENTREZID, rank_gene,dataset,log2FoldChange,is_DE) %>%
  filter(is_DE =="FDR<0.1") %>%
  arrange(rank_gene)  %>%
  mutate(UP_DOWN = ifelse(log2FoldChange >0, "up","down"),
         UP_DOWN = ifelse(is_DE == "FDR<0.1" , toupper(UP_DOWN),UP_DOWN)) %>%
  mutate(dataset="DEG Total\nFDR<0.1",
         ONTOLOGY = "DEG")


  write.table(total_oa_GO_ALL_X@result %>% filter(p.adjust<0.1) , file = "GO_clusterProfiler_PELICAN_total.tsv",
              sep = "\t", quote = FALSE)
  
  
total_oa_GO_ALL_X_ok <-  total_oa_GO_ALL_X@result %>% filter(p.adjust<0.05) 
GOID2NAME_list <- total_oa_GO_ALL_X@geneSets[total_oa_GO_ALL_X_ok %>% pull(ID)]

df_GOID2NAME <- bind_rows(lapply(names(GOID2NAME_list), function(go_id) {
  data.frame(ID = go_id, SYMBOL = GOID2NAME_list[[go_id]], stringsAsFactors = FALSE)
}))  %>%
  left_join(total_oa_GO_ALL_X_ok %>% dplyr::select(ID,ONTOLOGY,Description))



total_oa_gene_DE_score_symporter_df_simplify <- total_results_pelican_MGI_SYMBOL_ENTREZ_ok %>%
  left_join(gene_DE_score_df %>% dplyr::select(SYMBOL,dataset,is_DE,log2FoldChange) %>% filter(dataset == "Total") ) %>%
  dplyr::select(SYMBOL,ENTREZID, rank_gene,dataset,log2FoldChange,is_DE) %>%
  filter(SYMBOL %in% (df_GOID2NAME %>% filter(ONTOLOGY == "MF") %>% pull(SYMBOL))) %>%
  left_join(df_GOID2NAME %>% filter(ONTOLOGY == "MF") %>% dplyr::select(SYMBOL,Description,ONTOLOGY)) %>%
  mutate(UP_DOWN = ifelse(log2FoldChange >0, "up","down"),
         UP_DOWN = ifelse(is_DE == "FDR<0.1" , toupper(UP_DOWN),"ns"),
          dataset = str_wrap(Description, width = 20)) %>%
  arrange(rank_gene) %>%
  group_by(SYMBOL) %>%
  mutate(n=n(),
         n = ifelse(SYMBOL %in% total_oa_gene_DE_score_df_simplify$SYMBOL, n+1,n))


total_oa_gene_SLC_DE_score_df_l = bind_rows(total_oa_gene_DE_score_df_simplify,total_oa_gene_DE_score_symporter_df_simplify) %>%
 mutate(dataset = factor(dataset),
       dataset = fct_reorder2(dataset,.x=dataset,.y=ONTOLOGY),
       dataset = fct_relevel(dataset,"DEG Total\nFDR<0.1", after =Inf))



total_nb_gene_inf_th <- total_nb_gene_per_cat_th %>% filter(log10_aagtr_pval_cut_th=="inf_th") %>% pull(cumul_nb_gene) 

total_plot_DE_cross_pelican_oa <- ggplot(total_oa_gene_SLC_DE_score_df_l, aes(x=rank_gene, y = dataset)) +
  geom_line(aes(group = SYMBOL),  color = "grey80", linetype = "dashed")+
  geom_point(aes(shape = UP_DOWN, color = UP_DOWN, fill=UP_DOWN), size = 2) +
  scale_color_manual(values = c("UP"="orange","DOWN"="darkgreen","up"="#f9da78","down"="#a7d4a6","ns" = "grey80"), name = "DE gene") +
   scale_fill_manual(values = c("UP"="orange","DOWN"="darkgreen","up"="#f9da78","down"="#a7d4a6","ns" = "grey80"), name = "DE gene") +
  scale_shape_manual(values = c("UP"=24,"DOWN"=25,"up"=24,"down"=25, "ns" = 20), name = "DE gene") +
  theme_bw() +
  ggrepel::geom_text_repel(data = total_oa_gene_DE_score_df_simplify %>% filter(rank_gene <= total_nb_gene_inf_th), aes(label = SYMBOL, color = UP_DOWN ),
                           min.segment.length = 0,  max.overlaps = Inf,
                           nudge_y = ifelse(total_oa_gene_DE_score_df_simplify %>% filter(rank_gene <= total_nb_gene_inf_th) %>% pull(UP_DOWN) == "UP", 0.3, -0.3), 
                           segment.alpha = 0.4, size = 4, show.legend = F) +
  
    ggrepel::geom_text_repel(data = total_oa_gene_DE_score_symporter_df_simplify  %>% filter(rank_gene <= total_nb_gene_inf_th, n > 1), aes(label = SYMBOL, color = UP_DOWN ),
                           min.segment.length = 0, force = 20,   max.overlaps = Inf,
                           nudge_y = ifelse(total_oa_gene_DE_score_symporter_df_simplify %>% filter(rank_gene < total_nb_gene_inf_th) %>% pull(UP_DOWN) %in% c("UP","up"), 0.3, -0.3), 
                           segment.alpha = 0.4, size = 4, show.legend = F) +
  
  theme(axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 10),
        axis.ticks.length = unit(0.3, "cm"),
        axis.title.x = element_text(size = 10),
        legend.position = "top") +
  ylab("") +
  xlab("Pelican gene rank (Total_seq)")  +
    geom_vline(xintercept = total_nb_gene_inf_th, color = "#145db6" )

total_plot_DE_cross_pelican_oa
```



```{r, fig.width=16, fig.height=10}
total_plot_gs_oa <- cowplot::plot_grid(
  cowplot::plot_grid(
  cowplot::plot_grid(total_nb_signi_GS_oa_plot + theme(legend.position = "none"),  
                     total_p_best_DEG_oa + scale_pval + scale_size, ncol = 1),
  total_p_best_go_all_oa+ scale_pval + scale_size,ncol = 2, rel_widths = c(2,3)),
  total_plot_DE_cross_pelican_oa,
  ncol = 1, rel_heights = c(2,2,1)
)

total_plot_gs_oa

ggsave(plot=total_plot_gs_oa, filename = str_c(output_figures,"E3_pelican_oa_total.pdf"), width = 16, height = 10)

```

### Murinae OA th

```{r}
murinae_oa_GO_BP_X <- enrichGO(gene     = murinae_results_pelican_MGI_SYMBOL_ENTREZ_ok %>% filter(for_oa_th) %>% pull(SYMBOL),
                       universe = murinae_results_pelican_MGI_SYMBOL_ENTREZ_ok %>%  pull(SYMBOL),
                       keyType = "SYMBOL",
                       OrgDb         = org.Mm.eg.db,
                       ont          = "BP",
                       minGSSize    = 10,
                       maxGSSize    = 500,
                       pvalueCutoff = 1,
                       qvalueCutoff = 1,
                       readable = T)

murinae_oa_GO_MF_X <-enrichGO(gene     = murinae_results_pelican_MGI_SYMBOL_ENTREZ_ok %>% filter(for_oa_th) %>% pull(SYMBOL),
                       universe = murinae_results_pelican_MGI_SYMBOL_ENTREZ_ok %>%  pull(SYMBOL),
                      keyType = "SYMBOL",
                      OrgDb         = org.Mm.eg.db,
                      ont          = "MF",
                      minGSSize    = 10,
                      maxGSSize    = 500,
                      pvalueCutoff = 1,
                      qvalueCutoff = 1,
                      readable = T)

murinae_oa_GO_CC_X <-enrichGO(gene     = murinae_results_pelican_MGI_SYMBOL_ENTREZ_ok %>% filter(for_oa_th) %>% pull(SYMBOL),
                       universe = murinae_results_pelican_MGI_SYMBOL_ENTREZ_ok %>%  pull(SYMBOL),
                      keyType = "SYMBOL",
                      OrgDb         = org.Mm.eg.db,
                      ont          = "CC",
                      minGSSize    = 10,
                      maxGSSize    = 500,
                      pvalueCutoff = 1,
                      qvalueCutoff = 1,
                      readable = T)
    
murinae_oa_GO_ALL_X <-enrichGO(gene     = murinae_results_pelican_MGI_SYMBOL_ENTREZ_ok %>% filter(for_oa_th) %>% pull(SYMBOL),
                       universe = murinae_results_pelican_MGI_SYMBOL_ENTREZ_ok %>%  pull(SYMBOL),
                      keyType = "SYMBOL",
                      OrgDb         = org.Mm.eg.db,
                      ont          = "ALL",
                      minGSSize    = 10,
                      maxGSSize    = 500,
                      pvalueCutoff = 1,
                      qvalueCutoff = 1,
                      readable = T)
  


murinae_oa_pathway_X <- enrichPathway(gene     = as.character(murinae_results_pelican_MGI_SYMBOL_ENTREZ_ok %>% filter(for_oa_th) %>% pull(ENTREZID)),
                              universe = as.character(murinae_results_pelican_MGI_SYMBOL_ENTREZ_ok %>% pull(ENTREZID)),
                              organism = "mouse",
                              minGSSize    = 10,
                              maxGSSize    = 500,
                              pAdjustMethod = "BH",
                              pvalueCutoff = 1,
                              qvalueCutoff = 1,
                              readable=T)

murinae_oa_DEG_custom_X <- enricher(gene =   murinae_results_pelican_MGI_SYMBOL_ENTREZ_ok %>% filter(for_oa_th) %>% pull(SYMBOL),
                          universe = murinae_results_pelican_MGI_SYMBOL_ENTREZ_ok %>%  pull(SYMBOL),
                          TERM2GENE = TERM2GENE_DEG %>%  filter(grepl("Murinae", TERM)),
                          minGSSize    = 10,
                          maxGSSize    = 500,
                          pvalueCutoff = 1,
                          qvalueCutoff = 1)

```



```{r}
murinae_nb_GO_BP_sign_oa <- murinae_oa_GO_BP_X %>% as_tibble() %>% dplyr::filter(p.adjust < 0.1) %>% nrow()
murinae_nb_GO_MF_sign_oa <- murinae_oa_GO_MF_X %>% as_tibble() %>% dplyr::filter(p.adjust < 0.1) %>% nrow()
murinae_nb_GO_CC_sign_oa <- murinae_oa_GO_CC_X %>% as_tibble() %>% dplyr::filter(p.adjust < 0.1) %>% nrow()
murinae_nb_GO_ALL_sign_oa <- murinae_oa_GO_ALL_X %>% as_tibble() %>% dplyr::filter(p.adjust < 0.1) %>% nrow()
murinae_nb_GO_RA_sign_oa <- murinae_oa_pathway_X %>% as_tibble() %>% dplyr::filter(p.adjust < 0.1) %>% nrow()
murinae_nb_GS_DEG_oa <-  murinae_oa_DEG_custom_X %>% as_tibble() %>% dplyr::filter(p.adjust < 0.1) %>% nrow()


murinae_nb_signi_GS_oa_df <- tibble(DB = c("ALL GO","Reactome","DEG"), #"BP","MF","CC",
       nb_GO_inf0.1 = c(#murinae_nb_GO_BP_sign,
                        #murinae_nb_GO_MF_sign,
                        #murinae_nb_GO_CC_sign,
                        murinae_nb_GO_ALL_sign_oa,
                        murinae_nb_GO_RA_sign_oa,
                        murinae_nb_GS_DEG_oa)) %>%
  mutate(dataset = "Murinae")


murinae_nb_signi_GS_oa_plot <- murinae_nb_signi_GS_oa_df %>%
  ggplot(aes(x=DB, y = nb_GO_inf0.1) ) +
  geom_col() +
  geom_text(aes(label =nb_GO_inf0.1,size = 4 ),  vjust = -1, colour = "black") + 
  theme_bw() +
  ggtitle("OA : Nb geneset signi padj < 0.1") +
  xlab("") +
  ylab("nb_inf0.1") +
  ylim(c(0,12)) +
  facet_wrap(~dataset)

murinae_nb_signi_GS_oa_plot
```


```{r}
scale_pval <- scale_fill_gradientn( limits = c(0,1), breaks = c(0,0.1,0.2,0.3,0.4,0.5,1), colors = c("deepskyblue4","deepskyblue2","grey","grey","grey","grey","grey","grey","grey","grey","grey"))
scale_size <- scale_size_continuous( range = c(3,10), breaks = c(10,30, 50,100,200))

# best  GS

if(murinae_oa_GO_BP_X %>% filter(p.adjust < 0.1) %>% nrow() > 0) {
  murinae_p_best_go_bp_oa <-
    dotplot( murinae_oa_GO_BP_X %>%filter(p.adjust < 0.1), orderBy = "pvalue", decreasing = T, label_format = 100, showCategory = 160) +
    scale_y_discrete(limits=rev) + facet_wrap(~ "Murinae_seq - (Pelican)"+ "GO BP") + ggtitle("OA")
} else {
  murinae_p_best_go_bp_oa <- NULL
}
murinae_p_best_go_bp_oa + scale_pval + scale_size

if(murinae_oa_GO_MF_X %>% filter(p.adjust < 0.1) %>% nrow() > 0) {
  murinae_p_best_go_mf_oa <-
    dotplot( murinae_oa_GO_MF_X %>%filter(p.adjust < 0.1), orderBy = "pvalue", decreasing = T, label_format = 100, showCategory = 160) +
    scale_y_discrete(limits=rev) + facet_wrap(~ "Murinae_seq - (Pelican)"+ "GO MF") + ggtitle("OA")
} else {
  murinae_p_best_go_mf_oa <- NULL
}
murinae_p_best_go_mf_oa + scale_pval + scale_size


if(murinae_oa_GO_CC_X %>% filter(p.adjust < 0.1) %>% nrow() > 0) {
  murinae_p_best_go_cc_oa <-
    dotplot( murinae_oa_GO_CC_X %>%filter(p.adjust < 0.1), orderBy = "pvalue", decreasing = T, label_format = 100, showCategory = 160) +
    scale_y_discrete(limits=rev) + facet_wrap(~ "Murinae_seq - (Pelican)"+ "GO CC") + ggtitle("OA")
} else {
  murinae_p_best_go_cc_oa <- NULL
}
murinae_p_best_go_cc_oa + scale_pval + scale_size



if(murinae_oa_GO_ALL_X %>% filter(p.adjust < 0.1) %>% nrow() > 0) {
  murinae_p_best_go_all_oa <- dotplot( murinae_oa_GO_ALL_X %>% filter(p.adjust < 0.1),orderBy = "pvalue", decreasing = T, label_format = 100, showCategory = 160) +
    scale_y_discrete(limits=rev) + facet_wrap(~ "Murinae_seq - (Pelican) "+ ONTOLOGY, scale = "free", ncol = 1, strip.position = "right" ) + ggtitle("OA")
} else {
  murinae_p_best_go_all_oa <- NULL
}

murinae_p_best_go_all_oa+ scale_pval + scale_size

if(murinae_oa_pathway_X %>% filter(p.adjust < 0.1) %>% nrow() > 0) {
murinae_p_best_go_ra_oa <- dotplot( murinae_oa_pathway_X %>% filter(p.adjust < 0.1),orderBy = "pvalue", decreasing = T, label_format = 100, showCategory = 160,) +
  scale_y_discrete(limits=rev) + facet_wrap(~ "Murinae_seq - (Pelican)"+ "REACTOME") + ggtitle("OA")
} else {
  murinae_p_best_go_ra_oa <- NULL
}
murinae_p_best_go_ra_oa + scale_pval + scale_size

if(murinae_oa_DEG_custom_X %>% filter(p.adjust < 0.1) %>% nrow() > 0) {
murinae_p_best_DEG_oa <- dotplot( murinae_oa_DEG_custom_X%>% filter(p.adjust < 0.1), orderBy = "pvalue", decreasing = T, label_format = 100, showCategory = 160,) +
  scale_y_discrete(limits=rev) + facet_wrap(~ "Murinae_seq - (Pelican)"+ "DEG")  + ggtitle("OA")
} else {
  murinae_p_best_DEG_oa <- NULL
}
murinae_p_best_DEG_oa + scale_pval + scale_size


```






```{r,fig.width=12, fig.height=5}
murinae_oa_gene_DE_score_df_simplify <-murinae_results_pelican_MGI_SYMBOL_ENTREZ_ok %>%
  left_join(gene_DE_score_df %>% dplyr::select(SYMBOL,dataset,is_DE,log2FoldChange) %>% filter(dataset == "Murinae")) %>%
  dplyr::select(SYMBOL,ENTREZID, rank_gene,dataset,log2FoldChange,is_DE) %>%
  filter(is_DE =="FDR<0.1") %>%
  arrange(rank_gene)  %>%
  mutate(UP_DOWN = ifelse(log2FoldChange >0, "up","down"),
         UP_DOWN = ifelse(is_DE == "FDR<0.1" , toupper(UP_DOWN),UP_DOWN)) %>%
  mutate(dataset="DEG Murinae\nFDR<0.1",
         ONTOLOGY = "DEG") 


write.table(murinae_oa_GO_ALL_X@result %>% filter(p.adjust<0.1) , file = "GO_clusterProfiler_PELICAN_murinae.tsv",
              sep = "\t", quote = FALSE)
  
murinae_oa_GO_ALL_X_ok <-  murinae_oa_GO_ALL_X@result %>% filter(p.adjust<0.05)  %>% slice_min(pvalue, n = 1)
murinae_oa_GOID2NAME_list <- murinae_oa_GO_ALL_X@geneSets[murinae_oa_GO_ALL_X_ok %>% pull(ID)]

murinae_oa_df_GOID2NAME <- bind_rows(lapply(names(murinae_oa_GOID2NAME_list), function(go_id) {
  data.frame(ID = go_id, SYMBOL = murinae_oa_GOID2NAME_list[[go_id]], stringsAsFactors = FALSE)
}))  %>%
  left_join(murinae_oa_GO_ALL_X_ok %>% dplyr::select(ID,ONTOLOGY,Description))



murinae_oa_gene_DE_score_symporter_df_simplify <- murinae_results_pelican_MGI_SYMBOL_ENTREZ_ok %>%
  left_join(gene_DE_score_df %>% dplyr::select(SYMBOL,dataset,is_DE,log2FoldChange) %>% filter(dataset == "Murinae") ) %>%
  dplyr::select(SYMBOL,ENTREZID, rank_gene,dataset,log2FoldChange,is_DE) %>%
  filter(SYMBOL %in% (murinae_oa_df_GOID2NAME %>% filter(ONTOLOGY %in%  c("CC","MF","BP")) %>% pull(SYMBOL))) %>%
  left_join(murinae_oa_df_GOID2NAME %>% filter(ONTOLOGY %in%  c("CC","MF","BP")) %>% dplyr::select(SYMBOL,Description,ONTOLOGY)) %>%
  mutate(UP_DOWN = ifelse(log2FoldChange >0, "up","down"),
         UP_DOWN = ifelse(is_DE == "FDR<0.1" , toupper(UP_DOWN),"ns"),
         dataset = str_wrap(Description, width = 20)) %>%
  arrange(rank_gene) %>%
  group_by(SYMBOL) %>%
  mutate(n=n(),
         n = ifelse(SYMBOL %in% murinae_gene_DE_score_df_simplify$SYMBOL, n+1,n))


murinae_oa_gene_DE_score_df_simplify <- murinae_oa_gene_DE_score_df_simplify %>%
  left_join(murinae_oa_gene_DE_score_symporter_df_simplify %>% dplyr::select(SYMBOL,n) %>% unique())

murinae_oa_gene_SLC_DE_score_df_l = bind_rows(murinae_oa_gene_DE_score_df_simplify,murinae_oa_gene_DE_score_symporter_df_simplify) %>%
  mutate(dataset = factor(dataset),
          dataset = fct_reorder2(dataset,.x=dataset,.y=ONTOLOGY),
          dataset = fct_relevel(dataset,"DEG Murinae\nFDR<0.1", after =Inf))



murinae_nb_gene_inf_th <- murinae_nb_gene_per_cat_th %>% filter(log10_aagtr_pval_cut_th=="inf_th") %>% pull(cumul_nb_gene) 

murinae_plot_DE_cross_pelican_oa <- ggplot(murinae_oa_gene_SLC_DE_score_df_l, aes(x=rank_gene, y = dataset)) +
  geom_line(aes(group = SYMBOL),  color = "grey80", linetype = "dashed")+
  geom_point(aes(shape = UP_DOWN, color = UP_DOWN, fill=UP_DOWN), size = 2) +
  scale_color_manual(values = c("UP"="orange","DOWN"="darkgreen","up"="#f9da78","down"="#a7d4a6", "ns" = "grey80"), name = "DE gene") +
   scale_fill_manual(values = c("UP"="orange","DOWN"="darkgreen","up"="#f9da78","down"="#a7d4a6", "ns" = "grey80"), name = "DE gene") +
  scale_shape_manual(values = c("UP"=24,"DOWN"=25,"up"=24,"down"=25, "ns"= 20), name = "DE gene") +
   theme_bw() +
  ggrepel::geom_text_repel(data = murinae_oa_gene_DE_score_df_simplify %>% filter(rank_gene <= murinae_nb_gene_inf_th,ONTOLOGY %in%  c("DEG","MF","BP"), n >1), aes(label = SYMBOL, color = UP_DOWN ),
                           min.segment.length = 0, force = 30,  max.overlaps = Inf,
                           nudge_y = ifelse(murinae_oa_gene_DE_score_df_simplify %>% filter(rank_gene <= murinae_nb_gene_inf_th,ONTOLOGY %in%  c("DEG","MF","BP"), n >1) %>% pull(UP_DOWN) == "UP", 0.5, -0.5), 
                           segment.alpha = 0.4, size = 5, show.legend = F) +
  
  theme(axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 10),
        axis.ticks.length = unit(0.3, "cm"),
        axis.title.x = element_text(size = 10),
        legend.position = "top") +
  ylab("") +
  xlab("Pelican gene rank (Murinae_seq)")  +
    geom_vline(xintercept = murinae_nb_gene_inf_th, color = "#145db6" ) 

murinae_plot_DE_cross_pelican_oa



murinae_plot_DE_cross_pelican_oa_small <- ggplot(murinae_oa_gene_SLC_DE_score_df_l %>% filter(ONTOLOGY %in%  c("DEG","MF","BP")), aes(x=rank_gene, y = dataset)) +
  geom_line(aes(group = SYMBOL),  color = "grey80", linetype = "dashed")+
  geom_point(aes(shape = UP_DOWN, color = UP_DOWN, fill=UP_DOWN), size = 2) +
  scale_color_manual(values = c("UP"="orange","DOWN"="darkgreen","up"="#f9da78","down"="#a7d4a6"), name = "DE gene") +
   scale_fill_manual(values = c("UP"="orange","DOWN"="darkgreen","up"="#f9da78","down"="#a7d4a6"), name = "DE gene") +
  scale_shape_manual(values = c("UP"=24,"DOWN"=25,"up"=24,"down"=25, "ns"= 1), name = "DE gene") +
  theme_bw() +
  ggrepel::geom_text_repel(data = murinae_oa_gene_DE_score_df_simplify %>% filter(rank_gene <= murinae_nb_gene_inf_th,ONTOLOGY %in%  c("DEG","MF","BP"), n >1), aes(label = SYMBOL, color = UP_DOWN ),
                           min.segment.length = 0, force = 20, max.overlaps = 10,
                           nudge_y = ifelse(murinae_gene_DE_score_df_simplify %>% filter(rank_gene <= murinae_nb_gene_inf_th) %>% pull(UP_DOWN) == "UP", 0.3, -0.3), 
                           segment.alpha = 0.4, size = 2, show.legend = F) +
  
    ggrepel::geom_text_repel(data = murinae_oa_gene_DE_score_symporter_df_simplify  %>% filter(rank_gene <= murinae_nb_gene_inf_th,ONTOLOGY %in%  c("DEG","MF","BP"), n >1), aes(label = SYMBOL, color = UP_DOWN ),
                           min.segment.length = 0, force = 20, max.overlaps = 10,
                           nudge_y = ifelse(murinae_oa_gene_DE_score_symporter_df_simplify %>% filter(rank_gene < murinae_nb_gene_inf_th) %>% pull(UP_DOWN) %in% c("UP","up"), 0.3, -0.3), 
                           segment.alpha = 0.4, size = 2, show.legend = F) +
  
  theme(axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 12),
        axis.ticks.length = unit(0.3, "cm"),
        axis.title.x = element_text(size = 12),
        legend.position = "top") +
  ylab("") +
  xlab("Pelican gene rank (Murinae dataset)")  +
    geom_vline(xintercept = murinae_nb_gene_inf_th, color = "#145db6" , linetype = "dashed") 

murinae_plot_DE_cross_pelican_oa_small


```



```{r, fig.width=16, fig.height=10}
murinae_plot_gs_oa <- cowplot::plot_grid(
  cowplot::plot_grid(
  cowplot::plot_grid(murinae_nb_signi_GS_oa_plot + theme(legend.position = "none"),  
                     murinae_p_best_DEG_oa + scale_pval + scale_size, ncol = 1),
  murinae_p_best_go_all_oa+ scale_pval + scale_size,ncol = 2, rel_widths = c(2,3)),
    murinae_plot_DE_cross_pelican_oa_small,
  murinae_plot_DE_cross_pelican_oa,
  murinae_plot_DE_cross_pelican_line,
  ncol = 1, rel_heights = c(4,2,4,2)
)

murinae_plot_gs_oa

ggsave(plot=murinae_plot_gs_oa, filename = str_c(output_figures,"E4_pelican_oa_murinae.pdf"), width = 16, height = 15)

```



```{r}
total_plot_DE_cross_pelican_oa / murinae_plot_DE_cross_pelican_oa + geom_point(data = tibble(dataset ="", rank_gene = 1, color = "white"))
ggsave(plot= plot_grid(total_plot_DE_cross_pelican_oa,
                       murinae_plot_DE_cross_pelican_oa+ geom_point(data = tibble(dataset ="", rank_gene = 1, fill = "white")), ncol = 1, rel_heights = c(5,5)), filename = str_c(output_figures,"E6_pelican_oa_murinae_total.pdf"), width = 10, height = 6)

```




```{r}
total_p_best_go_all_oa + scale_pval + scale_size



ggsave(plot=  plot_grid(total_random_spike_violin,
                        plot_grid(total_p_best_go_all_oa + scale_pval + scale_size + ggtitle("Total_seq"),  murinae_p_best_go_all_oa+ scale_pval + scale_size + ggtitle("Murinae_seq"), ncol = 2, rel_heights = c(5,5), labels = c("B","C")),
                        ncol = 1, rel_heights = c(5,5), labels = c("A","")),
       filename = str_c(output_figures,"E7_pelican_oa_murinae_total.pdf"), width = 14, height = 8)


```
