14  Bonus for microarray

#pkg
# install.packages(here::here("data/microarray/org.Psativum1c.eg.db"), repos = NULL, type = "source")
library(org.Psativum1c.eg.db) # if not working install it
library(readxl)
library(tidyverse)
library(patchwork)
library(igraph)
library(factoextra)
library(pheatmap)
library(tidyverse)
library(clusterProfiler)  # BiocManager::install("clusterProfiler")
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/Evaluate_contrast.R"))

# function 
clean_session <- function() {
  # Get all objects in the global environment
  objects <- ls(envir = .GlobalEnv)
  
  # Check if each object is a data frame or a matrix
  is_data <- vapply(objects, function(x) {
    obj <- get(x, envir = .GlobalEnv)
    is.data.frame(obj) || is.matrix(obj)
  }, logical(1))  # Ensure logical vector output
  
  # Remove the identified objects
  rm(list = objects[is_data], envir = .GlobalEnv)
  
  # Reset all plots
  while (dev.cur() > 1) dev.off()
  
  cat_col("Session cleaned: data and plots removed.", "green")
}


cor_visualisation <- function(gene_1, gene_2){
  # gene_1 = "PsCam054590"
  # gene_2 = "PsCam049495"
  
  px = df_compile %>% 
    pivot_longer(cols = colnames(df_compile[,6:length(df_compile)]), names_to = "gene_id", values_to = "value") %>% 
    filter(gene_id %in% c(gene_1, gene_2)) %>% 
    pivot_wider(names_from = "gene_id", values_from = "value") %>% 
    ggplot(., aes_string(x = gene_1 , y =  gene_2, col = "genotype"))+
    geom_point()+
    geom_smooth(method = "lm", se = F)+
    facet_grid(.~sulfur_condition)
  return(px)
}

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

14.1 Co-regulation network (Clustering)

14.1.1 Clustering on deregulator genes (present in venn diagram)

First I have to select the input data. Do they come from everything, or do I select only what’s in the venn diagram?
Here I’ve chosen the venn diagram.

# data importation
load(file = here::here("data/microarray/output/venn/tmp0_peasulf_for_LMM_venn.RData"))
df_origin <- df_compile
rm(df_compile)

variable_v <- colnames(df_origin[6:length(df_origin)]) #peut être pas utile et a supprimer

matrix_scale_select <- df_origin %>%
  dplyr::select(sample_id, all_of(variable_v)) %>% 
  column_to_rownames("sample_id") %>% 
  as.matrix()

matrix_standardized <- scale(matrix_scale_select, center = TRUE, scale = TRUE)

# modifier pour avoir des variables uniquement entre 0 et 1 
normalize_min_max <- function(x) {
  (x - min(x, na.rm = TRUE)) / (max(x, na.rm = TRUE) - min(x, na.rm = TRUE)) * 2 - 1
}

matrix_standardized<-as.data.frame (matrix_scale_select) %>% 
  dplyr::mutate(across(everything(), normalize_min_max))

###### # Finding the optimal number of clusters (begin) #######
hc.cols <- hclust(dist(t(matrix_standardized)), method = "ward.D2")
hc.cols <- hclust(dist(t(matrix_standardized), method = "euclidean"), method = "complete")

plot(hc.cols, sub = "", xlab = "", main = "Clusters des paramètres")

# Calculate distance matrix
distance_matrix <- dist(t(matrix_standardized))
# fviz_nbclust(as.data.frame(distance_matrix), FUN = hcut, method = "wss")
# fviz_nbclust(as.data.frame(distance_matrix), FUN = hcut, method = "silhouette")
# fviz_nbclust(as.data.frame(distance_matrix), FUN = hcut, method = "gap_stat", nboot = 50)

# Select a random subset of columns
set.seed(100)  # For the reproductibility
sample_columns <- sample(ncol(matrix_standardized), 40)
matrix_sample <- matrix_standardized[, sample_columns]

# Calculate the distance matrix on the subset
distance_matrix_sample <- dist(t(matrix_sample))

# Use fviz_nbclust for the subset
fviz_nbclust(as.data.frame(distance_matrix_sample), FUN = hcut, method = "wss", k.max = 20)
fviz_nbclust(as.data.frame(distance_matrix_sample),FUN = hcut, method = "silhouette", k.max = 20)
fviz_nbclust(as.data.frame(distance_matrix_sample),FUN = hcut, method = "gap_stat", k.max = 20)

# or with PCA
# Dimension reduction with PCA (keeping the first 20 components)
pca_result <- prcomp(matrix_standardized, center = TRUE, scale. = TRUE)
matrix_reduced <- pca_result$x[, 1:32]

# Calculate distance matrix for principal components
distance_matrix_pca <- dist(matrix_reduced)

# Use fviz_nbclust to determine the number of clusters
fviz_nbclust(as.data.frame(distance_matrix_pca), FUN = hcut, method = "wss")

###### # Finding the optimal number of clusters (end) #######

# parameter
num_clusters <- 11  # Number of cluster

cluster_assignments <- cutree(hc.cols, k = num_clusters)

# Associate each variable with its cluster
clustered_variables <- data.frame(
  variable = colnames(matrix_standardized),
  cluster = cluster_assignments
)

annotation_row <- df_origin %>%
  distinct(sample_id, genotype, sulfur_condition) %>%
  column_to_rownames("sample_id") %>% 
  dplyr::rename(Genotype = genotype, 
         `Sulfur condition`= sulfur_condition
  )


##################
annotation_col <- data.frame(variable = colnames(matrix_standardized)) %>%
  full_join(., clustered_variables, by="variable") %>%
  mutate(cluster = as.character(cluster)) %>% 
  column_to_rownames("variable") %>%
  dplyr::rename(Cluster = cluster)# rename col cluster to be the same as annotation_colors

# levels(as.factor(annotation_col$variable_type))
my_color_palette <- read_excel(here::here("data/color_palette.xlsm"))
color_cluster_v=c(my_color_palette %>%
             filter(set=="cluster") %>%
             pull(color), "gray20", "gray70"
           )

# palette <- colorRampPalette(c("#01a2a4", "#FFD275", "#e3420e"))(500)
palette <- viridis::viridis(500)

annotation_colors <- list(
  `Sulfur condition` = c(sulfate_pallet),
  Genotype = c("WT1" = "#003049", "W78*" =  "#780000" , "WT2" =  "#7FACC7" , "E568K" = "#EC323E"),
  Cluster = setNames(color_cluster_v[1:11], as.character(1:11))
)

#heatmap(matrix_standardized,Colv=T, scale='none', col=palette,cexCol=0.6)
heatmap_all <- pheatmap(matrix_standardized,
         scale = "none", 
         cluster_rows = T, 
         cluster_cols = TRUE, 
         #main = paste("Heatmap global"),
         color = palette, 
         annotation_col = annotation_col,
         annotation_row = annotation_row,  
         annotation_colors = annotation_colors,
         cutree_cols = num_clusters, 
         show_rownames = T,
         show_colnames = F
         #cex = 0.8
         ) ; heatmap_all

# export
fig_export(path = "report/microarray/plot/heatmap/heatmap_all", plot_x = heatmap_all, height_i = 9, width_i = 18, res = 600)
write_csv(x = clustered_variables, file = here::here("data/microarray/output/clustered_genes_all.csv"))

####### tcheck proximity between cluster  using mean of expression profile of the genes

centroids <- aggregate(t(matrix_standardized), by = list(clustered_variables %>% deframe()), FUN = mean)  %>%
  rename(cluster = Group.1)
dist_matrix <- dist(centroids[,-1], method = "euclidean")
cor_matrix <- cor(t(centroids %>% column_to_rownames("cluster")), method = "pearson")
dist_matrix <- as.matrix(dist_matrix)  # convert to a matrix for pheatmap

# Create an annotation data frame for the clusters using the row names (which are the cluster labels)
cluster_annotation <- data.frame(Cluster = rownames(dist_matrix))
# Ensure the row names of the annotation match the distance matrix rows
rownames(cluster_annotation) <- rownames(dist_matrix)

mat_list <- list(dist_matrix, cor_matrix)
title_list <- c("Euclidean Distance Matrix", "Pearson Correlation Matrix")

# Initialize an empty list to store the ggplot objects
plot_list <- list()

# Loop over the matrices
for(i in seq_along(mat_list)) {
  # Create the pheatmap plot silently
  p <- pheatmap(mat_list[[i]],
                color = palette,
                annotation_row = cluster_annotation,
                annotation_col = cluster_annotation,
                annotation_colors = annotation_colors,
                silent = TRUE)
  
  # Convert the pheatmap grob to a ggplot object and add a title
  plot_list[[i]] <- ggplotify::as.ggplot(p$gtable) + ggtitle(title_list[i])
}

# Combine the plots side by side using patchwork
combined_plot <- plot_list[[1]] + plot_list[[2]]+ plot_annotation(title= "Assessing Cluster Similarity for Potential Merging")

fig_export(path = "report/microarray/plot/heatmap/similarity_between_cluster_all", plot_x = combined_plot, height_i = 8, width_i = 10, res = 300)

# Create a heatmap for each cluster with annotations

