2  Microarray FABER (2014)

Code 
#pkg
library(readxl)
library(affy)   # Affymetrix pre-processing
library(limma)  # two-color pre-processing; differential
library(tidyverse)
library(dplyr)
library(patchwork)
library(progress)
library(ggrepel)
library(plotly)
library(htmlwidgets)
library("FactoMineR")
library("factoextra")
library("corrplot")
library(missMDA)
library(ggnewscale)

# src
source(here::here("src/function/stat_function/stat_analysis_main.R")) # for make plot 
source(here::here("src/function/fig_export.R")) # This function saves a given plot (plot_x) as both a PDF and a high-resolution PNG file at specified dimensions.
source(here::here("src/function/upsetplot_condition_merge_sign.R")) # function "create_presence_matrix" and "upsetplot_condition_merge_sign"

# cosmetics
water_pallet=read_excel(here::here("data/color_palette.xlsm")) %>%
      filter(set == "water_condition") %>%
      dplyr::select(color, treatment) %>%
      pull(color) %>%
      setNames(read_excel(here::here("data/color_palette.xlsm")) %>%
                 filter(set == "water_condition") %>%
                 pull(treatment)
               )
Code 
# for the list of all comparison
## For nodule
comparison=read_excel(here::here("data/microarray/Resultats_12plex-PEA-2013-03_NODULES_mod_281014.xlsx"), sheet = "Complete", range = "L10:DD11", col_names = FALSE) %>% 
  t() %>% 
  as.data.frame() %>% 
  rownames_to_column("to_del") %>%
  dplyr::select(-"to_del") %>%
  mutate(comparison = paste0(V1, ".", V2)) %>% 
  distinct(comparison) %>% 
  tidyr::separate(comparison, into = c("Green", "Red"),sep = "\\.")

write_csv(x = comparison, file = here::here("data/microarray/output/comparison_microarray_nodule.csv"))

## For root
comparison=read_excel(here::here("data/microarray/Resultats_12plex-PEA-2013-03_RACINES_mod_281014.xlsx"), sheet = "Complete", range = "L10:DD11", col_names = FALSE) %>% 
  t() %>% 
  as.data.frame() %>% 
  rownames_to_column("to_del") %>%
  dplyr::select(-"to_del") %>%
  mutate(comparison = paste0(V1, ".", V2)) %>% 
  distinct(comparison) %>% 
  tidyr::separate(comparison, into = c("Green", "Red"),sep = "\\.")

write_csv(x = comparison, file = here::here("data/microarray/output/comparison_microarray_root.csv"))
Code 
read_csv(here::here("data/microarray/output/comparison_microarray_nodule.csv")) %>% 
knitr::kable(., caption = "List of all comparison")
List of all comparison
Green Red
NOD_WW-D1_1 NOD_WW-C1_1
NOD_WW-D2_2 NOD_WW-C2_2
NOD_WW-D3_3 NOD_WW-C3_3
NOD_WW-D4_4 NOD_WW-C4_4
NOD_WS-C1_1 NOD_WW-C1_1
NOD_WS-C2_2 NOD_WW-C2_2
NOD_WS-C3_3 NOD_WW-C3_3
NOD_WS-C4_4 NOD_WW-C4_4
NOD_WS-D1_1 NOD_WS-C1_1
NOD_WS-D2_2 NOD_WS-C2_2
NOD_WS-D3_3 NOD_WS-C3_3
NOD_WS-D4_4 NOD_WS-C4_4
NOD_WS-D1_1 NOD_WW-D1_1
NOD_WS-D2_2 NOD_WW-D2_2
NOD_WS-D3_3 NOD_WW-D3_3
NOD_WS-D4_4 NOD_WW-D4_4
Code 
read_csv(here::here("data/microarray/output/comparison_microarray_root.csv")) %>% 
knitr::kable(., caption = "List of all comparison")
List of all comparison
Green Red
RAC_WW-D1_1 RAC_WW-C1_1
RAC_WW-D2_2 RAC_WW-C2_2
RAC_WW-D3_3 RAC_WW-C3_3
RAC_WW-D4_4 RAC_WW-C4_4
RAC_WS-C1_1 RAC_WW-C1_1
RAC_WS-C2_2 RAC_WW-C2_2
RAC_WS-C3_3 RAC_WW-C3_3
RAC_WS-C4_4 RAC_WW-C4_4
RAC_WS-D1_1 RAC_WS-C1_1
RAC_WS-D2_2 RAC_WS-C2_2
RAC_WS-D3_3 RAC_WS-C3_3
RAC_WS-D4_4 RAC_WS-C4_4
RAC_WS-D1_1 RAC_WW-D1_1
RAC_WS-D2_2 RAC_WW-D2_2
RAC_WS-D3_3 RAC_WW-D3_3
RAC_WS-D4_4 RAC_WW-D4_4

