13 Exploratory analysis of Raw Data from microarray

13.1 Material and methods

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

# install.packages(here::here("data/microarray/org.Psativum1c.eg.db"), repos = NULL, type = "source")

#pkg
library(readxl)
library(affy) #BiocManager::install("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(ComplexHeatmap)
library(knitr)
library(kableExtra)
library(ggtext)
library(ggnewscale)
library(org.Psativum1c.eg.db) # if not working install it
library(clusterProfiler)
library(readxl)
library(GO.db)
library(AnnotationDbi) # BiocManager::install("AnnotationDbi")
library(ggnewscale) # to have two scale_fill

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

13.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. And average data and stats corresponding to one of the paired student tests, then corrected between the four biological repeats. This data has been previously analized here.)

Import raw data (green and red)

Code

# Function that import all data (take a lot of time)
source (here::here("src/function/microarray/import_raw_data_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_raw_data_microarray(path_i, path_info_range_i, sheet_i, begin_end_line)
write_csv(execut[["global_df"]], here::here("data/microarray/output/raw_data_microarray_leaf_PeaSulf.csv"))
write_csv(execut[["df_info_sample_clean_compile"]], here::here("data/microarray/output/list_microarray_microarray_leaf_PeaSulf.csv"))

On microarray chips, it’s quite common to find probes that do not belong to the species under study. These probes are typically placed there for various reasons, such as serving as exogenous controls (often called alien probes, spike-ins, or external controls) or to check hybridization specificity. Here is from human or bacteria. In this step i will import it.

Code

path = here::here("data/microarray/Resultats_4plex-POIS-2014_01_230415_modKG.xlsx")
path_info_range = here::here("data/microarray/info_range_Resultats_4plex-POIS-2014_01_230415_modKG.xlsx")
begin_end_line <- c("41123", "41159")
sheet_i = "Complete"

external_controls = import_raw_data_microarray(path,path_info_range, sheet_i, begin_end_line)
write_csv(external_controls[["global_df"]], here::here("data/microarray/output/external_controls_data_microarray_leaf_PeaSulf.csv"))

Code

external_controls = read_csv(here::here("data/microarray/output/external_controls_data_microarray_leaf_PeaSulf.csv"), show_col_types = FALSE)

df_mean = external_controls %>%  dplyr::summarise(
    global_mean = mean(unlist(.[, 9:ncol(.)]), na.rm = TRUE),
     global_max = max(unlist(.[, 9:ncol(.)]), na.rm = TRUE)
  )

2^df_mean

Code

# import raw data red green
# data1 = read_csv(here::here("data/microarray/output/raw_data_microarray_leaf_PeaSulf.csv"), show_col_types = FALSE) %>% 
#   as.data.frame() %>% 
#   dplyr::select(c("id_probe", "Mut_SD_5.1_Rep1_Red_1", "WT_SD_1.1_Rep1_Green_1"))
# 
# p_non_norm = data1 %>% 
#   pivot_longer(-"id_probe") %>% 
#   ggplot(aes(x = value, col = name)) +
#     geom_density() +
#     labs(
#       title = "Density of values",
#       x = "Values",
#       y = "Density"
#     ) +
#   theme_minimal()
# 
# data_longer_test = data1 %>% column_to_rownames("id_probe")
# 
# norm_data <- normalizeBetweenArrays(data_longer_test, method = "cyclicloess")
# 
# p_norm = norm_data %>% 
#   as.data.frame() %>% 
#   rownames_to_column("id_probe") %>% 
#  pivot_longer(-"id_probe") %>% 
#   ggplot(aes(x = value, col = name)) +
#     geom_density() +
#     labs(
#       title = "Density of normalized values",
#       x = "Values",
#       y = "Density"
#     ) +
#   theme_minimal()
# p_non_norm / p_norm

Creation of the RGlist file (with absolute value (befor log2))

Code

# Import raw data red green 
df_raw_RG_log2 = read_csv(here::here("data/microarray/output/raw_data_microarray_leaf_PeaSulf.csv"), show_col_types = FALSE) %>% 
  as.data.frame()

# Creation of the RG list for background stemps ####
# Converted into input data (exp of value with log2)
df_raw_RG<- df_raw_RG_log2
cols_to_transform <- grep("Red|Green", names(df_raw_RG))
df_raw_RG[cols_to_transform] <- lapply(df_raw_RG[cols_to_transform], function(col) 2^col)

# Verification before creation of the RGList
red_cols   <- grep("Red_\\d+$",   colnames(df_raw_RG), value = TRUE)
green_cols <- grep("Green_\\d+$", colnames(df_raw_RG), value = TRUE)

extract_index <- function(x, color = c("Red","Green")) {
  pattern <- paste0(".*", color, "_(\\d+)$") 
  sub(pattern, "\\1", x)
}

numbers_red   <- extract_index(red_cols,   "Red")
numbers_green <- extract_index(green_cols, "Green")

common_indices <- intersect(numbers_red, numbers_green)
cat("Same number find :", common_indices, "\n")

# Parameter of the array
n_spots  <- nrow(df_raw_RG)
n_arrays <- length(common_indices) # 20 = 4 rep * 5 comparison

R_matrix <- matrix(NA, nrow = n_spots, ncol = n_arrays) %>% as.data.frame()
G_matrix <- matrix(NA, nrow = n_spots, ncol = n_arrays) %>% as.data.frame()

k <- 1
for (i in common_indices) {
  
  # Trouver la/les colonnes Red_... qui finissent par Red_i
  # (en théorie, une seule par index i, mais grep renvoie un vecteur - on prend le premier)
  col_red   <- grep(paste0("Red_", i, "$"),   red_cols,   value = TRUE)
  col_green <- grep(paste0("Green_", i, "$"), green_cols, value = TRUE)
  
  # Remplir la matrice
  R_matrix[, k] <- df_raw_RG[[col_red]]
  G_matrix[, k] <- df_raw_RG[[col_green]]
  
  # Nommer la colonne pour s’y retrouver (ex : "i" ou le nom exact "Sample1")
  # Ici on va prendre "Red_xxxx" comme titre, ou "i" tout simplement
  colnames(R_matrix)[k] <- col_red
  colnames(G_matrix)[k] <- col_green
  
  k <- k + 1
}

# Creation of the RG list end ####

RG <- new("RGList")
RG$R <- R_matrix
RG$G <- G_matrix
RG$genes <- data.frame(ID = df_raw_RG$id_probe)

MA_raw <- normalizeWithinArrays(RG, method = "none") # without normalisation

# export RG
save(RG,MA_raw, file =here::here("data/microarray/output/raw_data_microarray_leaf_PeaSulf.RData"))

13.3 Standardization of the data cross-chip data

Bonus: If i apply double normalization ((loess and quartil) not necessary for this experiment))

In a two-channel microarray, we often apply within-array normalization (such as LOESS) to correct for individual slide biases (systematic variations linked to position on the chip, channel, etc.). Next, a so-called between-array normalization (such as quantile normalization) can be applied to make distributions between replicates comparable. Normalization functions are available in R/Bioconductor (limma package, etc.).

See also this paper Bolstad et al. (2003)

Code

# Load data
load(here::here("data/microarray/output/raw_data_microarray_leaf_PeaSulf.RData"))

# Raw data visualization ####
source (here::here("src/function/microarray/plot_rg_density.R")) # Load function for see density of the differentent normalization

p_raw <- plot_rg_density(RG)+labs(title = "plotDensities without normalization") # execut function

# Remove or adjust background values (noise) to obtain intensities closer to the true spot intensity ####
# Helps model and remove technical noise at lower intensities.
# The normexp method is commonly recommended because it avoids negative or zero intensities by adding an offset. In input is absolute value

RG_background <- backgroundCorrect(RG, method = "normexp", offset = 20) # 20 test multiple value. Perfect for delet noise 
rownames(RG_background$R) <- RG_background[["genes"]][["ID"]]
rownames(RG_background$G) <- RG_background[["genes"]][["ID"]]

# visualisation 
p_background <- plot_rg_density(RG_background)+labs(title = "plotDensities without normalization but with background correction")

log2(min(RG_background$G)) # need to be positive

MA_raw <- normalizeWithinArrays(RG, method = "none") # without normalisation
plotMA(MA_raw, array = 1) # interesting to se this plot

MA_background <- normalizeWithinArrays(RG_background, method = "none") # with normalisation
plotMA(MA_background, array = 1)

# Normalization loess ####
MA.p <-normalizeWithinArrays(RG_background, method = "loess")
p_loess <- plot_rg_density(MA.p, normalized = T)+
  labs(title = "plotDensities with normalization between the two chanel", subtitle = "Loess normalization")

# Normalisation Between Arrays (Aquantile)
MA.pAq <- normalizeBetweenArrays(MA.p, method="Aquantile")

p_quantile <- plot_rg_density(MA.pAq, normalized = T)+
  labs(title = "plotDensities with normalization between the two chanel and between sample", subtitle = "Quantile normalization")

# Filter on low value ? ####
RG.MA(MA.pAq)$R 

keep <- rowMeans(RG_background$R) > 1119.909 & rowMeans(RG_background$G) > 1119.909
keep <- rowMeans(RG_background$R) > 1119.909 & rowMeans(RG_background$G) > 1119.909

RG_filtered <- RG[keep, ]

# Calculate how many genes were removed
n_removed <- n_original - n_filtered

# Calculate the percentage of genes removed
percent_removed <- (n_removed / n_original) * 100

# Display results
cat("Number of genes removed:", n_removed, "\n")
cat("Percentage of genes removed:", round(percent_removed, 2), "%\n")

p_filtered <- plot_rg_density(RG_filtered)+labs(title = "plotDensities without normalization but with filtration on low intensity")
fig_export(here::here("report/microarray/plot/normalization/without_normalization_with_low_intensity"), p_filtered, height_i = 5, width_i = 12, res_i = 300)

# Export ####
fig_export(here::here("report/microarray/plot/normalization/without_normalization"), p_raw, height_i = 5, width_i = 12, res_i = 300)
fig_export(here::here("report/microarray/plot/normalization/without_normalization_with_background_substraction"), p_background, height_i = 5, width_i = 12, res_i = 300)
fig_export(here::here("report/microarray/plot/normalization/with_loess_normalization"), p_loess, height_i = 5, width_i = 12, res_i = 300)
fig_export(here::here("report/microarray/plot/normalization/with_quantile_normalization"), p_quantile, height_i = 5, width_i = 12, res_i = 300)

# all_fig <- p_raw/p_background/p_filtered/p_loess/p_quantile
all_fig <- p_raw/p_background/p_loess/p_quantile
fig_export(here::here("report/microarray/plot/normalization/all_step_normalization"), all_fig, height_i = 18, width_i = 12, res_i = 300)

# RG_filtered
save(RG,MA_raw, MA.pAq, file =here::here("data/microarray/output/normalized_data_microarray_leaf_PeaSulf.RData"))


# Delet of low intensity ####
n_original <- nrow(RG)   # Total genes before filtering
n_filtered <- nrow(RG_filtered)   # Total genes after filtering
red_row_means <- rowMeans(RG_background$R)

# 2. Combine with Gene IDs into a data frame
df_red_means <- data.frame(
  GeneID = RG_background$genes,  # Adjust column name if needed
  MeanRed = red_row_means
)

# 3. Subset for the specific gene of interest
df_red_means[df_red_means$ID == "PsCam042688", ]
df_red_means[df_red_means$ID == "PsCam036750", ]

row_means_red <- rowMeans(RG$R)

# 2. Convert to a data frame
df_row_means <- data.frame(mean_intensity = row_means_red)

# 3. Plot with ggplot2
ggplot(df_row_means, aes(x = log2(mean_intensity))) +
  geom_histogram(bins = 50, color = "black", fill = "skyblue") +
  theme_minimal() +
  labs(
    title = "Distribution of Row Means (Red Channel)",
    x = "Mean Intensity (Red Channel)",
    y = "Count"
)

ggplot(df_row_means, aes(x = log2(mean_intensity))) +
  geom_histogram(bins = 50, color = "black", fill = "skyblue") +
  scale_x_log10() +
  theme_minimal() +
  labs(
    title = "Distribution of Row Means (Red Channel) [Log Scale]",
    x = "Log2 (Mean Intensity) (Red Channel) [log10]",
    y = "Count"
  )
# RG is your RGList with RG$R and RG$G (non-log)
# Suppose you want to keep probes with average raw intensity above 100 in *both* channels


# Combination of all plots

Expression of PSULT4 with normalized data

Code

load(file =here::here("data/microarray/output/normalized_data_microarray_leaf_PeaSulf.RData"))
list_microarray <-read_csv(here::here("data/microarray/output/list_microarray_microarray_leaf_PeaSulf.csv"), show_col_types = FALSE) 

ID <- MA.pAq$genes[[1]]
A <- MA.pAq$A  # Average log intensity
M <- MA.pAq$M  # Log ratio

R <- 2^(A + (M / 2))  # Calculate Red intensities
G <- 2^(A - (M / 2))  # Calculate Green intensities

R <- as.data.frame(R)
G <- as.data.frame(G)

rownames(R) <- ID
rownames(G) <- ID
# Red had good name of column but not green data. i will change that
colnames(G) <- colnames(RG$G)


RG_norm_filtr <- cbind(R, G)
RG_norm_filtr_log2<-log2(RG_norm_filtr)

RG_norm_filtr_log2_h <- RG_norm_filtr_log2 %>% 
  rownames_to_column("ID") %>% 
  pivot_longer(-ID, values_to = "value", names_to = "sample_id") %>% 
  left_join(.,list_microarray, by = "sample_id")

#export data
write_csv(RG_norm_filtr_log2_h, here::here("data/microarray/output/normalized_data_microarray_leaf_PeaSulf.csv"))

Comparison of the two chanel

Code

RG_norm_filtr_log2_h = read_csv(here::here("data/microarray/output/normalized_data_microarray_leaf_PeaSulf.csv"), show_col_types = FALSE)

p_sult4 <-  RG_norm_filtr_log2_h %>% 
    mutate(sulfur_condition = forcats::fct_relevel(sulfur_condition, "SS", "SD"),
           genotype = forcats::fct_relevel(genotype, "WT1", "W78*", "WT2", "E568K")) %>% 
    filter(ID == "PsCam042688") %>%
  ggplot(aes(x = genotype, y = value, col= color, fill = sulfur_condition, group = interaction(genotype, sulfur_condition, color))) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(position = position_jitterdodge(jitter.width = 0.6)) +
  scale_fill_manual(
    values = c("SS" = "#FFD74C", "SD" = "gray"), # Adjust based on your `color` levels
    name = "Sulfure Condition"  # Custom title for color legend
  )+
    scale_color_manual(
    values = c("Red" = "red", "Green" = "forestgreen"), # Adjust based on your `color` levels
    name = "Channel"  # Custom title for color legend
  )+ 
  theme_minimal()+
    theme(
  # Hide panel borders and remove grid lines
  panel.border = element_blank(),
  panel.grid.major = element_blank(),
  panel.grid.minor = element_blank(),
  # Change axis line
  axis.line = element_line(colour = "black")
  )+
    labs(x = "Genotype", y = "Microarray gene expression of SULTR4")
  
fig_export(here::here("report/microarray/plot/exression_of_genes/SULTR4"), p_sult4, height_i = 4, width_i = 6, res_i = 300)

With Stats

Code

RG_norm_filtr_log2_h = read_csv(here::here("data/microarray/output/normalized_data_microarray_leaf_PeaSulf.csv"), show_col_types = FALSE)

df_sult4 <-  RG_norm_filtr_log2_h %>% 
  as.data.frame() %>% 
  mutate(sulfur_condition = forcats::fct_relevel(sulfur_condition, "SS", "SD"),
         genotype = forcats::fct_relevel(genotype, "WT1", "W78*", "WT2", "E568K")) %>% 
  filter(ID == "PsCam042688")

l_sult4 <- stat_analyse(
    data=df_sult4 %>% 
      #mutate(climat_condition=paste0(water_condition,"_",heat_condition)) %>% 
      #mutate(condition=factor(condition,levels=c("Sto_WW_OT","Stoc_WS_OT","Sto_WW_HS","Sto_WS_HS","Wen_WW_OT","Wen_WS_OT","Wen_WW_HS","Wen_WS_HS"))) %>% 
      #mutate(climat_condition=factor(climat_condition,levels=c("WW_OT","WS_OT","WW_HS","WS_HS"))) %>% 
      drop_na(value),
    column_value = "value",
    category_variables = c("sulfur_condition"),
    grp_var = "genotype",
    show_plot = T,
    outlier_show = F, 
    label_outlier = "sample_id",
    biologist_stats = T,
    # Ylab_i = expression(paste("Angle ", widehat(ABC2)," (degrees) for order 1")),
    Ylab_i = paste("Microarray gene expression of SULTR4"),
  control_conditions = c("SS"),
  strip_normale = F,
  hex_pallet = sulfate_pallet
)

p_sult4 <-l_sult4[["plot"]]+labs(color="Treatment",fill="Treatment")

fig_export(here::here("report/microarray/plot/exression_of_genes/SULTR4_stats"), p_sult4, height_i = 3.5, width_i = 5, res_i = 300)

13.3.1 Verification of the importance of data standardization

Code

# import no normalize data
load(here::here("data/microarray/output/raw_data_microarray_leaf_PeaSulf.RData"))

# creation of the quantile normalization dataset
MA_aquantil <- normalizeBetweenArrays(MA_raw, method="Aquantile")

RG_aquantil <- RG.MA(MA_aquantil)
colnames(RG_aquantil$G) = colnames(RG$G)

# creation of df global with all data
RG_global <- bind_cols(RG$genes, RG$R, RG$G) %>% 
  column_to_rownames("ID")

RG_aquantil_global <- bind_cols(RG_aquantil$genes, RG_aquantil$R, RG_aquantil$G) %>% 
  column_to_rownames("ID")

# transfor into log2
RG_global_log2 <-  RG_global %>% mutate(across(everything(), log2))
RG_aquantil_global_log2 <-  RG_aquantil_global %>% mutate(across(everything(), log2))


# mean for each column
mean_add_info <- function(data) {
  data %>%
    dplyr::summarise(across(everything(), mean, na.rm = TRUE)) %>%
    pivot_longer(cols = everything(), names_to = "sample_name", values_to = "Mean") %>%
    mutate(sulfur_condition = case_when(
      grepl("SS", sample_name) ~ "SS",
      grepl("SD", sample_name) ~ "SD",
      TRUE ~ "Other"
    )) %>%
    mutate(color = case_when(
      grepl("Red", sample_name) ~ "Red",
      grepl("Green", sample_name) ~ "Green",
      TRUE ~ "Other"
    )) %>%
    mutate(num_combination = str_extract(sample_name, "\\d+$")) %>%
    mutate(
      sample_num = str_extract(sample_name, "\\d+\\.\\d+"),
      sample_num = str_replace_all(sample_num, c("\\_" = "\\.")), 
      simplify_condition = str_extract(sample_name, "WT|Mut"),
      genotype = case_when(
        sample_num %in% c("1.1", "2.2", "3.3", "4.4", "17.1", "18.2", "19.3", "20.4") ~ "WT2",
        sample_num %in% c("5.1", "6.2", "7.3", "8.4", "21.1", "22.2", "23.3", "24.4") ~ "E568K",
        sample_num %in% c("9.5", "10.6", "11.7", "12.8", "25.5", "26.6", "27.7", "28.8") ~ "WT1",
        sample_num %in% c("13.5", "14.6", "15.7", "16.8", "29.5", "30.6", "31.7", "32.8") ~ "W78*",
        TRUE ~ sample_num
      ),
      simplify_condition = paste0(genotype, "_", sulfur_condition)
    ) %>% 
    mutate(simplify_condition = forcats::fct_relevel(simplify_condition, "WT1_SS", "W78*_SS","WT1_SD", "W78*_SD", "WT2_SS", "E568K_SS", "WT2_SD", "E568K_SD"))
}

RG_global_mean <- mean_add_info(RG_global_log2)
RG_aquantil_global_mean <- mean_add_info(RG_aquantil_global_log2)


# Visualisation of the different boxplot (for mean density)
p_raw <- ggplot(RG_global_mean, aes(x = simplify_condition, y = Mean, fill = color, colour = color)) +
  geom_point(position = position_jitterdodge()) +
  geom_boxplot(alpha = 0.3, outlier.shape = NA, col = "black") +
  theme_minimal() +
  labs(
    title = "Without normalisation",
    y = "Mean of log2", 
    fill = "Channel", 
    colour = "Channel"
  ) +
  scale_color_manual(values = c("Red" = "red", "Green" = "forestgreen")) +
  scale_fill_manual(values = c("Red" = "red", "Green" = "forestgreen"))

p_norm <- ggplot(RG_aquantil_global_mean, aes(x = simplify_condition, y = Mean, fill = color, colour = color)) +
  geom_point(position = position_jitterdodge()) +
  geom_boxplot(alpha = 0.3, outlier.shape = NA, col = "black") +
  theme_minimal() +
  labs(
    title = "With normalisation",
    y = "Mean of log2", 
    fill = "Channel", 
    colour = "Channel"
  ) +
  scale_color_manual(values = c("Red" = "red", "Green" = "forestgreen")) +
  scale_fill_manual(values = c("Red" = "red", "Green" = "forestgreen"))

p_comparison <- p_raw + p_norm + plot_layout(guides = "collect")

fig_export(here::here(paste0("report/microarray/plot/normalization/test_normalisation/boxplot_with_witout_normalization")), p_comparison <- p_raw + p_norm + plot_layout(guides = "collect"), height_i = 6, width_i = 16, res_i = 600)

# For PCA
## transform into matrix
make_matrix_add_info <- function(data) {
  data %>%
  t() %>% 
  as.data.frame() %>% 
  rownames_to_column("sample_id") %>% 
  left_join(.,read_csv(here::here("data/microarray/output/list_microarray_microarray_leaf_PeaSulf.csv"), show_col_types = FALSE) %>% 
  mutate(sample_id_simplify = paste(sep = "_", simplify_condition, genotype, num_combination), 
         condition = paste(sep = "_", simplify_condition, genotype)) %>% 
        dplyr::select(sample_id, condition, sulfur_condition, genotype), by = "sample_id") %>% 
   mutate(
    splitted = str_split_fixed(sample_id, "_", 6)
  ) %>%
  mutate(
    sample_id = paste(splitted[, 3], splitted[, 5], splitted[, 6], sep = "_"),
    genotype_simplify = splitted[,1]
  ) %>%
  dplyr::select(-splitted) %>% 
  column_to_rownames("sample_id") %>% 
  mutate(across(c(condition, genotype, sulfur_condition,genotype_simplify), as.factor)) %>% 
  mutate(condition = fct_relevel(condition, "WT_SS_WT1", "WT_SD_WT1", "Mut_SS_W78*", "Mut_SD_W78*",
                                 "WT_SS_WT2", "WT_SD_WT2", "Mut_SS_E568K", "Mut_SD_E568K"),
         genotype = fct_relevel(genotype, "WT1", "W78*", "WT2", "E568K"), 
         sulfur_condition = fct_relevel(sulfur_condition, "SS", "SD"),
         genotype_simplify = fct_relevel(genotype_simplify, "WT", "Mut")
         )  
}

RG_log2_t_info <- make_matrix_add_info(RG_global_log2)
RG_aquantil_log2_t_info <- make_matrix_add_info(RG_aquantil_global_log2)

## transform into pca for sulfur condition
RG_log2_t <- RG_log2_t_info %>% select(-c("condition", "genotype", "sulfur_condition","genotype_simplify"))
RG_aquantil_log2_t <- RG_aquantil_log2_t_info %>% select(-c("condition", "genotype", "sulfur_condition","genotype_simplify"))

res_pca <- PCA(RG_log2_t, quali.sup = 1:4, graph = F)
res_pca_aquantil <- PCA(RG_aquantil_log2_t, quali.sup = 1:4, graph = F)

fviz_eig(res_pca, addlabels = TRUE, ylim = c(0, 50))
fviz_eig(res_pca_aquantil, addlabels = TRUE, ylim = c(0, 50))
var <- get_pca_var(res_pca)
var_aquantil <- get_pca_var(res_pca_aquantil)

# By sulfur condition ####
PCA_indiv<-fviz_pca_ind(res_pca,
 #            geom.ind = "point", # Shows points only (not text)
             pointshape = 20,pointsize = 2,
             col.ind = RG_log2_t_info %>% pull(sulfur_condition), # colorer by groups
             palett=sulfate_pallet,
             addEllipses = T, # Concentration ellipses
             legend.title = "Sulfur condition", 
             ellipse.level = 0.95,
             ellipse.type = c("norm"),
             
)+ ggtitle("Visualizing Individual PCA wthout normalization \nby sulfur condition")
PCA_indiv

PCA_indiv_aquantil<-fviz_pca_ind(res_pca_aquantil,
 #            geom.ind = "point", # Shows points only (not text)
             pointshape = 20,pointsize = 2,
             col.ind = RG_log2_t_info %>% pull(sulfur_condition), # colorer by groups
             palett=sulfate_pallet,
             addEllipses = T, # Concentration ellipses
             legend.title = "Sulfur condition", 
             ellipse.level = 0.95,
             ellipse.type = c("norm"),
             
)+ ggtitle("Visualizing Individual PCA normalize with aquantil \nby sulfur condition")
PCA_indiv_aquantil

p_comparison <- PCA_indiv+PCA_indiv_aquantil+plot_layout(guides = "collect")
fig_export(here::here(paste0("report/microarray/plot/PCA/PCA_sulfur_condition_with_without_normalization")), p_comparison, height_i = 10, width_i = 22, res_i = 600)

# by genotype for each condition ####
make_PCA_genotype <- function(df_info, norm = T){
  df <- df_info %>% select(-c("condition", "genotype", "sulfur_condition","genotype_simplify"))
  res_pca <- PCA(df, quali.sup = 1:4, graph = F)
  
  #fviz_eig(res_pca, addlabels = TRUE, ylim = c(0, 50))
  var <- get_pca_var(res_pca)
  PCA_indiv<-fviz_pca_ind(res_pca,
   #            geom.ind = "point", # Shows points only (not text)
               #pointshape = 20,pointsize = 2,
               col.ind = df_info %>% pull(genotype), # colorer by groups
               shape.ind   = df_info$genotype_simplify,  # Forme selon le génotype
               palett= c("#003049", "#780000", "#7FACC7", "#EC323E"),
               addEllipses = F, # Concentration ellipses
               legend.title = "Treatment", 
               ellipse.level = 0.95,
               ellipse.type = c("norm")
               
  )+ ggtitle(paste0("Visualizing Individual PCA ", ifelse(norm ==T, "with", "without")," normalization \nby genotype for ",df_info$sulfur_condition[1]))
  return(PCA_indiv)
}

p_SS_raw = RG_log2_t_info %>% 
  filter(sulfur_condition =="SS") %>% 
  make_PCA_genotype(., norm=F)

p_SD_raw = RG_log2_t_info %>% 
  filter(sulfur_condition =="SD") %>% 
  make_PCA_genotype(., norm=F)

p_SS_norm = RG_aquantil_log2_t_info %>% 
  filter(sulfur_condition =="SS") %>% 
  make_PCA_genotype(., norm=T)

p_SD_norm = RG_aquantil_log2_t_info %>% 
  filter(sulfur_condition =="SD") %>% 
  make_PCA_genotype(., norm=T)

p_compil <- p_SS_raw+p_SS_norm+p_SD_raw+p_SD_norm+plot_layout(guides = "collect", ncol = 2)

fig_export(here::here(paste0("report/microarray/plot/PCA/PCA_raw_norm_genotype")), p_compil, height_i = 10, width_i = 14, res_i = 600)

# export for Marie-Laure
matrix_info <- RG_aquantil_global_mean %>% 
  dplyr::select(-Mean)
save(matrix_info, RG, MA_aquantil,file =  here::here("data/microarray/output/targets_design.RData"))

13.4 Differential expression of normalized data

Bonus: If i apply a double normalization ( not necessary for this experiment), here the vulcanoplot with normalization)

