---
title: "Simulation Expression_master"

output:
  html_document:
    toc: true
    toc_float: true
    number_sections: true
  pdf_document: default
---


```{r setup, include=T}
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE, cache = F)
print(getwd())
source("Expression_master_lib_simu.R")
resampling= T
nRand = 1000
test = F
run_whole_dataset = T
nCPU = 18
```

# Import dataset
```{r data import}
Alltable_onl = read.csv('../datalocalandonline_cds_counts_2025_final.tsv', h=T, sep = "\t", row.names=1) #15790   125
dim(Alltable_onl)


colnames(Alltable_onl) = gsub("loc_cds_SHPC43_Mastomys_kollmanspergeri", "loc_cds_SHPC43_Mastomys_kollmannspergeri", colnames(Alltable_onl))
colnames(Alltable_onl) = gsub("loc_cds_SHPC74_Mastomys_kollmanspergeri", "loc_cds_SHPC74_Mastomys_kollmannspergeri", colnames(Alltable_onl))
colnames(Alltable_onl) = gsub("loc_cds_SHPC75_Mastomys_kollmanspergeri", "loc_cds_SHPC75_Mastomys_kollmannspergeri", colnames(Alltable_onl))
colnames(Alltable_onl) = gsub("loc_cds_SHPC76_Mastomys_kollmanspergeri", "loc_cds_SHPC76_Mastomys_kollmannspergeri", colnames(Alltable_onl))
colnames(Alltable_onl) = gsub("SHPC46_Arvicanthis_niloticus_lmo", "SHPC46_Mastomys_natalensis", colnames(Alltable_onl)) #change Ani_lmo to Mna even if not using in expression)
colnames(Alltable_onl) = gsub("onl5_cds_SRR17216323_Peromyscus_maniculatus", "onl_cds_SRR17216323_Peromyscus_maniculatus", colnames(Alltable_onl))
colnames(Alltable_onl) = gsub("onl5_cds_SRR17216333_Peromyscus_maniculatus", "onl_cds_SRR17216333_Peromyscus_maniculatus", colnames(Alltable_onl))
colnames(Alltable_onl) = gsub("onl5_cds_SRR17216352_Peromyscus_maniculatus", "onl_cds_SRR17216352_Peromyscus_maniculatus", colnames(Alltable_onl))


#rownames(Alltable_onl) = Alltable_onl$Row.names
#Alltable_tpm = read.csv('datalocalandonline_cds_abundance_2020.tsv', h=T, sep = "\t", row.names=1)

# coldata
coldata_2025 <- read.csv("../coldata_July2025.tsv", sep = "\t")
coldata_simu = read.csv('coldata_localandonline_info20220630.tsv', h=T, sep = "\t") # 123  34

coldata_ok <- coldata_2025 %>% left_join(coldata_simu %>%
                                           dplyr::select(ID, total_permutation, recent_trans_permutation, ancient_trans_permutation, murinae_permutation),  by = "ID") %>%
 mutate(total_permutation = ifelse(species == "Peromyscus_maniculatus","g3",total_permutation),
        recent_trans_permutation = ifelse(species == "Peromyscus_maniculatus", "g9",recent_trans_permutation)
                                          )



# coldata SET TOTAL
coldata_tot = coldata_ok[coldata_ok$total_set == "yes",] # 80 32
dim(coldata_tot)
all(coldata_tot$ID %in% colnames(Alltable_onl))

# coldata SET ANCIENT TRANSITION
coldata_anc = coldata_ok[coldata_ok$ancient_set == "yes",] # 31 32
dim(coldata_anc)
all(coldata_anc$ID %in% colnames(Alltable_onl))

# coldata SET RECENT TRANSITION
coldata_rec = coldata_ok[coldata_ok$recent_set == "yes",] # 32 32
dim(coldata_rec)
all(coldata_rec$ID %in% colnames(Alltable_onl))

# coldata SET MURINAE
coldata_murinae= coldata_ok[coldata_ok$muridae_set == "yes",] # 29 32
dim(coldata_murinae)
all(coldata_murinae$ID %in% colnames(Alltable_onl))

```

# overview of the dataset



# Dataset expression 2022, SET ANCIENT TRANSITION , design : ~ family + cond_season


## ~ family + env (family in sva)