2.1 Data importation

  • Import data but only RebBio column (ratio and stats)
Code 
source (here::here("src/function/microarray/import_ratio_stats_microarray.R"))

# For root
path_i = here::here("data/microarray/Resultats_12plex-PEA-2013-03_RACINES_mod_281014.xlsx")
path_info_range_i = here::here("data/microarray/info_range_Resultats_12plex-PEA-2013-03_RACINES_mod_281014.xlsx")
begin_end_line <- c("12", "40730")
sheet_i = "Complete"

execut_root = import_ratio_stats_microarray(path_i, path_info_range_i, sheet_i, begin_end_line)

# For nodule
path_i = here::here("data/microarray/Resultats_12plex-PEA-2013-03_NODULES_mod_281014.xlsx")
path_info_range_i = here::here("data/microarray/info_range_Resultats_12plex-PEA-2013-03_NODULES_mod_281014.xlsx")
begin_end_line <- c("12", "40730")
sheet_i = "Complete"

execut_nodule = import_ratio_stats_microarray(path_i, path_info_range_i, sheet_i, begin_end_line)

# names(execut_nodule[["list_ratio_stats"]])

# Merge both data file
list_ratio_stats_global = c(execut_root[["list_ratio_stats"]],execut_nodule[["list_ratio_stats"]])
df_info_sample_clean_compile_global = rbind(execut_root[["df_info_sample_clean_compile"]], execut_nodule[["df_info_sample_clean_compile"]])

save(list_ratio_stats_global, file =  here::here("data/microarray/output/ratio_stats_microarray_root_nodule_FABER.RData"))
write_csv(df_info_sample_clean_compile_global, here::here("data/microarray/output/list_microarray_root_nodule_FABER_RepBio.csv"))

2.2 Normalisation for PCA ?

2.3 Differential expression

2.3.1 Of non-normalized data (stats from the platform)

I had to transform the statistics equal to 0 for the vulcanoplot 0 was replaced by 1.10-14

Because -log10(0) = Inf

Code 
# /!\ vey slow
load(file =  here::here("data/microarray/output/ratio_stats_microarray_root_nodule_FABER.RData"))
source (here::here("src/function/microarray/vulcanoplot_microarray_ratio_stats.R"))

for(i in 1:length(names(list_ratio_stats_global))){
  # create ggplot
  p = vulcanoplot_microarray_ratio_stats(comparison = names(list_ratio_stats_global)[i])
  fig_export(here::here(paste0("report/microarray/plot/volcanoplot/volcano_interactive_ratio_stats_",i)), p, height_i = 3.5*1.5, width_i = 5*1.5, res_i = 300)
  # export ggplot
  # export html for interactif
  ggplotly(p) %>% 
    saveWidget(here::here(paste0("report/microarray/plot/volcanoplot/volcano_interactive_ratio_stats_",i,".html")))
}
Code 
load(file =  here::here("data/microarray/output/ratio_stats_microarray_root_nodule_FABER.RData"))

For an interactive version : Click here

For an interactive version : Click here

For an interactive version : Click here