Of normalized data

Code

# /!\ vey slow
source (here::here("src/function/microarray/vulcanoplot_microarray_RG.R"))

for(i in 1:5){
  for(normalized_i in c(T, F)){
  # create ggplot
  p = vulcanoplot_microarray_RG(normalized = normalized_i, list_comparison_num_i =  i)
  fig_export(here::here(paste0("report/microarray/plot/volcanoplot/volcano_interactive_RG_",ifelse(normalized_i ==T, "normalized_","raw_"),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_RG_",ifelse(normalized_i ==T, "normalized_","raw_"),i,".html")))
  }
}

Code

load(file =  here::here("data/microarray/output/ratio_stats_microarray_leaf_PeaSulf.RData")) # only to show name of vulcanoplot

From normalisez data
From Raw data

For an interactive version : Click here

13.4.1 Using limma contrasts (method recommended by Marie Laure Martin-Magniette)

In multi-factor or multi-group microarray analyses, we often want to compare specific conditions such as different genotypes, treatments, or their combinations.

Why use contrasts? A linear model created with lmFit describes the contribution of each factor in the experiment, but it doesn’t automatically test the specific comparisons we are interested in. A contrast is a linear combination of the model’s coefficients that defines a particular comparison, for example, comparing genotype A to genotype B, or comparing interaction effects. Using a contrast matrix allows us to explicitly define and test these comparisons, ensuring that limma computes the appropriate log-fold changes and p-values.

How it works: First, we create a design matrix using model.matrix(…), which captures the main effects and possibly the interactions between experimental factors. Then we define a contrast matrix using a custom function (like Contrasts) or limma’s built-in tools. Each row of the contrast matrix corresponds to a specific biological question or comparison.

After this, the lmFit(…) function estimates gene-wise model coefficients. We apply our contrasts using contrasts.fit(…), which re-expresses the model in terms of the contrasts of interest. The eBayes(…) function then moderates the variance estimates and calculates statistical significance (p-values and adjusted p-values).

In the rest of the code, we use p.adjust(…) to control for multiple testing (FDR), identify significantly up- or down-regulated genes based on a fold-change threshold, and generate summary outputs such as barplots and volcano plots. This provides a clear overview of differential expression across the comparisons that matter most in the study.

In summary, contrasts are essential in differential analysis using limma, as they allow us to tailor the comparisons to the specific hypotheses of the experiment, especially when multiple factors and interactions are involved.

Code

#parameter
lfc_lim_i = 0
pval_i =0.05

load(file =  here::here("data/microarray/output/targets_design.RData"))
source(here::here("src/function/microarray/Contrasts.R"))

# creation of a complete dataframe
RG_aquantil <- limma::RG.MA(MA_aquantil)
colnames(RG_aquantil$G) <- colnames(RG$G)
df_RG_aquantil_log2 <- cbind(RG_aquantil$R, RG_aquantil$G, RG_aquantil$genes) %>% column_to_rownames("ID") %>%  mutate(across(everything(), log2))

save(df_RG_aquantil_log2, file = here::here("data/microarray/output/df_RG_aquantil_log2.RData"))

Target <- matrix_info %>% dplyr::select(sulfur_condition, genotype, num_combination)
Target=as.data.frame(Target)
Target$sulfur_condition = as.factor(Target$sulfur_condition)
Target$genotype = as.factor(Target$genotype)

model <- model.matrix(as.formula(paste("~", paste(colnames(Target)[1:2], collapse = " + "), paste(colnames(Target[,1:2]),collapse = ":"), sep = " + ")), data=Target)
contrast.matrix = Contrasts(model,Target,FALSE,TRUE)

mycontrast1 <- c(0, 0, +1, -1, +1, +1, -1, +1)
mycontrast2 <- c(0, 0, +1, -1, +1, 0, 0, 0)

contrast.matrix  # votre matrice après Contrasts(...)

# On l'augmente en ajoutant nos 2 lignes :
contrast.matrix <- rbind(
  contrast.matrix,
  mycontrast1,
  mycontrast2
)

# On renomme les deux nouvelles lignes :
n <- nrow(contrast.matrix)
rownames(contrast.matrix)[(n-1):n] <- c(
  "[SS_W78*-SS_WT1]-[SS_E568K-SS_WT2]",
  "[SD_W78*-SD_WT1]-[SD_E568K-SD_WT2]"
)

contrast.matrix


fit <- lmFit(df_RG_aquantil_log2, model)

#contrast.matrix <- makeContrasts(SpecialDiet-Control, SpecialDietDrug1-SpecialDiet, SpecialDietDrug1-Control, levels=design) 

fit2 <- contrasts.fit(fit, t(contrast.matrix))

fit2 <- eBayes(fit2)

# topTable(fit2, coef=1, adjust="fdr", sort.by="B", number=10) 
# creation of a plot to count the number of difference between eatch contraste
p_adj <- apply(fit2$p.value, 2, p.adjust, method = "BH")

is_significant <- p_adj < 0.05  # Matrice TRUE/FALSE pour FDR < 0.05
is_upregulated   <- is_significant & (fit2$coefficients > lfc_lim_i)  # logFC > 0
is_downregulated <- is_significant & (fit2$coefficients < -lfc_lim_i)  # logFC < 0

vector_with_interest <- c("[SD-SS]","[WT1-WT2]","[E568K-WT2]", "[W78*-WT1]", "[WT1_SD-WT1_SS]","[W78*_SD-W78*_SS]", "[WT2_SD-WT2_SS]", "[E568K_SD-E568K_SS]",  "[SS_W78*-SS_WT1]",  "[SS_E568K-SS_WT2]",  "[SD_W78*-SD_WT1]", "[SD_E568K-SD_WT2]", "[SD_W78*-SD_WT1]-[SS_W78*-SS_WT1]",   "[SD_E568K-SD_WT2]-[SS_E568K-SS_WT2]", "[SS_W78*-SS_WT1]-[SS_E568K-SS_WT2]", "[SD_W78*-SD_WT1]-[SD_E568K-SD_WT2]" )

summary_DEG <- tibble(
  Contrast = colnames(fit2$p.value),
  Upregulated = colSums(is_upregulated),
  Downregulated = colSums(is_downregulated)
) %>% 
  pivot_longer(cols = -c("Contrast"), names_to = "sign", values_to = "length") %>% 
  mutate(sign = recode(sign,
                         "Upregulated" = "Up",
                         "Downregulated" = "Down")) %>% 
  dplyr::group_by(Contrast) %>%
  mutate(
    total   = sum(length),
    percent = round(100 * length / total)
  ) %>%
  ungroup() %>%
  filter(Contrast %in% c(vector_with_interest)) %>% 
  mutate(sign = fct_relevel(sign, c("Up", "Down")), 
         Contrast = fct_relevel(Contrast, rev(vector_with_interest)))

# Show results 
px <- ggplot(summary_DEG, aes(x = Contrast, y = length, fill = sign)) +
  geom_col() +
  
  # 1) Texte pour les barres "assez grandes" (>= 5)
  geom_text(
    data = subset(summary_DEG, length >= ifelse(lfc_lim_i == 0, 600,30)),
    aes(
      label = paste0(length, "\n(", percent, "%)")
    ),
    position = position_stack(vjust = 0.5),
    color = "white",  # texte en blanc
    size  = 3,
    # hjust = 0.5 (par défaut) - pas indispensable
  ) +
  
  # 2) Texte pour les barres "trop petites" (< 5)
  geom_text(
    data = subset(summary_DEG, length < ifelse(lfc_lim_i == 0, 600,30)),
    aes(
      label = paste0(length, " (", percent, "%)"),
      color = sign,   # couleur du texte = même code couleur que "sign"
      vjust = ifelse(sign == "Up", -0.5, 1.5),
      y=total
    ),
    #position = position_stack(vjust = 0.5),
    hjust= -0.2,
    
    # On décale à droite (hors de la barre). En coord_flip(), 
    # hjust < 0 place le texte davantage à droite.
    #hjust = -0.2,
    size  = 3
  ) +
  theme_minimal() +
  labs(
    x     = "Comparison",
    y     = "Number of genes deregulated",
    title = paste0("Number of genes deregulated for each contrast (FDR=0.05 ; LogFC=", lfc_lim_i, ")")
  ) +
  coord_flip() +
  
  # Couleurs de remplissage (barres)
  scale_fill_manual(
    values = c('Down' = "#1d4877", 'Up' = "#ee3e32"),
    name   = "Sign"
  ) +
  
  # Couleurs du texte (identique aux barres), sans nouvelle légende
  scale_color_manual(
    values = c('Down' = "#1d4877", 'Up' = "#ee3e32"),
    guide = FALSE
  ) ; px

fig_export(path = here::here(paste0("report/microarray/plot/DEG/barplot_differential_expression_lfc_",lfc_lim_i)), px, height = 6.5, width = 12, res_i = 600)

# export eatch contrast
res_list <- list()

# Loop through each contrast in your vector of interest
for (contrast in vector_with_interest) {
  
  # Construct vector names, e.g., "[SD-SS](Up)" and "[SD-SS](Down)"
  up_name   <- paste0("(Up)", contrast)
  down_name <- paste0("(Down)", contrast)
  
  # Extract the rownames of all Upregulated genes for this contrast
  genes_up <- rownames(is_upregulated)[is_upregulated[, contrast]]
  
  # Extract the rownames of all Downregulated genes for this contrast
  genes_down <- rownames(is_downregulated)[is_downregulated[, contrast]]
  
  # Store them into the result list with meaningful names
  res_list[[up_name]]   <- genes_up
  res_list[[down_name]] <- genes_down
}

save(res_list, file = here::here(paste0("data/microarray/output/contraste_lfc_",lfc_lim_i,".RData")))

# make vulcanoplot for each contraste
for (contrast in vector_with_interest) {
  cat(paste0(contrast, "\n"))
  volcano_data <- tibble(
      Contrast = contrast,
      Gene = rownames(fit2$coefficients),
      logFC = fit2$coefficients[, contrast],
      pvalue = fit2$p.value[, contrast],
      adj_pvalue = p_adj[, contrast]
    ) %>%
      mutate(
        Significance = case_when(
          adj_pvalue < 0.05 & logFC > lfc_lim_i  ~ "Upregulated",
          adj_pvalue < 0.05 & logFC < -lfc_lim_i ~ "Downregulated",
          TRUE ~ "Not significant"
        )
      ) %>% 
    dplyr::rename(PsCam = Gene) %>% 
    # add info of genes
    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") %>% 
    left_join(.,read_csv(here::here("data/microarray/output/raw_data_microarray_leaf_PeaSulf.csv"), show_col_types = FALSE) %>%
                as.data.frame() %>%
        dplyr::select(c("id_probe", "CAT")) %>%
        dplyr::rename(PsCam = id_probe),
      by = "PsCam"
    )
  
    num_upregulated <- sum(volcano_data$Significance == "Upregulated")
    num_downregulated <- sum(volcano_data$Significance == "Downregulated")
    
    # Création du titre avec couleurs HTML
    plot_title <- paste0(
      "Volcano Plot - ", contrast, "<br>",
      " | <span style='color:#B02428;'>Upregulated: ", num_upregulated, "</span> | ",  
      "<span style='color:#6697EA;'>Downregulated: ", num_downregulated, "</span>"
    )

    p1 <- ggplot(volcano_data, aes(x = logFC, y = -log10(adj_pvalue), color = Significance)) +
    geom_point(alpha = 0.7) + 
    scale_color_manual(values = c("Upregulated" = "#B02428", "Downregulated" = "#6697EA", "Not significant" = "grey")) +
    geom_vline(xintercept = c(-lfc_lim_i, lfc_lim_i), linetype = "dashed") +
    geom_hline(yintercept = -log10(pval_i), linetype = "dashed") +
      geom_text_repel(
      data = volcano_data %>% filter(Significance %in% c("Upregulated", "Downregulated")) %>% 
        slice_max(abs(logFC),n = 20),    
      aes(label = sub("^symbols:", "", CAT), color = Significance),
      size = 3,
      max.overlaps = 0
    ) +
      
    
  theme_minimal() +
  labs(
    title = plot_title,
    x = expression(log[2]~Fold~Change),
    y = expression(-log[10](p-value))
  ) +
  
  theme(
    plot.title = element_markdown()  # Permet d'afficher du texte formaté en HTML
  )
    fig_export(path = here::here(paste0("report/microarray/plot/volcanoplot/volacano_contrast_", gsub("\\*", "", contrast), "_lfc_",lfc_lim_i)), p1, height = 6, width = 12, res_i = 300)
        
    # for interactive 
    p1_interactive <- ggplot(volcano_data, aes(x = logFC, y = -log10(adj_pvalue), color = Significance)) +
    geom_point(aes(text = paste("PsCam:", PsCam, "<br>Psat:", Psat, "<br>CAT:", CAT, "<br>logFC:", logFC, "<br>FDR:", adj_pvalue)), size = 1, alpha = .6) +
    scale_color_manual(
      values = c(
        "Not significant" = "gray",                           
        "Downregulated" = "#6697EA",                                  
        "Upregulated" = "#B02428"                                   
      )
    ) +
    theme_minimal() +
      labs(
    title = plot_title,
     x = "log2 Fold Change",  # Replace expression with plain text
      y = "-log10(p-value)"
  ) +
  
  theme(
    plot.title = element_markdown()  # Permet d'afficher du texte formaté en HTML
  )
  ggplotly(p1_interactive) %>% 
    saveWidget(here::here(paste0("report/microarray/plot/volcanoplot/volcano_interactive_contraste_",gsub("\\*", "", contrast),"_lfc_",lfc_lim_i,".html")))
}

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

For an interactive version : Click here

13.5 Gene onthologie on contraste

verification and simplification of the differents fonctions
add treemap plot

Code

# parameter
lfc_lim_i = 0
type_calcul_i = "hypergeometric" #clusterProfiler # hypergeometric # fisher # in reality its exactly the same

load(file = here::here(paste0("data/microarray/output/contraste_lfc_",lfc_lim_i,".RData")))

gene_table <- tibble(
  contrast = rep(names(res_list), lengths(res_list)),
  PsCam = unlist(res_list)
) %>% 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") %>% 
  drop_na(Psat)

nrow(gene_table %>% filter(contrast == "(Down)[SD-SS]"))

universe_vector <- unique(gene_table$Psat)

df_GO = GO_on_different_group(functional_roles = "BP",
                      group_info = gene_table,
                      group = "contrast",
                      ID = "Psat",
                      top = 10,
                      universe_i = universe_vector, 
                      type_calcul = type_calcul_i
                                            )

# combined_results = GO_on_different_group(functional_roles = "BP",
#                       group_info = gene_table,
#                       group = "contrast",
#                       ID = "Psat",
#                       top = FALSE,
#                       universe_i = universe_vector, 
#                       type_calcul = type_calcul_i
#                                             )



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, "(?<=\\().*?(?=\\))")
  )
# 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 = paste0("GO terme for all gene deregulated with lfc:",lfc_lim_i)) +
  #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

fig_export(path = paste0("report/microarray/plot/GO/GO_BP_contrast_",type_calcul_i,"_",lfc_lim_i), plot_x = px, height_i = 20, width_i = 12, res = 600)

Test

Code

################## test WT1 vs WT2 in contraste 
# test WT1 vs WT2 realy different also for sulfure



df_test = test %>% filter(Description_GO %in% c("sulfate transport", "sulfur compound transport"),
                group == "(Down)[WT1-WT2]") %>% pull(geneID) %>% str_split("/") %>%unlist() %>%unique() %>% tibble(Psat = .) %>% 
  left_join(., 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 = "Psat")) %>% 
  left_join(., left_join(.,read_excel(here::here("data/microarray/Resultats_4plex-POIS-2014_01_230415_modKG.xlsx"), sheet = "Complete", col_names = T, skip = 11)%>% 
              dplyr::rename(PsCam = "N°") %>% dplyr::select(PsCam, CAT ), by = "PsCam") )