```{r, eval = TRUE}
output_dir="2025_simu_dataset_ancient_season/"
system(paste0("mkdir -p " , output_dir))

data_expr_anc_season = list(var1_design = "ancient_group_batch", var2_design = "cond_season", Alltable = Alltable_onl, coldata = coldata_anc, output_dir=output_dir, resampling = resampling , nRand = nRand, nCPU = nCPU, test = test, column_permutation = "ancient_trans_permutation", dataset_name="ancient_trans")
```

```{r, fig.path = output_dir, eval = TRUE, fig.height=20, fig.width=10}

data_expr_anc_season_res_simu = make_exp_analyses_only_simu(data_expr_anc_season, run_whole_dataset = run_whole_dataset) 
data_expr_anc_season_res_simu %>% subset(couple == "All_data") %>%  kbl() %>% kable_styling()

```

```{r, fig.path = output_dir, eval = TRUE, fig.height=20, fig.width=10}

data_expr_anc_season_res_simu_r = make_exp_analyses_only_simu(data_expr_anc_season, run_whole_dataset = F, replace_signi_genes = T)
data_expr_anc_season_res_simu_r %>% subset(couple == "All_data") %>%  kbl() %>% kable_styling()

```


# Dataset expression 2022, SET RECENT TRANSITION , design : ~ family + cond_season

## ~ family + env (family in sva)

```{r, eval = TRUE}
output_dir="2025_simu_dataset_recent_season/"
system(paste0("mkdir -p " , output_dir))

data_expr_rec_season = list(var1_design = "recent_group_batch", var2_design = "cond_season", Alltable = Alltable_onl, coldata = coldata_rec, output_dir=output_dir, resampling = resampling , nRand = nRand, nCPU = nCPU, test = test, column_permutation = "recent_trans_permutation",  dataset_name="recent_trans")

```

```{r, fig.path = output_dir, eval = TRUE, fig.height=24, fig.width=10}

data_expr_rec_season_res_simu = make_exp_analyses_only_simu(data_expr_rec_season, run_whole_dataset = run_whole_dataset) 
data_expr_rec_season_res_simu %>% subset(couple == "All_data") %>%  kbl() %>% kable_styling()

```

```{r, fig.path = output_dir, eval = TRUE, fig.height=24, fig.width=10}

data_expr_rec_season_res_simu_r = make_exp_analyses_only_simu(data_expr_rec_season,run_whole_dataset = F, replace_signi_genes = T )
data_expr_rec_season_res_simu_r %>% subset(couple == "All_data") %>%  kbl() %>% kable_styling()

```


# Dataset expression 2022, SET MURINAE INTRA-FAMILY , design : ~ cond_season


##  ~ env (no sva)

```{r, eval = T}
output_dir="2025_simu_dataset_murinae_season/"
system(paste0("mkdir -p " , output_dir))

data_expr_mur_season = list(var1_design = NULL, var2_design = "cond_season", Alltable = Alltable_onl, coldata = coldata_murinae, output_dir=output_dir, resampling = resampling , nRand = nRand, nCPU = nCPU, test = test, column_permutation = "murinae_permutation", dataset_name="murinae")

```

```{r, fig.path = output_dir, eval = T, fig.height=25, fig.width=10}

data_expr_mur_season_res_simu = make_exp_analyses_only_simu(data_expr_mur_season, run_whole_dataset = run_whole_dataset) 
data_expr_mur_season_res_simu %>% subset(couple == "All_data") %>%  kbl() %>% kable_styling()

```

```{r, fig.path = output_dir, eval = T, fig.height=25, fig.width=10}

data_expr_mur_season_res_simu_r = make_exp_analyses_only_simu(data_expr_mur_season, replace_signi_genes = T, run_whole_dataset = F)
data_expr_mur_season_res_simu_r %>% subset(couple == "All_data") %>%  kbl() %>% kable_styling()

```

# Dataset expression 2022, TOTAL , design : ~ family + cond_season 


## ~ family + env (family in sva)

```{r, eval = T}
output_dir="13102025_simu_dataset_total_season/"
system(paste0("mkdir -p " , output_dir))

data_expr_total_fam_season = list(var1_design = "total_group_batch", var2_design = "cond_season", Alltable = Alltable_onl, coldata = coldata_tot, output_dir=output_dir, resampling = resampling , nRand = nRand, nCPU = nCPU, column_permutation = "total_permutation", test = test, dataset_name = "total_design_fam_cond")

```

```{r, fig.path = output_dir, eval = T, fig.height=20, fig.width=8}

data_expr_total_fam_season_res_simu = make_exp_analyses_only_simu(data_expr_total_fam_season, run_whole_dataset = run_whole_dataset)
data_expr_total_fam_season_res_simu %>% subset(couple == "All_data") %>%  kbl() %>% kable_styling()

```