For an interactive version : Click here

For an interactive version : Click here

For an interactive version : Click here

For an interactive version : Click here

For an interactive version : Click here

2.4 UpsetPlot

2.4.1 All comparison

Code 
# import data
load(file =  here::here("data/microarray/output/ratio_stats_microarray_root_nodule_FABER.RData"))

# create function to have list of genes
get_genes<- function(df_comparison_i, lfc_lim = 0, padj_lim= 0.05, sign_i){
  v_gene = df_comparison_i %>% 
    mutate(sign= ifelse(logFC>0, "Up", "Down")) %>% 
    filter(
      sign == sign_i,
      abs(logFC)>abs(lfc_lim), 
      BH<= padj_lim 
           ) %>% 
    pull(id_probe)
  return(v_gene)
}

# get_genes(df_comparison_i = list_ratio_stats$`Mut_SS_W78*_RepBio_4 vs WT_SS_WT1_RepBio_4`, sign_i = "Up")

 # creation of the list
lt_up <- list(
  `Root_WW_D vs Root_WW_C` = get_genes(list_ratio_stats_global$`Root_WW_D_1.1_rep_bio_1 vs Root_WW_C_1.1_rep_bio_1`, sign_i = "Up"),
  `Root_WS_C vs Root_WW_C` = get_genes(list_ratio_stats_global$`Root_WS_C_1.1_rep_bio_2 vs Root_WW_C_1.1_rep_bio_2`, sign_i = "Up"),
  `Root_WS_D vs Root_WS_C` = get_genes(list_ratio_stats_global$`Root_WS_D_1.1_rep_bio_3 vs Root_WS_C_1.1_rep_bio_3`, sign_i = "Up"),
  `Root_WS_D vs Root_WW_D` = get_genes(list_ratio_stats_global$`Root_WS_D_1.1_rep_bio_4 vs Root_WW_D_1.1_rep_bio_4`, sign_i = "Up"),
  
  `Nodule_WW_D vs Nodule_WW_C` = get_genes(list_ratio_stats_global$`Nodule_WW_D_1.1_Rep_Bio_1 vs Nodule_WW_C_1.1_Rep_Bio_1`, sign_i = "Up"),
  `Nodule_WS_C vs Nodule_WW_C` = get_genes(list_ratio_stats_global$`Nodule_WS_C_1.1_Rep_Bio_2 vs Nodule_WW_C_1.1_Rep_Bio_2`, sign_i = "Up"),
  `Nodule_WS_D vs Nodule_WS_C` = get_genes(list_ratio_stats_global$`Nodule_WS_D_1.1_Rep_Bio_3 vs Nodule_WS_C_1.1_Rep_Bio_3`, sign_i = "Up"),
  `Nodule_WS_D vs Nodule_WW_D` = get_genes(list_ratio_stats_global$`Nodule_WS_D_1.1_Rep_Bio_4 vs Nodule_WW_D_1.1_Rep_Bio_4`, sign_i = "Up")
)

m_up <- ComplexHeatmap::make_comb_mat(lt_up, mode = "distinct")