load(file =  here::here("data/microarray/output/targets_design.RData"))
source(here::here("src/function/stat_function/stat_analysis_main.R")) # for make plot 

# create data frame
RG_aquantil <- limma::RG.MA(MA_aquantil)
colnames(RG_aquantil$G) <- colnames(RG$G)
df_RG_aquantil_log2 <- cbind(RG_aquantil$R, RG_aquantil$G, RG_aquantil$genes) %>% column_to_rownames("ID") %>%  mutate(across(everything(), log2)) %>% 
  rownames_to_column("PsCam") %>% 
  pivot_longer(-PsCam, names_to = "sample_name") %>% 
   dplyr::mutate(
        sample_num =  str_extract(sample_name, "(?<=_[A-Z]{2}_)[0-9\\.]+"),
        sample_num = str_replace_all(sample_num, c("\\_" = "\\.")), 
        simplify_condition = str_extract(sample_name, "WT|Mut"),
        simplify_condition = str_extract(sample_name, "SS|SD"),
        sulfur_condition = str_sub(simplify_condition, -2, -1),
        genotype = ifelse(sample_num %in% c(1.1, 2.2, 3.3, 4.4, 17.1, 18.2, 19.3, 20.4), "WT2", sample_num),# Add genotype
        genotype = ifelse(sample_num %in% c(5.1, 6.2, 7.3, 8.4, 21.1, 22.2, 23.3, 24.4), "E568K", genotype),# Add genotype
        genotype = ifelse(sample_num %in% c(9.5, 10.6, 11.7, 12.8, 25.5, 26.6, 27.7, 28.8), "WT1", genotype),# Add genotype
        genotype = ifelse(sample_num %in% c(13.5, 14.6, 15.7, 16.8, 29.5, 30.6, 31.7, 32.8), "W78*", genotype),# Add genotype
        sample_id = paste(sep = "_", simplify_condition, genotype)
      )

df_prepared <- df_test_value %>%
  as.data.frame() %>%
  mutate(
    sulfur_condition = fct_relevel(sulfur_condition, "SS", "SD"),
    genotype         = fct_relevel(genotype, "WT1", "W78*", "WT2", "E568K")
  ) %>%
  drop_na(value)

# 2. Extraire les niveaux uniques de "psat_pscam" (plus rapide que as.factor(...) puis levels(...))
vector_sulfur <- unique(df_prepared$psat_pscam)

# 3. Boucler sur chaque valeur de psat_pscam pour générer la liste de graphiques
plots <- lapply(vector_sulfur, function(psat) {
  # Afficher le nom en vert (si cat_col est une fonction perso pour colorer la console)
  cat_col(paste0(psat, "\n"), color = "green")

  # Sous-ensemble des données correspondant à psat
  df_x <- df_prepared %>% 
    filter(psat_pscam == psat)

  # Analyse statistique
  l_sult4 <- stat_analyse(
    data               = df_x,
    column_value       = "value",
    category_variables = c("sulfur_condition"),
    grp_var            = "genotype",
    show_plot          = TRUE,
    outlier_show       = FALSE, 
    label_outlier      = "sample_name",
    biologist_stats    = TRUE,
    Ylab_i             = paste("Microarray gene expression of \n", df_x$CAT[1]),
    control_conditions = c("SS"),
    strip_normale      = FALSE,
    hex_pallet         = sulfate_pallet
  )

  # Ajout des labels et récupération du plot
  p_plot <- l_sult4[["plot"]] + labs(color = "Treatment", fill = "Treatment")

  # On retourne l’objet graphique
  return(p_plot)
})

final_plot <- wrap_plots(plots, ncol = 7) +
  plot_layout(guides = "collect") +
  plot_annotation() & theme(legend.position = 'bottom')

fig_export(here::here("report/microarray/plot/variation_gene_sulfur"), final_plot, height_i = 21, width_i = 29.7, res_i = 600)

With clusterprofiler
With fisher test or hypergeometric test

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

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

13.6 UpsetPlot

Code

# Define the logFC and p-value thresholds
lfc_lim_i <- 0
pval_i = 0.05
  
load(file =  here::here("data/microarray/output/targets_design.RData"))
source(here::here("src/function/microarray/Contrasts.R"))

# creation of a complete dataframe
RG_aquantil <- limma::RG.MA(MA_aquantil)
colnames(RG_aquantil$G) <- colnames(RG$G)
df_RG_aquantil_log2 <- cbind(RG_aquantil$R, RG_aquantil$G, RG_aquantil$genes) %>% column_to_rownames("ID") %>%  mutate(across(everything(), log2))

Target <- matrix_info %>% dplyr::select(sulfur_condition, genotype, num_combination )
Target=as.data.frame(Target)
Target$sulfur_condition = as.factor(Target$sulfur_condition)
Target$genotype = as.factor(Target$genotype)

model <- model.matrix(as.formula(paste("~", paste(colnames(Target)[1:2], collapse = " + "), paste(colnames(Target[,1:2]),collapse = ":"), sep = " + ")), data=Target)
contrast.matrix = Contrasts(model,Target,FALSE,TRUE)

fit <- lmFit(df_RG_aquantil_log2, model)

#contrast.matrix <- makeContrasts(SpecialDiet-Control, SpecialDietDrug1-SpecialDiet, SpecialDietDrug1-Control, levels=design) 

fit2 <- contrasts.fit(fit, t(contrast.matrix))

fit2 <- eBayes(fit2)

# Create a list of differentially expressed genes per contrast
list_genes_by_contrast <- lapply(colnames(fit2$p.value), function(contrast) {
  
  # Upregulated genes
  upregulated <- rownames(fit2)[is_upregulated[, contrast]]
  
  # Downregulated genes
  downregulated <- rownames(fit2)[is_downregulated[, contrast]]
  
  list(
    Up = upregulated,
    Down = downregulated
  )
})

# Name the list using the contrast names
names(list_genes_by_contrast) <- colnames(fit2$p.value)

# ➡️ Filter only the contrasts of interest
list_genes_by_contrast_interest <- list_genes_by_contrast[vector_with_interest]

# Check the result (for example, for the first contrast of interest)
# 
#  # 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` = list_genes_by_contrast_interest[[vector_with_interest[1]]]$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")))

Code

## Output DEG list only in EK mutant  
lfc_lim_i <- 0
pval_i = 0.05
  
load(file =  here::here("data/microarray/output/targets_design.RData"))
source(here::here("src/function/microarray/Contrasts.R"))

vector_with_interest <- c("[SD-SS]","[E568K-WT2]", "[W78*-WT1]", "[WT1_SD-WT1_SS]","[W78*_SD-W78*_SS]", "[WT2_SD-WT2_SS]", "[E568K_SD-E568K_SS]",  "[SS_W78*-SS_WT1]",  "[SS_E568K-SS_WT2]",  "[SD_W78*-SD_WT1]", "[SD_E568K-SD_WT2]", "[SD_W78*-SD_WT1]-[SS_W78*-SS_WT1]",   "[SD_E568K-SD_WT2]-[SS_E568K-SS_WT2]")

# creation of a complete dataframe
RG_aquantil <- limma::RG.MA(MA_aquantil)
colnames(RG_aquantil$G) <- colnames(RG$G)
df_RG_aquantil_log2 <- cbind(RG_aquantil$R, RG_aquantil$G, RG_aquantil$genes) %>% column_to_rownames("ID") %>%  mutate(across(everything(), log2))

Target <- matrix_info %>% dplyr::select(sulfur_condition, genotype, num_combination )
Target=as.data.frame(Target)
Target$sulfur_condition = as.factor(Target$sulfur_condition)
Target$genotype = as.factor(Target$genotype)