```{r, fig.path = output_dir, eval = T, fig.height=20, fig.width=8}

data_expr_total_fam_season_res_simu_r = make_exp_analyses_only_simu(data_expr_total_fam_season, replace_signi_genes = T, run_whole_dataset = F)
data_expr_total_fam_season_res_simu_r %>% subset(couple == "All_data") %>%  kbl() %>% kable_styling()

```


# Plot Final Bilan plot

```{r}
all_datasets_res = rbind(
  data_expr_rec_season_res_simu,
  data_expr_rec_season_res_simu_r,
  data_expr_anc_season_res_simu,
  data_expr_anc_season_res_simu_r,
  data_expr_mur_season_res_simu,
  data_expr_mur_season_res_simu_r,
  data_expr_total_fam_season_res_simu,
  data_expr_total_fam_season_res_simu_r
)

write_tsv(
  x = all_datasets_res %>%
    subset(couple == "All_data") %>%
    dplyr::select(-couple) %>%
    relocate("dataset"),
  file = "bilan_simulation_papier_raw_data.tsv")
```


```{r}
all_datasets_res_tmp = read_tsv("bilan_simulation_papier_raw_data.tsv")

library(ggsignif)

cat_order = c("Total","Up","Down","Conserved","Diversification","LFC.XM.P","LFC.XM.N")
count_table_order = c("True","Shuffled_all_genes","Shuffled_replace_DE")
dataset_type_order = c("replace_DE", "all_genes")

dataset_order = unique(all_datasets_res_tmp$dataset)
dataset_prefix_order = unique(gsub("_replace_signi","",dataset_order))

all_datasets_res_df = 
  all_datasets_res_tmp %>%
  pivot_longer(-dataset) %>%
  separate("name",sep = "\\.",into = c("variable","cat"), extra = "merge") %>% # split by last .
  mutate(count_table = ifelse(grepl("true",variable),"True","Other"),
         count_table = ifelse(grepl("shuffling",variable),"Shuffled_all_genes",count_table),
         count_table = ifelse(grepl("replace_signi",dataset) & count_table == "Shuffled_all_genes","Shuffled_replace_DE",count_table))

all_datasets_res_df_mean =
  subset(all_datasets_res_df, grepl("mean", all_datasets_res_df$variable)) %>%
  dplyr::rename(mean = "value")

all_datasets_res_df_mean_sd =
  subset(all_datasets_res_df, grepl("sd_[st]", all_datasets_res_df$variable)) %>%
  left_join(all_datasets_res_df_mean %>% dplyr::select(-variable), by = c("cat",  "dataset", "count_table")) %>%
  dplyr::rename(sd = "value") %>%
  mutate(cat = factor(cat, levels = cat_order),
         count_table = factor(count_table, levels = count_table_order),
         dataset = factor(dataset, levels = dataset_order),
         dataset_type = ifelse(grepl("replace_signi",dataset),"replace_DE","all_genes"),
         dataset_prefix = gsub("_replace_signi","",dataset),
         dataset_prefix = factor(dataset_prefix, levels = dataset_prefix_order),
         dataset_type = factor(dataset_type, levels = dataset_type_order))
         

all_datasets_res_df_quantile =
  subset(all_datasets_res_df, grepl("q[257]", all_datasets_res_df$variable)) %>%
  separate("variable",sep = "_",into = c("quantile","t_or_s"), extra = "merge") %>%
  dplyr::select(-t_or_s) %>%
  pivot_wider(names_from = quantile) %>%
  mutate(cat = factor(cat, levels = cat_order),
         count_table = factor(count_table, levels = count_table_order),
         dataset = factor(dataset, levels = dataset_order),
         dataset_type = ifelse(grepl("replace_signi",dataset),"replace_DE","all_genes"),
         dataset_prefix = gsub("_replace_signi","",dataset),
         dataset_prefix = factor(dataset_prefix, levels = dataset_prefix_order),
         dataset_type = factor(dataset_type, levels = dataset_type_order))

all_datasets_res_df_pvalues =
  all_datasets_res_df %>%
  subset(grepl("pvalues", variable)) %>%
  dplyr::select(-count_table) %>%
  pivot_wider(names_from = variable)  %>%
  mutate(cat = factor(cat, levels = cat_order))


marge_signi_bar = all_datasets_res_df_mean_sd %>%
  group_by(dataset_prefix, cat) %>%
  summarise(ymax = max(sd + mean)) %>%
  group_by(dataset_prefix) %>%
  summarise(
    marge_signi_bar_val = max(ymax) * 0.1
  ) 

makeStars <- function(x){
  stars <- c("****", "***", "**", "*", "ns")
  vec <- c(0, 0.0001, 0.001, 0.01, 0.05, 1)
  i <- findInterval(x, vec)
  stars[i]
}

all_datasets_res_df_increase = 
  all_datasets_res_df_mean  %>%
  dplyr::select(-variable) %>%
  pivot_wider(names_from = c(count_table), values_from = mean)  %>%
  mutate(increase = True / (ifelse(is.na(Shuffled_all_genes), Shuffled_replace_DE, Shuffled_all_genes))) %>%
  dplyr::select(dataset,cat,increase)


all_datasets_res_df_increase_med = 
  all_datasets_res_df_quantile  %>%
  dplyr::select(-q25,-q75) %>%
  pivot_wider(names_from = c(count_table), values_from = q50)  %>%
  mutate(increase = True / (ifelse(is.na(Shuffled_all_genes), Shuffled_replace_DE, Shuffled_all_genes))) %>%
  dplyr::select(dataset,cat,increase)


all_datasets_res_df_signig = 
  all_datasets_res_df_mean_sd %>%
  left_join(marge_signi_bar) %>%
  group_by(dataset, cat, dataset_type) %>%
  summarise(ymax = max(sd + mean) + marge_signi_bar_val) %>%
  ungroup()  %>% 
  mutate(cat = factor(cat, levels = cat_order),
         dataset = factor(dataset, levels = dataset_order),
         dataset_type = factor(dataset_type, levels = dataset_type_order)) %>%
  mutate(xmin = as.numeric(dataset_type) - 0.2,
         xmax = as.numeric(dataset_type) + 0.2,
         xmin_papier = as.numeric(cat) - 0.2,
         xmax_papier = as.numeric(cat) + 0.2) %>%
  left_join(all_datasets_res_df_pvalues) %>%
  left_join(all_datasets_res_df_increase_med) %>%
  mutate(cat = factor(cat, levels = cat_order),
         dataset = factor(dataset, levels = dataset_order),
         dataset_prefix = gsub("_replace_signi","",dataset),
         dataset_prefix = factor(dataset_prefix, levels = dataset_prefix_order)
  ) %>%
  mutate(pvalues_str = as.character(round(pvalues,4)),
         pvalues_w_str = as.character(round(pvalues_wilc_paired,4)),
         pvalues_w_star = ifelse(pvalues_wilc_paired > 0.05, makeStars(pvalues_wilc_paired),
                                 ifelse(increase > 1,
                                        paste0(makeStars(pvalues_wilc_paired), "(x", round(increase,1),")"),
                                        paste0(makeStars(pvalues_wilc_paired), "(/", round(1/increase,1),")")
                                 ))
  )


fake_line =  all_datasets_res_df_quantile %>%
  group_by(cat, dataset_prefix) %>%
  summarise(ymax = max(q75) * 1.5)
```