for (i in 1:num_clusters) {
  # Select columns belonging to cluster i
  cluster_columns <- clustered_variables %>% 
    filter(cluster == i) %>%
    pull(variable)
  
  n_var=length(cluster_columns)
  
  # Sub-matrix for the cluster
  matrix_cluster <- matrix_standardized[, cluster_columns, drop = FALSE]
  
  # Display heatmap for cluster with annotations
  if (ncol(matrix_cluster) > 1) {  # Checks if there is more than one column in the cluster
    heatmap_part<- pheatmap(matrix_cluster, 
            scale = "none", 
            cluster_rows = T, 
            cluster_cols = TRUE, 
            main = paste("Heatmap of cluster", i, " n_var:", n_var),
            color = palette, 
            annotation_col = annotation_col,
            annotation_row = annotation_row,  
            annotation_colors = annotation_colors,
            show_rownames = T,
            show_colnames = F
    )
  } else {
    message(paste("Cluster", i, "contains only one variable :", cluster_columns))
  }
  png(here::here(paste0("report/microarray/plot/heatmap/heatmap_part/heatmap_part_",i,".png")), height = 8*300, width = 16*300, res=300)
  print(heatmap_part)
  dev.off()
}

df_long <- as.data.frame(matrix_standardized) %>%
  #dplyr::mutate(across(everything(), normalize_min_max)) %>% 
  rownames_to_column("sample_id") %>%
  pivot_longer(-sample_id, names_to = "variable", values_to = "value") %>%
  left_join(clustered_variables, by = "variable") %>%       # Ajouter le cluster
  left_join(df_origin %>% dplyr::select(sample_id,genotype, sulfur_condition) %>% distinct(sample_id,genotype, sulfur_condition), by = c("sample_id")) %>%  # Ajouter les conditions
  select(sample_id, variable, value, cluster, genotype, sulfur_condition) %>% 
  mutate(condition = paste0(genotype, "_", sulfur_condition))

# Data summarization
df_summary <- df_long %>%
  dplyr::group_by(cluster, condition, sulfur_condition, genotype) %>%
  dplyr::summarise(
    mean_value = mean(value, na.rm = TRUE),
    sd_value = sd(value, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  dplyr::mutate(
    genotype = factor(genotype, levels = c("WT1", "W78*", "WT2", "E568K")),
    condition = factor(condition, levels = c(
      "WT1_SS", "W78*_SS", "WT2_SS", "E568K_SS", "WT1_SD", "W78*_SD", "WT2_SD", "E568K_SD"
    )),
    cluster = as.factor(cluster),
    cluster_col = color_cluster_v[as.numeric(cluster)],
    cluster = paste0("Cluster: ", cluster), 
    cluster = factor(cluster, levels = paste0("Cluster: ", 1:num_clusters)),
  )

# Prepare colors for facet strips
df_summary_color <- df_summary %>%
  distinct(cluster, cluster_col) %>%
  arrange(cluster) %>%
  mutate(couleur_texte_perso = sapply(cluster_col, evaluate_contrast))

# Create the plot
plot_x <- ggplot(df_summary, aes(x = condition, y = mean_value, fill = genotype)) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.8), color = "black") +
  geom_errorbar(
    aes(ymin = mean_value - sd_value, ymax = mean_value + sd_value),
    width = 0.2, position = position_dodge(width = 0.8)
  ) +
  labs(x = "Condition",   y = "Mean and standard deviation") +
  theme_bw() +
  scale_fill_manual(name = "Genotype", values = c("WT1" = "#003049", "W78*" =  "#780000" , "WT2" =  "#7FACC7" , "E568K" = "#EC323E")) +
  theme(
    panel.spacing = unit(0.5, "lines"),
    axis.ticks.y = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1),
    legend.text = element_text(size = 12),
    legend.title = element_text(size = 12),
    legend.key.size = unit(10, "pt"),
    plot.caption = element_text(size = 8),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
  ) +
  ggh4x::facet_wrap2(
    ~cluster, nrow = 5,
    strip = ggh4x::strip_themed(
      background_x = ggh4x::elem_list_rect(fill = df_summary_color$cluster_col),
      text_x = ggh4x::elem_list_text(colour = df_summary_color$couleur_texte_perso)
    )
  )

# Display the plot
print(plot_x)

fig_export(path = "report/microarray/plot/heatmap/barplot_cluster", plot_x = plot_x, height_i = 12, width_i = 10, res = 600)

# 
clustered_variables %>% filter(variable == "PsCam042688") %>% pull(cluster)

Note

In which cluster is psult4? Is in cluster 6

Bonnus: Heatmap by cluster

1

2

3

4

5

6

7

8

9

10

11

14.1.2 Clusterisitation only on genes deregulated under sulfur stress.

# data importation
load(file = here::here("data/microarray/output/venn/tmp0_peasulf_for_LMM_venn.RData"))
load(file = here::here("data/microarray/output/upset_result_condition_sign.RData"))

df_origin <- df_compile %>% filter(sulfur_condition == "SD")
rm(df_compile)

variable_v <- unique(c(
  lt_up$`SD_WT1 vs SD_W78*`,
  lt_up$`SD_WT2 vs SD_E568K`,

  lt_down$`SD_WT1 vs SD_W78*`,
  lt_down$`SD_WT2 vs SD_E568K`
))

length(variable_v)

#variable_v <- colnames(df_origin[6:length(df_origin)]) #peut être pas utile et a supprimer

matrix_scale_select <- df_origin %>%
  dplyr::select(sample_id, all_of(variable_v)) %>% 
  column_to_rownames("sample_id") %>% 
  as.matrix()

matrix_standardized <- scale(matrix_scale_select, center = TRUE, scale = TRUE)

# modifier pour avoir des variables uniquement entre 0 et 1 
normalize_min_max <- function(x) {
  (x - min(x, na.rm = TRUE)) / (max(x, na.rm = TRUE) - min(x, na.rm = TRUE)) * 2 - 1
}

matrix_standardized<-as.data.frame (matrix_scale_select) %>% 
  dplyr::mutate(across(everything(), normalize_min_max))

###### # Finding the optimal number of clusters (begin) #######
hc.cols <- hclust(dist(t(matrix_standardized)), method = "ward.D2")
hc.cols <- hclust(dist(t(matrix_standardized), method = "euclidean"), method = "complete")

plot(hc.cols, sub = "", xlab = "", main = "Clusters des paramètres")

# Calculate distance matrix
distance_matrix <- dist(t(matrix_standardized))
# fviz_nbclust(as.data.frame(distance_matrix), FUN = hcut, method = "wss")
# fviz_nbclust(as.data.frame(distance_matrix), FUN = hcut, method = "silhouette")
# fviz_nbclust(as.data.frame(distance_matrix), FUN = hcut, method = "gap_stat", nboot = 50)

# Select a random subset of columns
set.seed(100)  # For the reproductibility
sample_columns <- sample(ncol(matrix_standardized), 40)
matrix_sample <- matrix_standardized[, sample_columns]

# Calculate the distance matrix on the subset
distance_matrix_sample <- dist(t(matrix_sample))

# Use fviz_nbclust for the subset
fviz_nbclust(as.data.frame(distance_matrix_sample), FUN = hcut, method = "wss", k.max = 20)
fviz_nbclust(as.data.frame(distance_matrix_sample),FUN = hcut, method = "silhouette", k.max = 20)
fviz_nbclust(as.data.frame(distance_matrix_sample),FUN = hcut, method = "gap_stat", k.max = 20)

# or with PCA
# Dimension reduction with PCA (keeping the first 20 components)
pca_result <- prcomp(matrix_standardized, center = TRUE, scale. = TRUE)
matrix_reduced <- pca_result$x[, 1:20]

# Calculate distance matrix for principal components
distance_matrix_pca <- dist(matrix_reduced)

# Use fviz_nbclust to determine the number of clusters
fviz_nbclust(as.data.frame(distance_matrix_pca), FUN = hcut, method = "wss")

###### # Finding the optimal number of clusters (end) #######

# parameter
num_clusters <- 8  # Number of cluster

cluster_assignments <- cutree(hc.cols, k = num_clusters)

# Associate each variable with its cluster
clustered_variables <- data.frame(
  variable = colnames(matrix_standardized),
  cluster = cluster_assignments
)

annotation_row <- df_origin %>%
  distinct(sample_id, genotype, sulfur_condition) %>%
  column_to_rownames("sample_id") %>% 
  dplyr::rename(Genotype = genotype, 
         `Sulfur condition`= sulfur_condition
  )


##################
annotation_col <- data.frame(variable = colnames(matrix_standardized)) %>%
  full_join(., clustered_variables, by="variable") %>%
  mutate(cluster = as.character(cluster)) %>% 
  column_to_rownames("variable") %>%
  dplyr::rename(Cluster = cluster)# rename col cluster to be the same as annotation_colors

# levels(as.factor(annotation_col$variable_type))
my_color_palette <- read_excel(here::here("data/color_palette.xlsm"))
color_cluster_v=c(my_color_palette %>%
             filter(set=="cluster") %>%
             pull(color), "gray20", "gray70"
           )

color_cluster_v= color_cluster_v[1:num_clusters]

#palette <- colorRampPalette(c("#01a2a4", "#FFD275", "#e3420e"))(500)
palette <- viridis::viridis(500)

annotation_colors <- list(
  `Sulfur condition` = c(sulfate_pallet),
  Genotype = c("WT1" = "#003049", "W78*" =  "#780000" , "WT2" =  "#7FACC7" , "E568K" = "#EC323E"),
  Cluster = setNames(color_cluster_v[1:num_clusters], as.character(1:num_clusters))
)