model <- model.matrix(as.formula(paste("~", paste(colnames(Target)[1:2], collapse = " + "), paste(colnames(Target[,1:2]),collapse = ":"), sep = " + ")), data=Target)
contrast.matrix = Contrasts(model,Target,FALSE,TRUE)

fit <- lmFit(df_RG_aquantil_log2, model)

#contrast.matrix <- makeContrasts(SpecialDiet-Control, SpecialDietDrug1-SpecialDiet, SpecialDietDrug1-Control, levels=design) 

fit2 <- contrasts.fit(fit, t(contrast.matrix))

fit2 <- eBayes(fit2)

# Adjust p-values using BH (Benjamini–Hochberg)
p_adj <- apply(fit2$p.value, 2, p.adjust, method = "BH")

# Determine significance based on 'pval_i'
is_significant <- p_adj < pval_i

# Define up-/downregulated based on 'lfc_lim_i' and significance
is_upregulated   <- is_significant & (fit2$coefficients >  lfc_lim_i)
is_downregulated <- is_significant & (fit2$coefficients < -lfc_lim_i)

# Now create the list of genes per contrast
list_genes_by_contrast <- lapply(colnames(fit2$p.value), function(contrast) {
  
  # Upregulated genes based on your chosen thresholds
  upregulated <- rownames(fit2)[is_upregulated[, contrast]]
  
  # Downregulated genes based on your chosen thresholds
  downregulated <- rownames(fit2)[is_downregulated[, contrast]]
  
  list(
    Up   = upregulated,
    Down = downregulated
  )
})

names(list_genes_by_contrast) <- colnames(fit2$p.value)

# Subset contrasts of interest
list_genes_by_contrast_interest <- list_genes_by_contrast[vector_with_interest]

# For Up ####
lt_up <- list(
  `SS_WT1 vs SS_W78*` = list_genes_by_contrast_interest[["[SS_W78*-SS_WT1]"]]$Up, # only PsCam042688 down 
  `SS_WT2 vs SS_E568K` = list_genes_by_contrast_interest[["[SS_E568K-SS_WT2]"]]$Up,
  `SD_WT1 vs SD_W78*` = list_genes_by_contrast_interest[["[SD_W78*-SD_WT1]"]]$Up,
  `SD_WT2 vs SD_E568K` = list_genes_by_contrast_interest[["[SD_E568K-SD_WT2]"]]$Up
)

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

# Down-regulated genes
lt_down <- list(
  `SS_WT1 vs SS_W78*` = list_genes_by_contrast_interest[["[SS_W78*-SS_WT1]"]]$Down, # only PsCam042688 down 
  `SS_WT2 vs SS_E568K` = list_genes_by_contrast_interest[["[SS_E568K-SS_WT2]"]]$Down,
  `SD_WT1 vs SD_W78*` = list_genes_by_contrast_interest[["[SD_W78*-SD_WT1]"]]$Down,
  `SD_WT2 vs SD_E568K` = list_genes_by_contrast_interest[["[SD_E568K-SD_WT2]"]]$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_contrast(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_contrast",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)),
     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)),
     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_contrast",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_contrast",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_contrast_",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_contrast_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*

In Up Regulated

Which group does the psult4 gene belong to?

Code

load(file = here::here("data/microarray/output/upset_result_condition_sign_lfc_contrast_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

13.7 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_contrast_",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"))

#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_contrast_",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_contrast_",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_lfc0 <- list(
  I1 = 'SS_WT2 vs SS_E568K', 
  I2 = 'SD_WT2 vs SD_E568K',
  I3  = c('SS_WT2 vs SS_E568K', 'SD_WT2 vs SD_E568K'),
  I4  = c('SD_WT1 vs SD_W78*', 'SD_WT2 vs SD_E568K'),
  I5 = 'SD_WT1 vs SD_W78*',
  I6 = c('SS_WT2 vs SS_E568K', 'SD_WT2 vs SD_E568K','SD_WT1 vs SD_W78*'),
  I7 = c('SS_WT2 vs SS_E568K','SD_WT1 vs SD_W78*'),
  I8 = c('SS_WT2 vs SS_E568K','SD_WT1 vs SD_W78*','SS_WT1 vs SS_W78*')
)

intersection_list_lfc1 <- list(
  I1 = 'SD_WT2 vs SD_E568K',
  I2 = 'SS_WT2 vs SS_E568K', 
  I3  = c('SS_WT2 vs SS_E568K', 'SD_WT2 vs SD_E568K'),
  I4  = c('SD_WT1 vs SD_W78*', 'SD_WT2 vs SD_E568K'),
  I5 = 'SD_WT1 vs SD_W78*',
  I6 = c('SD_WT1 vs SD_W78*','SS_WT1 vs SS_W78*')
)

if(lfc_lim_i == 1){
  intersection_list = intersection_list_lfc1
}else if(lfc_lim_i == 0){
  intersection_list = intersection_list_lfc0
}

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_contrast_", 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

13.8 Plot for publication

Bonus: Details are provided here.

For this plot I only took the standardised data and the comparisons or similarities that I found interesting. For the comparisons I generated a contrast matrix for a linear model which I then used to carry out the differential analysis in limma. For the similarity, I looked at which genes were similar to the two groups. I’m only going to do this for a log fold change equal to 1.

13.8.1 DEG

Code

lfc_lim_i = 0
load(file = here::here(paste0("data/microarray/output/contraste_lfc_",lfc_lim_i,".RData")))


## ─── 1. Contrastes de base ──────────────────────────────────────────────────────
c1    <- "[W78*-WT1]"
c2    <- "[E568K-WT2]"
c3_1  <- "[SS_W78*-SS_WT1]"
c3_2  <- "[SS_E568K-SS_WT2]"
c4_1  <- "[SD_W78*-SD_WT1]"
c4_2  <- "[SD_E568K-SD_WT2]"


## ─── 2. Commons (même clé partout) ──────────────────────────────────────────────
common1_key <- paste0(c1 , " common ", c2 )
common2_key <- paste0(c3_1, " common ", c3_2)
common3_key <- paste0(c4_1, " common ", c4_2)
common4_key <- paste0(c3_1, " common ", c4_1)
common5_key <- paste0(c3_2, " common ", c4_2)


for (pair in list(
  list(up  = paste0("(Up)"  , c1),   down  = paste0("(Down)"  , c1),
       up2 = paste0("(Up)"  , c2),   down2 = paste0("(Down)"  , c2),
       common_key = common1_key),
  list(up  = paste0("(Up)"  , c3_1), down  = paste0("(Down)" , c3_1),
       up2 = paste0("(Up)"  , c3_2), down2 = paste0("(Down)" , c3_2),
       common_key = common2_key),
  list(up  = paste0("(Up)"  , c4_1), down  = paste0("(Down)" , c4_1),
       up2 = paste0("(Up)"  , c4_2), down2 = paste0("(Down)" , c4_2),
       common_key = common3_key),
  
    list(up  = paste0("(Up)"  , c3_1), down  = paste0("(Down)" , c3_1),
       up2 = paste0("(Up)"  , c4_1), down2 = paste0("(Down)" , c4_1),
       common_key = common4_key),
  list(up  = paste0("(Up)"  , c3_2), down  = paste0("(Down)" , c3_2),
       up2 = paste0("(Up)"  , c4_2), down2 = paste0("(Down)" , c4_2),
       common_key = common5_key))
) {
  ## 2.1  Intersection Up / Down
  up_genes   <- intersect(res_list[[pair$up ]], res_list[[pair$up2 ]])
  down_genes <- intersect(res_list[[pair$down]], res_list[[pair$down2]])

  ## 2.2  Stockage avec un nom cohérent                ↓↓↓↓↓
  res_list[[paste0("(Up_common)"  , pair$common_key)]] <- up_genes
  res_list[[paste0("(Down_common)", pair$common_key)]] <- down_genes
}

## ─── 3. Ordre d’affichage désiré ───────────────────────────────────────────────
vector_contrast0 <- c(
  "[SD-SS]",
  common1_key,
  common2_key,
  common3_key,
  "[SS_W78*-SS_WT1]-[SS_E568K-SS_WT2]",
  "[SD_W78*-SD_WT1]-[SD_E568K-SD_WT2]",
  "[SD_W78*-SD_WT1]-[SS_W78*-SS_WT1]",
  "[SD_E568K-SD_WT2]-[SS_E568K-SS_WT2]",
  common4_key,
  common5_key#,
  # "[E568K-WT2]",
  # "[W78*-WT1]",
  # "[WT1_SD-WT1_SS]",
  # "[W78*_SD-W78*_SS]",
  # "[WT2_SD-WT2_SS]",
  # "[E568K_SD-E568K_SS]",
  # "[SS_W78*-SS_WT1]",
  # "[SS_E568K-SS_WT2]",
  # "[SD_W78*-SD_WT1]",
  # "[SD_E568K-SD_WT2]"
  
)

## ─── 4. Étiquette “common” plus lisible ────────────────────────────────────────
#common_label <- stringr::str_replace_all(Contrast0, "_vs_", " inter ")

## ─── 5. Reconstruction de summary_DEG ───────────────────────────────────────────
summary_DEG <- enframe(res_list, name = "name", value = "genes") %>%
  mutate(
    sign = case_when(
      startsWith(name, "(Up_common)"  ) ~ "Up",
      startsWith(name, "(Down_common)") ~ "Down",
      startsWith(name, "(Up)"         ) ~ "Up",
      TRUE                               ~ "Down"
    ),
    Contrast0 = sub(
      "^\\(Up_common\\)|^\\(Down_common\\)|^\\(Up\\)|^\\(Down\\)",
      "",
      name
    ),
    length = lengths(genes)
  ) %>%
  filter(Contrast0 %in% vector_contrast0) %>%
  group_by(Contrast0) %>%
  mutate(
    total   = sum(length),
    percent = round(100 * length / total, 1)
  ) %>%
  ungroup() %>%
  mutate(
    Contrast =
      Contrast0
    ,
    sign     = factor(sign, levels = c("Up", "Down")),
    Contrast = factor(
      Contrast,
      levels = rev(vector_contrast0)
        )
  )%>%
  select(Contrast, sign, length, percent, total)
# ─── 6. Replot ────────────────────────────────────────────────────────────────────
px <- ggplot(summary_DEG, aes(x = Contrast, y = length, fill = sign)) +
  geom_col() +
  geom_text(
    data     = filter(summary_DEG, length >= ifelse(lfc_lim_i == 0, 600, 30)),
    aes(label = paste0(length, "\n(", percent, "%)")),
    position = position_stack(vjust = 0.5),
    color    = "white", size = 3
  ) +
  geom_text(
    data = filter(summary_DEG, length < ifelse(lfc_lim_i == 0, 600, 30)),
    aes(
      label = paste0(length, " (", percent, "%)"),
      color = sign, vjust = ifelse(sign=="Up", -0.5, 1.5), y = total
    ),
    hjust = -0.2, size = 3
  ) +
  coord_flip() +
  theme_minimal() +
  labs(
    x     = "Comparison",
    y     = "Number of genes deregulated",
    title = paste0("Number of genes deregulated (including commons)\n",
                   "(FDR = 0.05 ; LogFC = ", lfc_lim_i, ")")
  ) +
  scale_fill_manual(values = c(Down = "#1d4877", Up = "#ee3e32"), name = "Sign") +
  scale_color_manual(values = c(Down = "#1d4877", Up = "#ee3e32"), guide = FALSE)

print(px)

# export 
fig_export(path = here::here(paste0("report/microarray/plot/DEG/barplot_publication_differential_expression_lfc_",lfc_lim_i)), px, height = 6.5, width = 12, res_i = 600)

save(res_list, file = here::here(paste0("data/microarray/output/contraste_publication_lfc_",lfc_lim_i,".RData")))

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

13.8.2 GO

Code

lfc_lim_i = 0
type_calcul_i = "hypergeometric" #clusterProfiler # hypergeometric # fisher # in reality its exactly the same

load(file = here::here(paste0("data/microarray/output/contraste_publication_lfc_",lfc_lim_i,".RData")))

vector_contrast0 <- c(
  "[SD-SS]",
  common1_key,
  common2_key,
  common3_key,
  "[SS_W78*-SS_WT1]-[SS_E568K-SS_WT2]",
  "[SD_W78*-SD_WT1]-[SD_E568K-SD_WT2]",
  "[SD_W78*-SD_WT1]-[SS_W78*-SS_WT1]",
  "[SD_E568K-SD_WT2]-[SS_E568K-SS_WT2]",
  common4_key,
  common5_key#,
  # "[E568K-WT2]",
  # "[W78*-WT1]",
  # "[WT1_SD-WT1_SS]",
  # "[W78*_SD-W78*_SS]",
  # "[WT2_SD-WT2_SS]",
  # "[E568K_SD-E568K_SS]",
  # "[SS_W78*-SS_WT1]",
  # "[SS_E568K-SS_WT2]",
  # "[SD_W78*-SD_WT1]",
  # "[SD_E568K-SD_WT2]"
  
)

common_keys <- c(common1_key, common2_key, common3_key, common4_key, common5_key)

# 2. Pour chaque contraste brut, on colle les bons préfixes
vector_contrast0_sign <- unlist(
  lapply(vector_contrast0, function(ct) {
    # si c'est un common → on utilise les préfixes "_common"
    if (ct %in% common_keys) {
      paste0(c("(Up_common)", "(Down_common)"), ct)
    } else {
      paste0(c("(Up)", "(Down)"), ct)
    }
  })
)

gene_table <- tibble(
  contrast = rep(names(res_list), lengths(res_list)),
  PsCam = unlist(res_list)
) %>% 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") %>% 
  drop_na(Psat) %>% 
  filter(contrast %in% vector_contrast0_sign)

nrow(gene_table %>% filter(contrast == "(Down)[SD-SS]"))

universe_vector <- unique(gene_table$Psat)

df_GO = GO_on_different_group(functional_roles = "BP",
                      group_info = gene_table,
                      group = "contrast",
                      ID = "Psat",
                      top = 10,
                      universe_i = universe_vector, 
                      type_calcul = type_calcul_i
                                            )

# combined_results = GO_on_different_group(functional_roles = "BP",
#                       group_info = gene_table,
#                       group = "contrast",
#                       ID = "Psat",
#                       top = FALSE,
#                       universe_i = universe_vector, 
#                       type_calcul = type_calcul_i
#                                             )



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, "(?<=\\().*?(?=\\))")
  )
# 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 = paste0("GO terme for all gene deregulated with lfc:",lfc_lim_i)) +
  #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

fig_export(path = paste0("report/microarray/plot/GO/GO_BP_contrast_",type_calcul_i,"_",lfc_lim_i), plot_x = px, height_i = 20, width_i = 12, res = 600)