# Down-regulated genes
lt_down <- list(
  `Root_WW_D vs Root_WW_C` = get_genes(list_ratio_stats_global$`Root_WW_D_1.1_rep_bio_1 vs Root_WW_C_1.1_rep_bio_1`, sign_i = "Down"),
  `Root_WS_C vs Root_WW_C` = get_genes(list_ratio_stats_global$`Root_WS_C_1.1_rep_bio_2 vs Root_WW_C_1.1_rep_bio_2`, sign_i = "Down"),
  `Root_WS_D vs Root_WS_C` = get_genes(list_ratio_stats_global$`Root_WS_D_1.1_rep_bio_3 vs Root_WS_C_1.1_rep_bio_3`, sign_i = "Down"),
  `Root_WS_D vs Root_WW_D` = get_genes(list_ratio_stats_global$`Root_WS_D_1.1_rep_bio_4 vs Root_WW_D_1.1_rep_bio_4`, sign_i = "Down"),
  
  `Nodule_WW_D vs Nodule_WW_C` = get_genes(list_ratio_stats_global$`Nodule_WW_D_1.1_Rep_Bio_1 vs Nodule_WW_C_1.1_Rep_Bio_1`, sign_i = "Down"),
  `Nodule_WS_C vs Nodule_WW_C` = get_genes(list_ratio_stats_global$`Nodule_WS_C_1.1_Rep_Bio_2 vs Nodule_WW_C_1.1_Rep_Bio_2`, sign_i = "Down"),
  `Nodule_WS_D vs Nodule_WS_C` = get_genes(list_ratio_stats_global$`Nodule_WS_D_1.1_Rep_Bio_3 vs Nodule_WS_C_1.1_Rep_Bio_3`, sign_i = "Down"),
  `Nodule_WS_D vs Nodule_WW_D` = get_genes(list_ratio_stats_global$`Nodule_WS_D_1.1_Rep_Bio_4 vs Nodule_WW_D_1.1_Rep_Bio_4`, sign_i = "Down")
)


m_down <- ComplexHeatmap::make_comb_mat(lt_down, mode = "distinct")

df_DEG_condition <- rbind(create_presence_matrix(lt_up, "Up"),create_presence_matrix(lt_down, "Down")) %>% 
  mutate(Sign=as.factor(Sign),
         Sign=relevel(Sign,ref="Up")
         )

upsetplot_condition_all <-upsetplot_condition_merge_sign_all(df_deregulate_condition = df_DEG_condition, 
                                                     title_i = paste0("\n Deregulated genes"),
                                                     color_palette = c('Down' = "#1d4877", 'Up' = "#ee3e32")
                                                     )

# Export figure ####
fig_export(path = here::here("report/microarray/plot/upsetplot/upset_up_down_condition_all"), upsetplot_condition_all, height = 4, width = 14.5*10, res_i = 600)

set.seed(123)
p_down <- plot(eulerr::euler(lt_down, shape = "ellipse"),
     quantities = TRUE, fills = c("#6f1926", "#de324c",  "#f4895f",  "#f8e16f", "#95cf92",  "#369acc",  "#9656a2", "#cbabd1" ),
     legend = list(side = "right"),
     main = list(label = "Euler diagram of down regulated genes", cex = 1.3))
p_up <- plot(eulerr::euler(lt_up, shape = "ellipse"),
     quantities = TRUE, fills = c("#6f1926", "#de324c",  "#f4895f",  "#f8e16f", "#95cf92",  "#369acc",  "#9656a2", "#cbabd1" ),
     legend = list(side = "right"),
     main = list(label = "Euler diagram of upregulated genes", cex = 1.3))

fig_export(path = here::here("report/microarray/plot/venn_diagram/euler_down_all"), p_down, height = 21/3, width = 29.7/3, res_i = 600)
fig_export(path = here::here("report/microarray/plot/venn_diagram/euler_up_all"), p_up, height = 21/3, width = 29.7/3, res_i = 600)

## simple plot (right of the upsetplot)

df_down <- tibble(
  vector = names(lt_down),
  taille  = sapply(lt_down, length),
  sign  = "Down"
)

df_up <- tibble(
  vector = names(lt_up),
  taille  = sapply(lt_up, length),
  sign  = "Up"
)



# Création d'un data frame pour lt_up
df_x <- bind_rows(df_down, df_up) %>%
  group_by(vector) %>%
  mutate(
    total   = sum(taille),
    percent = round(100 * taille / total)
  ) %>%
  ungroup() %>%
  mutate(sign = fct_relevel(sign, c("Up", "Down")),
         vector =  forcats::fct_relevel(vector, "Root_WW_D vs Root_WW_C",  "Root_WS_C vs Root_WW_C", "Root_WS_D vs Root_WS_C", "Root_WS_D vs Root_WW_D",
                                            "Nodule_WW_D vs Nodule_WW_C", "Nodule_WS_C vs Nodule_WW_C", "Nodule_WS_D vs Nodule_WS_C",  "Nodule_WS_D vs Nodule_WW_D")) 