#heatmap(matrix_standardized,Colv=T, scale='none', col=palette,cexCol=0.6)
heatmap_all <- pheatmap(matrix_standardized,
         scale = "none", 
         cluster_rows = T, 
         cluster_cols = TRUE, 
         #main = paste("Heatmap global"),
         color = palette, 
         annotation_col = annotation_col,
         annotation_row = annotation_row,  
         annotation_colors = annotation_colors,
         cutree_cols = num_clusters, 
         show_rownames = T,
         show_colnames = F
         #cex = 0.8
         ) ; heatmap_all

# export
fig_export(path = "report/microarray/plot/heatmap/heatmap_all_SD", plot_x = heatmap_all, height_i = 9, width_i = 18, res = 600)


####### tcheck proximity between cluster  using mean of expression profile of the genes

centroids <- aggregate(t(matrix_standardized), by = list(clustered_variables %>% deframe()), FUN = mean)  %>%
  rename(cluster = Group.1)
dist_matrix <- dist(centroids[,-1], method = "euclidean")
cor_matrix <- cor(t(centroids %>% column_to_rownames("cluster")), method = "pearson")
dist_matrix <- as.matrix(dist_matrix)  # convert to a matrix for pheatmap

# Create an annotation data frame for the clusters using the row names (which are the cluster labels)
cluster_annotation <- data.frame(Cluster = rownames(dist_matrix))
# Ensure the row names of the annotation match the distance matrix rows
rownames(cluster_annotation) <- rownames(dist_matrix)

mat_list <- list(dist_matrix, cor_matrix)
title_list <- c("Euclidean Distance Matrix", "Pearson Correlation Matrix")

# Initialize an empty list to store the ggplot objects
plot_list <- list()

# Loop over the matrices
for(i in seq_along(mat_list)) {
  # Create the pheatmap plot silently
  p <- pheatmap(mat_list[[i]],
                color = palette,
                annotation_row = cluster_annotation,
                annotation_col = cluster_annotation,
                annotation_colors = annotation_colors,
                silent = TRUE)
  
  # Convert the pheatmap grob to a ggplot object and add a title
  plot_list[[i]] <- ggplotify::as.ggplot(p$gtable) + ggtitle(title_list[i])
}

# Combine the plots side by side using patchwork
combined_plot <- plot_list[[1]] + plot_list[[2]]+ plot_annotation(title= "Assessing Cluster Similarity for Potential Merging")

fig_export(path = "report/microarray/plot/heatmap/similarity_between_cluster_SD", plot_x = combined_plot, height_i = 8, width_i = 10, res = 300)
write_csv(x = clustered_variables, file = here::here("data/microarray/output/clustered_genes_SD.csv"))

# Create a heatmap for each cluster with annotations

for (i in 1:num_clusters) {
  # Select columns belonging to cluster i
  cluster_columns <- clustered_variables %>% 
    filter(cluster == i) %>%
    pull(variable)
  
  n_var=length(cluster_columns)
  
  # Sub-matrix for the cluster
  matrix_cluster <- matrix_standardized[, cluster_columns, drop = FALSE]
  
  # Display heatmap for cluster with annotations
  if (ncol(matrix_cluster) > 1) {  # Checks if there is more than one column in the cluster
    heatmap_part<- pheatmap(matrix_cluster, 
            scale = "none", 
            cluster_rows = T, 
            cluster_cols = TRUE, 
            main = paste("Heatmap of cluster", i, " n_var:", n_var),
            color = palette, 
            annotation_col = annotation_col,
            annotation_row = annotation_row,  
            annotation_colors = annotation_colors,
            show_rownames = T,
            show_colnames = F
    )
  } else {
    message(paste("Cluster", i, "contains only one variable :", cluster_columns))
  }
  png(here::here(paste0("report/microarray/plot/heatmap/heatmap_part/heatmap_part_",i,"_SD.png")), height = 8*300, width = 16*300, res=300)
  print(heatmap_part)
  dev.off()
}

df_long <- as.data.frame(matrix_standardized) %>%
  #dplyr::mutate(across(everything(), normalize_min_max)) %>% 
  rownames_to_column("sample_id") %>%
  pivot_longer(-sample_id, names_to = "variable", values_to = "value") %>%
  left_join(clustered_variables, by = "variable") %>%       # Ajouter le cluster
  left_join(df_origin %>% dplyr::select(sample_id,genotype, sulfur_condition) %>% distinct(sample_id,genotype, sulfur_condition), by = c("sample_id")) %>%  # Ajouter les conditions
  select(sample_id, variable, value, cluster, genotype, sulfur_condition) %>% 
  mutate(condition = paste0(genotype, "_", sulfur_condition))

# Data summarization
df_summary <- df_long %>%
  dplyr::group_by(cluster, condition, sulfur_condition, genotype) %>%
  dplyr::summarise(
    mean_value = mean(value, na.rm = TRUE),
    sd_value = sd(value, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  dplyr::mutate(
    genotype = factor(genotype, levels = c("WT1", "W78*", "WT2", "E568K")),
    condition = factor(condition, levels = c(
      "WT1_SS", "W78*_SS", "WT2_SS", "E568K_SS", "WT1_SD", "W78*_SD", "WT2_SD", "E568K_SD"
    )),
    cluster = as.factor(cluster),
    cluster_col = color_cluster_v[as.numeric(cluster)],
    cluster = paste0("Cluster: ", cluster), 
    cluster = factor(cluster, levels = paste0("Cluster: ", 1:num_clusters)),
  )

# Prepare colors for facet strips
df_summary_color <- df_summary %>%
  distinct(cluster, cluster_col) %>%
  arrange(cluster) %>%
  mutate(couleur_texte_perso = sapply(cluster_col, evaluate_contrast))

# Create the plot
plot_x <- ggplot(df_summary, aes(x = condition, y = mean_value, fill = genotype)) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.8), color = "black") +
  geom_errorbar(
    aes(ymin = mean_value - sd_value, ymax = mean_value + sd_value),
    width = 0.2, position = position_dodge(width = 0.8)
  ) +
  labs(x = "Condition",   y = "Mean and standard deviation") +
  theme_bw() +
  scale_fill_manual(name = "Genotype", values = c("WT1" = "#003049", "W78*" =  "#780000" , "WT2" =  "#7FACC7" , "E568K" = "#EC323E")) +
  theme(
    panel.spacing = unit(0.5, "lines"),
    axis.ticks.y = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1),
    legend.text = element_text(size = 12),
    legend.title = element_text(size = 12),
    legend.key.size = unit(10, "pt"),
    plot.caption = element_text(size = 8),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
  ) +
  ggh4x::facet_wrap2(
    ~cluster, nrow = 2,
    strip = ggh4x::strip_themed(
      background_x = ggh4x::elem_list_rect(fill = df_summary_color$cluster_col),
      text_x = ggh4x::elem_list_text(colour = df_summary_color$couleur_texte_perso)
    )
  )

# Display the plot
print(plot_x)

fig_export(path = "report/microarray/plot/heatmap/barplot_cluster_SD", plot_x = plot_x, height_i = 6, width_i = 6, res = 600)

Note

In which cluster is Psult4? Is in cluster 6

Bonnus: Heatmap by cluster for SD condition

1

2

3

4

5

6

7

8

14.1.3 GO enrichment

# parameter
# "All" for All gene and the different cluster or "SD" for gene only deregulated for SD condition
data_type <- "All"

# import data
load(here::here("data/microarray/output/raw_data_microarray_leaf_PeaSulf.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")
# all_GO_BP <- read_csv(file = here::here("data/microarray/output/all_GO_BP.csv"))


# Import result from clustering and replace PsCam by Psat
cluster_info <- read_csv(file = here::here(paste0("data/microarray/output/clustered_genes_",data_type,".csv")),show_col_types = FALSE) %>% 
  dplyr::rename(PsCam = variable) %>% 
  left_join(., df_info_gene, by = "PsCam") %>% 
  mutate(ID = Psat) %>% 
  drop_na(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)))

# 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() +
  scale_x_discrete(breaks = seq(1,length(unique(cluster_info$cluster)),1))+
  theme(axis.text = element_text(size = 8, colour = "black"),
        panel.grid.major = element_blank(),
        plot.background = element_rect(fill = "white", colour = "white")) +
  labs(fill = "-log10(FDR)",
       x= "Biological process",
       y = "Cluster", 
       title = ifelse(data_type == "All", "GO terme for all gene deregulated", "GO terme for all gene deregulated in \nSD condition")) +
  #scale_size_manual(values = c(dot = 2, no_dot = NA), guide = "none")+
  new_scale_fill() +
 scale_fill_manual(
    values = my_color_palette %>%
      filter(set == "cluster") %>%
      pull(color),
    name = "Cluster"
  ) +
geom_tile(
    aes(
      x = cluster,
      y = -0.025,     # If you truly want it at a negative y-value, 
                      # ensure your y-scale is continuous or can handle this
      fill = cluster,
      width = 0.95,
      height = 0.60
    ),
    data = Cluster_selected_GO_filled,
    color = "black",
    alpha = 1,
    inherit.aes = FALSE
  )

# export 
fig_export(path = paste0("report/microarray/plot/GO/GO_",data_type), plot_x = px, height_i = 13, width_i = 7.5, res = 600)
#test_to_verif = Cluster_selected_GO_filled %>% filter(str_detect(geneID, "Psat3g185920"))