--- output: html_document editor_options: chunk_output_type: console editor: markdown: wrap: 72 --- ## Exploratory analysis of Raw Data from microarray {#analyse_data_normalized} ## Material and methods <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} # install.packages(here::here("data/microarray/org.Psativum1c.eg.db"), repos = NULL, type = "source") #pkg library(readxl) library(affy) #BiocManager::install("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(ComplexHeatmap) library(knitr) library(kableExtra) library(ggtext) library(ggnewscale) library(org.Psativum1c.eg.db) # if not working install it library(clusterProfiler) library(readxl) library(GO.db) library(AnnotationDbi) # BiocManager::install("AnnotationDbi") library(ggnewscale) # to have two scale_fill # 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 ```{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. And average data and stats corresponding to one of the paired student tests, then corrected between the four biological repeats. [This data has been previously analized here.](A_preprocessed_data_analysis.qmd#analyse_data_pretrated)) **Import raw data (green and red)** ```{r, eval = F} # Function that import all data (take a lot of time) source (here::here("src/function/microarray/import_raw_data_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_raw_data_microarray(path_i, path_info_range_i, sheet_i, begin_end_line) write_csv(execut[["global_df"]], here::here("data/microarray/output/raw_data_microarray_leaf_PeaSulf.csv")) write_csv(execut[["df_info_sample_clean_compile"]], here::here("data/microarray/output/list_microarray_microarray_leaf_PeaSulf.csv")) ``` On microarray chips, it's quite common to find probes that do not belong to the species under study. These probes are typically placed there for various reasons, such as serving as exogenous controls (often called alien probes, spike-ins, or external controls) or to check hybridization specificity. Here is from human or bacteria. In this step i will import it. ```{r, eval = F} path = here::here("data/microarray/Resultats_4plex-POIS-2014_01_230415_modKG.xlsx") path_info_range = here::here("data/microarray/info_range_Resultats_4plex-POIS-2014_01_230415_modKG.xlsx") begin_end_line <- c("41123", "41159") sheet_i = "Complete" external_controls = import_raw_data_microarray(path,path_info_range, sheet_i, begin_end_line) write_csv(external_controls[["global_df"]], here::here("data/microarray/output/external_controls_data_microarray_leaf_PeaSulf.csv")) ``` ```{r, eval =FALSE} external_controls = read_csv(here::here("data/microarray/output/external_controls_data_microarray_leaf_PeaSulf.csv"), show_col_types = FALSE) df_mean = external_controls %>% dplyr::summarise( global_mean = mean(unlist(.[, 9:ncol(.)]), na.rm = TRUE), global_max = max(unlist(.[, 9:ncol(.)]), na.rm = TRUE) ) 2^df_mean ``` ```{r verify normalisation, eval=F} # import raw data red green # data1 = read_csv(here::here("data/microarray/output/raw_data_microarray_leaf_PeaSulf.csv"), show_col_types = FALSE) %>% # as.data.frame() %>% # dplyr::select(c("id_probe", "Mut_SD_5.1_Rep1_Red_1", "WT_SD_1.1_Rep1_Green_1")) # # p_non_norm = data1 %>% # pivot_longer(-"id_probe") %>% # ggplot(aes(x = value, col = name)) + # geom_density() + # labs( # title = "Density of values", # x = "Values", # y = "Density" # ) + # theme_minimal() # # data_longer_test = data1 %>% column_to_rownames("id_probe") # # norm_data <- normalizeBetweenArrays(data_longer_test, method = "cyclicloess") # # p_norm = norm_data %>% # as.data.frame() %>% # rownames_to_column("id_probe") %>% # pivot_longer(-"id_probe") %>% # ggplot(aes(x = value, col = name)) + # geom_density() + # labs( # title = "Density of normalized values", # x = "Values", # y = "Density" # ) + # theme_minimal() # p_non_norm / p_norm ``` Creation of the RGlist file (with absolute value (befor log2)) ```{r, eval = F} # Import raw data red green df_raw_RG_log2 = read_csv(here::here("data/microarray/output/raw_data_microarray_leaf_PeaSulf.csv"), show_col_types = FALSE) %>% as.data.frame() # Creation of the RG list for background stemps #### # Converted into input data (exp of value with log2) df_raw_RG<- df_raw_RG_log2 cols_to_transform <- grep("Red|Green", names(df_raw_RG)) df_raw_RG[cols_to_transform] <- lapply(df_raw_RG[cols_to_transform], function(col) 2^col) # Verification before creation of the RGList red_cols <- grep("Red_\\d+$", colnames(df_raw_RG), value = TRUE) green_cols <- grep("Green_\\d+$", colnames(df_raw_RG), value = TRUE) extract_index <- function(x, color = c("Red","Green")) { pattern <- paste0(".*", color, "_(\\d+)$") sub(pattern, "\\1", x) } numbers_red <- extract_index(red_cols, "Red") numbers_green <- extract_index(green_cols, "Green") common_indices <- intersect(numbers_red, numbers_green) cat("Same number find :", common_indices, "\n") # Parameter of the array n_spots <- nrow(df_raw_RG) n_arrays <- length(common_indices) # 20 = 4 rep * 5 comparison R_matrix <- matrix(NA, nrow = n_spots, ncol = n_arrays) %>% as.data.frame() G_matrix <- matrix(NA, nrow = n_spots, ncol = n_arrays) %>% as.data.frame() k <- 1 for (i in common_indices) { # Trouver la/les colonnes Red_... qui finissent par Red_i # (en théorie, une seule par index i, mais grep renvoie un vecteur - on prend le premier) col_red <- grep(paste0("Red_", i, "$"), red_cols, value = TRUE) col_green <- grep(paste0("Green_", i, "$"), green_cols, value = TRUE) # Remplir la matrice R_matrix[, k] <- df_raw_RG[[col_red]] G_matrix[, k] <- df_raw_RG[[col_green]] # Nommer la colonne pour s’y retrouver (ex : "i" ou le nom exact "Sample1") # Ici on va prendre "Red_xxxx" comme titre, ou "i" tout simplement colnames(R_matrix)[k] <- col_red colnames(G_matrix)[k] <- col_green k <- k + 1 } # Creation of the RG list end #### RG <- new("RGList") RG$R <- R_matrix RG$G <- G_matrix RG$genes <- data.frame(ID = df_raw_RG$id_probe) MA_raw <- normalizeWithinArrays(RG, method = "none") # without normalisation # export RG save(RG,MA_raw, file =here::here("data/microarray/output/raw_data_microarray_leaf_PeaSulf.RData")) ``` ## Standardization of the data cross-chip data <details> <summary style="color: #007bff;"> Bonus: If i apply double normalization ((loess and quartil) not necessary for this experiment)) </summary> In a two-channel microarray, we often apply within-array normalization (such as LOESS) to correct for individual slide biases (systematic variations linked to position on the chip, channel, etc.). Next, a so-called between-array normalization (such as quantile normalization) can be applied to make distributions between replicates comparable. Normalization functions are available in R/Bioconductor (limma package, etc.). See also this paper @bolstad_comparison_2003 ```{r, eval = F} # Load data load(here::here("data/microarray/output/raw_data_microarray_leaf_PeaSulf.RData")) # Raw data visualization #### source (here::here("src/function/microarray/plot_rg_density.R")) # Load function for see density of the differentent normalization p_raw <- plot_rg_density(RG)+labs(title = "plotDensities without normalization") # execut function # Remove or adjust background values (noise) to obtain intensities closer to the true spot intensity #### # Helps model and remove technical noise at lower intensities. # The normexp method is commonly recommended because it avoids negative or zero intensities by adding an offset. In input is absolute value RG_background <- backgroundCorrect(RG, method = "normexp", offset = 20) # 20 test multiple value. Perfect for delet noise rownames(RG_background$R) <- RG_background[["genes"]][["ID"]] rownames(RG_background$G) <- RG_background[["genes"]][["ID"]] # visualisation p_background <- plot_rg_density(RG_background)+labs(title = "plotDensities without normalization but with background correction") log2(min(RG_background$G)) # need to be positive MA_raw <- normalizeWithinArrays(RG, method = "none") # without normalisation plotMA(MA_raw, array = 1) # interesting to se this plot MA_background <- normalizeWithinArrays(RG_background, method = "none") # with normalisation plotMA(MA_background, array = 1) # Normalization loess #### MA.p <-normalizeWithinArrays(RG_background, method = "loess") p_loess <- plot_rg_density(MA.p, normalized = T)+ labs(title = "plotDensities with normalization between the two chanel", subtitle = "Loess normalization") # Normalisation Between Arrays (Aquantile) MA.pAq <- normalizeBetweenArrays(MA.p, method="Aquantile") p_quantile <- plot_rg_density(MA.pAq, normalized = T)+ labs(title = "plotDensities with normalization between the two chanel and between sample", subtitle = "Quantile normalization") # Filter on low value ? #### RG.MA(MA.pAq)$R keep <- rowMeans(RG_background$R) > 1119.909 & rowMeans(RG_background$G) > 1119.909 keep <- rowMeans(RG_background$R) > 1119.909 & rowMeans(RG_background$G) > 1119.909 RG_filtered <- RG[keep, ] # Calculate how many genes were removed n_removed <- n_original - n_filtered # Calculate the percentage of genes removed percent_removed <- (n_removed / n_original) * 100 # Display results cat("Number of genes removed:", n_removed, "\n") cat("Percentage of genes removed:", round(percent_removed, 2), "%\n") p_filtered <- plot_rg_density(RG_filtered)+labs(title = "plotDensities without normalization but with filtration on low intensity") fig_export(here::here("report/microarray/plot/normalization/without_normalization_with_low_intensity"), p_filtered, height_i = 5, width_i = 12, res_i = 300) # Export #### fig_export(here::here("report/microarray/plot/normalization/without_normalization"), p_raw, height_i = 5, width_i = 12, res_i = 300) fig_export(here::here("report/microarray/plot/normalization/without_normalization_with_background_substraction"), p_background, height_i = 5, width_i = 12, res_i = 300) fig_export(here::here("report/microarray/plot/normalization/with_loess_normalization"), p_loess, height_i = 5, width_i = 12, res_i = 300) fig_export(here::here("report/microarray/plot/normalization/with_quantile_normalization"), p_quantile, height_i = 5, width_i = 12, res_i = 300) # all_fig <- p_raw/p_background/p_filtered/p_loess/p_quantile all_fig <- p_raw/p_background/p_loess/p_quantile fig_export(here::here("report/microarray/plot/normalization/all_step_normalization"), all_fig, height_i = 18, width_i = 12, res_i = 300) # RG_filtered save(RG,MA_raw, MA.pAq, file =here::here("data/microarray/output/normalized_data_microarray_leaf_PeaSulf.RData")) # Delet of low intensity #### n_original <- nrow(RG) # Total genes before filtering n_filtered <- nrow(RG_filtered) # Total genes after filtering red_row_means <- rowMeans(RG_background$R) # 2. Combine with Gene IDs into a data frame df_red_means <- data.frame( GeneID = RG_background$genes, # Adjust column name if needed MeanRed = red_row_means ) # 3. Subset for the specific gene of interest df_red_means[df_red_means$ID == "PsCam042688", ] df_red_means[df_red_means$ID == "PsCam036750", ] row_means_red <- rowMeans(RG$R) # 2. Convert to a data frame df_row_means <- data.frame(mean_intensity = row_means_red) # 3. Plot with ggplot2 ggplot(df_row_means, aes(x = log2(mean_intensity))) + geom_histogram(bins = 50, color = "black", fill = "skyblue") + theme_minimal() + labs( title = "Distribution of Row Means (Red Channel)", x = "Mean Intensity (Red Channel)", y = "Count" ) ggplot(df_row_means, aes(x = log2(mean_intensity))) + geom_histogram(bins = 50, color = "black", fill = "skyblue") + scale_x_log10() + theme_minimal() + labs( title = "Distribution of Row Means (Red Channel) [Log Scale]", x = "Log2 (Mean Intensity) (Red Channel) [log10]", y = "Count" ) # RG is your RGList with RG$R and RG$G (non-log) # Suppose you want to keep probes with average raw intensity above 100 in *both* channels # Combination of all plots ``` ![](plot/normalization/all_step_normalization.png) **Expression of PSULT4 with normalized data** ```{r, eval=F} load(file =here::here("data/microarray/output/normalized_data_microarray_leaf_PeaSulf.RData")) list_microarray <-read_csv(here::here("data/microarray/output/list_microarray_microarray_leaf_PeaSulf.csv"), show_col_types = FALSE) ID <- MA.pAq$genes[[1]] A <- MA.pAq$A # Average log intensity M <- MA.pAq$M # Log ratio R <- 2^(A + (M / 2)) # Calculate Red intensities G <- 2^(A - (M / 2)) # Calculate Green intensities R <- as.data.frame(R) G <- as.data.frame(G) rownames(R) <- ID rownames(G) <- ID # Red had good name of column but not green data. i will change that colnames(G) <- colnames(RG$G) RG_norm_filtr <- cbind(R, G) RG_norm_filtr_log2<-log2(RG_norm_filtr) RG_norm_filtr_log2_h <- RG_norm_filtr_log2 %>% rownames_to_column("ID") %>% pivot_longer(-ID, values_to = "value", names_to = "sample_id") %>% left_join(.,list_microarray, by = "sample_id") #export data write_csv(RG_norm_filtr_log2_h, here::here("data/microarray/output/normalized_data_microarray_leaf_PeaSulf.csv")) ``` Comparison of the two chanel ```{r, eval = F} RG_norm_filtr_log2_h = read_csv(here::here("data/microarray/output/normalized_data_microarray_leaf_PeaSulf.csv"), show_col_types = FALSE) p_sult4 <- RG_norm_filtr_log2_h %>% mutate(sulfur_condition = forcats::fct_relevel(sulfur_condition, "SS", "SD"), genotype = forcats::fct_relevel(genotype, "WT1", "W78*", "WT2", "E568K")) %>% filter(ID == "PsCam042688") %>% ggplot(aes(x = genotype, y = value, col= color, fill = sulfur_condition, group = interaction(genotype, sulfur_condition, color))) + geom_boxplot(outlier.shape = NA) + geom_jitter(position = position_jitterdodge(jitter.width = 0.6)) + scale_fill_manual( values = c("SS" = "#FFD74C", "SD" = "gray"), # Adjust based on your `color` levels name = "Sulfure Condition" # Custom title for color legend )+ scale_color_manual( values = c("Red" = "red", "Green" = "forestgreen"), # Adjust based on your `color` levels name = "Channel" # Custom title for color legend )+ theme_minimal()+ theme( # Hide panel borders and remove grid lines panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # Change axis line axis.line = element_line(colour = "black") )+ labs(x = "Genotype", y = "Microarray gene expression of SULTR4") fig_export(here::here("report/microarray/plot/exression_of_genes/SULTR4"), p_sult4, height_i = 4, width_i = 6, res_i = 300) ``` With Stats ```{r, eval = F} RG_norm_filtr_log2_h = read_csv(here::here("data/microarray/output/normalized_data_microarray_leaf_PeaSulf.csv"), show_col_types = FALSE) df_sult4 <- RG_norm_filtr_log2_h %>% as.data.frame() %>% mutate(sulfur_condition = forcats::fct_relevel(sulfur_condition, "SS", "SD"), genotype = forcats::fct_relevel(genotype, "WT1", "W78*", "WT2", "E568K")) %>% filter(ID == "PsCam042688") l_sult4 <- stat_analyse( data=df_sult4 %>% #mutate(climat_condition=paste0(water_condition,"_",heat_condition)) %>% #mutate(condition=factor(condition,levels=c("Sto_WW_OT","Stoc_WS_OT","Sto_WW_HS","Sto_WS_HS","Wen_WW_OT","Wen_WS_OT","Wen_WW_HS","Wen_WS_HS"))) %>% #mutate(climat_condition=factor(climat_condition,levels=c("WW_OT","WS_OT","WW_HS","WS_HS"))) %>% drop_na(value), column_value = "value", category_variables = c("sulfur_condition"), grp_var = "genotype", show_plot = T, outlier_show = F, label_outlier = "sample_id", biologist_stats = T, # Ylab_i = expression(paste("Angle ", widehat(ABC2)," (degrees) for order 1")), Ylab_i = paste("Microarray gene expression of SULTR4"), control_conditions = c("SS"), strip_normale = F, hex_pallet = sulfate_pallet ) p_sult4 <-l_sult4[["plot"]]+labs(color="Treatment",fill="Treatment") fig_export(here::here("report/microarray/plot/exression_of_genes/SULTR4_stats"), p_sult4, height_i = 3.5, width_i = 5, res_i = 300) ``` ::: {layout-ncol="2"} ![](plot/exression_of_genes/SULTR4.png) ![](plot/exression_of_genes/SULTR4_stats.png) ::: </details> ### Verification of the importance of data standardization ```{r, eval = F} # import no normalize data load(here::here("data/microarray/output/raw_data_microarray_leaf_PeaSulf.RData")) # creation of the quantile normalization dataset MA_aquantil <- normalizeBetweenArrays(MA_raw, method="Aquantile") RG_aquantil <- RG.MA(MA_aquantil) colnames(RG_aquantil$G) = colnames(RG$G) # creation of df global with all data RG_global <- bind_cols(RG$genes, RG$R, RG$G) %>% column_to_rownames("ID") RG_aquantil_global <- bind_cols(RG_aquantil$genes, RG_aquantil$R, RG_aquantil$G) %>% column_to_rownames("ID") # transfor into log2 RG_global_log2 <- RG_global %>% mutate(across(everything(), log2)) RG_aquantil_global_log2 <- RG_aquantil_global %>% mutate(across(everything(), log2)) # mean for each column mean_add_info <- function(data) { data %>% dplyr::summarise(across(everything(), mean, na.rm = TRUE)) %>% pivot_longer(cols = everything(), names_to = "sample_name", values_to = "Mean") %>% mutate(sulfur_condition = case_when( grepl("SS", sample_name) ~ "SS", grepl("SD", sample_name) ~ "SD", TRUE ~ "Other" )) %>% mutate(color = case_when( grepl("Red", sample_name) ~ "Red", grepl("Green", sample_name) ~ "Green", TRUE ~ "Other" )) %>% mutate(num_combination = str_extract(sample_name, "\\d+$")) %>% mutate( sample_num = str_extract(sample_name, "\\d+\\.\\d+"), sample_num = str_replace_all(sample_num, c("\\_" = "\\.")), simplify_condition = str_extract(sample_name, "WT|Mut"), genotype = case_when( sample_num %in% c("1.1", "2.2", "3.3", "4.4", "17.1", "18.2", "19.3", "20.4") ~ "WT2", sample_num %in% c("5.1", "6.2", "7.3", "8.4", "21.1", "22.2", "23.3", "24.4") ~ "E568K", sample_num %in% c("9.5", "10.6", "11.7", "12.8", "25.5", "26.6", "27.7", "28.8") ~ "WT1", sample_num %in% c("13.5", "14.6", "15.7", "16.8", "29.5", "30.6", "31.7", "32.8") ~ "W78*", TRUE ~ sample_num ), simplify_condition = paste0(genotype, "_", sulfur_condition) ) %>% mutate(simplify_condition = forcats::fct_relevel(simplify_condition, "WT1_SS", "W78*_SS","WT1_SD", "W78*_SD", "WT2_SS", "E568K_SS", "WT2_SD", "E568K_SD")) } RG_global_mean <- mean_add_info(RG_global_log2) RG_aquantil_global_mean <- mean_add_info(RG_aquantil_global_log2) # Visualisation of the different boxplot (for mean density) p_raw <- ggplot(RG_global_mean, aes(x = simplify_condition, y = Mean, fill = color, colour = color)) + geom_point(position = position_jitterdodge()) + geom_boxplot(alpha = 0.3, outlier.shape = NA, col = "black") + theme_minimal() + labs( title = "Without normalisation", y = "Mean of log2", fill = "Channel", colour = "Channel" ) + scale_color_manual(values = c("Red" = "red", "Green" = "forestgreen")) + scale_fill_manual(values = c("Red" = "red", "Green" = "forestgreen")) p_norm <- ggplot(RG_aquantil_global_mean, aes(x = simplify_condition, y = Mean, fill = color, colour = color)) + geom_point(position = position_jitterdodge()) + geom_boxplot(alpha = 0.3, outlier.shape = NA, col = "black") + theme_minimal() + labs( title = "With normalisation", y = "Mean of log2", fill = "Channel", colour = "Channel" ) + scale_color_manual(values = c("Red" = "red", "Green" = "forestgreen")) + scale_fill_manual(values = c("Red" = "red", "Green" = "forestgreen")) p_comparison <- p_raw + p_norm + plot_layout(guides = "collect") fig_export(here::here(paste0("report/microarray/plot/normalization/test_normalisation/boxplot_with_witout_normalization")), p_comparison <- p_raw + p_norm + plot_layout(guides = "collect"), height_i = 6, width_i = 16, res_i = 600) # For PCA ## transform into matrix make_matrix_add_info <- function(data) { data %>% t() %>% as.data.frame() %>% rownames_to_column("sample_id") %>% left_join(.,read_csv(here::here("data/microarray/output/list_microarray_microarray_leaf_PeaSulf.csv"), show_col_types = FALSE) %>% mutate(sample_id_simplify = paste(sep = "_", simplify_condition, genotype, num_combination), condition = paste(sep = "_", simplify_condition, genotype)) %>% dplyr::select(sample_id, condition, sulfur_condition, genotype), by = "sample_id") %>% mutate( splitted = str_split_fixed(sample_id, "_", 6) ) %>% mutate( sample_id = paste(splitted[, 3], splitted[, 5], splitted[, 6], sep = "_"), genotype_simplify = splitted[,1] ) %>% dplyr::select(-splitted) %>% column_to_rownames("sample_id") %>% mutate(across(c(condition, genotype, sulfur_condition,genotype_simplify), as.factor)) %>% mutate(condition = fct_relevel(condition, "WT_SS_WT1", "WT_SD_WT1", "Mut_SS_W78*", "Mut_SD_W78*", "WT_SS_WT2", "WT_SD_WT2", "Mut_SS_E568K", "Mut_SD_E568K"), genotype = fct_relevel(genotype, "WT1", "W78*", "WT2", "E568K"), sulfur_condition = fct_relevel(sulfur_condition, "SS", "SD"), genotype_simplify = fct_relevel(genotype_simplify, "WT", "Mut") ) } RG_log2_t_info <- make_matrix_add_info(RG_global_log2) RG_aquantil_log2_t_info <- make_matrix_add_info(RG_aquantil_global_log2) ## transform into pca for sulfur condition RG_log2_t <- RG_log2_t_info %>% select(-c("condition", "genotype", "sulfur_condition","genotype_simplify")) RG_aquantil_log2_t <- RG_aquantil_log2_t_info %>% select(-c("condition", "genotype", "sulfur_condition","genotype_simplify")) res_pca <- PCA(RG_log2_t, quali.sup = 1:4, graph = F) res_pca_aquantil <- PCA(RG_aquantil_log2_t, quali.sup = 1:4, graph = F) fviz_eig(res_pca, addlabels = TRUE, ylim = c(0, 50)) fviz_eig(res_pca_aquantil, addlabels = TRUE, ylim = c(0, 50)) var <- get_pca_var(res_pca) var_aquantil <- get_pca_var(res_pca_aquantil) # By sulfur condition #### PCA_indiv<-fviz_pca_ind(res_pca, # geom.ind = "point", # Shows points only (not text) pointshape = 20,pointsize = 2, col.ind = RG_log2_t_info %>% pull(sulfur_condition), # colorer by groups palett=sulfate_pallet, addEllipses = T, # Concentration ellipses legend.title = "Sulfur condition", ellipse.level = 0.95, ellipse.type = c("norm"), )+ ggtitle("Visualizing Individual PCA wthout normalization \nby sulfur condition") PCA_indiv PCA_indiv_aquantil<-fviz_pca_ind(res_pca_aquantil, # geom.ind = "point", # Shows points only (not text) pointshape = 20,pointsize = 2, col.ind = RG_log2_t_info %>% pull(sulfur_condition), # colorer by groups palett=sulfate_pallet, addEllipses = T, # Concentration ellipses legend.title = "Sulfur condition", ellipse.level = 0.95, ellipse.type = c("norm"), )+ ggtitle("Visualizing Individual PCA normalize with aquantil \nby sulfur condition") PCA_indiv_aquantil p_comparison <- PCA_indiv+PCA_indiv_aquantil+plot_layout(guides = "collect") fig_export(here::here(paste0("report/microarray/plot/PCA/PCA_sulfur_condition_with_without_normalization")), p_comparison, height_i = 10, width_i = 22, res_i = 600) # by genotype for each condition #### make_PCA_genotype <- function(df_info, norm = T){ df <- df_info %>% select(-c("condition", "genotype", "sulfur_condition","genotype_simplify")) res_pca <- PCA(df, quali.sup = 1:4, graph = F) #fviz_eig(res_pca, addlabels = TRUE, ylim = c(0, 50)) var <- get_pca_var(res_pca) PCA_indiv<-fviz_pca_ind(res_pca, # geom.ind = "point", # Shows points only (not text) #pointshape = 20,pointsize = 2, col.ind = df_info %>% pull(genotype), # colorer by groups shape.ind = df_info$genotype_simplify, # Forme selon le génotype palett= c("#003049", "#780000", "#7FACC7", "#EC323E"), addEllipses = F, # Concentration ellipses legend.title = "Treatment", ellipse.level = 0.95, ellipse.type = c("norm") )+ ggtitle(paste0("Visualizing Individual PCA ", ifelse(norm ==T, "with", "without")," normalization \nby genotype for ",df_info$sulfur_condition[1])) return(PCA_indiv) } p_SS_raw = RG_log2_t_info %>% filter(sulfur_condition =="SS") %>% make_PCA_genotype(., norm=F) p_SD_raw = RG_log2_t_info %>% filter(sulfur_condition =="SD") %>% make_PCA_genotype(., norm=F) p_SS_norm = RG_aquantil_log2_t_info %>% filter(sulfur_condition =="SS") %>% make_PCA_genotype(., norm=T) p_SD_norm = RG_aquantil_log2_t_info %>% filter(sulfur_condition =="SD") %>% make_PCA_genotype(., norm=T) p_compil <- p_SS_raw+p_SS_norm+p_SD_raw+p_SD_norm+plot_layout(guides = "collect", ncol = 2) fig_export(here::here(paste0("report/microarray/plot/PCA/PCA_raw_norm_genotype")), p_compil, height_i = 10, width_i = 14, res_i = 600) # export for Marie-Laure matrix_info <- RG_aquantil_global_mean %>% dplyr::select(-Mean) save(matrix_info, RG, MA_aquantil,file = here::here("data/microarray/output/targets_design.RData")) ``` ![](plot/normalization/test_normalisation/boxplot_with_witout_normalization.png) ![](plot/PCA/PCA_sulfur_condition_with_without_normalization.png) ![](plot/PCA/PCA_raw_norm_genotype.png) ## Differential expression of normalized data <details> <summary style="color: #007bff;"> Bonus: If i apply a double normalization ( not necessary for this experiment), here the vulcanoplot with normalization) </summary> **Of normalized data** ```{r, eval = F} # /!\ vey slow source (here::here("src/function/microarray/vulcanoplot_microarray_RG.R")) for(i in 1:5){ for(normalized_i in c(T, F)){ # create ggplot p = vulcanoplot_microarray_RG(normalized = normalized_i, list_comparison_num_i = i) fig_export(here::here(paste0("report/microarray/plot/volcanoplot/volcano_interactive_RG_",ifelse(normalized_i ==T, "normalized_","raw_"),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_RG_",ifelse(normalized_i ==T, "normalized_","raw_"),i,".html"))) } } ``` ```{r} load(file = here::here("data/microarray/output/ratio_stats_microarray_leaf_PeaSulf.RData")) # only to show name of vulcanoplot ``` ::: panel-tabset ## From normalisez data ::: panel-tabset ### `r gsub("_RepBio_\\d+", "", names(list_ratio_stats)[2])` normalized For an interactive version : [Click here](plot/volcanoplot/volcano_interactive_RG_normalized_2.html) ![](plot/volcanoplot/volcano_interactive_RG_normalized_2.png) ### `r gsub("_RepBio_\\d+", "", names(list_ratio_stats)[1])` normalized For an interactive version : [Click here](plot/volcanoplot/volcano_interactive_RG_normalized_1.html) ![](plot/volcanoplot/volcano_interactive_RG_normalized_1.png) ### `r gsub("_RepBio_\\d+", "", names(list_ratio_stats)[3])` normalized For an interactive version : [Click here](plot/volcanoplot/volcano_interactive_RG_normalized_3.html) ![](plot/volcanoplot/volcano_interactive_RG_normalized_3.png) ### `r gsub("_RepBio_\\d+", "", names(list_ratio_stats)[4])` normalized For an interactive version : [Click here](plot/volcanoplot/volcano_interactive_RG_normalized_4.html) ![](plot/volcanoplot/volcano_interactive_RG_normalized_4.png) ### `r gsub("_RepBio_\\d+", "", names(list_ratio_stats)[5])` normalized For an interactive version : [Click here](plot/volcanoplot/volcano_interactive_RG_normalized_5.html) ![](plot/volcanoplot/volcano_interactive_RG_normalized_5.png) ::: ## From Raw data ::: panel-tabset ### `r gsub("_RepBio_\\d+", "", names(list_ratio_stats)[2])` For an interactive version : [Click here](plot/volcanoplot/volcano_interactive_RG_raw_2.html) ![](plot/volcanoplot/volcano_interactive_RG_raw_2.png) ### `r gsub("_RepBio_\\d+", "", names(list_ratio_stats)[1])` For an interactive version : [Click here](plot/volcanoplot/volcano_interactive_RG_raw_1.html) ![](plot/volcanoplot/volcano_interactive_RG_raw_1.png) ### `r gsub("_RepBio_\\d+", "", names(list_ratio_stats)[3])` For an interactive version : [Click here](plot/volcanoplot/volcano_interactive_RG_raw_3.html) ![](plot/volcanoplot/volcano_interactive_RG_raw_3.png) ### `r gsub("_RepBio_\\d+", "", names(list_ratio_stats)[4])` For an interactive version : [Click here](plot/volcanoplot/volcano_interactive_RG_raw_4.html) ![](plot/volcanoplot/volcano_interactive_RG_raw_4.png) ### `r gsub("_RepBio_\\d+", "", names(list_ratio_stats)[5])` For an interactive version : [Click here](plot/volcanoplot/volcano_interactive_RG_raw_5.html) ![](plot/volcanoplot/volcano_interactive_RG_raw_5.png) ::: ::: </details> ### Using limma contrasts (method recommended by Marie Laure Martin-Magniette) In multi-factor or multi-group microarray analyses, we often want to compare specific conditions such as different genotypes, treatments, or their combinations. Why use contrasts? A linear model created with lmFit describes the contribution of each factor in the experiment, but it doesn’t automatically test the specific comparisons we are interested in. A contrast is a linear combination of the model’s coefficients that defines a particular comparison, for example, comparing genotype A to genotype B, or comparing interaction effects. Using a contrast matrix allows us to explicitly define and test these comparisons, ensuring that limma computes the appropriate log-fold changes and p-values. How it works: First, we create a design matrix using model.matrix(...), which captures the main effects and possibly the interactions between experimental factors. Then we define a contrast matrix using a custom function (like Contrasts) or limma’s built-in tools. Each row of the contrast matrix corresponds to a specific biological question or comparison. After this, the lmFit(...) function estimates gene-wise model coefficients. We apply our contrasts using contrasts.fit(...), which re-expresses the model in terms of the contrasts of interest. The eBayes(...) function then moderates the variance estimates and calculates statistical significance (p-values and adjusted p-values). In the rest of the code, we use p.adjust(...) to control for multiple testing (FDR), identify significantly up- or down-regulated genes based on a fold-change threshold, and generate summary outputs such as barplots and volcano plots. This provides a clear overview of differential expression across the comparisons that matter most in the study. In summary, contrasts are essential in differential analysis using limma, as they allow us to tailor the comparisons to the specific hypotheses of the experiment, especially when multiple factors and interactions are involved. ```{r, eval = F} #parameter lfc_lim_i = 0 pval_i =0.05 load(file = here::here("data/microarray/output/targets_design.RData")) source(here::here("src/function/microarray/Contrasts.R")) # creation of a complete dataframe RG_aquantil <- limma::RG.MA(MA_aquantil) colnames(RG_aquantil$G) <- colnames(RG$G) df_RG_aquantil_log2 <- cbind(RG_aquantil$R, RG_aquantil$G, RG_aquantil$genes) %>% column_to_rownames("ID") %>% mutate(across(everything(), log2)) save(df_RG_aquantil_log2, file = here::here("data/microarray/output/df_RG_aquantil_log2.RData")) Target <- matrix_info %>% dplyr::select(sulfur_condition, genotype, num_combination) Target=as.data.frame(Target) Target$sulfur_condition = as.factor(Target$sulfur_condition) Target$genotype = as.factor(Target$genotype) model <- model.matrix(as.formula(paste("~", paste(colnames(Target)[1:2], collapse = " + "), paste(colnames(Target[,1:2]),collapse = ":"), sep = " + ")), data=Target) contrast.matrix = Contrasts(model,Target,FALSE,TRUE) mycontrast1 <- c(0, 0, +1, -1, +1, +1, -1, +1) mycontrast2 <- c(0, 0, +1, -1, +1, 0, 0, 0) contrast.matrix # votre matrice après Contrasts(...) # On l'augmente en ajoutant nos 2 lignes : contrast.matrix <- rbind( contrast.matrix, mycontrast1, mycontrast2 ) # On renomme les deux nouvelles lignes : n <- nrow(contrast.matrix) rownames(contrast.matrix)[(n-1):n] <- c( "[SS_W78*-SS_WT1]-[SS_E568K-SS_WT2]", "[SD_W78*-SD_WT1]-[SD_E568K-SD_WT2]" ) contrast.matrix fit <- lmFit(df_RG_aquantil_log2, model) #contrast.matrix <- makeContrasts(SpecialDiet-Control, SpecialDietDrug1-SpecialDiet, SpecialDietDrug1-Control, levels=design) fit2 <- contrasts.fit(fit, t(contrast.matrix)) fit2 <- eBayes(fit2) # topTable(fit2, coef=1, adjust="fdr", sort.by="B", number=10) # creation of a plot to count the number of difference between eatch contraste p_adj <- apply(fit2$p.value, 2, p.adjust, method = "BH") is_significant <- p_adj < 0.05 # Matrice TRUE/FALSE pour FDR < 0.05 is_upregulated <- is_significant & (fit2$coefficients > lfc_lim_i) # logFC > 0 is_downregulated <- is_significant & (fit2$coefficients < -lfc_lim_i) # logFC < 0 vector_with_interest <- c("[SD-SS]","[WT1-WT2]","[E568K-WT2]", "[W78*-WT1]", "[WT1_SD-WT1_SS]","[W78*_SD-W78*_SS]", "[WT2_SD-WT2_SS]", "[E568K_SD-E568K_SS]", "[SS_W78*-SS_WT1]", "[SS_E568K-SS_WT2]", "[SD_W78*-SD_WT1]", "[SD_E568K-SD_WT2]", "[SD_W78*-SD_WT1]-[SS_W78*-SS_WT1]", "[SD_E568K-SD_WT2]-[SS_E568K-SS_WT2]", "[SS_W78*-SS_WT1]-[SS_E568K-SS_WT2]", "[SD_W78*-SD_WT1]-[SD_E568K-SD_WT2]" ) summary_DEG <- tibble( Contrast = colnames(fit2$p.value), Upregulated = colSums(is_upregulated), Downregulated = colSums(is_downregulated) ) %>% pivot_longer(cols = -c("Contrast"), names_to = "sign", values_to = "length") %>% mutate(sign = recode(sign, "Upregulated" = "Up", "Downregulated" = "Down")) %>% dplyr::group_by(Contrast) %>% mutate( total = sum(length), percent = round(100 * length / total) ) %>% ungroup() %>% filter(Contrast %in% c(vector_with_interest)) %>% mutate(sign = fct_relevel(sign, c("Up", "Down")), Contrast = fct_relevel(Contrast, rev(vector_with_interest))) # Show results px <- ggplot(summary_DEG, aes(x = Contrast, y = length, fill = sign)) + geom_col() + # 1) Texte pour les barres "assez grandes" (>= 5) geom_text( data = subset(summary_DEG, length >= ifelse(lfc_lim_i == 0, 600,30)), aes( label = paste0(length, "\n(", percent, "%)") ), position = position_stack(vjust = 0.5), color = "white", # texte en blanc size = 3, # hjust = 0.5 (par défaut) - pas indispensable ) + # 2) Texte pour les barres "trop petites" (< 5) geom_text( data = subset(summary_DEG, length < ifelse(lfc_lim_i == 0, 600,30)), aes( label = paste0(length, " (", percent, "%)"), color = sign, # couleur du texte = même code couleur que "sign" vjust = ifelse(sign == "Up", -0.5, 1.5), y=total ), #position = position_stack(vjust = 0.5), hjust= -0.2, # On décale à droite (hors de la barre). En coord_flip(), # hjust < 0 place le texte davantage à droite. #hjust = -0.2, size = 3 ) + theme_minimal() + labs( x = "Comparison", y = "Number of genes deregulated", title = paste0("Number of genes deregulated for each contrast (FDR=0.05 ; LogFC=", lfc_lim_i, ")") ) + coord_flip() + # Couleurs de remplissage (barres) scale_fill_manual( values = c('Down' = "#1d4877", 'Up' = "#ee3e32"), name = "Sign" ) + # Couleurs du texte (identique aux barres), sans nouvelle légende scale_color_manual( values = c('Down' = "#1d4877", 'Up' = "#ee3e32"), guide = FALSE ) ; px fig_export(path = here::here(paste0("report/microarray/plot/DEG/barplot_differential_expression_lfc_",lfc_lim_i)), px, height = 6.5, width = 12, res_i = 600) # export eatch contrast res_list <- list() # Loop through each contrast in your vector of interest for (contrast in vector_with_interest) { # Construct vector names, e.g., "[SD-SS](Up)" and "[SD-SS](Down)" up_name <- paste0("(Up)", contrast) down_name <- paste0("(Down)", contrast) # Extract the rownames of all Upregulated genes for this contrast genes_up <- rownames(is_upregulated)[is_upregulated[, contrast]] # Extract the rownames of all Downregulated genes for this contrast genes_down <- rownames(is_downregulated)[is_downregulated[, contrast]] # Store them into the result list with meaningful names res_list[[up_name]] <- genes_up res_list[[down_name]] <- genes_down } save(res_list, file = here::here(paste0("data/microarray/output/contraste_lfc_",lfc_lim_i,".RData"))) # make vulcanoplot for each contraste for (contrast in vector_with_interest) { cat(paste0(contrast, "\n")) volcano_data <- tibble( Contrast = contrast, Gene = rownames(fit2$coefficients), logFC = fit2$coefficients[, contrast], pvalue = fit2$p.value[, contrast], adj_pvalue = p_adj[, contrast] ) %>% mutate( Significance = case_when( adj_pvalue < 0.05 & logFC > lfc_lim_i ~ "Upregulated", adj_pvalue < 0.05 & logFC < -lfc_lim_i ~ "Downregulated", TRUE ~ "Not significant" ) ) %>% dplyr::rename(PsCam = Gene) %>% # add info of genes 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") %>% left_join(.,read_csv(here::here("data/microarray/output/raw_data_microarray_leaf_PeaSulf.csv"), show_col_types = FALSE) %>% as.data.frame() %>% dplyr::select(c("id_probe", "CAT")) %>% dplyr::rename(PsCam = id_probe), by = "PsCam" ) num_upregulated <- sum(volcano_data$Significance == "Upregulated") num_downregulated <- sum(volcano_data$Significance == "Downregulated") # Création du titre avec couleurs HTML plot_title <- paste0( "Volcano Plot - ", contrast, "<br>", " | <span style='color:#B02428;'>Upregulated: ", num_upregulated, "</span> | ", "<span style='color:#6697EA;'>Downregulated: ", num_downregulated, "</span>" ) p1 <- ggplot(volcano_data, aes(x = logFC, y = -log10(adj_pvalue), color = Significance)) + geom_point(alpha = 0.7) + scale_color_manual(values = c("Upregulated" = "#B02428", "Downregulated" = "#6697EA", "Not significant" = "grey")) + geom_vline(xintercept = c(-lfc_lim_i, lfc_lim_i), linetype = "dashed") + geom_hline(yintercept = -log10(pval_i), linetype = "dashed") + geom_text_repel( data = volcano_data %>% filter(Significance %in% c("Upregulated", "Downregulated")) %>% slice_max(abs(logFC),n = 20), aes(label = sub("^symbols:", "", CAT), color = Significance), size = 3, max.overlaps = 0 ) + theme_minimal() + labs( title = plot_title, x = expression(log[2]~Fold~Change), y = expression(-log[10](p-value)) ) + theme( plot.title = element_markdown() # Permet d'afficher du texte formaté en HTML ) fig_export(path = here::here(paste0("report/microarray/plot/volcanoplot/volacano_contrast_", gsub("\\*", "", contrast), "_lfc_",lfc_lim_i)), p1, height = 6, width = 12, res_i = 300) # for interactive p1_interactive <- ggplot(volcano_data, aes(x = logFC, y = -log10(adj_pvalue), color = Significance)) + geom_point(aes(text = paste("PsCam:", PsCam, "<br>Psat:", Psat, "<br>CAT:", CAT, "<br>logFC:", logFC, "<br>FDR:", adj_pvalue)), size = 1, alpha = .6) + scale_color_manual( values = c( "Not significant" = "gray", "Downregulated" = "#6697EA", "Upregulated" = "#B02428" ) ) + theme_minimal() + labs( title = plot_title, x = "log2 Fold Change", # Replace expression with plain text y = "-log10(p-value)" ) + theme( plot.title = element_markdown() # Permet d'afficher du texte formaté en HTML ) ggplotly(p1_interactive) %>% saveWidget(here::here(paste0("report/microarray/plot/volcanoplot/volcano_interactive_contraste_",gsub("\\*", "", contrast),"_lfc_",lfc_lim_i,".html"))) } ``` ::: panel-tabset ## With log fold change of 0 ![](plot/DEG/barplot_differential_expression_lfc_0.png) ## With log fold change of 1 ![](plot/DEG/barplot_differential_expression_lfc_1.png) ::: ### Volcanoplot ```{r, echo=FALSE, results = "asis"} # On liste uniquement les fichiers débutant par "volacano_contrast_" et finissant par ".png" image_files <- list.files( path = "plot/volcanoplot", pattern = "^volacano_contrast_.*\\.png$", full.names = TRUE ) # Fonction pour extraire contraste et LFC via regex extract_info <- function(filepath) { # Nom de fichier sans extension fname_noext <- tools::file_path_sans_ext(basename(filepath)) # Ex : "volacano_contrast_[WT1_SD-WT1_SS]_lfc_1" # On veut : contrast = "[WT1_SD-WT1_SS]", lfc = "1" # Première regex : capture tout ce qu’il y a entre "volacano_contrast_" et "_lfc_" # Deuxième capture : la valeur numérique (0 ou 1) après "_lfc_" # ^ : début de chaîne # volacano_contrast_ : littéral # (.*) : capture tout, de manière gourmande # _lfc_ : littéral # ([0-9]+) : capture 1 ou plusieurs chiffres # $ : fin de chaîne pattern <- "^volacano_contrast_(.*)_lfc_([0-9]+)$" contrast <- sub(pattern, "\\1", fname_noext) lfc_val <- sub(pattern, "\\2", fname_noext) # S’il n’y a pas correspondance, on peut gérer un cas « vide » if (contrast == fname_noext) { # cas où la regex ne correspond pas contrast <- NA lfc_val <- NA } return(list(contrast = contrast, lfc = lfc_val)) } # Extraire les infos pour chaque fichier infos <- lapply(image_files, extract_info) # Construire un data.frame df <- data.frame( file = image_files, contrast = sapply(infos, `[[`, "contrast"), lfc = sapply(infos, `[[`, "lfc"), stringsAsFactors = FALSE ) # Filtrer uniquement les lignes valides (i.e. pas de NA, et lfc dans 0 ou 1) df <- subset(df, !is.na(contrast) & !is.na(lfc) & lfc %in% c("0","1")) # Génération du code en panneaux tabset tabset_code <- ":::: {.panel-tabset}\n\n" # On boucle sur lfc = 0, puis lfc = 1 (ou récupérez directement unique(df$lfc)) for(lfc_val in c("0", "1")) { sub_df <- subset(df, lfc == lfc_val) if(nrow(sub_df) == 0) next # Panel principal : with log fold change 0 ou 1 tabset_code <- paste0( tabset_code, "## with log fold change ", lfc_val, "\n\n", ":::: {.panel-tabset}\n\n" ) # Boucle sur les contrastes contrasts <- unique(sub_df$contrast) for(con in contrasts){ tmp <- subset(sub_df, contrast == con) tabset_code <- paste0( tabset_code, "### ", con, "\n\n" ) # Si plusieurs images pour le même contraste, on boucle for(i in seq_len(nrow(tmp))) { # Chemin vers la version interactive # Ex : volcano_interactive_contraste_[WT1_SD-WT1_SS]_lfc_1.html interactive_path <- gsub( pattern = "^volacano_contrast_", replacement = "volcano_interactive_contraste_", x = basename(tmp$file[i]) ) interactive_path <- sub("\\.png$", ".html", interactive_path) # Lien vers la version interactive tabset_code <- paste0( tabset_code, "For an interactive version : [Click here](plot/volcanoplot/", interactive_path, ")\n\n" ) # Image statique tabset_code <- paste0( tabset_code, "![](", tmp$file[i], ")\n\n" ) } } tabset_code <- paste0(tabset_code, "::::\n\n") } tabset_code <- paste0(tabset_code, "::::\n") cat(tabset_code) ``` ## Gene onthologie on contraste - verification and simplification of the differents fonctions - add treemap plot ```{r, eval = F} # parameter lfc_lim_i = 0 type_calcul_i = "hypergeometric" #clusterProfiler # hypergeometric # fisher # in reality its exactly the same load(file = here::here(paste0("data/microarray/output/contraste_lfc_",lfc_lim_i,".RData"))) gene_table <- tibble( contrast = rep(names(res_list), lengths(res_list)), PsCam = unlist(res_list) ) %>% 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") %>% drop_na(Psat) nrow(gene_table %>% filter(contrast == "(Down)[SD-SS]")) universe_vector <- unique(gene_table$Psat) df_GO = GO_on_different_group(functional_roles = "BP", group_info = gene_table, group = "contrast", ID = "Psat", top = 10, universe_i = universe_vector, type_calcul = type_calcul_i ) # combined_results = GO_on_different_group(functional_roles = "BP", # group_info = gene_table, # group = "contrast", # ID = "Psat", # top = FALSE, # universe_i = universe_vector, # type_calcul = type_calcul_i # ) 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, "(?<=\$).*?(?=\$)") ) # 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 = paste0("GO terme for all gene deregulated with lfc:",lfc_lim_i)) + #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 fig_export(path = paste0("report/microarray/plot/GO/GO_BP_contrast_",type_calcul_i,"_",lfc_lim_i), plot_x = px, height_i = 20, width_i = 12, res = 600) ``` Test ```{r, eval = F} ################## test WT1 vs WT2 in contraste # test WT1 vs WT2 realy different also for sulfure df_test = test %>% filter(Description_GO %in% c("sulfate transport", "sulfur compound transport"), group == "(Down)[WT1-WT2]") %>% pull(geneID) %>% str_split("/") %>%unlist() %>%unique() %>% tibble(Psat = .) %>% left_join(., 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 = "Psat")) %>% left_join(., left_join(.,read_excel(here::here("data/microarray/Resultats_4plex-POIS-2014_01_230415_modKG.xlsx"), sheet = "Complete", col_names = T, skip = 11)%>% dplyr::rename(PsCam = "N°") %>% dplyr::select(PsCam, CAT ), by = "PsCam") ) load(file = here::here("data/microarray/output/targets_design.RData")) source(here::here("src/function/stat_function/stat_analysis_main.R")) # for make plot # create data frame RG_aquantil <- limma::RG.MA(MA_aquantil) colnames(RG_aquantil$G) <- colnames(RG$G) df_RG_aquantil_log2 <- cbind(RG_aquantil$R, RG_aquantil$G, RG_aquantil$genes) %>% column_to_rownames("ID") %>% mutate(across(everything(), log2)) %>% rownames_to_column("PsCam") %>% pivot_longer(-PsCam, names_to = "sample_name") %>% dplyr::mutate( sample_num = str_extract(sample_name, "(?<=_[A-Z]{2}_)[0-9\\.]+"), sample_num = str_replace_all(sample_num, c("\\_" = "\\.")), simplify_condition = str_extract(sample_name, "WT|Mut"), simplify_condition = str_extract(sample_name, "SS|SD"), sulfur_condition = str_sub(simplify_condition, -2, -1), genotype = ifelse(sample_num %in% c(1.1, 2.2, 3.3, 4.4, 17.1, 18.2, 19.3, 20.4), "WT2", sample_num),# Add genotype genotype = ifelse(sample_num %in% c(5.1, 6.2, 7.3, 8.4, 21.1, 22.2, 23.3, 24.4), "E568K", genotype),# Add genotype genotype = ifelse(sample_num %in% c(9.5, 10.6, 11.7, 12.8, 25.5, 26.6, 27.7, 28.8), "WT1", genotype),# Add genotype genotype = ifelse(sample_num %in% c(13.5, 14.6, 15.7, 16.8, 29.5, 30.6, 31.7, 32.8), "W78*", genotype),# Add genotype sample_id = paste(sep = "_", simplify_condition, genotype) ) df_prepared <- df_test_value %>% as.data.frame() %>% mutate( sulfur_condition = fct_relevel(sulfur_condition, "SS", "SD"), genotype = fct_relevel(genotype, "WT1", "W78*", "WT2", "E568K") ) %>% drop_na(value) # 2. Extraire les niveaux uniques de "psat_pscam" (plus rapide que as.factor(...) puis levels(...)) vector_sulfur <- unique(df_prepared$psat_pscam) # 3. Boucler sur chaque valeur de psat_pscam pour générer la liste de graphiques plots <- lapply(vector_sulfur, function(psat) { # Afficher le nom en vert (si cat_col est une fonction perso pour colorer la console) cat_col(paste0(psat, "\n"), color = "green") # Sous-ensemble des données correspondant à psat df_x <- df_prepared %>% filter(psat_pscam == psat) # Analyse statistique l_sult4 <- stat_analyse( data = df_x, column_value = "value", category_variables = c("sulfur_condition"), grp_var = "genotype", show_plot = TRUE, outlier_show = FALSE, label_outlier = "sample_name", biologist_stats = TRUE, Ylab_i = paste("Microarray gene expression of \n", df_x$CAT[1]), control_conditions = c("SS"), strip_normale = FALSE, hex_pallet = sulfate_pallet ) # Ajout des labels et récupération du plot p_plot <- l_sult4[["plot"]] + labs(color = "Treatment", fill = "Treatment") # On retourne l’objet graphique return(p_plot) }) final_plot <- wrap_plots(plots, ncol = 7) + plot_layout(guides = "collect") + plot_annotation() & theme(legend.position = 'bottom') fig_export(here::here("report/microarray/plot/variation_gene_sulfur"), final_plot, height_i = 21, width_i = 29.7, res_i = 600) ``` :::: panel-tabset # With clusterprofiler ::: panel-tabset ## With log fold change of 0 ![](plot/GO/GO_BP_contrast_clusterProfiler_0.png) ## With log fold change of 1 ![](plot/GO/GO_BP_contrast_clusterProfiler_1.png) ::: # With fisher test or hypergeometric test ::: panel-tabset ## With log fold change of 0 ![](plot/GO/GO_BP_contrast_hypergeometric_0.png) ## With log fold change of 1 ![](plot/GO/GO_BP_contrast_hypergeometric_1.png) ::: :::: ## UpsetPlot ```{r, eval = F} # Define the logFC and p-value thresholds lfc_lim_i <- 0 pval_i = 0.05 load(file = here::here("data/microarray/output/targets_design.RData")) source(here::here("src/function/microarray/Contrasts.R")) # creation of a complete dataframe RG_aquantil <- limma::RG.MA(MA_aquantil) colnames(RG_aquantil$G) <- colnames(RG$G) df_RG_aquantil_log2 <- cbind(RG_aquantil$R, RG_aquantil$G, RG_aquantil$genes) %>% column_to_rownames("ID") %>% mutate(across(everything(), log2)) Target <- matrix_info %>% dplyr::select(sulfur_condition, genotype, num_combination ) Target=as.data.frame(Target) Target$sulfur_condition = as.factor(Target$sulfur_condition) Target$genotype = as.factor(Target$genotype) model <- model.matrix(as.formula(paste("~", paste(colnames(Target)[1:2], collapse = " + "), paste(colnames(Target[,1:2]),collapse = ":"), sep = " + ")), data=Target) contrast.matrix = Contrasts(model,Target,FALSE,TRUE) fit <- lmFit(df_RG_aquantil_log2, model) #contrast.matrix <- makeContrasts(SpecialDiet-Control, SpecialDietDrug1-SpecialDiet, SpecialDietDrug1-Control, levels=design) fit2 <- contrasts.fit(fit, t(contrast.matrix)) fit2 <- eBayes(fit2) # Create a list of differentially expressed genes per contrast list_genes_by_contrast <- lapply(colnames(fit2$p.value), function(contrast) { # Upregulated genes upregulated <- rownames(fit2)[is_upregulated[, contrast]] # Downregulated genes downregulated <- rownames(fit2)[is_downregulated[, contrast]] list( Up = upregulated, Down = downregulated ) }) # Name the list using the contrast names names(list_genes_by_contrast) <- colnames(fit2$p.value) # ➡️ Filter only the contrasts of interest list_genes_by_contrast_interest <- list_genes_by_contrast[vector_with_interest] # Check the result (for example, for the first contrast of interest) # # # 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` = list_genes_by_contrast_interest[[vector_with_interest[1]]]$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"))) ``` ```{r, eval = F} ## Output DEG list only in EK mutant lfc_lim_i <- 0 pval_i = 0.05 load(file = here::here("data/microarray/output/targets_design.RData")) source(here::here("src/function/microarray/Contrasts.R")) vector_with_interest <- c("[SD-SS]","[E568K-WT2]", "[W78*-WT1]", "[WT1_SD-WT1_SS]","[W78*_SD-W78*_SS]", "[WT2_SD-WT2_SS]", "[E568K_SD-E568K_SS]", "[SS_W78*-SS_WT1]", "[SS_E568K-SS_WT2]", "[SD_W78*-SD_WT1]", "[SD_E568K-SD_WT2]", "[SD_W78*-SD_WT1]-[SS_W78*-SS_WT1]", "[SD_E568K-SD_WT2]-[SS_E568K-SS_WT2]") # creation of a complete dataframe RG_aquantil <- limma::RG.MA(MA_aquantil) colnames(RG_aquantil$G) <- colnames(RG$G) df_RG_aquantil_log2 <- cbind(RG_aquantil$R, RG_aquantil$G, RG_aquantil$genes) %>% column_to_rownames("ID") %>% mutate(across(everything(), log2)) Target <- matrix_info %>% dplyr::select(sulfur_condition, genotype, num_combination ) Target=as.data.frame(Target) Target$sulfur_condition = as.factor(Target$sulfur_condition) Target$genotype = as.factor(Target$genotype) model <- model.matrix(as.formula(paste("~", paste(colnames(Target)[1:2], collapse = " + "), paste(colnames(Target[,1:2]),collapse = ":"), sep = " + ")), data=Target) contrast.matrix = Contrasts(model,Target,FALSE,TRUE) fit <- lmFit(df_RG_aquantil_log2, model) #contrast.matrix <- makeContrasts(SpecialDiet-Control, SpecialDietDrug1-SpecialDiet, SpecialDietDrug1-Control, levels=design) fit2 <- contrasts.fit(fit, t(contrast.matrix)) fit2 <- eBayes(fit2) # Adjust p-values using BH (Benjamini–Hochberg) p_adj <- apply(fit2$p.value, 2, p.adjust, method = "BH") # Determine significance based on 'pval_i' is_significant <- p_adj < pval_i # Define up-/downregulated based on 'lfc_lim_i' and significance is_upregulated <- is_significant & (fit2$coefficients > lfc_lim_i) is_downregulated <- is_significant & (fit2$coefficients < -lfc_lim_i) # Now create the list of genes per contrast list_genes_by_contrast <- lapply(colnames(fit2$p.value), function(contrast) { # Upregulated genes based on your chosen thresholds upregulated <- rownames(fit2)[is_upregulated[, contrast]] # Downregulated genes based on your chosen thresholds downregulated <- rownames(fit2)[is_downregulated[, contrast]] list( Up = upregulated, Down = downregulated ) }) names(list_genes_by_contrast) <- colnames(fit2$p.value) # Subset contrasts of interest list_genes_by_contrast_interest <- list_genes_by_contrast[vector_with_interest] # For Up #### lt_up <- list( `SS_WT1 vs SS_W78*` = list_genes_by_contrast_interest[["[SS_W78*-SS_WT1]"]]$Up, # only PsCam042688 down `SS_WT2 vs SS_E568K` = list_genes_by_contrast_interest[["[SS_E568K-SS_WT2]"]]$Up, `SD_WT1 vs SD_W78*` = list_genes_by_contrast_interest[["[SD_W78*-SD_WT1]"]]$Up, `SD_WT2 vs SD_E568K` = list_genes_by_contrast_interest[["[SD_E568K-SD_WT2]"]]$Up ) m_up <- ComplexHeatmap::make_comb_mat(lt_up, mode = "distinct") # Down-regulated genes lt_down <- list( `SS_WT1 vs SS_W78*` = list_genes_by_contrast_interest[["[SS_W78*-SS_WT1]"]]$Down, # only PsCam042688 down `SS_WT2 vs SS_E568K` = list_genes_by_contrast_interest[["[SS_E568K-SS_WT2]"]]$Down, `SD_WT1 vs SD_W78*` = list_genes_by_contrast_interest[["[SD_W78*-SD_WT1]"]]$Down, `SD_WT2 vs SD_E568K` = list_genes_by_contrast_interest[["[SD_E568K-SD_WT2]"]]$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_contrast(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_contrast",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)), 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)), 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_contrast",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_contrast",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_contrast_",lfc_lim_i,".RData"))) ``` :::: panel-tabset ### With log fold change of 0 ![](plot/upsetplot/upset_up_down_condition_lfc_contrast0.png) ::: {layout-ncol="2"} ![](plot/venn_diagram/euler_down_lfc_contrast0.png) ![](plot/venn_diagram/euler_up_lfc_contrast0.png) ::: **Which group does the psult4 gene belong to?** ```{r} load(file = here::here("data/microarray/output/upset_result_condition_sign_lfc_contrast_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_contrast1.png) ::: {layout-ncol="2"} ![](plot/venn_diagram/euler_down_lfc_contrast1.png) ![](plot/venn_diagram/euler_up_lfc_contrast1.png) ::: **Which group does the psult4 gene belong to?** ```{r} load(file = here::here("data/microarray/output/upset_result_condition_sign_lfc_contrast_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_contrast_",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")) #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_contrast_",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_contrast_",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_lfc0 <- list( I1 = 'SS_WT2 vs SS_E568K', I2 = 'SD_WT2 vs SD_E568K', I3 = c('SS_WT2 vs SS_E568K', 'SD_WT2 vs SD_E568K'), I4 = c('SD_WT1 vs SD_W78*', 'SD_WT2 vs SD_E568K'), I5 = 'SD_WT1 vs SD_W78*', I6 = c('SS_WT2 vs SS_E568K', 'SD_WT2 vs SD_E568K','SD_WT1 vs SD_W78*'), I7 = c('SS_WT2 vs SS_E568K','SD_WT1 vs SD_W78*'), I8 = c('SS_WT2 vs SS_E568K','SD_WT1 vs SD_W78*','SS_WT1 vs SS_W78*') ) intersection_list_lfc1 <- list( I1 = 'SD_WT2 vs SD_E568K', I2 = 'SS_WT2 vs SS_E568K', I3 = c('SS_WT2 vs SS_E568K', 'SD_WT2 vs SD_E568K'), I4 = c('SD_WT1 vs SD_W78*', 'SD_WT2 vs SD_E568K'), I5 = 'SD_WT1 vs SD_W78*', I6 = c('SD_WT1 vs SD_W78*','SS_WT1 vs SS_W78*') ) if(lfc_lim_i == 1){ intersection_list = intersection_list_lfc1 }else if(lfc_lim_i == 0){ intersection_list = intersection_list_lfc0 } 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_contrast_", 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_contrast_0.png) ##### With log fold change of 1 ![](plot/GO/GO_BP_all_absolute_lfc_contrast_1.png) ::: #### On Euler diagram (or upset plot) ::: panel-tabset ##### With log fold change of 0 ![](plot/GO/GO_BP_euler_lfc_contrast_0.png) ##### With log fold change of 1 ![](plot/GO/GO_BP_euler_lfc_contrast_1.png) ::: :::: ## Plot for publication <details> <summary style="color: #007bff;"> Bonus: Details are provided here. </summary> For this plot I only took the standardised data and the comparisons or similarities that I found interesting. For the comparisons I generated a contrast matrix for a linear model which I then used to carry out the differential analysis in limma. For the similarity, I looked at which genes were similar to the two groups. I'm only going to do this for a log fold change equal to 1. ### DEG ```{r, eval = F} lfc_lim_i = 0 load(file = here::here(paste0("data/microarray/output/contraste_lfc_",lfc_lim_i,".RData"))) ## ─── 1. Contrastes de base ────────────────────────────────────────────────────── c1 <- "[W78*-WT1]" c2 <- "[E568K-WT2]" c3_1 <- "[SS_W78*-SS_WT1]" c3_2 <- "[SS_E568K-SS_WT2]" c4_1 <- "[SD_W78*-SD_WT1]" c4_2 <- "[SD_E568K-SD_WT2]" ## ─── 2. Commons (même clé partout) ────────────────────────────────────────────── common1_key <- paste0(c1 , " common ", c2 ) common2_key <- paste0(c3_1, " common ", c3_2) common3_key <- paste0(c4_1, " common ", c4_2) common4_key <- paste0(c3_1, " common ", c4_1) common5_key <- paste0(c3_2, " common ", c4_2) for (pair in list( list(up = paste0("(Up)" , c1), down = paste0("(Down)" , c1), up2 = paste0("(Up)" , c2), down2 = paste0("(Down)" , c2), common_key = common1_key), list(up = paste0("(Up)" , c3_1), down = paste0("(Down)" , c3_1), up2 = paste0("(Up)" , c3_2), down2 = paste0("(Down)" , c3_2), common_key = common2_key), list(up = paste0("(Up)" , c4_1), down = paste0("(Down)" , c4_1), up2 = paste0("(Up)" , c4_2), down2 = paste0("(Down)" , c4_2), common_key = common3_key), list(up = paste0("(Up)" , c3_1), down = paste0("(Down)" , c3_1), up2 = paste0("(Up)" , c4_1), down2 = paste0("(Down)" , c4_1), common_key = common4_key), list(up = paste0("(Up)" , c3_2), down = paste0("(Down)" , c3_2), up2 = paste0("(Up)" , c4_2), down2 = paste0("(Down)" , c4_2), common_key = common5_key)) ) { ## 2.1 Intersection Up / Down up_genes <- intersect(res_list[[pair$up ]], res_list[[pair$up2 ]]) down_genes <- intersect(res_list[[pair$down]], res_list[[pair$down2]]) ## 2.2 Stockage avec un nom cohérent ↓↓↓↓↓ res_list[[paste0("(Up_common)" , pair$common_key)]] <- up_genes res_list[[paste0("(Down_common)", pair$common_key)]] <- down_genes } ## ─── 3. Ordre d’affichage désiré ─────────────────────────────────────────────── vector_contrast0 <- c( "[SD-SS]", common1_key, common2_key, common3_key, "[SS_W78*-SS_WT1]-[SS_E568K-SS_WT2]", "[SD_W78*-SD_WT1]-[SD_E568K-SD_WT2]", "[SD_W78*-SD_WT1]-[SS_W78*-SS_WT1]", "[SD_E568K-SD_WT2]-[SS_E568K-SS_WT2]", common4_key, common5_key#, # "[E568K-WT2]", # "[W78*-WT1]", # "[WT1_SD-WT1_SS]", # "[W78*_SD-W78*_SS]", # "[WT2_SD-WT2_SS]", # "[E568K_SD-E568K_SS]", # "[SS_W78*-SS_WT1]", # "[SS_E568K-SS_WT2]", # "[SD_W78*-SD_WT1]", # "[SD_E568K-SD_WT2]" ) ## ─── 4. Étiquette “common” plus lisible ──────────────────────────────────────── #common_label <- stringr::str_replace_all(Contrast0, "_vs_", " inter ") ## ─── 5. Reconstruction de summary_DEG ─────────────────────────────────────────── summary_DEG <- enframe(res_list, name = "name", value = "genes") %>% mutate( sign = case_when( startsWith(name, "(Up_common)" ) ~ "Up", startsWith(name, "(Down_common)") ~ "Down", startsWith(name, "(Up)" ) ~ "Up", TRUE ~ "Down" ), Contrast0 = sub( "^\$Up_common\$|^\$Down_common\$|^\$Up\$|^\$Down\$", "", name ), length = lengths(genes) ) %>% filter(Contrast0 %in% vector_contrast0) %>% group_by(Contrast0) %>% mutate( total = sum(length), percent = round(100 * length / total, 1) ) %>% ungroup() %>% mutate( Contrast = Contrast0 , sign = factor(sign, levels = c("Up", "Down")), Contrast = factor( Contrast, levels = rev(vector_contrast0) ) )%>% select(Contrast, sign, length, percent, total) # ─── 6. Replot ──────────────────────────────────────────────────────────────────── px <- ggplot(summary_DEG, aes(x = Contrast, y = length, fill = sign)) + geom_col() + geom_text( data = filter(summary_DEG, length >= ifelse(lfc_lim_i == 0, 600, 30)), aes(label = paste0(length, "\n(", percent, "%)")), position = position_stack(vjust = 0.5), color = "white", size = 3 ) + geom_text( data = filter(summary_DEG, length < ifelse(lfc_lim_i == 0, 600, 30)), aes( label = paste0(length, " (", percent, "%)"), color = sign, vjust = ifelse(sign=="Up", -0.5, 1.5), y = total ), hjust = -0.2, size = 3 ) + coord_flip() + theme_minimal() + labs( x = "Comparison", y = "Number of genes deregulated", title = paste0("Number of genes deregulated (including commons)\n", "(FDR = 0.05 ; LogFC = ", lfc_lim_i, ")") ) + scale_fill_manual(values = c(Down = "#1d4877", Up = "#ee3e32"), name = "Sign") + scale_color_manual(values = c(Down = "#1d4877", Up = "#ee3e32"), guide = FALSE) print(px) # export fig_export(path = here::here(paste0("report/microarray/plot/DEG/barplot_publication_differential_expression_lfc_",lfc_lim_i)), px, height = 6.5, width = 12, res_i = 600) save(res_list, file = here::here(paste0("data/microarray/output/contraste_publication_lfc_",lfc_lim_i,".RData"))) ``` ::: panel-tabset ## With log fold change of 0 ![](plot/DEG/barplot_publication_differential_expression_lfc_0.png) ## With log fold change of 1 ![](plot/DEG/barplot_publication_differential_expression_lfc_1.png) ::: ### GO ```{r, eval = F} lfc_lim_i = 0 type_calcul_i = "hypergeometric" #clusterProfiler # hypergeometric # fisher # in reality its exactly the same load(file = here::here(paste0("data/microarray/output/contraste_publication_lfc_",lfc_lim_i,".RData"))) vector_contrast0 <- c( "[SD-SS]", common1_key, common2_key, common3_key, "[SS_W78*-SS_WT1]-[SS_E568K-SS_WT2]", "[SD_W78*-SD_WT1]-[SD_E568K-SD_WT2]", "[SD_W78*-SD_WT1]-[SS_W78*-SS_WT1]", "[SD_E568K-SD_WT2]-[SS_E568K-SS_WT2]", common4_key, common5_key#, # "[E568K-WT2]", # "[W78*-WT1]", # "[WT1_SD-WT1_SS]", # "[W78*_SD-W78*_SS]", # "[WT2_SD-WT2_SS]", # "[E568K_SD-E568K_SS]", # "[SS_W78*-SS_WT1]", # "[SS_E568K-SS_WT2]", # "[SD_W78*-SD_WT1]", # "[SD_E568K-SD_WT2]" ) common_keys <- c(common1_key, common2_key, common3_key, common4_key, common5_key) # 2. Pour chaque contraste brut, on colle les bons préfixes vector_contrast0_sign <- unlist( lapply(vector_contrast0, function(ct) { # si c'est un common → on utilise les préfixes "_common" if (ct %in% common_keys) { paste0(c("(Up_common)", "(Down_common)"), ct) } else { paste0(c("(Up)", "(Down)"), ct) } }) ) gene_table <- tibble( contrast = rep(names(res_list), lengths(res_list)), PsCam = unlist(res_list) ) %>% 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") %>% drop_na(Psat) %>% filter(contrast %in% vector_contrast0_sign) nrow(gene_table %>% filter(contrast == "(Down)[SD-SS]")) universe_vector <- unique(gene_table$Psat) df_GO = GO_on_different_group(functional_roles = "BP", group_info = gene_table, group = "contrast", ID = "Psat", top = 10, universe_i = universe_vector, type_calcul = type_calcul_i ) # combined_results = GO_on_different_group(functional_roles = "BP", # group_info = gene_table, # group = "contrast", # ID = "Psat", # top = FALSE, # universe_i = universe_vector, # type_calcul = type_calcul_i # ) 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, "(?<=\$).*?(?=\$)") ) # 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 = paste0("GO terme for all gene deregulated with lfc:",lfc_lim_i)) + #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 fig_export(path = paste0("report/microarray/plot/GO/GO_BP_contrast_",type_calcul_i,"_",lfc_lim_i), plot_x = px, height_i = 20, width_i = 12, res = 600) ``` </details>