# Plot
px <- ggplot(df_x, aes(x = vector, y = taille, fill = sign)) +
  # First layer: stacked columns by Up/Down
  geom_col() +
  geom_text(
    aes(label = paste0(taille, "\n(", percent, "%)")),
    position = position_stack(vjust = 0.5),
    color = "white", 
    size  = 3
  ) +
  theme_minimal() +
  labs(
    x     = "Comparison",
    y     = "Number of genes deregulated",
    title = "Number of genes deregulated"
  ) +
   
  coord_flip() +
  # First scale: fill for Up/Down
  scale_fill_manual(
    values = c('Down' = "#1d4877", 'Up' = "#ee3e32"),
    name   = "Sign"
  ) +
   scale_x_discrete(limits = rev(levels(df_x$vector)))+ 
  # Reset the fill scale so we can define a new one for the tile
  new_scale_fill() +

  # Second layer: a tile for each comparison at y = -1
  geom_tile(
    aes(
      x     = vector,
      y     = -600,  # negative y just for demonstration
      fill  = vector,
      width = 0.88,
      height= 400
    ),
    data          = df_x,
    color         = "black",
    alpha         = 1,
    inherit.aes   = FALSE
  ) +
  # Second scale: fill for the 'comparison'
  scale_fill_manual(
    values = c("#6f1926", "#de324c", "#f4895f", "#f8e16f",
               "#95cf92", "#369acc", "#9656a2", "#cbabd1"),
    name   = "Comparison"
  ) ; px
  # Optionally, expand the y-limits to see the tile better
  # (since we placed it at y=-1)
  #scale_y_continuous(expand = expansion(mult = c(0.1, 0.1)))
fig_export(path = here::here("report/microarray/plot/barplot_differential_expression"), px, height = 3.5, width = 12, res_i = 600)

# Export results for Gene Enrichment ####
save(m_up, m_down, lt_down, lt_up, file = here::here("data/microarray/output/upset_result_condition_sign_all.RData"))

2.4.2 Distinction between the two compartment

2.4.2.1 Root

Code 
# import data
load(file =  here::here("data/microarray/output/ratio_stats_microarray_root_nodule_FABER.RData"))

# create function to have list of genes
get_genes<- function(df_comparison_i, lfc_lim = 0, padj_lim= 0.05, sign_i){
  v_gene = df_comparison_i %>% 
    mutate(sign= ifelse(logFC>0, "Up", "Down")) %>% 
    filter(
      sign == sign_i,
      abs(logFC)>abs(lfc_lim), 
      BH<= padj_lim 
           ) %>% 
    pull(id_probe)
  return(v_gene)
}

 # creation of the list
lt_up <- list(
  `Root_WW_D vs Root_WW_C` = get_genes(list_ratio_stats_global$`Root_WW_D_1.1_rep_bio_1 vs Root_WW_C_1.1_rep_bio_1`, sign_i = "Up"),
  `Root_WS_C vs Root_WW_C` = get_genes(list_ratio_stats_global$`Root_WS_C_1.1_rep_bio_2 vs Root_WW_C_1.1_rep_bio_2`, sign_i = "Up"),
  `Root_WS_D vs Root_WS_C` = get_genes(list_ratio_stats_global$`Root_WS_D_1.1_rep_bio_3 vs Root_WS_C_1.1_rep_bio_3`, sign_i = "Up"),
  `Root_WS_D vs Root_WW_D` = get_genes(list_ratio_stats_global$`Root_WS_D_1.1_rep_bio_4 vs Root_WW_D_1.1_rep_bio_4`, sign_i = "Up")
)

m_up <- ComplexHeatmap::make_comb_mat(lt_up, mode = "distinct")