PSULT4 is annoted as monoatomic anion transport

14.1.3.1 List of genes deregulated only in the EK mutant

# Take data from Upsetplot. Gene DEG only in mutant EK compare to WT2 in SS and SD
load(here::here("data/microarray/output/raw_data_microarray_leaf_PeaSulf.RData"))
load(file = here::here("data/microarray/output/upset_result_condition_sign.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)

levels(as.factor(cluster_info$comparison))

in_mutant_EK <- cluster_info %>% filter(comparison %in% c("SS_WT2 vs SS_E568K", "SD_WT2 vs SD_E568K")) %>% pull(Psat)

in_other <- cluster_info %>% filter(comparison %in% c("SD_WT1 vs SD_W78*", "SS_WT1 vs SS_W78*" )) %>% pull(Psat)

length(only_in_EK)

only_in_EK <- setdiff(in_mutant_EK, in_other) %>% tibble(ID = .) %>% mutate(group = "nos_group")

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

write.csv(df_GO, file = "data/microarray/output/gene_only_in_mutant_EK.csv")

14.2 Co-regulation network (LMM)

#pkg
library(readxl)
library(tidyverse)
library(patchwork)
library(igraph)
library(factoextra)
library(pheatmap)

# 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/stat_function/stat_analysis_main.R")) # for make plot 
source(here::here("src/function/Evaluate_contrast.R"))

# function 
clean_session <- function() {
  # Get all objects in the global environment
  objects <- ls(envir = .GlobalEnv)
  
  # Check if each object is a data frame or a matrix
  is_data <- vapply(objects, function(x) {
    obj <- get(x, envir = .GlobalEnv)
    is.data.frame(obj) || is.matrix(obj)
  }, logical(1))  # Ensure logical vector output
  
  # Remove the identified objects
  rm(list = objects[is_data], envir = .GlobalEnv)
  
  # Reset all plots
  while (dev.cur() > 1) dev.off()
  
  cat_col("Session cleaned: data and plots removed.", "green")
}


cor_visualisation <- function(gene_1, gene_2){
  # gene_1 = "PsCam054590"
  # gene_2 = "PsCam049495"
  
  px = df_compile %>% 
    pivot_longer(cols = colnames(df_compile[,6:length(df_compile)]), names_to = "gene_id", values_to = "value") %>% 
    filter(gene_id %in% c(gene_1, gene_2)) %>% 
    pivot_wider(names_from = "gene_id", values_from = "value") %>% 
    ggplot(., aes_string(x = gene_1 , y =  gene_2, col = "genotype"))+
    geom_point()+
    geom_smooth(method = "lm", se = F)+
    facet_grid(.~sulfur_condition)
  return(px)
}

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

14.2.1 Data formatting

# Data importation
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_aquantil_global <- bind_cols(RG_aquantil$genes, RG_aquantil$R, RG_aquantil$G) %>% 
  column_to_rownames("ID")

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

# transform into vertical df
RG_aquantil_global_log2_good_format <- RG_aquantil_global_log2 %>% 
  rownames_to_column("variable") %>% 
  pivot_longer(-variable, names_to = "sample_id", values_to = "value") %>% 
  left_join(.,read_csv(here::here("data/microarray/output/list_microarray_microarray_leaf_PeaSulf.csv"), show_col_types = FALSE) %>% 
              dplyr::select(genotype, sulfur_condition,simplify_condition, sample_id, color),
            by = "sample_id") %>% 
  pivot_wider(names_from = "variable", values_from = "value")

df_compile <- RG_aquantil_global_log2_good_format

save(df_compile, file = here::here("data/microarray/output/cor_08/tmp0_peasulf_for_LMM_all.RData"))

##################### filter to take only gene from venn diagram ################

# data from venn diragram
load(file = here::here("data/microarray/output/upset_result_condition_sign.RData"))

all_list <- c(lt_up, lt_down)

all_genes <- unlist(all_list, use.names = FALSE)

unique_genes <- unique(all_genes)

n <-length(unique_genes)

# nb_combination
nb_combination <- n * (n - 1) / 2
nb_combination
timing <-  nb_combination/100000*10/60/24

# 10 for 10 min 
# donc a peut près 17 jours


df_compile <- df_compile %>% 
  dplyr::select(sample_id, genotype, sulfur_condition,simplify_condition,color, all_of(unique_genes))

save(df_compile, file = here::here("data/microarray/output/venn/tmp0_peasulf_for_LMM_venn.RData"))

14.2.2 Calcul of all comparisons with MLM

See code for the computing cluster.

Code for analysing create co-regulation network

Download the R script here

14.2.2.1 Merge of results from the calculation server

# reassemble to recreate tmp1
reassemble_export_tmp1 <- function(end_with="75"){
  list_files <- list.files(path = here::here("data/microarray/output/cor_08/LMM_cut_df/"), pattern = paste0(end_with,"\\.csv$"), full.names = TRUE)
  combined_data <- do.call(rbind, lapply(list_files, function(f) {
  data <- read.csv(f)
  data <- data[, -1]
  return(data)
}))
  df_result=combined_data
  save(df_result, file = here::here("data/microarray/output/cor_08/tmp1_LMM_peasulf_microarray.RData"))
}

reassemble_export_tmp1()

14.2.3 Selection based on pvalue and R²

R2m (Marginal R²): This represents the variance explained by the fixed effects in your model. It is similar in interpretation to the R² value in traditional linear regression models. R2c (Conditional R²): This represents the variance explained by both the fixed and random effects in your model. It considers both fixed and random effects, providing a more comprehensive measure of goodness-of-fit for mixed-effects models. Both Rm2 and R2c can be useful in understanding the proportion of variance explained by the predictors in your model.

clean_session() #  Remove data frames and matrices and  Reset plots

# data importation
load(file = here::here("data/microarray/output/cor_08/tmp1_LMM_peasulf_microarray.RData"))

# parameter
r2c_i = 0.85
r2m_i = 0.5
fdr_lim = 0.05

df_result = df_result %>%  
  mutate(fdr = p.adjust(pval, method = "fdr"),
         log10fdr = -log10(fdr)) %>% 
  mutate(color=ifelse(r2c*1.5 - r2m*1.5 > 0.5, "red",
                      ifelse(r2m<r2m_i, ifelse(r2c<r2c_i,"#DD6342","orange"),
                             ifelse(fdr<fdr_lim, "green","forestgreen")
                             )
                      )
         )
# Visualization 
## Use of randomized data to avoid displaying all points
random_result <- df_result %>%
  sample_n(1000)

random_result%>% 
  ggplot(aes(x = r2m, y = -log10(pval), size = cor, col = -log10(pval))) +
  geom_point(alpha = 0.7) +  
  scale_color_viridis_c() + 
  labs(title = "Dispersion diagramme", 
       x = "R2m", 
       y = "-log10(fdr)", 
       color = "-log10(fdr)") +
  theme_minimal()+

random_result %>% 
  ggplot(aes(x = r2c, y = -log10(pval), size = cor, col = -log10(pval))) +
  geom_point(alpha = 0.7) +  
  scale_color_viridis_c() + 
  labs(title = "Dispersion diagramme", 
       x = "R2c", 
       y = "-log10(fdr)", 
       color = "-log10(fdr)") +
  theme_minimal()

random_result %>% 
  ggplot(aes(x = r2m, y = r2c, col = color, size = abs(cor))) +
  geom_point(aes(shape = cor < 0), na.rm = TRUE) +
  scale_shape_manual(values = c(`TRUE` = 17, `FALSE` = 16)) +
  scale_color_manual(values = c("red" = "red", "#DD6342"="#DD6342","orange" = "orange", "green" = "green", "forestgreen" = "forestgreen")) +  # Utilisation de couleurs fixes
  geom_vline(xintercept = r2m_i, linetype = "dashed", color = "red", linewidth = 1) +
  geom_hline(yintercept = r2c_i, linetype = "dashed", color = "red", linewidth = 1)

# filter the result to have variable
df_result_select <- df_result %>% 
  filter(color=="green") %>% 
  arrange(pval)

cat_col(paste("Number of variable:",length(unique(c(df_result_select$V1,df_result_select$V2))), "\n"), "green")

# export
save(df_result_select, file = here::here("data/microarray/output/cor_08/tmp2_LMM_peasulf_microarray.RData"))

round(dim(df_result_select)[1]*100/dim(df_result)[1], 2)

14.2.4 Test with a prefiltre based on matrix correlation with filter of 0.8

14.2.4.1 Network plot with Igraph

clean_session() #  Remove data frames and matrices and  Reset plots
load(file = here::here("data/microarray/output/cor_08/tmp2_LMM_peasulf_microarray.RData"))
load(file = here::here("data/microarray/output/cor_08/tmp0_peasulf_for_LMM_all.RData"))

# pb <- progress::progress_bar$new(
#   format = "  Calcul [:bar] :percent | Étape :current/:total | Temps écoulé :elapsed | Temps restant :eta",
#   total = nrow(df_result_select),
#   clear = FALSE,
#   width = 100
# )

# Fonction pour calculer la corrélation avec mise à jour de la barre de progression
# calculate_corr <- function(var1, var2) {
#   pb$tick()
#   
#   # Vérifier si les variables existent dans df_compile
#   if (all(c(var1, var2) %in% colnames(df_compile))) {
#     # Calculer la corrélation
#     cor_value <- cor(df_compile[[var1]], df_compile[[var2]], use = "pairwise.complete.obs")
#     return(cor_value)
#   } else {
#     return(NA) # Retourne NA si une variable manque
#   }
# }
# 
# # Utiliser mapply pour vectoriser le calcul
# df_result_select$corr <- mapply(calculate_corr, df_result_select$V1, df_result_select$V2)


# filter the data previously generate


df_select=df_result_select %>% 
  #filter(select=="yes") %>% 
  filter(color=="green") %>% 
  filter(between_IC01=="no")

#hist(df_select$slope_modif)
# creation of the node 

nodes<-data.frame(variable = unique(c(df_select$V1,df_select$V2))) %>% 
  mutate(id = paste0("s", 1:length(unique(c(df_select$V1,df_select$V2))))) %>% 
  relocate(id, .before=variable)
  
# links
links<-df_select %>% 
  dplyr::select(V1,V2, cor) %>% 
    dplyr::mutate(weight = abs(cor)) %>%
  dplyr::mutate(weight_repulstion = 1-cor) %>%
  inner_join(nodes %>% dplyr::rename(V1=variable) %>% dplyr::rename(X1=id),.,by="V1") %>% 
  inner_join(nodes %>% dplyr::rename(V2=variable) %>% dplyr::rename(X2=id),.,by="V2") %>% 
  relocate(X1, .before=X2) %>% 
  relocate(V1, .before=V2) 

# Converting the data to an igraph object:
net <- graph.data.frame(links, nodes, directed=T) 

save(nodes, links, net,file = here::here("data/microarray/output/cor_08/tmp3_LMM_peasulf_microarray.RData"))
load(file = here::here("data/microarray/output/cor_08/tmp3_LMM_peasulf_microarray.RData"))


# Compute node degree (#links) and use it to set node size:
deg <- degree(net, mode="all")
tableau_connexions <- data.frame(id = names(deg), nb_connexions = deg) %>% 
  left_join(., nodes, by = "id") %>% 
  arrange(desc(nb_connexions))

tableau_connexions %>% head(10)

# if(fc_or_pval=="pval"){
#   V(net)$size <- V(net)$log10_pval_test_t*6 #ixi modifie la taille des points
# }else if(fc_or_pval=="fc"){
#   V(net)$size <- V(net)$abslog2fc*10 #ixi modifie la taille des points
#   }

myrhombus <- function(coords, v = NULL, params) {
  vertex.color <- params("vertex", "color")
  if (length(vertex.color) != 1 && !is.null(v)) {
    vertex.color <- vertex.color[v]
  }
  vertex.size <- 1/200 * params("vertex", "size")
  if (length(vertex.size) != 1 && !is.null(v)) {
    vertex.size <- vertex.size[v]
  }

  symbols(x = coords[, 1], y = coords[, 2], bg = vertex.color,
          stars = cbind(1.2*vertex.size, vertex.size, 1.2*vertex.size, vertex.size),
          add = TRUE, inches = FALSE)
}
# clips as a circle
add_shape("rhombus", clip = shapes("circle")$clip,
          plot = myrhombus)

## Function for plotting an elliptical node
myellipse <- function(coords, v=NULL, params) {
  vertex.color <- params("vertex", "color")
  if (length(vertex.color) != 1 && !is.null(v)) {
    vertex.color <- vertex.color[v]
  }
  vertex.size <- 1/70 * params("vertex", "size") # largeur des elipse
  if (length(vertex.size) != 1 && !is.null(v)) {
    vertex.size <- vertex.size[v]
  }

  draw.ellipse(x=coords[,1], y=coords[,2],
    a = vertex.size, b=0.012, col=vertex.color) #b = hauteur
}

## Register the shape with igraph
add_shape("ellipse", clip=shapes("circle")$clip,
                 plot=myellipse)

# Fonction pour dessiner un hexagone
myhexagon <- function(coords, v = NULL, params) {
  vertex.color <- params("vertex", "color")
  if(length(vertex.color) != 1 && !is.null(v)){
    vertex.color <- vertex.color[v]
  }
  vertex.size <- 1/200 * params("vertex", "size")
  if(length(vertex.size) !=1 && !is.null(v)){
    vertex.size <- vertex.size[v]
  }

  for(i in seq_len(nrow(coords))){
    x <- coords[i, 1]
    y <- coords[i, 2]
    size <- vertex.size[i]
    theta <- seq(0, 2*pi, length.out = 7)[-7]  # Angles pour les sommets de l'hexagone
    xs <- x + size * cos(theta)
    ys <- y + size * sin(theta)
    polygon(x = xs, y = ys, col = vertex.color[i], border = "black")
  }
}

hexagon_clip <- function(coords, el, params) {
  vertex.size <- 1/200 * params("vertex", "size")
  angles <- seq(0, 2 * pi, length.out = 7)[-7]

  clip_coords <- sapply(1:nrow(coords), function(i) {
    size <- vertex.size[i]
    angle <- atan2(el[2, 2] - el[1, 2], el[2, 1] - el[1, 1])  # Edge direction
    idx <- which.min(abs(angles - angle))  # Find closest angle in hexagon
    x <- coords[i, 1] + size * cos(angles[idx])
    y <- coords[i, 2] + size * sin(angles[idx])
    c(x, y)
  })

  return(t(clip_coords))
}

# Enregistrer la nouvelle forme avec igraph
add_shape("hexagon", clip = shapes("circle")$clip, plot = myhexagon)

# Set edge width based on weight:
#E(net)$width <- -log10(1-abs(E(net)$slope_modif))*1.5# avt /6
E(net)$width <- (abs(E(net)$cor)-min(abs(E(net)$cor))+0.01)*5# avt /6

#change arrow size and edge color:
E(net)$arrow.size <- .2
E(net)$edge.color <- "gray80"

# Let's color the edges of the graph based on their source node color.
# We'll get the starting node for each edge with "get.edges"
edge.start <- get.edges(net, 1:ecount(net))[,1]
edge.col <- V(net)$color[edge.start]

vec_color=ifelse(E(net)$cor>0,"red","blue")
edge.col <-vec_color #mettre red if is upper green if is lower 0

l=layout_in_circle
l=layout_with_kk # Kamada Kawai #Like Fruchterman Reingold, it attempts to minimize the energy in a spring system.
l=layout_on_grid
l=layout_with_fr #Fruchterman-Reingold
l <-layout_with_fr(net, niter = 10, grid = "nogrid")

for (sign in c("sign", "unsign")){
  if (sign =="unsign"){
  scaling_factor_attraction_force  <- 3 # A higher value for this scaling factor amplifies the effect of the weights, making heavily weighted edges pull nodes closer together more strongly.
  l <- igraph::layout_with_fr(net, weights = scaling_factor_attraction_force*(igraph::E(net)$weight-min(igraph::E(net)$weight)+0.001), niter = 50, grid = "nogrid") # without sign effect
  }else if (sign == "sign"){
    scaling_factor_attraction_force  <- 1.1 # A higher value for this scaling factor amplifies the effect of the weights, making heavily weighted edges pull nodes closer together more strongly.
    l <- igraph::layout_with_fr(net, weights = scaling_factor_attraction_force*(igraph::E(net)$cor+abs(min(igraph::E(net)$cor))+0.001), niter = 50, grid = "nogrid") # withe sign effect
  }
  # export
  # svg(width=15, height=15,filename = here::here(paste0("report/all_",sign,".svg")))
  # set.seed(1)
  # plot(net, layout = l, 
  #      edge.lty = 1,
  #      edge.arrow.size = 0,
  #      vertex.label.cex = 0.5,
  #             vertex.label = NA,         # Suppression des labels
  #      vertex.size = .1,           # Taille des sommets fixée très petite (proche de 0)
  #      edge.color = edge.col,
  #      edge.curved = .15)
  # title(main = "Network plot based on MLM results show correlation between variables")
  # mtext(paste0("Positive correlation: ",length(df_select %>% filter(cor>0) %>% pull(slope))),
  #       side = 1, line = 2, cex = 0.8,col = "red")
  # mtext(paste0("Negative correlation: ",length(df_select %>% filter(cor<0) %>% pull(slope))),
  #       side = 1, line = 1, cex = 0.8,col = "blue")
  # 
  # dev.off()
  
  png(filename = here::here(paste0("report/microarray/", sign, ".png")), width = 4500, height = 4500, res = 600)
  set.seed(1)
  plot(net, layout = l, 
       edge.lty = 1,
       edge.arrow.size = 0,
       vertex.label.cex = 0.5,
       vertex.label = NA,         # Suppression des labels
       vertex.size = 0.1,         # Taille des sommets très petite
       edge.color = edge.col,
       edge.curved = 0.15)
  title(main = "Network plot based on MLM results show correlation between variables")
  mtext(paste0("Positive correlation: ", length(df_select %>% filter(cor > 0) %>% pull(slope))),
        side = 1, line = 2, cex = 0.8, col = "red")
  mtext(paste0("Negative correlation: ", length(df_select %>% filter(cor < 0) %>% pull(slope))),
        side = 1, line = 1, cex = 0.8, col = "blue")
  dev.off()
}

# save plot to analyse in cytoskape
write_graph(net, file = here::here("data/multi_omics/output/B/network.graphml"), format = "graphml")

## Degree Distribution Histogram

df_deg=tibble(nodes=names(deg),deg=as.vector(deg)) %>%
  full_join(nodes %>% as.data.frame() %>% dplyr::rename(nodes=id),.,by="nodes") %>%
  arrange(desc(deg)) %>% 
  #mutate(variable_cleaned = ifelse(variable_type==1|variable_type==4,gsub("[0-9_]", "", variable),variable)) %>%
  # mutate(variable_cleaned = gsub("[_]","",variable_cleaned)) %>%
  # mutate(variable_cleaned = gsub("(concentration)", "", variable_cleaned)) %>% 
  # mutate(variable_cleaned = gsub("(stem)", "", variable_cleaned)) %>% 
  # mutate(variable_cleaned = gsub("(root)", " ", variable_cleaned)) %>% 
  # mutate(variable_cleaned = gsub("(leaf)", "  ", variable_cleaned)) %>% 
  mutate(organe=factor(organe, levels = c("leaf", "stem", "root","autre", "sEUpE"))) %>% 
  #mutate(variable_cleaned = gsub("sEUpE", "sEUE ", variable_cleaned)) %>% 
  mutate(variable_cleaned = paste0(variable_cleaned, " ", str_to_title(organe)))
  # mutate(variable_cleaned = gsub("(stem|leaf|root|concentration)", "", variable_cleaned))

px=df_deg %>% 
  filter(deg > 4) %>% 
  ggplot(aes(x = deg, y = reorder(variable_cleaned, -deg), fill = organe, label = deg)) +
  geom_col(position = 'dodge') +  
  geom_text(hjust = -0.4) +
  scale_fill_manual(values = colrs2, breaks = c("leaf","stem","root", "autre", "sEUpE")) +
  # scale_x_continuous(limits = c(0, 15),breaks = seq(0, 15, by = 1))+
  theme_minimal() +
  theme(axis.title.y = element_blank(),
        axis.text.x = element_blank(),
        axis.ticks.y = element_blank(),
        panel.grid = element_blank(),
        legend.position = "none")+
  labs(x="Degree", fill="Organe")

ggsave(here::here("report/multi_omics/plot/LMM_metabo_iono/degree_all.svg"), px, height = 15,width = 10)

Sign coregulation network

Unsign coregulation network

14.2.4.2 List of genes strongly connected to psult4

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

df_result_select_psult4 <- df_result_select %>% 
  filter(if_any(c(V1, V2), ~ . == "PsCam042688"))

write_csv(x = df_result_select_psult4, file = here::here("data/microarray/output/cor_08/genes_conected_to_sult4.csv"))

df_result_select %>% 
  arrange(AICc) %>% 
  filter(cor<.82, r2m>.95) %>% 
 head(10) 

# # For Sult 3 
# df_result_select_test <- df_result_select %>% 
#   filter(if_any(c(V1, V2), ~ . == "PsCam009589")) %>% 
#   arrange(abs(cor))

# cor_visualisation(gene_1 = "PsCam009589", gene_2 = "PsCam044134")
# cor_visualisation(gene_1 = "PsCam009589", gene_2 = "PsCam034863") # brest r2c
# cor_visualisation(gene_1 = "PsCam009589", gene_2 = "PsCam052646")

Unfortunately, there are only 9 genes using a filter of 0.8 on the correlation before the mixed model.

read_csv(here::here("data/microarray/output/cor_08/genes_conected_to_sult4.csv")) %>% 
knitr::kable(., caption = "List of all genes connected to SULT4 in the co-regulation network")

List of all genes connected to SULT4 in the co-regulation network

V1	V2	AICc	intercept	slope	pval	residuals	CI_sup05	CI_inf05	CI_sup01	CI_inf01	CI_sup001	CI_inf001	r2m	r2c	cor	between_IC05	between_IC01	between_IC001	fdr	log10fdr	color
PsCam053023	PsCam042688	57.96303	-0.0396126	-1.1395659	0.0000000	0.0753240	-0.9919336	-1.2871982	-0.9455441	-1.3335877	-0.8917103	-1.3874215	0.6912317	0.9338230	-0.8298061	no	no	no	0.0000000	14.559175	green
PsCam006543	PsCam042688	60.64571	-0.0776213	-1.1387179	0.0000000	0.0779787	-0.9858825	-1.2915532	-0.9378582	-1.3395776	-0.8821270	-1.3953087	0.6717523	0.9310115	-0.8013309	no	no	no	0.0000000	14.189943	green
PsCam033648	PsCam042688	61.64654	0.0106030	1.0609829	0.0000000	0.0795595	1.2169167	0.9050492	1.2659146	0.8560513	1.3227756	0.7991903	0.6766409	0.9160831	0.8222350	no	no	no	0.0000000	13.318080	green
PsCam051429	PsCam042688	81.94363	0.0147925	-0.8247786	0.0000000	0.1025982	-0.6236898	-1.0258673	-0.5605031	-1.0890540	-0.4871765	-1.1623806	0.6568484	0.7308065	-0.8217331	no	no	no	0.0000000	8.070045	green
PsCam044134	PsCam042688	54.92919	0.0339039	0.9045165	0.0000312	0.1525473	1.2035036	0.6055293	1.2974522	0.5115807	1.4064774	0.4025556	0.6448819	0.9087736	0.8408067	no	no	no	0.0000407	4.390304	green
PsCam024809	PsCam042688	65.31686	-0.0405056	-0.8707506	0.0003232	0.1678012	-0.5418662	-1.1996350	-0.4385232	-1.3029780	-0.3185961	-1.4229051	0.6033741	0.8759641	-0.8115685	no	no	no	0.0003839	3.415820	green
PsCam034922	PsCam042688	81.38056	0.0071993	-0.7331428	0.0030776	0.1746697	-0.3907965	-1.0754890	-0.2832236	-1.1830620	-0.1583876	-1.3078979	0.5238252	0.7478385	-0.8011649	no	no	no	0.0034137	2.466778	green
PsCam054904	PsCam042688	68.98067	-0.0001044	0.8620707	0.0037183	0.1674086	1.1901856	0.5339557	1.2932867	0.4308546	1.4129333	0.3112080	0.5918757	0.8623411	0.8082885	no	no	no	0.0041038	2.386815	green
PsCam006995	PsCam042688	65.80364	0.0223099	0.9185042	0.0405996	0.1475960	1.2077870	0.6292213	1.2986863	0.5383220	1.4041728	0.4328355	0.6444955	0.8783158	0.8173133	no	no	no	0.0424485	1.372138	green

14.2.4.3 Manual verification

The genes most connected with psult4 are shown below. We can see that they are in fact weakly correlated. This may be due to my pre-filter based on correlations.

# load data
load(file = here::here("data/microarray/output/cor_08/tmp0_peasulf_for_LMM_all.RData"))

px <- cor_visualisation(gene_1 = "PsCam053023", gene_2 = "PsCam042688") # best SULT4 

fig_export(here::here(paste0("report/microarray/plot/network_verification/psult4/best_combination")), px, height_i = 4, width_i = 8, res_i = 300)

# two genes are well conected if they are close to that
  
p_exemple_1 <- cor_visualisation(gene_1 = "PsCam070023", gene_2 = "PsCam070031")
fig_export(here::here(paste0("report/microarray/plot/network_verification/exemple/best_combination_ex1")), p_exemple_1, height_i = 4, width_i = 8, res_i = 300)

p_exemple_2 <- cor_visualisation(gene_1 = "PsCam035544", gene_2 = "PsCam053495")
fig_export(here::here(paste0("report/microarray/plot/network_verification/exemple/best_combination_ex2")), p_exemple_2, height_i = 4, width_i = 8, res_i = 300)

Exemple 1

Exemple 2

14.2.5 Test without prefiltered on correlation but by using only genes from venn diagram

There is 13980 genes in the venn diagram in down and up. So is 84233710 posibility. That will take 15 days !!! In the end there were a total of 614 genes co-regulated with Psult4.1.

clean_session() #  Remove data frames and matrices and  Reset plots

# data importation
load(file = here::here("data/microarray/output/venn/tmp2_LMM_peasulf_microarray.RData"))

# parameter
r2c_i = 0.90
r2m_i = 0.6
fdr_lim = 0.05

df_result = df_result_select %>%  
  mutate(fdr = p.adjust(pval, method = "fdr"),
         log10fdr = -log10(fdr)) %>% 
  mutate(color=ifelse(r2c*1.5 - r2m*1.5 > 0.5, "red",
                      ifelse(r2m<r2m_i, ifelse(r2c<r2c_i,"#DD6342","orange"),
                             ifelse(fdr<fdr_lim, "green","forestgreen")
                             )
                      )
         )
# Visualization 
## Use of randomized data to avoid displaying all points
random_result <- df_result %>%
  sample_n(1000)

random_result%>% 
  ggplot(aes(x = r2m, y = -log10(pval), size = cor, col = -log10(pval))) +
  geom_point(alpha = 0.7) +  
  scale_color_viridis_c() + 
  labs(title = "Dispersion diagramme", 
       x = "R2m", 
       y = "-log10(fdr)", 
       color = "-log10(fdr)") +
  theme_minimal()+

random_result %>% 
  ggplot(aes(x = r2c, y = -log10(pval), size = cor, col = -log10(pval))) +
  geom_point(alpha = 0.7) +  
  scale_color_viridis_c() + 
  labs(title = "Dispersion diagramme", 
       x = "R2c", 
       y = "-log10(fdr)", 
       color = "-log10(fdr)") +
  theme_minimal()

random_result %>% 
  ggplot(aes(x = r2m, y = r2c, col = color, size = abs(cor))) +
  geom_point(aes(shape = cor < 0), na.rm = TRUE) +
  scale_shape_manual(values = c(`TRUE` = 17, `FALSE` = 16)) +
  scale_color_manual(values = c("red" = "red", "#DD6342"="#DD6342","orange" = "orange", "green" = "green", "forestgreen" = "forestgreen")) +  # Utilisation de couleurs fixes
  geom_vline(xintercept = r2m_i, linetype = "dashed", color = "red", linewidth = 1) +
  geom_hline(yintercept = r2c_i, linetype = "dashed", color = "red", linewidth = 1)

# filter the result to have variable
df_result_select <- df_result %>% 
  filter(color=="green") %>% 
  arrange(AICc)

cat_col(paste("Number of variable:",length(unique(c(df_result_select$V1,df_result_select$V2))), "\n"), "green")

# export
save(df_result_select, file = here::here("data/microarray/output/venn/tmp2_LMM_peasulf_microarray_bis.RData"))

round(dim(df_result_select)[1]*100/dim(df_result)[1], 2)

14.2.5.1 Analyse results

clean_session() #  Remove data frames and matrices and  Reset plots
load(file = here::here("data/microarray/output/venn/tmp2_LMM_peasulf_microarray_bis.RData"))
load(file = here::here("data/microarray/output/venn/tmp0_peasulf_for_LMM_all.RData"))

# pb <- progress::progress_bar$new(
#   format = "  Calcul [:bar] :percent | Étape :current/:total | Temps écoulé :elapsed | Temps restant :eta",
#   total = nrow(df_result_select),
#   clear = FALSE,
#   width = 100
# )

# Fonction pour calculer la corrélation avec mise à jour de la barre de progression
# calculate_corr <- function(var1, var2) {
#   pb$tick()
#   
#   # Vérifier si les variables existent dans df_compile
#   if (all(c(var1, var2) %in% colnames(df_compile))) {
#     # Calculer la corrélation
#     cor_value <- cor(df_compile[[var1]], df_compile[[var2]], use = "pairwise.complete.obs")
#     return(cor_value)
#   } else {
#     return(NA) # Retourne NA si une variable manque
#   }
# }
# 
# # Utiliser mapply pour vectoriser le calcul
# df_result_select$corr <- mapply(calculate_corr, df_result_select$V1, df_result_select$V2)


# filter the data previously generate

df_result_select_psult4 <- df_result_select %>% 
  filter(if_any(c(V1, V2), ~ . == "PsCam042688")) %>% 
  arrange(AICc)# %>% 
  #filter(AICc<300)

length(df_result_select_psult4$V1)
  
# write_csv(x = df_result_select_psult4, file = here::here("data/microarray/output/venn/genes_conected_to_sult4.csv"))

df_result_select_psult4_info1 <- add_info_onto(unique(c(df_result_select_psult4 %>% pull(V2), df_result_select_psult4 %>% pull(V1)))) %>% 
  #left_join(., df_result_select_psult4 %>% dplyr::rename(PsCam = V1), by = "PsCam") %>% 
  mutate(cluster = "first-order_neighbors")

genes_of_interest_1 <- unique(c(df_result_select_psult4 %>% pull(V1), df_result_select_psult4 %>% pull(V2)))

df_result_select_psult4_2 <- df_result_select %>% 
  filter(if_any(c(V1, V2), ~ . %in% genes_of_interest_1)) %>% 
  filter(AICc <20 & r2c >0.90 & r2m > 0.6) #%>% 

#ggplot(., aes(x = log10fdr, y = AICc))+geom_point()

length(df_result_select_psult4_2$V1)

df_result_select_psult4_info2 <- add_info_onto(unique(c(df_result_select_psult4_2 %>% pull(V1), df_result_select_psult4_2 %>% pull(V2)))) %>% 
  #left_join(., df_result_select_psult4_2 %>% dplyr::rename(PsCam = V1), by = "PsCam") %>% 
  mutate(cluster = "second-order_neighbors")

test = df_result_select_psult4_info2 %>% unnest(cols = term_BP)

#convert_to_PsCam = add_info_onto(type_input = "Psat", vector_i = "Psat2g155480")
# verification 
#px <- cor_visualisation(gene_1 = "PsCam004281", gene_2 = "PsCam042688") ; px # best SULT4 

df_combined <- bind_rows(df_result_select_psult4_info1, df_result_select_psult4_info2)

df_for_color_igraph <- df_combined %>% 
  dplyr::mutate(
    col_perso = case_when(
      PsCam == "PsCam042688" ~ "#ffc90a",
      is.na(term_BP) ~ "white",  # Si term_BP est NA, la couleur est noire
      map_lgl(term_BP, ~ any(. %in% c(""))) ~ "white",
      map_lgl(term_BP, ~ length(.) == 0 | is.null(.)) ~ "white",  # Détecte character(0) ou NULL
      map_lgl(term_BP, ~ any(str_detect(., "photosynthesis|chlorophyll"))) ~ "forestgreen",  
      map_lgl(term_BP, ~ any(str_detect(., "sulphur|sulfur"))) ~ "#ff9c00",
      map_lgl(GENENAME, ~ any(str_detect(., "sulphur|sulfur"))) ~ "#ff9c00",
      map_lgl(term_BP, ~ any(str_detect(., "flavonoid"))) ~ "#8338EC",  
      map_lgl(term_BP, ~ any(str_detect(., "seed"))) ~ "#9B2226",  
      map_lgl(term_BP, ~ any(str_detect(., "nitrogen"))) ~ "#55aaff",  
      map_lgl(term_BP, ~ any(str_detect(., "auxin"))) ~ "#94D2BD",  
      #map_lgl(term_BP, ~ any(str_detect(., "transport"))) ~ "#3A86FF",  
      map_lgl(term_BP, ~ any(str_detect(., "pollen"))) ~ "#d2e347",  
      TRUE ~ "gray50"
    )
  ) %>% 
  dplyr::rename(variable = PsCam) %>% 
  distinct(variable, col_perso, .keep_all = TRUE)

df_combine_igraph = bind_rows(df_result_select_psult4, df_result_select_psult4_2)

# df_combine_igraph %>% filter(V1 == "PsCam004281" & V2 == "PsCam042688" )

# GO terme
source(here::here("src/function/microarray/GO_on_different_group.R")) # function that find GO terme
df_GO = GO_on_different_group(functional_roles = "BP",
                      group_info = df_combined,
                      group = "cluster",pvalueCutoff_i = 0.05,
                      ID = "Psat",
                      top = 20
                      ) ; df_GO

px <- ggplot(df_GO, aes(x = group, 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 \nderegulated in LMM connected \nwith Sult4.1")  ; px

# export 
fig_export(path = "report/microarray/plot/GO/GO_sult4_venn", plot_x = px, height_i = 4, width_i = 4.6, res = 200)

df_select=df_combine_igraph %>% 
  #filter(select=="yes") %>% 
  filter(color=="green") %>% 
  filter(between_IC01=="no")

#hist(df_select$slope_modif)
# creation of the node 

nodes<-data.frame(variable = unique(c(df_select$V1,df_select$V2))) %>% 
  distinct(variable, .keep_all = T) %>% 
  mutate(id = paste0("s", 1:length(unique(c(df_select$V1,df_select$V2))))) %>% 
  relocate(id, .before=variable) %>% 
  left_join(., df_for_color_igraph %>% dplyr::select(-c(id, cluster)), by = "variable") %>% 
  distinct(variable, .keep_all = T)
  
# links
links<-df_select %>% 
  dplyr::select(V1,V2, cor,AICc) %>% 
    dplyr::mutate(weight = abs(cor)) %>%
  dplyr::mutate(weight_repulstion = 1-cor) %>%
  inner_join(nodes %>% dplyr::rename(V1=variable) %>% dplyr::rename(X1=id),.,by="V1") %>% 
  inner_join(nodes %>% dplyr::rename(V2=variable) %>% dplyr::rename(X2=id),.,by="V2") %>% 
  relocate(X1, .before=X2) %>% 
  relocate(V1, .before=V2) 

# Converting the data to an igraph object:
net <- graph_from_data_frame(links, nodes, directed=T) 

save(nodes, links, net,file = here::here("data/microarray/output/venn/tmp3_LMM_peasulf_microarray.RData"))
load(file = here::here("data/microarray/output/venn/tmp3_LMM_peasulf_microarray.RData"))


# Compute node degree (#links) and use it to set node size:
deg <- degree(net, mode="all")
tableau_connexions <- data.frame(id = names(deg), nb_connexions = deg) %>% 
  left_join(., nodes, by = "id") %>% 
  arrange(desc(nb_connexions))

tableau_connexions %>% head(10)

tableau_connexions_info = add_info_onto(vector_i = tableau_connexions %>% head(20) %>% pull("variable"))


# if(fc_or_pval=="pval"){
#   V(net)$size <- V(net)$log10_pval_test_t*6 #ixi modifie la taille des points
# }else if(fc_or_pval=="fc"){
#   V(net)$size <- V(net)$abslog2fc*10 #ixi modifie la taille des points
#   }
# 
# myrhombus <- function(coords, v = NULL, params) {
#   vertex.color <- params("vertex", "color")
#   if (length(vertex.color) != 1 && !is.null(v)) {
#     vertex.color <- vertex.color[v]
#   }
#   vertex.size <- 1/200 * params("vertex", "size")
#   if (length(vertex.size) != 1 && !is.null(v)) {
#     vertex.size <- vertex.size[v]
#   }
# 
#   symbols(x = coords[, 1], y = coords[, 2], bg = vertex.color,
#           stars = cbind(1.2*vertex.size, vertex.size, 1.2*vertex.size, vertex.size),
#           add = TRUE, inches = FALSE)
# }
# # clips as a circle
# add_shape("rhombus", clip = shapes("circle")$clip,
#           plot = myrhombus)
# 
# ## Function for plotting an elliptical node
# myellipse <- function(coords, v=NULL, params) {
#   vertex.color <- params("vertex", "color")
#   if (length(vertex.color) != 1 && !is.null(v)) {
#     vertex.color <- vertex.color[v]
#   }
#   vertex.size <- 1/70 * params("vertex", "size") # largeur des elipse
#   if (length(vertex.size) != 1 && !is.null(v)) {
#     vertex.size <- vertex.size[v]
#   }
# 
#   draw.ellipse(x=coords[,1], y=coords[,2],
#     a = vertex.size, b=0.012, col=vertex.color) #b = hauteur
# }
# 
# ## Register the shape with igraph
# add_shape("ellipse", clip=shapes("circle")$clip,
#                  plot=myellipse)
# 
# # Fonction pour dessiner un hexagone
# myhexagon <- function(coords, v = NULL, params) {
#   vertex.color <- params("vertex", "color")
#   if(length(vertex.color) != 1 && !is.null(v)){
#     vertex.color <- vertex.color[v]
#   }
#   vertex.size <- 1/200 * params("vertex", "size")
#   if(length(vertex.size) !=1 && !is.null(v)){
#     vertex.size <- vertex.size[v]
#   }
# 
#   for(i in seq_len(nrow(coords))){
#     x <- coords[i, 1]
#     y <- coords[i, 2]
#     size <- vertex.size[i]
#     theta <- seq(0, 2*pi, length.out = 7)[-7]  # Angles pour les sommets de l'hexagone
#     xs <- x + size * cos(theta)
#     ys <- y + size * sin(theta)
#     polygon(x = xs, y = ys, col = vertex.color[i], border = "black")
#   }
# }
# 
# hexagon_clip <- function(coords, el, params) {
#   vertex.size <- 1/200 * params("vertex", "size")
#   angles <- seq(0, 2 * pi, length.out = 7)[-7]
# 
#   clip_coords <- sapply(1:nrow(coords), function(i) {
#     size <- vertex.size[i]
#     angle <- atan2(el[2, 2] - el[1, 2], el[2, 1] - el[1, 1])  # Edge direction
#     idx <- which.min(abs(angles - angle))  # Find closest angle in hexagon
#     x <- coords[i, 1] + size * cos(angles[idx])
#     y <- coords[i, 2] + size * sin(angles[idx])
#     c(x, y)
#   })
# 
#   return(t(clip_coords))
# }
# 
# # Enregistrer la nouvelle forme avec igraph
# add_shape("hexagon", clip = shapes("circle")$clip, plot = myhexagon)

# Set edge width based on weight:
#E(net)$width <- -log10(1-abs(E(net)$slope_modif))*1.5# avt /6
E(net)$width <- (abs(E(net)$cor)-min(abs(E(net)$cor))+0.01)*3# avt /6
#E(net)$weight <- E(net)$AICc

#change arrow size and edge color:
E(net)$arrow.size <- .2
E(net)$edge.color <- "gray80"

# Let's color the edges of the graph based on their source node color.
# We'll get the starting node for each edge with "get.edges"
edge.start <- get.edges(net, 1:ecount(net))[,1]
edge.col <- V(net)$color[edge.start]

vec_color=ifelse(E(net)$cor>0,"red","blue")
edge.col <-vec_color #mettre red if is upper green if is lower 0

l=layout_in_circle
l=layout_with_kk # Kamada Kawai #Like Fruchterman Reingold, it attempts to minimize the energy in a spring system.
l=layout_on_grid
#l=layout_with_fr #Fruchterman-Reingold
l <-layout_with_fr(net, niter = 10, grid = "nogrid")

for (sign in c("sign", "unsign")){
  if (sign =="unsign"){
  scaling_factor_attraction_force  <- 1.1 # A higher value for this scaling factor amplifies the effect of the weights, making heavily weighted edges pull nodes closer together more strongly.
  l <- igraph::layout_with_fr(net, weights = scaling_factor_attraction_force*(igraph::E(net)$weight-min(igraph::E(net)$weight)+0.001), niter = 50, grid = "nogrid") # without sign effect
  }else if (sign == "sign"){
    scaling_factor_attraction_force  <- 1.1 # A higher value for this scaling factor amplifies the effect of the weights, making heavily weighted edges pull nodes closer together more strongly.
    l <- igraph::layout_with_fr(net, weights = scaling_factor_attraction_force*(igraph::E(net)$cor+abs(min(igraph::E(net)$cor))+0.001), niter = 50, grid = "nogrid") # withe sign effect
  }
  # export
  # svg(width=15, height=15,filename = here::here(paste0("report/all_",sign,".svg")))
  # set.seed(1)
  # plot(net, layout = l, 
  #      edge.lty = 1,
  #      edge.arrow.size = 0,
  #      vertex.label.cex = 0.5,
  #             vertex.label = NA,         # Suppression des labels
  #      vertex.size = .1,           # Taille des sommets fixée très petite (proche de 0)
  #      edge.color = edge.col,
  #      edge.curved = .15)
  # title(main = "Network plot based on MLM results show correlation between variables")
  # mtext(paste0("Positive correlation: ",length(df_select %>% filter(cor>0) %>% pull(slope))),
  #       side = 1, line = 2, cex = 0.8,col = "red")
  # mtext(paste0("Negative correlation: ",length(df_select %>% filter(cor<0) %>% pull(slope))),
  #       side = 1, line = 1, cex = 0.8,col = "blue")
  # 
  # dev.off()
  
  png(filename = here::here(paste0("report/microarray/plot/coregulation_network/", sign, ".png")), width = 4500*1.5, height = 4500, res = 800)
  set.seed(1)
  plot(net, layout = l, 
       edge.lty = 1,
       vertex.frame.width = 0.5,
       edge.arrow.size = 0,
       vertex.label.cex = 0.4,
       vertex.label = NA,         # Suppression des labels
       vertex.size = ifelse(V(net)$col_perso == "white",1.5,2.8),         # Taille des sommets très petite
       edge.color = edge.col,
       vertex.color = V(net)$col_perso,
       edge.curved = 0.15)
  title(main = "Network plot based on MLM results show correlation between variables")
  mtext(paste0("Positive correlation: ", length(df_select %>% filter(cor > 0) %>% pull(slope))),
        side = 1, line = 2, cex = 0.8, col = "red")
  mtext(paste0("Negative correlation: ", length(df_select %>% filter(cor < 0) %>% pull(slope))),
        side = 1, line = 1, cex = 0.8, col = "blue")
  dev.off()
}

# save plot to analyse in cytoskape
write_graph(net, file = here::here("data/multi_omics/output/B/network.graphml"), format = "graphml")

## Degree Distribution Histogram

df_deg=tibble(nodes=names(deg),deg=as.vector(deg)) %>%
  full_join(nodes %>% as.data.frame() %>% dplyr::rename(nodes=id),.,by="nodes") %>%
  arrange(desc(deg)) %>% 
  #mutate(variable_cleaned = ifelse(variable_type==1|variable_type==4,gsub("[0-9_]", "", variable),variable)) %>%
  # mutate(variable_cleaned = gsub("[_]","",variable_cleaned)) %>%
  # mutate(variable_cleaned = gsub("(concentration)", "", variable_cleaned)) %>% 
  # mutate(variable_cleaned = gsub("(stem)", "", variable_cleaned)) %>% 
  # mutate(variable_cleaned = gsub("(root)", " ", variable_cleaned)) %>% 
  # mutate(variable_cleaned = gsub("(leaf)", "  ", variable_cleaned)) %>% 
  mutate(organe=factor(organe, levels = c("leaf", "stem", "root","autre", "sEUpE"))) %>% 
  #mutate(variable_cleaned = gsub("sEUpE", "sEUE ", variable_cleaned)) %>% 
  mutate(variable_cleaned = paste0(variable_cleaned, " ", str_to_title(organe)))
  # mutate(variable_cleaned = gsub("(stem|leaf|root|concentration)", "", variable_cleaned))

px=df_deg %>% 
  filter(deg > 4) %>% 
  ggplot(aes(x = deg, y = reorder(variable_cleaned, -deg), fill = organe, label = deg)) +
  geom_col(position = 'dodge') +  
  geom_text(hjust = -0.4) +
  scale_fill_manual(values = colrs2, breaks = c("leaf","stem","root", "autre", "sEUpE")) +
  # scale_x_continuous(limits = c(0, 15),breaks = seq(0, 15, by = 1))+
  theme_minimal() +
  theme(axis.title.y = element_blank(),
        axis.text.x = element_blank(),
        axis.ticks.y = element_blank(),
        panel.grid = element_blank(),
        legend.position = "none")+
  labs(x="Degree", fill="Organe")

ggsave(here::here("report/multi_omics/plot/LMM_metabo_iono/degree_all.svg"), px, height = 15,width = 10)

Unsign coregulation network