12 Exploratory analysis of preprocessed microarray data

12.1 Material and methods

Data acquisition was performed on experiment 1.

Bonus: Here you’ll find leaf transcriptomic protocols for this experiment

Leaves from four biological replicates of plants (i.e. independent plants) were used. Each biological replicate corresponds to the three last leaves formed before S deficiency imposition (i.e. vegetative leaves), which were collected 25 Days after Flowering (DAF). The leaves were immediately frozen in liquid nitrogen, then stored at -80°C. RNA was extracted from 100 mg of frozen powder using the RNeasy Plant Mini Kit according to manufacturer’s protocol (Qiagen, Courtaboeuf, France). RNA quality was checked on agarose gel 1.5%, then using the Agilent 2100 Bioanalyzer. The subsequent steps were performed at IPS2 (Institute of Plant Sciences Paris-Saclay). The Pea NimbleGen-microarrays were developed to profile expression of 40795 sequences: 40454 mRNA originating from the PsCameor_Uni_Lowcopy set (Alves‐Carvalho et al. 2015), 323 putative precursors of miRNA predicted in the “Test assembly multiple k-mer” contig set (Alves‐Carvalho et al. 2015), and 18 controls. The Ambion MessageAmpTM II aRNA Amplification Kit was used to amplify sufficient amounts of copy RNA extracted from leaves of the four biological replicates. The Double stranded cDNA synthesis was realized using T7-oligo-dT and the antisense RNA (aRNA) was created by in vitro transcription according to manufacturer’s protocol (Life technologies SAS, Saint Aubin, France). The labeling with Cy3 or Cy5 was performed by reverse transcription of aRNA using labeled nucleotides (Cy3-dUTP or Cy5-dUTP, Perkin-Elmer-NEN Life Science Products). For each sample, the following co-hybridizations were performed: (1) W78* mutant vs. wildtype 1 under +S; (2) E568K mutant vs. wild-type 2 under +S; (3) W78* mutant vs. wild-type 1 under –S; (4) E568K mutant vs. wild-type 2 under –S; (5) wild-types under –S vs. wild-type under +S (for this comparison, the biological replicates under each condition were made of two wildtype 1 plants and two wild-type 2 plants). For each comparison, a dye swap was realized. These probes were spotted in triplicates on the GENOPEA array. The hybridization of labeled samples on the slides, scanning and data normalization were performed as previously described (Lurin et al. 2004).

Differential analysis was based on the log₂ ratios averaged on the dye-swap: the technical replicates were averaged to get one log₂ ratio per biological replicate and these values were used to perform a paired t-test. The raw P-values were adjusted by the Benjamini Hochberg method, which controls the family wise error rate, and probes were considered as differentially expressed when the Benjamini Hochberg Pvalue was <0.05. Transcriptome datasets were deposited in the NCBI Gene Expression Omnibus database with the accession numbers GSE121967. All pea sequences with “PsCam” accession numbers could be retrieved from the pea RNAseq gene atlas at http://bios.dijon. inra.fr/ (PsUniLowCopy data set).

Code

#pkg
library(tidyverse)
library(readxl)
library(affy) #BiocManager::install("affy")  # Affymetrix pre-processing
library(limma)  # two-color pre-processing; differential
library(patchwork)
library(progress)
library(ggrepel)
library(plotly)
library(htmlwidgets)
library("FactoMineR")
library("factoextra")
library("corrplot")
library(missMDA)
library(ComplexHeatmap)
library(knitr)
library(kableExtra)
library(ggnewscale) # to have two scale_fill

# for GO terme
library(clusterProfiler)
library(readxl)
library(GO.db)
library(AnnotationDbi) # BiocManager::install("AnnotationDbi")

# 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"
source(here::here("src/function/microarray/GO_on_different_group.R")) # function that find GO terme

# cosmetics
sulfate_pallet=read_excel(here::here("data/color_palette.xlsm")) %>%
      filter(set == "sulfure_condition") %>%
      dplyr::select(color, treatment) %>%
      pull(color) %>%
      setNames(read_excel(here::here("data/color_palette.xlsm")) %>%
                 filter(set == "sulfure_condition") %>%
                 pull(treatment)
               )

my_color_palette <- read_excel(here::here("data/color_palette.xlsm"))

12.2 Data importation

Code