```{r fig.height=15, fig.width=12, eval = F}
fig_bilan = ggplot() + theme_bw() + theme(legend.position = "top") +  
  geom_bar(data = all_datasets_res_df_mean_sd, 
           aes(x=dataset_type,y=mean,fill = count_table), col = "black",
           stat = "identity", position = position_dodge(.9)) +
  facet_grid( dataset_prefix ~ cat, scale = "free_y") +
  geom_errorbar(data = all_datasets_res_df_mean_sd, 
                aes(x=dataset_type,ymin=mean-sd, ymax=mean+sd,fill = count_table), width=.2,
                 position=position_dodge(.9)) +
  geom_signif( data = all_datasets_res_df_signig,
               aes(y_position=ymax, xmin=xmin, xmax=xmax,annotations=pvalues_w_star, group = dataset),
                tip_length = 0.01, manual = T, vjust = -0.5 , map_signif_level=F) +
  geom_hline(data = fake_line, aes(yintercept = ymax), alpha = 0) + ggtitle("Figure Bilan Simulation", subtitle = "pvalue wilcoxon test paired :\n0 < **** <0.0001 < *** < 0.001 < ** < 0.01 < * < 0.05 < ns < 1") +
  xlab ("Used genes") + ylab("mean number of DE genes")

fig_bilan
ggsave("bilan_simulation_tous_les_datasets.pdf", plot = fig_bilan, width = 12, height = 15)

```