# Down-regulated genes
lt_down <- list(
  `Root_WW_D vs Root_WW_C` = get_genes(list_ratio_stats_global$`Root_WW_D_1.1_rep_bio_1 vs Root_WW_C_1.1_rep_bio_1`, sign_i = "Down"),
  `Root_WS_C vs Root_WW_C` = get_genes(list_ratio_stats_global$`Root_WS_C_1.1_rep_bio_2 vs Root_WW_C_1.1_rep_bio_2`, sign_i = "Down"),
  `Root_WS_D vs Root_WS_C` = get_genes(list_ratio_stats_global$`Root_WS_D_1.1_rep_bio_3 vs Root_WS_C_1.1_rep_bio_3`, sign_i = "Down"),
  `Root_WS_D vs Root_WW_D` = get_genes(list_ratio_stats_global$`Root_WS_D_1.1_rep_bio_4 vs Root_WW_D_1.1_rep_bio_4`, sign_i = "Down")
)


m_down <- ComplexHeatmap::make_comb_mat(lt_down, mode = "distinct")

df_DEG_condition <- rbind(create_presence_matrix(lt_up, "Up"),create_presence_matrix(lt_down, "Down")) %>% 
  mutate(Sign=as.factor(Sign),
         Sign=relevel(Sign,ref="Up")
         )

upsetplot_condition_root <-upsetplot_condition_merge_sign_all(df_deregulate_condition = df_DEG_condition, 
                                                     title_i = paste0("\n Deregulated genes"),
                                                     color_palette = c('Down' = "#1d4877", 'Up' = "#ee3e32")
                                                     )

# Export figure ####
fig_export(path = here::here("report/microarray/plot/upsetplot/upset_up_down_condition_root"), upsetplot_condition_root, height = 4, width = 14.5, res_i = 600)

set.seed(123)
p_down <- plot(eulerr::euler(lt_down, shape = "ellipse"),
     quantities = TRUE, fills = c("#6f1926", "#de324c",  "#f4895f",  "#f8e16f"),
     legend = list(side = "right"),
     main = list(label = "Euler diagram of down regulated genes for root", cex = 1.3))
p_up <- plot(eulerr::euler(lt_up, shape = "ellipse"),
     quantities = TRUE, fills = c("#6f1926", "#de324c",  "#f4895f",  "#f8e16f"),
     legend = list(side = "right"),
     main = list(label = "Euler diagram of upregulated genes for root", cex = 1.3))

fig_export(path = here::here("report/microarray/plot/venn_diagram/euler_down_root"), p_down, height = 21/3, width = 29.7/3, res_i = 600)
fig_export(path = here::here("report/microarray/plot/venn_diagram/euler_up_root"), p_up, height = 21/3, width = 29.7/3, res_i = 600)

# Export results for Gene enrichment ####
#save(m_up, m_down,file = here::here("data/rnaseq/output/genotype/root_m_clim_sign.RData"))

2.4.2.2 Nodule

Code 
# import data
load(file =  here::here("data/microarray/output/ratio_stats_microarray_root_nodule_FABER.RData"))

# create function to have list of genes
get_genes<- function(df_comparison_i, lfc_lim = 0, padj_lim= 0.05, sign_i){
  v_gene = df_comparison_i %>% 
    mutate(sign= ifelse(logFC>0, "Up", "Down")) %>% 
    filter(
      sign == sign_i,
      abs(logFC)>abs(lfc_lim), 
      BH<= padj_lim 
           ) %>% 
    pull(id_probe)
  return(v_gene)
}

# get_genes(df_comparison_i = list_ratio_stats$`Mut_SS_W78*_RepBio_4 vs WT_SS_WT1_RepBio_4`, sign_i = "Up")

 # creation of the list