# for the list of all comparison
comparison=read_excel(here::here("data/microarray/Resultats_4plex-POIS-2014_01_230415_modKG.xlsx"), sheet = "Complete", range = "L10:EF11", 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.csv"))

Code

# show results
read_csv(here::here("data/microarray/output/comparison_microarray.csv")) %>% 
knitr::kable(., caption = "List of all comparison")

List of all comparison
Green	Red
5_1Mut-	1_1WT-
6_2Mut-	2_2WT-
7_3Mut-	3_3WT-
8_4Mut-	4_4WT-
13_5Mut-	9_5WT-
14_6Mut-	10_6WT-
15_7Mut-	11_7WT-
16_8Mut-	12_8WT-
21_1Mut+	17_1WT+
22_2Mut+	18_2WT+
23_3Mut+	19_3WT+
24_4Mut+	20_4WT+
29_5Mut+	25_5WT+
30_6Mut+	26_6WT+
31_7Mut+	27_7WT+
32_8Mut+	28_8WT+
1_1WT-	17_1WT+
4_4WT-	19_3WT+
9_5WT-	25_5WT+
12_8WT-	28_8WT+

The data is very poorly organized in the data table, so it will take some time to import. There are two types of data. Columns that come directly from the chips and correspond to the log2 of the original values but normalized chip by chip between Red and Green (This data will be analyzed later, see this part of the document). And average data and stats corresponding to one of the paired student tests, then corrected between the four biological repeats.

Import RebBio data (only ratio and stats)

Code

source (here::here("src/function/microarray/import_ratio_stats_microarray.R"))

path_i = here::here("data/microarray/Resultats_4plex-POIS-2014_01_230415_modKG.xlsx")
path_info_range_i = here::here("data/microarray/info_range_Resultats_4plex-POIS-2014_01_230415_modKG.xlsx")
begin_end_line <- c("12", "41122")
sheet_i = "Complete"

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

# correct name of list
list_ratio_stats <-execut[["list_ratio_stats"]]
names(list_ratio_stats)<- c("Mut_SD_E568K_RepBio_1 vs WT_SD_WT2_RepBio_1", "Mut_SD_W78*_RepBio_2 vs WT_SD_WT1_RepBio_2", "Mut_SS_E568K_RepBio_3 vs WT_SS_WT2_RepBio_3", "Mut_SS_W78*_RepBio_4 vs WT_SS_WT1_RepBio_4", "WT_SD_RepBio_5 vs WT_SS_RepBio_5") # because there is a mix between the two reference genotypes
# export results ####

save(list_ratio_stats, file =  here::here("data/microarray/output/ratio_stats_microarray_leaf_PeaSulf.RData"))
write_csv(execut[["df_info_sample_clean_compile"]], here::here("data/microarray/output/list_microarray_leaf_PeaSulf_RepBio.csv"))

12.3 Differential expression 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_leaf_PeaSulf.RData"))
source (here::here("src/function/microarray/vulcanoplot_microarray_ratio_stats.R"))

for(i in 1:length(names(list_ratio_stats))){
  # create ggplot
  p = vulcanoplot_microarray_ratio_stats(comparison = names(list_ratio_stats)[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_leaf_PeaSulf.RData"))

For an interactive version : Click here

12.4 Simplified differential analysis

12.4.1 Venn diagram

Code

# import data
load(file =  here::here("data/microarray/output/ratio_stats_microarray_leaf_PeaSulf.RData"))
lfc_lim_i <- 0

# create function to have list of genes
get_genes<- function(df_comparison_i, lfc_lim = lfc_lim_i, padj_lim= 0.05, sign_i){ # change value here to change parameter
  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")

lt_up <- list(
  `SS_WT1 vs SS_W78*` = get_genes(list_ratio_stats$`Mut_SS_W78*_RepBio_4 vs WT_SS_WT1_RepBio_4`, sign_i = "Up",lfc_lim = lfc_lim_i),
  `SS_WT2 vs SS_E568K` = get_genes(list_ratio_stats$`Mut_SS_E568K_RepBio_3 vs WT_SS_WT2_RepBio_3`, sign_i = "Up",lfc_lim = lfc_lim_i),
  `SD_WT1 vs SD_W78*` = get_genes(list_ratio_stats$`Mut_SD_W78*_RepBio_2 vs WT_SD_WT1_RepBio_2`, sign_i = "Up",lfc_lim = lfc_lim_i),
  `SD_WT2 vs SD_E568K` = get_genes(list_ratio_stats$`Mut_SD_E568K_RepBio_1 vs WT_SD_WT2_RepBio_1`, sign_i = "Up",lfc_lim = lfc_lim_i), 
  `SS_WT vs SD_WT` = get_genes(list_ratio_stats$`WT_SD_RepBio_5 vs WT_SS_RepBio_5`, sign_i = "Up")
)

lt_down <- list(
  `SS_WT1 vs SS_W78*` = get_genes(list_ratio_stats$`Mut_SS_W78*_RepBio_4 vs WT_SS_WT1_RepBio_4`, sign_i = "Down",lfc_lim = lfc_lim_i),
  `SS_WT2 vs SS_E568K` = get_genes(list_ratio_stats$`Mut_SS_E568K_RepBio_3 vs WT_SS_WT2_RepBio_3`, sign_i = "Down",lfc_lim = lfc_lim_i),
  `SD_WT1 vs SD_W78*` = get_genes(list_ratio_stats$`Mut_SD_W78*_RepBio_2 vs WT_SD_WT1_RepBio_2`, sign_i = "Down",lfc_lim = lfc_lim_i),
  `SD_WT2 vs SD_E568K` = get_genes(list_ratio_stats$`Mut_SD_E568K_RepBio_1 vs WT_SD_WT2_RepBio_1`, sign_i = "Down",lfc_lim = lfc_lim_i), 
  `SS_WT vs SD_WT` = get_genes(list_ratio_stats$`WT_SD_RepBio_5 vs WT_SS_RepBio_5`, sign_i = "Down")
)

length(lt_down$`SS_WT vs SD_WT`)

lt_up_simplified <- list(
  `Mutant vs WT under SS` = intersect(lt_up$`SS_WT1 vs SS_W78*`, lt_up$`SS_WT2 vs SS_E568K`),
  `Mutant vs WT under SD` = intersect(lt_up$`SD_WT1 vs SD_W78*`, lt_up$`SD_WT2 vs SD_E568K`),
  `SS_WT vs SD_WT` = lt_up$`SS_WT vs SD_WT`
)

lt_down_simplified <- list(
  `Mutant vs WT under SS` = intersect(lt_down$`SS_WT1 vs SS_W78*`, lt_down$`SS_WT2 vs SS_E568K`),
  `Mutant vs WT under SD` = intersect(lt_down$`SD_WT1 vs SD_W78*`, lt_down$`SD_WT2 vs SD_E568K`),
  `SS_WT vs SD_WT` = lt_down$`SS_WT vs SD_WT`
)

p_down <- plot(eulerr::euler(lt_down_simplified, shape = "ellipse"),
     quantities = TRUE, fills = c(sulfate_pallet[1], sulfate_pallet[2], "#9381FF"),
     legend = list(side = "right"),
     main = list(label = paste0("Euler diagram of down regulated genes with LFC ", lfc_lim_i), cex = 1.3))

p_up <- plot(eulerr::euler(lt_up_simplified, shape = "ellipse"),
     quantities = TRUE, fills = c(sulfate_pallet[1], sulfate_pallet[2], "#9381FF"),
     legend = list(side = "right"),
     main = list(label = paste0("Euler diagram of up regulated genes with LFC ", lfc_lim_i), cex = 1.3))

fig_export(path = here::here(paste0("report/microarray/plot/venn_diagram/euler_simplified_down_lfc_",lfc_lim_i)), p_down, height = 21/3, width = 29.7/3, res_i = 600)
fig_export(path = here::here(paste0("report/microarray/plot/venn_diagram/euler_simplified_up_lfc_",lfc_lim_i)), p_up, height = 21/3, width = 29.7/3, res_i = 600)

# diagram venn 
fills <- c(
  as.character(sulfate_pallet[1]), as.character(sulfate_pallet[2]),"SS_WT vs SD_WT"        = "#9381FF"
)

names(lt_down_simplified) <- c("Mutant vs WT (SS)   ", "    Mutant vs WT (SD)", "SS_WT vs SD_WT")
names(lt_up_simplified) <- c("Mutant vs WT (SS)   ", "    Mutant vs WT (SD)", "SS_WT vs SD_WT")
ggvenn_plot_down = ggvenn::ggvenn(
  lt_down_simplified,
  fill_color = unname(fills),
  fill_alpha = 0.55,
  stroke_size = 0.4
) +
  ggtitle(paste0(
    "Down-regulated genes")
  ) +
  theme(plot.title = element_text(face = "bold", size = 15),
        legend.position = "none")

ggvenn_plot_up = ggvenn::ggvenn(
  lt_up_simplified,
  fill_color = unname(fills),
  fill_alpha = 0.55,
  stroke_size = 0.4
) +
  ggtitle(paste0(
    "Up-regulated genes")
  ) +
  theme(plot.title = element_text(face = "bold", size = 15),
        legend.position = "none")


print(ggvenn_plot_down) ; print(ggvenn_plot_up)

fig_export(path = here::here(paste0("report/microarray/plot/venn_diagram/venn_simplified_up_lfc_",lfc_lim_i)), ggvenn_plot_up, height = 21/3, width = 29.7/3, res_i = 600)
fig_export(path = here::here(paste0("report/microarray/plot/venn_diagram/venn_simplified_down_lfc_",lfc_lim_i)), ggvenn_plot_down, height = 21/3, width = 29.7/3, res_i = 600)

# save plot for article
venn_list <- list(
  up   = ggvenn_plot_up,
  down = ggvenn_plot_down
)

save(venn_list,file = here::here(paste0("report/microarray/plot/venn_diagram/venn_simplified_lfc_",lfc_lim_i,".RData")))

# export results ####
save(lt_up_simplified, lt_down_simplified,file = here::here(paste0("data/microarray/output/venn_result_simplified_condition_sign_lfc_",lfc_lim_i,".RData")))

12.4.2 GO enrichment

Code

#parameter
lfc_lim_i = 0
# collect data and add Psat to PsCAM
load(here::here("data/microarray/output/raw_data_microarray_leaf_PeaSulf.RData"))
load(file = here::here(paste0("data/microarray/output/venn_result_simplified_condition_sign_lfc_",lfc_lim_i,".RData")))

df_info_gene <- RG$genes %>% 
  dplyr::rename(PsCam = ID) %>% 
  left_join(.,read_excel(here::here("data/microarray/Resultats_4plex-POIS-2014_01_230415_modKG.xlsx"), sheet = "Pscam_to_Psat_v1c_mrna_besthit-", col_names = T), by = "PsCam")

cluster_info <- rbind(purrr::imap_dfr(lt_down_simplified, ~ tibble(PsCam = .x, comparison = .y)) %>% 
                        mutate(cluster = paste0("(Down) ", comparison)), 
                      purrr::imap_dfr(lt_up_simplified, ~ tibble(PsCam = .x, comparison = .y)) %>% 
                        mutate(cluster = paste0("(Up) ", comparison))
) %>% full_join(., df_info_gene, by = "PsCam") %>% 
  mutate(ID = Psat) %>% 
  drop_na(comparison, ID)

df_GO = GO_on_different_group(functional_roles = "BP", #"BP"
                      group_info = cluster_info,
                      group = "cluster",
                      ID = "ID",
                      top = 10
                      )

# Merge this with your original data to fill missing combinations
Cluster_selected_GO_filled <-df_GO %>% # on s'arrete ici pour la fonction. 
  dplyr::rename(cluster = group) %>% 
  mutate(cluster = factor(cluster,
  levels = unique(cluster)), 
  comparison =   str_trim(str_remove(cluster, "\\s*\\(.*?\\)")),
  sign = str_extract(cluster, "(?<=\\().*?(?=\\))")  
  ) %>% 
  mutate(comparison =  forcats::fct_relevel(comparison, "Mutant vs WT under SS",  "Mutant vs WT under SD", "SS_WT vs SD_WT"))

#vector_color <- c("#067B5B", "#6AD6AD", "#B87100", "#FFAC5C", "#9381FF", "#067B5B", "#6AD6AD", "#B87100", "#FFAC5C", "#9381FF")
# Plot the enrichment by GO
px <- ggplot(Cluster_selected_GO_filled, aes(x = cluster, y = Description_GO, fill = `|-log10(Pval)|`)) + 
  geom_tile(color = "black") + 
  scale_fill_gradient2(low = "white", high = "red", na.value = "white", limits = c(0, NA)) +
  theme_minimal() +
  theme(axis.text = element_text(size = 8, colour = "black"),
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
        panel.grid.major = element_blank(),
        plot.background = element_rect(fill = "white", colour = "white")) +
  labs(fill = "-log10(FDR)",
       x= "Biological process",
       y = "Cluster", 
       title = "GO terms for all gene deregulated") +
  #scale_size_manual(values = c(dot = 2, no_dot = NA), guide = "none")+
  new_scale_fill() +
 scale_fill_manual(
    values =  c("Mutant vs WT under SS" = as.character(sulfate_pallet[1]),
                "Mutant vs WT under SD" = as.character(sulfate_pallet[2]),
                "SS_WT vs SD_WT" = "#9381FF"),
    name = "Comparison"
  ) +
geom_tile(
    aes(
      x = cluster,
      y = -1,     # If you truly want it at a negative y-value, 
                      # ensure your y-scale is continuous or can handle this
      fill = comparison,
      width = 0.90,
      height = 0.80
    ),
    data = Cluster_selected_GO_filled,
    color = "black",
    alpha = 1,
    inherit.aes = FALSE
  )+
  new_scale_fill() +
  scale_fill_manual(
    values = c('Down' = "#1d4877", 'Up' = "#ee3e32"),
    name = "Sign"
  ) +
geom_tile(
    aes(
      x = cluster,
      y = -0.05,     # If you truly want it at a negative y-value,
                      # ensure your y-scale is continuous or can handle this
      fill = sign,
      width = 0.90,
      height = 0.80
    ),
    data = Cluster_selected_GO_filled,
    color = "black",
    alpha = 1,
    inherit.aes = FALSE
  ) ; px

# export 
fig_export(path =paste0("report/microarray/plot/GO/GO_BP_all_simplified_absolute_lfc_",lfc_lim_i), plot_x = px, height_i = 13, width_i = 9, res = 600)

# export for article
load(file = here::here(paste0("report/microarray/plot/venn_diagram/venn_simplified_lfc_",lfc_lim_i,".RData"))) # load the list of the two venne diagram

venn_list$up   <- venn_list$up   + labs(tag = "A")
venn_list$down <- venn_list$down + labs(tag = "B")
px             <- px             + labs(tag = "C")


venn_stack <- venn_list$up / venn_list$down 

combined_plot <- (venn_stack | px) +     # | = côte-à-côte
  plot_layout(widths = c(2, 1), guides = "collect") &
  theme(plot.margin = margin(5, 5, 5, 5))+
   theme(plot.margin = margin(5, 5, 5, 5),
        plot.tag      = element_text(size = 14))

# ─────────────────────────────────────────────────────────
# 3) APERÇU RAPIDE
# ─────────────────────────────────────────────────────────
print(combined_plot)

fig_export(path =paste0("report/microarray/plot/fig_lfc_",lfc_lim_i), plot_x = combined_plot, height_i = 13, width_i = 15, res = 600, format = "png")

12.5 Differentiation between the two genotypes

12.5.1 UpsetPlot

Code

# import data
load(file =  here::here("data/microarray/output/ratio_stats_microarray_leaf_PeaSulf.RData"))
lfc_lim_i <- 1
  
# create function to have list of genes
get_genes<- function(df_comparison_i, lfc_lim = lfc_lim_i, padj_lim= 0.05, sign_i){ # change value here to change parameter
  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(
  `SS_WT1 vs SS_W78*` = get_genes(list_ratio_stats$`Mut_SS_W78*_RepBio_4 vs WT_SS_WT1_RepBio_4`, sign_i = "Up"),
  `SS_WT2 vs SS_E568K` = get_genes(list_ratio_stats$`Mut_SS_E568K_RepBio_3 vs WT_SS_WT2_RepBio_3`, sign_i = "Up"),
  `SD_WT1 vs SD_W78*` = get_genes(list_ratio_stats$`Mut_SD_W78*_RepBio_2 vs WT_SD_WT1_RepBio_2`, sign_i = "Up"),
  `SD_WT2 vs SD_E568K` = get_genes(list_ratio_stats$`Mut_SD_E568K_RepBio_1 vs WT_SD_WT2_RepBio_1`, sign_i = "Up"), 
  `SS_WT vs SD_WT` = get_genes(list_ratio_stats$`WT_SD_RepBio_5 vs WT_SS_RepBio_5`, sign_i = "Up")
)

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

# Down-regulated genes
lt_down <- list(
  `SS_WT1 vs SS_W78*` = get_genes(list_ratio_stats$`Mut_SS_W78*_RepBio_4 vs WT_SS_WT1_RepBio_4`, sign_i = "Down"),
  `SS_WT2 vs SS_E568K` = get_genes(list_ratio_stats$`Mut_SS_E568K_RepBio_3 vs WT_SS_WT2_RepBio_3`, sign_i = "Down"),
  `SD_WT1 vs SD_W78*` = get_genes(list_ratio_stats$`Mut_SD_W78*_RepBio_2 vs WT_SD_WT1_RepBio_2`, sign_i = "Down"),
  `SD_WT2 vs SD_E568K` = get_genes(list_ratio_stats$`Mut_SD_E568K_RepBio_1 vs WT_SD_WT2_RepBio_1`, sign_i = "Down"), 
  `SS_WT vs SD_WT` = get_genes(list_ratio_stats$`WT_SD_RepBio_5 vs WT_SS_RepBio_5`, sign_i = "Down")
)

m_down <- ComplexHeatmap::make_comb_mat(lt_down, mode = "distinct")
# 
# create_upset_plot_rnaseq(m_up, selected_columns = c("100000", "010000", "001000","000100","000010","000001","100100","010010","001001","101101","011011","101101","111111"), 
#                   "Up", organ_i)
# 
# create_upset_plot_rnaseq(m_down, selected_columns = c("100000", "010000", "001000","000100","000010","000001","100100","010010","001001","101101","011011","101101","111111"), 
#                   "Down", organ_i)

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 <-upsetplot_condition_merge_sign(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(paste0("report/microarray/plot/upsetplot/upset_up_down_condition_lfc_",lfc_lim_i)), upsetplot_condition, height = 4, width = 14.5, res_i = 600)

set.seed(123)
p_down <- plot(eulerr::euler(lt_down, shape = "ellipse"),
     quantities = TRUE, fills = c(lighten(as.character(sulfate_pallet[1]), amount = -.2), lighten(as.character(sulfate_pallet[1]), amount = .5), lighten(as.character(sulfate_pallet[2]), amount = -.2), lighten(as.character(sulfate_pallet[2]), amount = .3), "#9381FF"),
     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(lighten(as.character(sulfate_pallet[1]), amount = -.2), lighten(as.character(sulfate_pallet[1]), amount = .5), lighten(as.character(sulfate_pallet[2]), amount = -.2), lighten(as.character(sulfate_pallet[2]), amount = .3), "#9381FF"),
     legend = list(side = "right"),
     main = list(label = "Euler diagram of upregulated genes", cex = 1.3))

fig_export(path = here::here(paste0("report/microarray/plot/venn_diagram/euler_down_lfc_",lfc_lim_i)), p_down, height = 21/3, width = 29.7/3, res_i = 600)
fig_export(path = here::here(paste0("report/microarray/plot/venn_diagram/euler_up_lfc_",lfc_lim_i)), p_up, height = 21/3, width = 29.7/3, res_i = 600)

# verification if i had the same results as in excel.
length(intersect(lt_up$`SD_WT1 vs SD_W78*`, lt_up$`SD_WT2 vs SD_E568K`))
length(intersect(lt_down$`SD_WT1 vs SD_W78*`, lt_down$`SD_WT2 vs SD_E568K`))
length(intersect(lt_up$`SS_WT1 vs SS_W78*`, lt_up$`SS_WT2 vs SS_E568K`))
length(intersect(lt_down$`SS_WT1 vs SS_W78*`, lt_down$`SS_WT2 vs SS_E568K`))

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

With log fold change of 0
With log fold change of 1

Which group does the psult4 gene belong to?

Code

load(file = here::here("data/microarray/output/upset_result_condition_sign_lfc_0.RData"))

my_gene_of_interest <- "PsCam042688"

in_up <- names(keep(lt_up, ~ my_gene_of_interest %in% .x))
in_down <- names(keep(lt_down, ~ my_gene_of_interest %in% .x))

df_up <- data.frame("In Up Regulated" = in_up, check.names = FALSE)
df_down <- data.frame("In Down Regulated" = in_down, check.names = FALSE)

# Generate HTML for each table with kable
table1_html <- kable(df_up, format = "html", table.attr = "class='table'") %>% 
  kable_styling(full_width = FALSE)
table2_html <- kable(df_down, format = "html", table.attr = "class='table'") %>% 
  kable_styling(full_width = FALSE)

# Combine the two panels side by side using a flexbox divider
html_output <- paste0(
  "<div style='display: flex; justify-content: center; gap: 20px;'>",
    "<div>", table2_html, "</div>",
    "<div>", table1_html, "</div>",
  "</div>"
)

# Display HTML without escapement
knitr::asis_output(html_output)

In Down Regulated
SS_WT1 vs SS_W78*
SS_WT2 vs SS_E568K
SD_WT1 vs SD_W78*
SD_WT2 vs SD_E568K

In Up Regulated
SS_WT vs SD_WT

Which group does the psult4 gene belong to?

Code

load(file = here::here("data/microarray/output/upset_result_condition_sign_lfc_1.RData"))

my_gene_of_interest <- "PsCam042688"

in_up <- names(keep(lt_up, ~ my_gene_of_interest %in% .x))
in_down <- names(keep(lt_down, ~ my_gene_of_interest %in% .x))

df_up <- data.frame("In Up Regulated" = in_up, check.names = FALSE)
df_down <- data.frame("In Down Regulated" = in_down, check.names = FALSE)

# Generate HTML for each table with kable
table1_html <- kable(df_up, format = "html", table.attr = "class='table'") %>% 
  kable_styling(full_width = FALSE)
table2_html <- kable(df_down, format = "html", table.attr = "class='table'") %>% 
  kable_styling(full_width = FALSE)

# Combine the two panels side by side using a flexbox divider
html_output <- paste0(
  "<div style='display: flex; justify-content: center; gap: 20px;'>",
    "<div>", table2_html, "</div>",
    "<div>", table1_html, "</div>",
  "</div>"
)

# Display HTML without escapement
knitr::asis_output(html_output)

In Down Regulated
SS_WT1 vs SS_W78*
SD_WT1 vs SD_W78*

In Up Regulated
SS_WT vs SD_WT

12.5.2 GO enrichment

Code

#parameter
lfc_lim_i = 1
# collect data and add Psat to PsCAM
load(here::here("data/microarray/output/raw_data_microarray_leaf_PeaSulf.RData"))
load(file = here::here(paste0("data/microarray/output/upset_result_condition_sign_lfc_",lfc_lim_i,".RData")))

df_info_gene <- RG$genes %>% 
  dplyr::rename(PsCam = ID) %>% 
  left_join(.,read_excel(here::here("data/microarray/Resultats_4plex-POIS-2014_01_230415_modKG.xlsx"), sheet = "Pscam_to_Psat_v1c_mrna_besthit-", col_names = T), by = "PsCam")

cluster_info <- rbind(purrr::imap_dfr(lt_down, ~ tibble(PsCam = .x, comparison = .y)) %>% 
                        mutate(cluster = paste0("(Down) ", comparison)), 
                      purrr::imap_dfr(lt_up, ~ tibble(PsCam = .x, comparison = .y)) %>% 
                        mutate(cluster = paste0("(Up) ", comparison))
) %>% full_join(., df_info_gene, by = "PsCam") %>% 
  mutate(ID = Psat) %>% 
  drop_na(comparison, ID)

df_GO = GO_on_different_group(functional_roles = "BP",
                      group_info = cluster_info,
                      group = "cluster",
                      ID = "ID",
                      top = 10
                      )

# Merge this with your original data to fill missing combinations
Cluster_selected_GO_filled <-df_GO %>% # on s'arrete ici pour la fonction. 
  dplyr::rename(cluster = group) %>% 
  mutate(cluster = factor(cluster,
  levels = unique(cluster)), 
  comparison =   str_trim(str_remove(cluster, "\\s*\\(.*?\\)")),
  sign = str_extract(cluster, "(?<=\\().*?(?=\\))")  
  ) %>% 
  mutate(comparison =  forcats::fct_relevel(comparison, "SS_WT1 vs SS_W78*",  "SS_WT2 vs SS_E568K", "SD_WT1 vs SD_W78*", "SD_WT2 vs SD_E568K", "SS_WT vs SD_WT"))

#vector_color <- c("#067B5B", "#6AD6AD", "#B87100", "#FFAC5C", "#9381FF", "#067B5B", "#6AD6AD", "#B87100", "#FFAC5C", "#9381FF")
# Plot the enrichment by GO
px <- ggplot(Cluster_selected_GO_filled, aes(x = cluster, y = Description_GO, fill = `|-log10(Pval)|`)) + 
  geom_tile(color = "black") + 
  scale_fill_gradient2(low = "white", high = "red", na.value = "white", limits = c(0, NA)) +
  theme_minimal() +
  theme(axis.text = element_text(size = 8, colour = "black"),
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
        panel.grid.major = element_blank(),
        plot.background = element_rect(fill = "white", colour = "white")) +
  labs(fill = "-log10(FDR)",
       x= "Biological process",
       y = "Cluster", 
       title = "GO terme for all gene deregulated") +
  #scale_size_manual(values = c(dot = 2, no_dot = NA), guide = "none")+
  new_scale_fill() +
 scale_fill_manual(
    values =  c(lighten(as.character(sulfate_pallet[1]), amount = -.2),
  lighten(as.character(sulfate_pallet[1]), amount = .5),
  lighten(as.character(sulfate_pallet[2]), amount = -.2),
  lighten(as.character(sulfate_pallet[2]), amount = .3),
  "#9381FF"
),
    name = "Comparison"
  ) +
geom_tile(
    aes(
      x = cluster,
      y = -1,     # If you truly want it at a negative y-value, 
                      # ensure your y-scale is continuous or can handle this
      fill = comparison,
      width = 0.90,
      height = 0.80
    ),
    data = Cluster_selected_GO_filled,
    color = "black",
    alpha = 1,
    inherit.aes = FALSE
  )+
  new_scale_fill() +
  scale_fill_manual(
    values = c('Down' = "#1d4877", 'Up' = "#ee3e32"),
    name = "Sign"
  ) +
geom_tile(
    aes(
      x = cluster,
      y = -0.05,     # If you truly want it at a negative y-value,
                      # ensure your y-scale is continuous or can handle this
      fill = sign,
      width = 0.90,
      height = 0.80
    ),
    data = Cluster_selected_GO_filled,
    color = "black",
    alpha = 1,
    inherit.aes = FALSE
  ) ; px

# export 
fig_export(path =paste0("report/microarray/plot/GO/GO_BP_all_absolute_lfc_",lfc_lim_i), plot_x = px, height_i = 13, width_i = 9, res = 600)

Code

# euler (or upset)
#parameter
lfc_lim_i = 1
# collect data and add Psat to PsCAM
load(here::here("data/microarray/output/raw_data_microarray_leaf_PeaSulf.RData"))
load(file = here::here(paste0("data/microarray/output/upset_result_condition_sign_lfc_",lfc_lim_i,".RData")))

df_info_gene <- RG$genes %>% 
  dplyr::rename(PsCam = ID) %>% 
  left_join(.,read_excel(here::here("data/microarray/Resultats_4plex-POIS-2014_01_230415_modKG.xlsx"), sheet = "Pscam_to_Psat_v1c_mrna_besthit-", col_names = T), by = "PsCam")

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")
  )

gene_id_col <- colnames(df_DEG_condition)[1]
v_conditions <- colnames(df_DEG_condition)[2:(ncol(df_DEG_condition) - 1)]

df_gene_intersections <- df_DEG_condition %>%
    rowwise() %>%
  filter(Sign =="Down") %>% 
  mutate(intersection = paste0(c_across(all_of(v_conditions)), collapse = "")) %>%
  ungroup()

intersection_list <- list(
  I1  = c('SS_WT vs SD_WT', 'SS_WT1 vs SS_W78*', 'SS_WT2 vs SS_E568K', 'SD_WT1 vs SD_W78*', 'SD_WT2 vs SD_E568K'),
  I2  = c('SS_WT1 vs SS_W78*', 'SS_WT2 vs SS_E568K', 'SD_WT1 vs SD_W78*', 'SD_WT2 vs SD_E568K'),
  I3  = c('SS_WT1 vs SS_W78*', 'SS_WT2 vs SS_E568K', 'SS_WT vs SD_WT'),
  I4  = c('SD_WT1 vs SD_W78*', 'SD_WT2 vs SD_E568K', 'SS_WT vs SD_WT'),
  I5  = c('SS_WT1 vs SS_W78*', 'SS_WT2 vs SS_E568K'),
  I6  = c('SD_WT1 vs SD_W78*', 'SD_WT2 vs SD_E568K'),
  I7  = c('SS_WT1 vs SS_W78*', 'SD_WT1 vs SD_W78*'),
  I8  = c('SS_WT2 vs SS_E568K', 'SD_WT2 vs SD_E568K'),
  I9  = c('SS_WT vs SD_WT', 'SS_WT1 vs SS_W78*'),
  I10 = c('SS_WT vs SD_WT', 'SS_WT2 vs SS_E568K'),
  I11 = c('SS_WT vs SD_WT', 'SD_WT1 vs SD_W78*'),
  I12 = c('SS_WT vs SD_WT', 'SD_WT2 vs SD_E568K'),
  I13 = 'SS_WT1 vs SS_W78*',
  I14 = 'SS_WT2 vs SS_E568K',
  I15 = 'SD_WT1 vs SD_W78*',
  I16 = 'SD_WT2 vs SD_E568K',
  I17 = 'SS_WT vs SD_WT'
)

intersection_patterns <- map(intersection_list, function(conds) {
  sapply(v_conditions, function(x) as.integer(x %in% conds)) %>% paste0(collapse = "")
})

signs <- c("Down", "Up")
df_gene_intersections_by_sign <- map(signs, function(s) {
  df_DEG_condition %>%
    filter(Sign == s) %>%
    rowwise() %>%
    mutate(intersection = paste0(c_across(all_of(v_conditions)), collapse = "")) %>%
    ungroup()
})
names(df_gene_intersections_by_sign) <- signs


intersection_genes_all <- list()
for(i in names(intersection_list)) {
  for(s in signs) {
    pattern <- intersection_patterns[[i]]
    genes <- df_gene_intersections_by_sign[[s]] %>%
      filter(intersection == pattern) %>%
      pull(!!sym(gene_id_col))
    # Nommer l'élément comme "I1 (Up)" ou "I1 (Down)", etc.
    intersection_genes_all[[paste0("(", s, ") ", i)]] <- genes
  }
}

# and then same hase before
cluster_info <- purrr::imap_dfr(intersection_genes_all, ~ tibble(PsCam = .x, intersection = .y)) %>% 
                        mutate(cluster = intersection) %>% full_join(., df_info_gene, by = "PsCam") %>% 
  mutate(ID = Psat) %>% 
  drop_na(intersection, ID)

df_GO = GO_on_different_group(functional_roles = "BP",
                      group_info = cluster_info,
                      group = "cluster",
                      ID = "ID",
                      top = 10
                      )

# Merge this with your original data to fill missing combinations
Cluster_selected_GO_filled <- df_GO %>% 
  dplyr::rename(cluster = group) %>% 
  mutate(cluster = factor(cluster,
  levels = unique(cluster)), 
  comparison =   str_trim(str_remove(cluster, "\\s*\\(.*?\\)")),
  sign = str_extract(cluster, "(?<=\\().*?(?=\\))")  
  ) 

#vector_color <- c("#067B5B", "#6AD6AD", "#B87100", "#FFAC5C", "#9381FF", "#067B5B", "#6AD6AD", "#B87100", "#FFAC5C", "#9381FF")
# Plot the enrichment by GO
px <- ggplot(Cluster_selected_GO_filled, aes(x = cluster, y = Description_GO, fill = `|-log10(Pval)|`)) + 
  geom_tile(color = "black") + 
  scale_fill_gradient2(low = "white", high = "red", na.value = "white", limits = c(0, NA)) +
  theme_minimal() +
  theme(axis.text = element_text(size = 8, colour = "black"),
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
        panel.grid.major = element_blank(),
        plot.background = element_rect(fill = "white", colour = "white")) +
  labs(fill = "-log10(FDR)",
       x= "Biological process",
       y = "Cluster", 
       title = "GO terme for all gene deregulated") +
  #scale_size_manual(values = c(dot = 2, no_dot = NA), guide = "none")+
  
  new_scale_fill() +
  scale_fill_manual(
    values = c('Down' = "#1d4877", 'Up' = "#ee3e32"),
    name = "Sign"
  ) +
geom_tile(
    aes(
      x = cluster,
      y = -0.05,     # If you truly want it at a negative y-value,
                      # ensure your y-scale is continuous or can handle this
      fill = sign,
      width = 0.90,
      height = 0.80
    ),
    data = Cluster_selected_GO_filled,
    color = "black",
    alpha = 1,
    inherit.aes = FALSE
  ) ; px

# export 
fig_export(path = paste0("report/microarray/plot/GO/GO_BP_euler_lfc_", lfc_lim_i), plot_x = px, height_i = 20, width_i = 12, res = 600)

GO enrichment on deregulated genes
On Euler diagram (or upset plot)

With log fold change of 0
With log fold change of 1

With log fold change of 0
With log fold change of 1

--- output: html_document editor_options: chunk_output_type: console editor: markdown: wrap: 72 --- ## Exploratory analysis of preprocessed microarray data {#transcriptomic} ## Material and methods Data acquisition was performed on [experiment 1](#experimental_design). <details> <summary style="color: #007bff;"> Bonus: Here you'll find leaf transcriptomic protocols for this experiment </summary> Leaves from four biological replicates of plants (i.e. independent plants) were used. Each biological replicate corresponds to the three last leaves formed before S deficiency imposition (i.e. vegetative leaves), which were collected **25 Days after Flowering (DAF)**. The leaves were immediately frozen in liquid nitrogen, then stored at -80°C. RNA was extracted from 100 mg of frozen powder using the RNeasy Plant Mini Kit according to manufacturer's protocol (Qiagen, Courtaboeuf, France). RNA quality was checked on agarose gel 1.5%, then using the Agilent 2100 Bioanalyzer. The subsequent steps were performed at IPS2 (Institute of Plant Sciences Paris-Saclay). The Pea NimbleGen-microarrays were developed to profile expression of 40795 sequences: 40454 mRNA originating from the PsCameor_Uni_Lowcopy set [@alvescarvalho_fulllength_2015], 323 putative precursors of miRNA predicted in the "Test assembly multiple k-mer" contig set [@alvescarvalho_fulllength_2015], and 18 controls. The Ambion MessageAmpTM II aRNA Amplification Kit was used to amplify sufficient amounts of copy RNA extracted from leaves of the four biological replicates. The Double stranded cDNA synthesis was realized using T7-oligo-dT and the antisense RNA (aRNA) was created by in vitro transcription according to manufacturer's protocol (Life technologies SAS, Saint Aubin, France). The labeling with Cy3 or Cy5 was performed by reverse transcription of aRNA using labeled nucleotides (Cy3-dUTP or Cy5-dUTP, Perkin-Elmer-NEN Life Science Products). For each sample, the following co-hybridizations were performed: (1) W78\* mutant vs. wildtype 1 under +S; (2) E568K mutant vs. wild-type 2 under +S; (3) W78\* mutant vs. wild-type 1 under --S; (4) E568K mutant vs. wild-type 2 under --S; (5) wild-types under --S vs. wild-type under +S (for this comparison, the biological replicates under each condition were made of two wildtype 1 plants and two wild-type 2 plants). For each comparison, a dye swap was realized. These probes were spotted in triplicates on the GENOPEA array. The hybridization of labeled samples on the slides, scanning and data normalization were performed as previously described [@lurin_genome-wide_2004]. Differential analysis was based on the log~2~ ratios averaged on the dye-swap: the technical replicates were averaged to get one log~2~ ratio per biological replicate and these values were used to perform a paired t-test. The raw P-values were adjusted by the Benjamini Hochberg method, which controls the family wise error rate, and probes were considered as differentially expressed when the Benjamini Hochberg Pvalue was \<0.05. Transcriptome datasets were deposited in the NCBI Gene Expression Omnibus database with the accession numbers [GSE121967](http://www.ncbi.nlm.nih.gov/projects/geo/query/acc.cgi?acc=GSE121967). All pea sequences with "PsCam" accession numbers could be retrieved from the pea RNAseq gene atlas at http://bios.dijon. inra.fr/ (PsUniLowCopy data set). </details> ```{r} #pkg library(tidyverse) library(readxl) library(affy) #BiocManager::install("affy") # Affymetrix pre-processing library(limma) # two-color pre-processing; differential library(patchwork) library(progress) library(ggrepel) library(plotly) library(htmlwidgets) library("FactoMineR") library("factoextra") library("corrplot") library(missMDA) library(ComplexHeatmap) library(knitr) library(kableExtra) library(ggnewscale) # to have two scale_fill # for GO terme library(clusterProfiler) library(readxl) library(GO.db) library(AnnotationDbi) # BiocManager::install("AnnotationDbi") # 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" source(here::here("src/function/microarray/GO_on_different_group.R")) # function that find GO terme # cosmetics sulfate_pallet=read_excel(here::here("data/color_palette.xlsm")) %>% filter(set == "sulfure_condition") %>% dplyr::select(color, treatment) %>% pull(color) %>% setNames(read_excel(here::here("data/color_palette.xlsm")) %>% filter(set == "sulfure_condition") %>% pull(treatment) ) my_color_palette <- read_excel(here::here("data/color_palette.xlsm")) ``` ## Data importation{#analyse_data_pretrated} ```{r, eval=FALSE} # for the list of all comparison comparison=read_excel(here::here("data/microarray/Resultats_4plex-POIS-2014_01_230415_modKG.xlsx"), sheet = "Complete", range = "L10:EF11", 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.csv")) ``` ```{r} # show results read_csv(here::here("data/microarray/output/comparison_microarray.csv")) %>% knitr::kable(., caption = "List of all comparison") ``` The data is very poorly organized in the data table, so it will take some time to import. There are two types of data. Columns that come directly from the chips and correspond to the log2 of the original values but normalized chip by chip between Red and Green (This data will be analyzed later, [see this part of the document](B_raw_data_analysis.qmd#analyse_data_normalized)). And average data and stats corresponding to one of the paired student tests, then corrected between the four biological repeats. **Import RebBio data (only ratio and stats)** ```{r, eval=FALSE} source (here::here("src/function/microarray/import_ratio_stats_microarray.R")) path_i = here::here("data/microarray/Resultats_4plex-POIS-2014_01_230415_modKG.xlsx") path_info_range_i = here::here("data/microarray/info_range_Resultats_4plex-POIS-2014_01_230415_modKG.xlsx") begin_end_line <- c("12", "41122") sheet_i = "Complete" execut = import_ratio_stats_microarray(path_i, path_info_range_i, sheet_i, begin_end_line) # correct name of list list_ratio_stats <-execut[["list_ratio_stats"]] names(list_ratio_stats)<- c("Mut_SD_E568K_RepBio_1 vs WT_SD_WT2_RepBio_1", "Mut_SD_W78*_RepBio_2 vs WT_SD_WT1_RepBio_2", "Mut_SS_E568K_RepBio_3 vs WT_SS_WT2_RepBio_3", "Mut_SS_W78*_RepBio_4 vs WT_SS_WT1_RepBio_4", "WT_SD_RepBio_5 vs WT_SS_RepBio_5") # because there is a mix between the two reference genotypes # export results #### save(list_ratio_stats, file = here::here("data/microarray/output/ratio_stats_microarray_leaf_PeaSulf.RData")) write_csv(execut[["df_info_sample_clean_compile"]], here::here("data/microarray/output/list_microarray_leaf_PeaSulf_RepBio.csv")) ``` </details> ## Differential expression of non-normalized data (stats from the platform) ::: callout-warning ## I had to transform the statistics equal to 0 for the vulcanoplot 0 was replaced by 1.10<sup>-14 Because `-log10(0) = Inf` ::: ```{r, eval = F} # /!\ vey slow load(file = here::here("data/microarray/output/ratio_stats_microarray_leaf_PeaSulf.RData")) source (here::here("src/function/microarray/vulcanoplot_microarray_ratio_stats.R")) for(i in 1:length(names(list_ratio_stats))){ # create ggplot p = vulcanoplot_microarray_ratio_stats(comparison = names(list_ratio_stats)[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"))) } ``` ```{r} load(file = here::here("data/microarray/output/ratio_stats_microarray_leaf_PeaSulf.RData")) ``` ::: panel-tabset ### `r gsub("_RepBio_\\d+", "", names(list_ratio_stats)[2])` For an interactive version : [Click here](plot/volcanoplot/volcano_interactive_ratio_stats_2.html) ![](plot/volcanoplot/volcano_interactive_ratio_stats_2.png) ### `r gsub("_RepBio_\\d+", "", names(list_ratio_stats)[1])` For an interactive version : [Click here](plot/volcanoplot/volcano_interactive_ratio_stats_1.html) ![](plot/volcanoplot/volcano_interactive_ratio_stats_1.png) ### `r gsub("_RepBio_\\d+", "", names(list_ratio_stats)[3])` For an interactive version : [Click here](plot/volcanoplot/volcano_interactive_ratio_stats_3.html) ![](plot/volcanoplot/volcano_interactive_ratio_stats_3.png) ### `r gsub("_RepBio_\\d+", "", names(list_ratio_stats)[4])` For an interactive version : [Click here](plot/volcanoplot/volcano_interactive_ratio_stats_4.html) ![](plot/volcanoplot/volcano_interactive_ratio_stats_4.png) ### `r gsub("_RepBio_\\d+", "", names(list_ratio_stats)[5])` For an interactive version : [Click here](plot/volcanoplot/volcano_interactive_ratio_stats_5.html) ![](plot/volcanoplot/volcano_interactive_ratio_stats_5.png) ::: ## Simplified differential analysis ### Venn diagram ```{r, eval = F} # import data load(file = here::here("data/microarray/output/ratio_stats_microarray_leaf_PeaSulf.RData")) lfc_lim_i <- 0 # create function to have list of genes get_genes<- function(df_comparison_i, lfc_lim = lfc_lim_i, padj_lim= 0.05, sign_i){ # change value here to change parameter 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") lt_up <- list( `SS_WT1 vs SS_W78*` = get_genes(list_ratio_stats$`Mut_SS_W78*_RepBio_4 vs WT_SS_WT1_RepBio_4`, sign_i = "Up",lfc_lim = lfc_lim_i), `SS_WT2 vs SS_E568K` = get_genes(list_ratio_stats$`Mut_SS_E568K_RepBio_3 vs WT_SS_WT2_RepBio_3`, sign_i = "Up",lfc_lim = lfc_lim_i), `SD_WT1 vs SD_W78*` = get_genes(list_ratio_stats$`Mut_SD_W78*_RepBio_2 vs WT_SD_WT1_RepBio_2`, sign_i = "Up",lfc_lim = lfc_lim_i), `SD_WT2 vs SD_E568K` = get_genes(list_ratio_stats$`Mut_SD_E568K_RepBio_1 vs WT_SD_WT2_RepBio_1`, sign_i = "Up",lfc_lim = lfc_lim_i), `SS_WT vs SD_WT` = get_genes(list_ratio_stats$`WT_SD_RepBio_5 vs WT_SS_RepBio_5`, sign_i = "Up") ) lt_down <- list( `SS_WT1 vs SS_W78*` = get_genes(list_ratio_stats$`Mut_SS_W78*_RepBio_4 vs WT_SS_WT1_RepBio_4`, sign_i = "Down",lfc_lim = lfc_lim_i), `SS_WT2 vs SS_E568K` = get_genes(list_ratio_stats$`Mut_SS_E568K_RepBio_3 vs WT_SS_WT2_RepBio_3`, sign_i = "Down",lfc_lim = lfc_lim_i), `SD_WT1 vs SD_W78*` = get_genes(list_ratio_stats$`Mut_SD_W78*_RepBio_2 vs WT_SD_WT1_RepBio_2`, sign_i = "Down",lfc_lim = lfc_lim_i), `SD_WT2 vs SD_E568K` = get_genes(list_ratio_stats$`Mut_SD_E568K_RepBio_1 vs WT_SD_WT2_RepBio_1`, sign_i = "Down",lfc_lim = lfc_lim_i), `SS_WT vs SD_WT` = get_genes(list_ratio_stats$`WT_SD_RepBio_5 vs WT_SS_RepBio_5`, sign_i = "Down") ) length(lt_down$`SS_WT vs SD_WT`) lt_up_simplified <- list( `Mutant vs WT under SS` = intersect(lt_up$`SS_WT1 vs SS_W78*`, lt_up$`SS_WT2 vs SS_E568K`), `Mutant vs WT under SD` = intersect(lt_up$`SD_WT1 vs SD_W78*`, lt_up$`SD_WT2 vs SD_E568K`), `SS_WT vs SD_WT` = lt_up$`SS_WT vs SD_WT` ) lt_down_simplified <- list( `Mutant vs WT under SS` = intersect(lt_down$`SS_WT1 vs SS_W78*`, lt_down$`SS_WT2 vs SS_E568K`), `Mutant vs WT under SD` = intersect(lt_down$`SD_WT1 vs SD_W78*`, lt_down$`SD_WT2 vs SD_E568K`), `SS_WT vs SD_WT` = lt_down$`SS_WT vs SD_WT` ) p_down <- plot(eulerr::euler(lt_down_simplified, shape = "ellipse"), quantities = TRUE, fills = c(sulfate_pallet[1], sulfate_pallet[2], "#9381FF"), legend = list(side = "right"), main = list(label = paste0("Euler diagram of down regulated genes with LFC ", lfc_lim_i), cex = 1.3)) p_up <- plot(eulerr::euler(lt_up_simplified, shape = "ellipse"), quantities = TRUE, fills = c(sulfate_pallet[1], sulfate_pallet[2], "#9381FF"), legend = list(side = "right"), main = list(label = paste0("Euler diagram of up regulated genes with LFC ", lfc_lim_i), cex = 1.3)) fig_export(path = here::here(paste0("report/microarray/plot/venn_diagram/euler_simplified_down_lfc_",lfc_lim_i)), p_down, height = 21/3, width = 29.7/3, res_i = 600) fig_export(path = here::here(paste0("report/microarray/plot/venn_diagram/euler_simplified_up_lfc_",lfc_lim_i)), p_up, height = 21/3, width = 29.7/3, res_i = 600) # diagram venn fills <- c( as.character(sulfate_pallet[1]), as.character(sulfate_pallet[2]),"SS_WT vs SD_WT" = "#9381FF" ) names(lt_down_simplified) <- c("Mutant vs WT (SS) ", " Mutant vs WT (SD)", "SS_WT vs SD_WT") names(lt_up_simplified) <- c("Mutant vs WT (SS) ", " Mutant vs WT (SD)", "SS_WT vs SD_WT") ggvenn_plot_down = ggvenn::ggvenn( lt_down_simplified, fill_color = unname(fills), fill_alpha = 0.55, stroke_size = 0.4 ) + ggtitle(paste0( "Down-regulated genes") ) + theme(plot.title = element_text(face = "bold", size = 15), legend.position = "none") ggvenn_plot_up = ggvenn::ggvenn( lt_up_simplified, fill_color = unname(fills), fill_alpha = 0.55, stroke_size = 0.4 ) + ggtitle(paste0( "Up-regulated genes") ) + theme(plot.title = element_text(face = "bold", size = 15), legend.position = "none") print(ggvenn_plot_down) ; print(ggvenn_plot_up) fig_export(path = here::here(paste0("report/microarray/plot/venn_diagram/venn_simplified_up_lfc_",lfc_lim_i)), ggvenn_plot_up, height = 21/3, width = 29.7/3, res_i = 600) fig_export(path = here::here(paste0("report/microarray/plot/venn_diagram/venn_simplified_down_lfc_",lfc_lim_i)), ggvenn_plot_down, height = 21/3, width = 29.7/3, res_i = 600) # save plot for article venn_list <- list( up = ggvenn_plot_up, down = ggvenn_plot_down ) save(venn_list,file = here::here(paste0("report/microarray/plot/venn_diagram/venn_simplified_lfc_",lfc_lim_i,".RData"))) # export results #### save(lt_up_simplified, lt_down_simplified,file = here::here(paste0("data/microarray/output/venn_result_simplified_condition_sign_lfc_",lfc_lim_i,".RData"))) ``` ### GO enrichment ```{r, eval = F} #parameter lfc_lim_i = 0 # collect data and add Psat to PsCAM load(here::here("data/microarray/output/raw_data_microarray_leaf_PeaSulf.RData")) load(file = here::here(paste0("data/microarray/output/venn_result_simplified_condition_sign_lfc_",lfc_lim_i,".RData"))) df_info_gene <- RG$genes %>% dplyr::rename(PsCam = ID) %>% left_join(.,read_excel(here::here("data/microarray/Resultats_4plex-POIS-2014_01_230415_modKG.xlsx"), sheet = "Pscam_to_Psat_v1c_mrna_besthit-", col_names = T), by = "PsCam") cluster_info <- rbind(purrr::imap_dfr(lt_down_simplified, ~ tibble(PsCam = .x, comparison = .y)) %>% mutate(cluster = paste0("(Down) ", comparison)), purrr::imap_dfr(lt_up_simplified, ~ tibble(PsCam = .x, comparison = .y)) %>% mutate(cluster = paste0("(Up) ", comparison)) ) %>% full_join(., df_info_gene, by = "PsCam") %>% mutate(ID = Psat) %>% drop_na(comparison, ID) df_GO = GO_on_different_group(functional_roles = "BP", #"BP" group_info = cluster_info, group = "cluster", ID = "ID", top = 10 ) # Merge this with your original data to fill missing combinations Cluster_selected_GO_filled <-df_GO %>% # on s'arrete ici pour la fonction. dplyr::rename(cluster = group) %>% mutate(cluster = factor(cluster, levels = unique(cluster)), comparison = str_trim(str_remove(cluster, "\\s*\$.*?\$")), sign = str_extract(cluster, "(?<=\$).*?(?=\$)") ) %>% mutate(comparison = forcats::fct_relevel(comparison, "Mutant vs WT under SS", "Mutant vs WT under SD", "SS_WT vs SD_WT")) #vector_color <- c("#067B5B", "#6AD6AD", "#B87100", "#FFAC5C", "#9381FF", "#067B5B", "#6AD6AD", "#B87100", "#FFAC5C", "#9381FF") # Plot the enrichment by GO px <- ggplot(Cluster_selected_GO_filled, aes(x = cluster, y = Description_GO, fill = `|-log10(Pval)|`)) + geom_tile(color = "black") + scale_fill_gradient2(low = "white", high = "red", na.value = "white", limits = c(0, NA)) + theme_minimal() + theme(axis.text = element_text(size = 8, colour = "black"), axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1), panel.grid.major = element_blank(), plot.background = element_rect(fill = "white", colour = "white")) + labs(fill = "-log10(FDR)", x= "Biological process", y = "Cluster", title = "GO terms for all gene deregulated") + #scale_size_manual(values = c(dot = 2, no_dot = NA), guide = "none")+ new_scale_fill() + scale_fill_manual( values = c("Mutant vs WT under SS" = as.character(sulfate_pallet[1]), "Mutant vs WT under SD" = as.character(sulfate_pallet[2]), "SS_WT vs SD_WT" = "#9381FF"), name = "Comparison" ) + geom_tile( aes( x = cluster, y = -1, # If you truly want it at a negative y-value, # ensure your y-scale is continuous or can handle this fill = comparison, width = 0.90, height = 0.80 ), data = Cluster_selected_GO_filled, color = "black", alpha = 1, inherit.aes = FALSE )+ new_scale_fill() + scale_fill_manual( values = c('Down' = "#1d4877", 'Up' = "#ee3e32"), name = "Sign" ) + geom_tile( aes( x = cluster, y = -0.05, # If you truly want it at a negative y-value, # ensure your y-scale is continuous or can handle this fill = sign, width = 0.90, height = 0.80 ), data = Cluster_selected_GO_filled, color = "black", alpha = 1, inherit.aes = FALSE ) ; px # export fig_export(path =paste0("report/microarray/plot/GO/GO_BP_all_simplified_absolute_lfc_",lfc_lim_i), plot_x = px, height_i = 13, width_i = 9, res = 600) # export for article load(file = here::here(paste0("report/microarray/plot/venn_diagram/venn_simplified_lfc_",lfc_lim_i,".RData"))) # load the list of the two venne diagram venn_list$up <- venn_list$up + labs(tag = "A") venn_list$down <- venn_list$down + labs(tag = "B") px <- px + labs(tag = "C") venn_stack <- venn_list$up / venn_list$down combined_plot <- (venn_stack | px) + # | = côte-à-côte plot_layout(widths = c(2, 1), guides = "collect") & theme(plot.margin = margin(5, 5, 5, 5))+ theme(plot.margin = margin(5, 5, 5, 5), plot.tag = element_text(size = 14)) # ───────────────────────────────────────────────────────── # 3) APERÇU RAPIDE # ───────────────────────────────────────────────────────── print(combined_plot) fig_export(path =paste0("report/microarray/plot/fig_lfc_",lfc_lim_i), plot_x = combined_plot, height_i = 13, width_i = 15, res = 600, format = "png") ``` ### Results :::: panel-tabset ### With log fold change of 0 ::: {layout-ncol="2"} ![](plot/venn_diagram/euler_simplified_down_lfc_0.png) ![](plot/venn_diagram/venn_simplified_down_lfc_0.png) ![](plot/venn_diagram/euler_simplified_up_lfc_0.png) ![](plot/venn_diagram/venn_simplified_up_lfc_0.png) ::: ![](plot/GO/GO_BP_all_simplified_absolute_lfc_0.png) ### With log fold change of 1 ::: {layout-ncol="2"} ![](plot/venn_diagram/euler_simplified_down_lfc_1.png) ![](plot/venn_diagram/euler_simplified_up_lfc_1.png) ::: ![](plot/GO/GO_BP_all_simplified_absolute_lfc_1.png) :::: ## Differentiation between the two genotypes ### UpsetPlot ```{r, eval = F} # import data load(file = here::here("data/microarray/output/ratio_stats_microarray_leaf_PeaSulf.RData")) lfc_lim_i <- 1 # create function to have list of genes get_genes<- function(df_comparison_i, lfc_lim = lfc_lim_i, padj_lim= 0.05, sign_i){ # change value here to change parameter 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( `SS_WT1 vs SS_W78*` = get_genes(list_ratio_stats$`Mut_SS_W78*_RepBio_4 vs WT_SS_WT1_RepBio_4`, sign_i = "Up"), `SS_WT2 vs SS_E568K` = get_genes(list_ratio_stats$`Mut_SS_E568K_RepBio_3 vs WT_SS_WT2_RepBio_3`, sign_i = "Up"), `SD_WT1 vs SD_W78*` = get_genes(list_ratio_stats$`Mut_SD_W78*_RepBio_2 vs WT_SD_WT1_RepBio_2`, sign_i = "Up"), `SD_WT2 vs SD_E568K` = get_genes(list_ratio_stats$`Mut_SD_E568K_RepBio_1 vs WT_SD_WT2_RepBio_1`, sign_i = "Up"), `SS_WT vs SD_WT` = get_genes(list_ratio_stats$`WT_SD_RepBio_5 vs WT_SS_RepBio_5`, sign_i = "Up") ) m_up <- ComplexHeatmap::make_comb_mat(lt_up, mode = "distinct") # Down-regulated genes lt_down <- list( `SS_WT1 vs SS_W78*` = get_genes(list_ratio_stats$`Mut_SS_W78*_RepBio_4 vs WT_SS_WT1_RepBio_4`, sign_i = "Down"), `SS_WT2 vs SS_E568K` = get_genes(list_ratio_stats$`Mut_SS_E568K_RepBio_3 vs WT_SS_WT2_RepBio_3`, sign_i = "Down"), `SD_WT1 vs SD_W78*` = get_genes(list_ratio_stats$`Mut_SD_W78*_RepBio_2 vs WT_SD_WT1_RepBio_2`, sign_i = "Down"), `SD_WT2 vs SD_E568K` = get_genes(list_ratio_stats$`Mut_SD_E568K_RepBio_1 vs WT_SD_WT2_RepBio_1`, sign_i = "Down"), `SS_WT vs SD_WT` = get_genes(list_ratio_stats$`WT_SD_RepBio_5 vs WT_SS_RepBio_5`, sign_i = "Down") ) m_down <- ComplexHeatmap::make_comb_mat(lt_down, mode = "distinct") # # create_upset_plot_rnaseq(m_up, selected_columns = c("100000", "010000", "001000","000100","000010","000001","100100","010010","001001","101101","011011","101101","111111"), # "Up", organ_i) # # create_upset_plot_rnaseq(m_down, selected_columns = c("100000", "010000", "001000","000100","000010","000001","100100","010010","001001","101101","011011","101101","111111"), # "Down", organ_i) 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 <-upsetplot_condition_merge_sign(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(paste0("report/microarray/plot/upsetplot/upset_up_down_condition_lfc_",lfc_lim_i)), upsetplot_condition, height = 4, width = 14.5, res_i = 600) set.seed(123) p_down <- plot(eulerr::euler(lt_down, shape = "ellipse"), quantities = TRUE, fills = c(lighten(as.character(sulfate_pallet[1]), amount = -.2), lighten(as.character(sulfate_pallet[1]), amount = .5), lighten(as.character(sulfate_pallet[2]), amount = -.2), lighten(as.character(sulfate_pallet[2]), amount = .3), "#9381FF"), 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(lighten(as.character(sulfate_pallet[1]), amount = -.2), lighten(as.character(sulfate_pallet[1]), amount = .5), lighten(as.character(sulfate_pallet[2]), amount = -.2), lighten(as.character(sulfate_pallet[2]), amount = .3), "#9381FF"), legend = list(side = "right"), main = list(label = "Euler diagram of upregulated genes", cex = 1.3)) fig_export(path = here::here(paste0("report/microarray/plot/venn_diagram/euler_down_lfc_",lfc_lim_i)), p_down, height = 21/3, width = 29.7/3, res_i = 600) fig_export(path = here::here(paste0("report/microarray/plot/venn_diagram/euler_up_lfc_",lfc_lim_i)), p_up, height = 21/3, width = 29.7/3, res_i = 600) # verification if i had the same results as in excel. length(intersect(lt_up$`SD_WT1 vs SD_W78*`, lt_up$`SD_WT2 vs SD_E568K`)) length(intersect(lt_down$`SD_WT1 vs SD_W78*`, lt_down$`SD_WT2 vs SD_E568K`)) length(intersect(lt_up$`SS_WT1 vs SS_W78*`, lt_up$`SS_WT2 vs SS_E568K`)) length(intersect(lt_down$`SS_WT1 vs SS_W78*`, lt_down$`SS_WT2 vs SS_E568K`)) # Export results for Gene enrichment #### save(lt_up, lt_down, m_up, m_down,file = here::here(paste0("data/microarray/output/upset_result_condition_sign_lfc_",lfc_lim_i,".RData"))) ``` :::: panel-tabset ### With log fold change of 0 ![](plot/upsetplot/upset_up_down_condition_lfc_0.png) ::: {layout-ncol="2"} ![](plot/venn_diagram/euler_down_lfc_0.png) ![](plot/venn_diagram/euler_up_lfc_0.png) ::: **Which group does the psult4 gene belong to?** ```{r} load(file = here::here("data/microarray/output/upset_result_condition_sign_lfc_0.RData")) my_gene_of_interest <- "PsCam042688" in_up <- names(keep(lt_up, ~ my_gene_of_interest %in% .x)) in_down <- names(keep(lt_down, ~ my_gene_of_interest %in% .x)) df_up <- data.frame("In Up Regulated" = in_up, check.names = FALSE) df_down <- data.frame("In Down Regulated" = in_down, check.names = FALSE) # Generate HTML for each table with kable table1_html <- kable(df_up, format = "html", table.attr = "class='table'") %>% kable_styling(full_width = FALSE) table2_html <- kable(df_down, format = "html", table.attr = "class='table'") %>% kable_styling(full_width = FALSE) # Combine the two panels side by side using a flexbox divider html_output <- paste0( "<div style='display: flex; justify-content: center; gap: 20px;'>", "<div>", table2_html, "</div>", "<div>", table1_html, "</div>", "</div>" ) # Display HTML without escapement knitr::asis_output(html_output) ``` ### With log fold change of 1 ![](plot/upsetplot/upset_up_down_condition_lfc_1.png) ::: {layout-ncol="2"} ![](plot/venn_diagram/euler_down_lfc_1.png) ![](plot/venn_diagram/euler_up_lfc_1.png) ::: **Which group does the psult4 gene belong to?** ```{r} load(file = here::here("data/microarray/output/upset_result_condition_sign_lfc_1.RData")) my_gene_of_interest <- "PsCam042688" in_up <- names(keep(lt_up, ~ my_gene_of_interest %in% .x)) in_down <- names(keep(lt_down, ~ my_gene_of_interest %in% .x)) df_up <- data.frame("In Up Regulated" = in_up, check.names = FALSE) df_down <- data.frame("In Down Regulated" = in_down, check.names = FALSE) # Generate HTML for each table with kable table1_html <- kable(df_up, format = "html", table.attr = "class='table'") %>% kable_styling(full_width = FALSE) table2_html <- kable(df_down, format = "html", table.attr = "class='table'") %>% kable_styling(full_width = FALSE) # Combine the two panels side by side using a flexbox divider html_output <- paste0( "<div style='display: flex; justify-content: center; gap: 20px;'>", "<div>", table2_html, "</div>", "<div>", table1_html, "</div>", "</div>" ) # Display HTML without escapement knitr::asis_output(html_output) ``` :::: ### GO enrichment ```{r, eval = F} #parameter lfc_lim_i = 1 # collect data and add Psat to PsCAM load(here::here("data/microarray/output/raw_data_microarray_leaf_PeaSulf.RData")) load(file = here::here(paste0("data/microarray/output/upset_result_condition_sign_lfc_",lfc_lim_i,".RData"))) df_info_gene <- RG$genes %>% dplyr::rename(PsCam = ID) %>% left_join(.,read_excel(here::here("data/microarray/Resultats_4plex-POIS-2014_01_230415_modKG.xlsx"), sheet = "Pscam_to_Psat_v1c_mrna_besthit-", col_names = T), by = "PsCam") cluster_info <- rbind(purrr::imap_dfr(lt_down, ~ tibble(PsCam = .x, comparison = .y)) %>% mutate(cluster = paste0("(Down) ", comparison)), purrr::imap_dfr(lt_up, ~ tibble(PsCam = .x, comparison = .y)) %>% mutate(cluster = paste0("(Up) ", comparison)) ) %>% full_join(., df_info_gene, by = "PsCam") %>% mutate(ID = Psat) %>% drop_na(comparison, ID) df_GO = GO_on_different_group(functional_roles = "BP", group_info = cluster_info, group = "cluster", ID = "ID", top = 10 ) # Merge this with your original data to fill missing combinations Cluster_selected_GO_filled <-df_GO %>% # on s'arrete ici pour la fonction. dplyr::rename(cluster = group) %>% mutate(cluster = factor(cluster, levels = unique(cluster)), comparison = str_trim(str_remove(cluster, "\\s*\$.*?\$")), sign = str_extract(cluster, "(?<=\$).*?(?=\$)") ) %>% mutate(comparison = forcats::fct_relevel(comparison, "SS_WT1 vs SS_W78*", "SS_WT2 vs SS_E568K", "SD_WT1 vs SD_W78*", "SD_WT2 vs SD_E568K", "SS_WT vs SD_WT")) #vector_color <- c("#067B5B", "#6AD6AD", "#B87100", "#FFAC5C", "#9381FF", "#067B5B", "#6AD6AD", "#B87100", "#FFAC5C", "#9381FF") # Plot the enrichment by GO px <- ggplot(Cluster_selected_GO_filled, aes(x = cluster, y = Description_GO, fill = `|-log10(Pval)|`)) + geom_tile(color = "black") + scale_fill_gradient2(low = "white", high = "red", na.value = "white", limits = c(0, NA)) + theme_minimal() + theme(axis.text = element_text(size = 8, colour = "black"), axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1), panel.grid.major = element_blank(), plot.background = element_rect(fill = "white", colour = "white")) + labs(fill = "-log10(FDR)", x= "Biological process", y = "Cluster", title = "GO terme for all gene deregulated") + #scale_size_manual(values = c(dot = 2, no_dot = NA), guide = "none")+ new_scale_fill() + scale_fill_manual( values = c(lighten(as.character(sulfate_pallet[1]), amount = -.2), lighten(as.character(sulfate_pallet[1]), amount = .5), lighten(as.character(sulfate_pallet[2]), amount = -.2), lighten(as.character(sulfate_pallet[2]), amount = .3), "#9381FF" ), name = "Comparison" ) + geom_tile( aes( x = cluster, y = -1, # If you truly want it at a negative y-value, # ensure your y-scale is continuous or can handle this fill = comparison, width = 0.90, height = 0.80 ), data = Cluster_selected_GO_filled, color = "black", alpha = 1, inherit.aes = FALSE )+ new_scale_fill() + scale_fill_manual( values = c('Down' = "#1d4877", 'Up' = "#ee3e32"), name = "Sign" ) + geom_tile( aes( x = cluster, y = -0.05, # If you truly want it at a negative y-value, # ensure your y-scale is continuous or can handle this fill = sign, width = 0.90, height = 0.80 ), data = Cluster_selected_GO_filled, color = "black", alpha = 1, inherit.aes = FALSE ) ; px # export fig_export(path =paste0("report/microarray/plot/GO/GO_BP_all_absolute_lfc_",lfc_lim_i), plot_x = px, height_i = 13, width_i = 9, res = 600) ``` ```{r, eval = F} # euler (or upset) #parameter lfc_lim_i = 1 # collect data and add Psat to PsCAM load(here::here("data/microarray/output/raw_data_microarray_leaf_PeaSulf.RData")) load(file = here::here(paste0("data/microarray/output/upset_result_condition_sign_lfc_",lfc_lim_i,".RData"))) df_info_gene <- RG$genes %>% dplyr::rename(PsCam = ID) %>% left_join(.,read_excel(here::here("data/microarray/Resultats_4plex-POIS-2014_01_230415_modKG.xlsx"), sheet = "Pscam_to_Psat_v1c_mrna_besthit-", col_names = T), by = "PsCam") 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") ) gene_id_col <- colnames(df_DEG_condition)[1] v_conditions <- colnames(df_DEG_condition)[2:(ncol(df_DEG_condition) - 1)] df_gene_intersections <- df_DEG_condition %>% rowwise() %>% filter(Sign =="Down") %>% mutate(intersection = paste0(c_across(all_of(v_conditions)), collapse = "")) %>% ungroup() intersection_list <- list( I1 = c('SS_WT vs SD_WT', 'SS_WT1 vs SS_W78*', 'SS_WT2 vs SS_E568K', 'SD_WT1 vs SD_W78*', 'SD_WT2 vs SD_E568K'), I2 = c('SS_WT1 vs SS_W78*', 'SS_WT2 vs SS_E568K', 'SD_WT1 vs SD_W78*', 'SD_WT2 vs SD_E568K'), I3 = c('SS_WT1 vs SS_W78*', 'SS_WT2 vs SS_E568K', 'SS_WT vs SD_WT'), I4 = c('SD_WT1 vs SD_W78*', 'SD_WT2 vs SD_E568K', 'SS_WT vs SD_WT'), I5 = c('SS_WT1 vs SS_W78*', 'SS_WT2 vs SS_E568K'), I6 = c('SD_WT1 vs SD_W78*', 'SD_WT2 vs SD_E568K'), I7 = c('SS_WT1 vs SS_W78*', 'SD_WT1 vs SD_W78*'), I8 = c('SS_WT2 vs SS_E568K', 'SD_WT2 vs SD_E568K'), I9 = c('SS_WT vs SD_WT', 'SS_WT1 vs SS_W78*'), I10 = c('SS_WT vs SD_WT', 'SS_WT2 vs SS_E568K'), I11 = c('SS_WT vs SD_WT', 'SD_WT1 vs SD_W78*'), I12 = c('SS_WT vs SD_WT', 'SD_WT2 vs SD_E568K'), I13 = 'SS_WT1 vs SS_W78*', I14 = 'SS_WT2 vs SS_E568K', I15 = 'SD_WT1 vs SD_W78*', I16 = 'SD_WT2 vs SD_E568K', I17 = 'SS_WT vs SD_WT' ) intersection_patterns <- map(intersection_list, function(conds) { sapply(v_conditions, function(x) as.integer(x %in% conds)) %>% paste0(collapse = "") }) signs <- c("Down", "Up") df_gene_intersections_by_sign <- map(signs, function(s) { df_DEG_condition %>% filter(Sign == s) %>% rowwise() %>% mutate(intersection = paste0(c_across(all_of(v_conditions)), collapse = "")) %>% ungroup() }) names(df_gene_intersections_by_sign) <- signs intersection_genes_all <- list() for(i in names(intersection_list)) { for(s in signs) { pattern <- intersection_patterns[[i]] genes <- df_gene_intersections_by_sign[[s]] %>% filter(intersection == pattern) %>% pull(!!sym(gene_id_col)) # Nommer l'élément comme "I1 (Up)" ou "I1 (Down)", etc. intersection_genes_all[[paste0("(", s, ") ", i)]] <- genes } } # and then same hase before cluster_info <- purrr::imap_dfr(intersection_genes_all, ~ tibble(PsCam = .x, intersection = .y)) %>% mutate(cluster = intersection) %>% full_join(., df_info_gene, by = "PsCam") %>% mutate(ID = Psat) %>% drop_na(intersection, ID) df_GO = GO_on_different_group(functional_roles = "BP", group_info = cluster_info, group = "cluster", ID = "ID", top = 10 ) # Merge this with your original data to fill missing combinations Cluster_selected_GO_filled <- df_GO %>% dplyr::rename(cluster = group) %>% mutate(cluster = factor(cluster, levels = unique(cluster)), comparison = str_trim(str_remove(cluster, "\\s*\$.*?\$")), sign = str_extract(cluster, "(?<=\$).*?(?=\$)") ) #vector_color <- c("#067B5B", "#6AD6AD", "#B87100", "#FFAC5C", "#9381FF", "#067B5B", "#6AD6AD", "#B87100", "#FFAC5C", "#9381FF") # Plot the enrichment by GO px <- ggplot(Cluster_selected_GO_filled, aes(x = cluster, y = Description_GO, fill = `|-log10(Pval)|`)) + geom_tile(color = "black") + scale_fill_gradient2(low = "white", high = "red", na.value = "white", limits = c(0, NA)) + theme_minimal() + theme(axis.text = element_text(size = 8, colour = "black"), axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1), panel.grid.major = element_blank(), plot.background = element_rect(fill = "white", colour = "white")) + labs(fill = "-log10(FDR)", x= "Biological process", y = "Cluster", title = "GO terme for all gene deregulated") + #scale_size_manual(values = c(dot = 2, no_dot = NA), guide = "none")+ new_scale_fill() + scale_fill_manual( values = c('Down' = "#1d4877", 'Up' = "#ee3e32"), name = "Sign" ) + geom_tile( aes( x = cluster, y = -0.05, # If you truly want it at a negative y-value, # ensure your y-scale is continuous or can handle this fill = sign, width = 0.90, height = 0.80 ), data = Cluster_selected_GO_filled, color = "black", alpha = 1, inherit.aes = FALSE ) ; px # export fig_export(path = paste0("report/microarray/plot/GO/GO_BP_euler_lfc_", lfc_lim_i), plot_x = px, height_i = 20, width_i = 12, res = 600) ``` :::: panel-tabset #### GO enrichment on deregulated genes ::: panel-tabset ##### With log fold change of 0 ![](plot/GO/GO_BP_all_absolute_lfc_0.png) ##### With log fold change of 1 ![](plot/GO/GO_BP_all_absolute_lfc_1.png) ::: #### On Euler diagram (or upset plot) ::: panel-tabset ##### With log fold change of 0 ![](plot/GO/GO_BP_euler_lfc_0.png) ##### With log fold change of 1 ![](plot/GO/GO_BP_euler_lfc_1.png) ::: ::::