```{r fig.height=15, fig.width=12}
fig_bilan = ggplot() + theme_bw() + theme(legend.position = "top") +  
  geom_bar(data = all_datasets_res_df_quantile, 
           aes(x=dataset_type,y=q50,fill = count_table), col = "black",
           stat = "identity", position = position_dodge(.9)) +
  facet_grid( dataset_prefix ~ cat, scale = "free_y") +
  geom_errorbar(data = all_datasets_res_df_quantile, 
                aes(x=dataset_type,ymin=q25, ymax=q75,fill = count_table), width=.2,
                 position=position_dodge(.9)) +
  geom_signif( data = all_datasets_res_df_signig,
               aes(y_position=ymax, xmin=xmin, xmax=xmax,annotations=pvalues_w_star, group = dataset),
                tip_length = 0.01, manual = T, vjust = -0.5 , map_signif_level=F) +
  geom_hline(data = fake_line, aes(yintercept = ymax), alpha = 0) + ggtitle("Figure Bilan Simulation", subtitle = "pvalue wilcoxon test paired :\n0 < **** <0.0001 < *** < 0.001 < ** < 0.01 < * < 0.05 < ns < 1\nQ25,Med,Q75") +
  xlab ("Used genes") + ylab("mean number of DE genes")

fig_bilan
ggsave("bilan_simulation_tous_les_datasets.pdf", plot = fig_bilan, width = 12, height = 15)

```