lt_up <- list(
  `Nodule_WW_D vs Nodule_WW_C` = get_genes(list_ratio_stats_global$`Nodule_WW_D_1.1_Rep_Bio_1 vs Nodule_WW_C_1.1_Rep_Bio_1`, sign_i = "Up"),
  `Nodule_WS_C vs Nodule_WW_C` = get_genes(list_ratio_stats_global$`Nodule_WS_C_1.1_Rep_Bio_2 vs Nodule_WW_C_1.1_Rep_Bio_2`, sign_i = "Up"),
  `Nodule_WS_D vs Nodule_WS_C` = get_genes(list_ratio_stats_global$`Nodule_WS_D_1.1_Rep_Bio_3 vs Nodule_WS_C_1.1_Rep_Bio_3`, sign_i = "Up"),
  `Nodule_WS_D vs Nodule_WW_D` = get_genes(list_ratio_stats_global$`Nodule_WS_D_1.1_Rep_Bio_4 vs Nodule_WW_D_1.1_Rep_Bio_4`, sign_i = "Up")
)

m_up <- ComplexHeatmap::make_comb_mat(lt_up, mode = "distinct")

# Down-regulated genes
lt_down <- list(
  `Nodule_WW_D vs Nodule_WW_C` = get_genes(list_ratio_stats_global$`Nodule_WW_D_1.1_Rep_Bio_1 vs Nodule_WW_C_1.1_Rep_Bio_1`, sign_i = "Down"),
  `Nodule_WS_C vs Nodule_WW_C` = get_genes(list_ratio_stats_global$`Nodule_WS_C_1.1_Rep_Bio_2 vs Nodule_WW_C_1.1_Rep_Bio_2`, sign_i = "Down"),
  `Nodule_WS_D vs Nodule_WS_C` = get_genes(list_ratio_stats_global$`Nodule_WS_D_1.1_Rep_Bio_3 vs Nodule_WS_C_1.1_Rep_Bio_3`, sign_i = "Down"),
  `Nodule_WS_D vs Nodule_WW_D` = get_genes(list_ratio_stats_global$`Nodule_WS_D_1.1_Rep_Bio_4 vs Nodule_WW_D_1.1_Rep_Bio_4`, sign_i = "Down")
)


m_down <- ComplexHeatmap::make_comb_mat(lt_down, mode = "distinct")

df_DEG_condition <- rbind(create_presence_matrix(lt_up, "Up"),create_presence_matrix(lt_down, "Down")) %>% 
  mutate(Sign=as.factor(Sign),
         Sign=relevel(Sign,ref="Up")
         )

upsetplot_condition_nodule <-upsetplot_condition_merge_sign_nodule(df_deregulate_condition = df_DEG_condition, 
                                                     title_i = paste0("\n Deregulated genes"),
                                                     color_palette = c('Down' = "#1d4877", 'Up' = "#ee3e32")
                                                     )

# Export figure ####
fig_export(path = here::here("report/microarray/plot/upsetplot/upset_up_down_condition_nodule"), upsetplot_condition_all, height = 4, width = 14.5*1.2, res_i = 600)

set.seed(123)
p_down <- plot(eulerr::euler(lt_down, shape = "ellipse"),
     quantities = TRUE, fills = c("#95cf92",  "#369acc",  "#9656a2", "#cbabd1" ),
     legend = list(side = "right"),
     main = list(label = "Euler diagram of down regulated genes for nodule", cex = 1.3))
p_up <- plot(eulerr::euler(lt_up, shape = "ellipse"),
     quantities = TRUE, fills = c("#95cf92",  "#369acc",  "#9656a2", "#cbabd1" ),
     legend = list(side = "right"),
     main = list(label = "Euler diagram of upregulated genes for Nodule", cex = 1.3))

fig_export(path = here::here("report/microarray/plot/venn_diagram/euler_down_nodule"), p_down, height = 21/3, width = 29.7/3, res_i = 600)
fig_export(path = here::here("report/microarray/plot/venn_diagram/euler_up_nodule"), p_up, height = 21/3, width = 29.7/3, res_i = 600)

# Export results for Gene enrichment ####
#save(m_up, m_down,file = here::here("data/rnaseq/output/genotype/root_m_clim_sign.RData"))