```{r fig.height=10, fig.width=6}
all_datasets_res_tmp = read_tsv("bilan_simulation_papier_raw_data.tsv")

new_dataset_name = c("Total" = "total_design_fam_cond",
                     "Ancient" = "ancient_trans",
                     "Recent" = "recent_trans",
                     "Murinae" = "murinae"
                     )
                     

all_datasets_res_tmp$dataset = plyr::mapvalues(all_datasets_res_tmp$dataset , unname(new_dataset_name),names(new_dataset_name))
all_datasets_res_tmp$dataset = plyr::mapvalues(all_datasets_res_tmp$dataset , paste0(unname(new_dataset_name),"_replace_signi"),paste0(names(new_dataset_name),"_replace_signi"))

cat_order = c("Total","LFC.XM.P","LFC.XM.N")
count_table_order = c("True","Shuffled_all_genes","Shuffled_replace_DE")
dataset_type_order = c("all_genes","replace_DE")

dataset_order = unique(all_datasets_res_tmp$dataset)
dataset_prefix_order = names(new_dataset_name)

all_datasets_res_df = 
  all_datasets_res_tmp %>%
  pivot_longer(-dataset) %>%
  separate("name",sep = "\\.",into = c("variable","cat"), extra = "merge") %>% # split by last .
  mutate(count_table = ifelse(grepl("true",variable),"True","Other"),
         count_table = ifelse(grepl("shuffling",variable),"Shuffled_all_genes",count_table),
         count_table = ifelse(grepl("replace_signi",dataset) & count_table == "Shuffled_all_genes","Shuffled_replace_DE",count_table)) %>%
  filter(cat %in% cat_order)

all_datasets_res_df_quantile =
  subset(all_datasets_res_df, grepl("q[257]", all_datasets_res_df$variable)) %>%
  separate("variable",sep = "_",into = c("quantile","t_or_s"), extra = "merge") %>%
  dplyr::select(-t_or_s) %>%
  pivot_wider(names_from = quantile) %>%
  mutate(cat = factor(cat, levels = cat_order),
         count_table = factor(count_table, levels = count_table_order),
         dataset = factor(dataset, levels = dataset_order),
         dataset_type = ifelse(grepl("replace_signi",dataset),"replace_DE","all_genes"),
         dataset_prefix = gsub("_replace_signi","",dataset),
         dataset_prefix = factor(dataset_prefix, levels = dataset_prefix_order),
         dataset_type = factor(dataset_type, levels = dataset_type_order))

all_datasets_res_df_pvalues =
  all_datasets_res_df %>%
  subset(grepl("pvalues", variable)) %>%
  dplyr::select(-count_table) %>%
  pivot_wider(names_from = variable)  %>%
  mutate(cat = factor(cat, levels = cat_order))



marge_signi_bar = all_datasets_res_df_quantile %>%
  group_by(dataset_prefix, cat) %>%
  summarise(ymax = max(q75)) %>%
  group_by(dataset_prefix) %>%
  summarise(
    marge_signi_bar_val = max(ymax) * 0.1
  ) 


makeStars <- function(x){
  stars <- c("****", "***", "**", "*", "ns")
  vec <- c(0, 0.0001, 0.001, 0.01, 0.05, 1)
  i <- findInterval(x, vec)
  stars[i]
}


all_datasets_res_df_increase_med = 
  all_datasets_res_df_quantile  %>%
  dplyr::select(-q25,-q75) %>%
  pivot_wider(names_from = c(count_table), values_from = q50)  %>%
  mutate(increase = True / (ifelse(is.na(Shuffled_all_genes), Shuffled_replace_DE, Shuffled_all_genes))) %>%
  dplyr::select(dataset,cat,increase)


all_datasets_res_df_signig = 
all_datasets_res_df_quantile  %>%
  left_join(marge_signi_bar) %>%
  group_by(dataset, cat, dataset_type) %>%
  summarise(ymax = max(q75) + marge_signi_bar_val)  %>%
  ungroup()  %>% 
  mutate(cat = factor(cat, levels = cat_order),
         dataset = factor(dataset, levels = dataset_order),
         dataset_type = factor(dataset_type, levels = dataset_type_order)) %>%
  mutate(xmin = as.numeric(dataset_type) - 0.2,
         xmax = as.numeric(dataset_type) + 0.2,
         xmin_papier = as.numeric(cat) - 0.2,
         xmax_papier = as.numeric(cat) + 0.2) %>%
  left_join(all_datasets_res_df_pvalues) %>%
  left_join(all_datasets_res_df_increase_med) %>%
  mutate(cat = factor(cat, levels = cat_order),
         dataset = factor(dataset, levels = dataset_order),
         dataset_prefix = gsub("_replace_signi","",dataset),
         dataset_prefix = factor(dataset_prefix, levels = dataset_prefix_order)
  ) %>%
  mutate(pvalues_str = as.character(round(pvalues,4)),
         pvalues_w_str = as.character(round(pvalues_wilc_paired,4)),
         pvalues_w_star = ifelse(pvalues_wilc_paired > 0.05, makeStars(pvalues_wilc_paired),
                                 ifelse(increase > 1,
                                        paste0(makeStars(pvalues_wilc_paired), "(x", round(increase,1),")"),
                                        paste0(makeStars(pvalues_wilc_paired), "(/", round(1/increase,1),")")
                                 ))
  )


fake_line =  all_datasets_res_df_quantile %>%
  group_by(cat, dataset_prefix) %>%
  summarise(ymax = max(q75) * 1.5)



removed_dataset = c("total_design_cond", "total_design_cond_replace_signi", "total_design_fam_cond", "total_design_fam_cond_replace_signi")
removed_dataset = unique(all_datasets_res_df_quantile$dataset)[! grepl("_replace_signi",unique(all_datasets_res_df_quantile$dataset))]


all_datasets_res_df_quantile_papier =
  filter(all_datasets_res_df_quantile,
         ! dataset  %in% removed_dataset &
           dataset_type == "replace_DE" &
           cat %in% c("Total","LFC.XM.P","LFC.XM.N")
         )

all_datasets_res_df_quantile_papier$cat = factor(all_datasets_res_df_quantile_papier$cat, levels = c("Total","LFC.XM.P","LFC.XM.N"))
all_datasets_res_df_signig_papier = 
  filter(all_datasets_res_df_signig, ! dataset  %in%  removed_dataset &
           dataset_type == "replace_DE" &
           cat %in% c("Total","LFC.XM.P","LFC.XM.N"))

all_datasets_res_df_signig_papier$cat = factor(all_datasets_res_df_signig_papier$cat, levels = c("Total","LFC.XM.P","LFC.XM.N"))

fake_line_papier =  
  filter(fake_line,
           cat %in% c("Total","LFC.XM.P","LFC.XM.N")
         )
fake_line_papier$cat = factor(fake_line_papier$cat, levels = c("Total","LFC.XM.P","LFC.XM.N"))

  #subset(fake_line, !  dataset_prefix  %in%  removed_dataset)

fig_papier1 = ggplot() + theme_bw() + theme(legend.position = "top") +  
  geom_bar(data = all_datasets_res_df_quantile_papier, 
           aes(x=cat,y=q50,fill = count_table), col = "black",
           stat = "identity", position = position_dodge(.9)) +
  ggforce::facet_col(~dataset_prefix, scales = "free_y") +
  geom_errorbar(data = all_datasets_res_df_quantile_papier, 
                aes(x=cat,ymin=q25, ymax=q75,fill = count_table), width=.2,
                 position=position_dodge(.9)) +
  geom_signif( data = all_datasets_res_df_signig_papier,
               aes(y_position=ymax, xmin=xmin_papier, xmax=xmax_papier,annotations=pvalues_w_star, group = cat),
                tip_length = 0.01, manual = T, vjust = -0.5 , map_signif_level=F) +
  geom_hline(data = fake_line_papier, aes(yintercept = ymax), alpha = 0) + ggtitle("Figure Papier à retravailler sous inkscape", subtitle = "pvalue wilcoxon test paired :\n0 < **** <0.0001 < *** < 0.001 < ** < 0.01 < * < 0.05 < ns < 1\nQ25,Med,Q75") +
  xlab ("") + ylab("median number of DE genes") +
  scale_fill_grey()

fig_papier1
ggsave("bilan_simulation_papier_à_retravailler1.pdf", plot = fig_papier1, width = 6, height = 8)
```
```{r}
all_datasets_res_tmp = read_tsv("bilan_simulation_papier_raw_data.tsv")

new_dataset_name = c("Total" = "total_design_fam_cond",
                     "Ancient" = "ancient_trans",
                     "Recent" = "recent_trans",
                     "Murinae" = "murinae"
                     )
                     

all_datasets_res_tmp$dataset = plyr::mapvalues(all_datasets_res_tmp$dataset , unname(new_dataset_name),names(new_dataset_name))
all_datasets_res_tmp$dataset = plyr::mapvalues(all_datasets_res_tmp$dataset , paste0(unname(new_dataset_name),"_replace_signi"),paste0(names(new_dataset_name),"_replace_signi"))

cat_order = c("Total")
count_table_order = c("True","Shuffled_all_genes","Shuffled_replace_DE")
dataset_type_order = c("all_genes","replace_DE")

dataset_order = unique(all_datasets_res_tmp$dataset)
dataset_prefix_order = names(new_dataset_name)

all_datasets_res_df = 
  all_datasets_res_tmp %>%
  pivot_longer(-dataset) %>%
  separate("name",sep = "\\.",into = c("variable","cat"), extra = "merge") %>% # split by last .
  mutate(count_table = ifelse(grepl("true",variable),"True","Other"),
         count_table = ifelse(grepl("shuffling",variable),"Shuffled_all_genes",count_table),
         count_table = ifelse(grepl("replace_signi",dataset) & count_table == "Shuffled_all_genes","Shuffled_replace_DE",count_table)) %>%
  filter(cat %in% cat_order)

all_datasets_res_df_quantile =
  subset(all_datasets_res_df, grepl("q[257]", all_datasets_res_df$variable)) %>%
  separate("variable",sep = "_",into = c("quantile","t_or_s"), extra = "merge") %>%
  dplyr::select(-t_or_s) %>%
  pivot_wider(names_from = quantile) %>%
  mutate(cat = factor(cat, levels = cat_order),
         count_table = factor(count_table, levels = count_table_order),
         dataset = factor(dataset, levels = dataset_order),
         dataset_type = ifelse(grepl("replace_signi",dataset),"replace_DE","all_genes"),
         dataset_prefix = gsub("_replace_signi","",dataset),
         dataset_prefix = factor(dataset_prefix, levels = dataset_prefix_order),
         dataset_type = factor(dataset_type, levels = dataset_type_order))

all_datasets_res_df_pvalues =
  all_datasets_res_df %>%
  subset(grepl("pvalues", variable)) %>%
  dplyr::select(-count_table) %>%
  pivot_wider(names_from = variable)  %>%
  mutate(cat = factor(cat, levels = cat_order))



marge_signi_bar = all_datasets_res_df_quantile %>%
  group_by(dataset_prefix, cat) %>%
  summarise(ymax = max(q75)) %>%
  group_by(dataset_prefix) %>%
  summarise(
    marge_signi_bar_val = max(ymax) * 0.1
  ) 


makeStars <- function(x){
  stars <- c("****", "***", "**", "*", "ns")
  vec <- c(0, 0.0001, 0.001, 0.01, 0.05, 1)
  i <- findInterval(x, vec)
  stars[i]
}


all_datasets_res_df_increase_med = 
  all_datasets_res_df_quantile  %>%
  dplyr::select(-q25,-q75) %>%
  pivot_wider(names_from = c(count_table), values_from = q50)  %>%
  mutate(increase = True / (ifelse(is.na(Shuffled_all_genes), Shuffled_replace_DE, Shuffled_all_genes))) %>%
  dplyr::select(dataset,cat,increase)


all_datasets_res_df_signig = 
all_datasets_res_df_quantile  %>%
  left_join(marge_signi_bar) %>%
  group_by(dataset, cat, dataset_type) %>%
  summarise(ymax = max(q75) + marge_signi_bar_val)  %>%
  ungroup()  %>% 
  mutate(cat = factor(cat, levels = cat_order),
         dataset = factor(dataset, levels = dataset_order),
         dataset_type = factor(dataset_type, levels = dataset_type_order)) %>%
  mutate(xmin = as.numeric(dataset_type) - 0.2,
         xmax = as.numeric(dataset_type) + 0.2,
         xmin_papier = as.numeric(cat) - 0.2,
         xmax_papier = as.numeric(cat) + 0.2) %>%
  left_join(all_datasets_res_df_pvalues) %>%
  left_join(all_datasets_res_df_increase_med) %>%
  mutate(cat = factor(cat, levels = cat_order),
         dataset = factor(dataset, levels = dataset_order),
         dataset_prefix = gsub("_replace_signi","",dataset),
         dataset_prefix = factor(dataset_prefix, levels = dataset_prefix_order)
  ) %>%
  mutate(pvalues_str = as.character(round(pvalues,4)),
         pvalues_w_str = as.character(round(pvalues_wilc_paired,4)),
         pvalues_w_star = ifelse(pvalues_wilc_paired > 0.05, makeStars(pvalues_wilc_paired),
                                 ifelse(increase > 1,
                                        paste0(makeStars(pvalues_wilc_paired), "(x", round(increase,1),")"),
                                        paste0(makeStars(pvalues_wilc_paired), "(/", round(1/increase,1),")")
                                 ))
  )


fake_line =  all_datasets_res_df_quantile %>%
  group_by(cat, dataset_prefix) %>%
  summarise(ymax = max(q75) * 1.5)



removed_dataset = c("total_design_cond", "total_design_cond_replace_signi", "total_design_fam_cond", "total_design_fam_cond_replace_signi")
removed_dataset = unique(all_datasets_res_df_quantile$dataset)[! grepl("_replace_signi",unique(all_datasets_res_df_quantile$dataset))]


all_datasets_res_df_quantile_papier =
  filter(all_datasets_res_df_quantile,
         ! dataset  %in% removed_dataset &
           dataset_type == "replace_DE"
         )

all_datasets_res_df_quantile_papier$cat = factor(all_datasets_res_df_quantile_papier$cat, levels = cat_order)
all_datasets_res_df_signig_papier = 
  filter(all_datasets_res_df_signig, ! dataset  %in%  removed_dataset &
           dataset_type == "replace_DE" )

all_datasets_res_df_signig_papier$cat = factor(all_datasets_res_df_signig_papier$cat, levels = cat_order)


fake_line_papier$cat = factor(fake_line_papier$cat, levels = cat_order)

  #subset(fake_line, !  dataset_prefix  %in%  removed_dataset)

fig_papier2 = ggplot() + theme_bw() + theme(legend.position = "top") +  
  geom_bar(data = all_datasets_res_df_quantile_papier, 
           aes(x=cat,y=q50,fill = count_table), col = "black",
           stat = "identity", position = position_dodge(.9)) +
  ggforce::facet_row(~dataset_prefix, scales = "free_y") +
  geom_errorbar(data = all_datasets_res_df_quantile_papier, 
                aes(x=cat,ymin=q25, ymax=q75,fill = count_table), width=.2,
                 position=position_dodge(.9)) +
  geom_signif( data = all_datasets_res_df_signig_papier,
               aes(y_position=ymax, xmin=xmin_papier, xmax=xmax_papier,annotations=pvalues_w_star, group = cat),
                tip_length = 0.01, manual = T, vjust = -0.5 , map_signif_level=F) +
  geom_hline(data = fake_line_papier, aes(yintercept = ymax), alpha = 0) + ggtitle("Figure Papier à retravailler sous inkscape", subtitle = "pvalue wilcoxon test paired :\n0 < **** <0.0001 < *** < 0.001 < ** < 0.01 < * < 0.05 < ns < 1\nQ25,Med,Q75") +
  xlab ("") + ylab("median number of DE genes") +
  scale_fill_grey()

fig_papier2
ggsave("bilan_simulation_papier_à_retravailler2.pdf", plot = fig_papier2, width = 8, height = 4)
```

