5 Visualization function

5.1 Visualisation of the density of the intensity in log2 for red and green channel

Code

plot_rg_density <- function(data_list, normalized = F) {
  
  if(normalized == T){
    # Convert to Red and Green matrices
    ID <- data_list[["genes"]][[1]]
    A <- data_list$A  # Average log intensity
    M <- data_list$M  # Log ratio
    
    R <- 2^(A + (M / 2))  # Calculate Red intensities
    G <- 2^(A - (M / 2))  # Calculate Green intensities
    
    rownames(R) <- ID
    rownames(G) <- ID
    # Convert into matrices
    R_matrix <- as.matrix(R)
    G_matrix <- as.matrix(G)
    
    # Combine into a list for the function
    data_list <- list(R = R_matrix, G = G_matrix)
  }
  
  # Load required libraries
  # 1) Extract R and G matrices from the input list.
  #    Each row = gene, each column = sample.
  Rmat <- data_list$R
  Gmat <- data_list$G
  
  # Make sure column names exist (these will be your sample IDs).
  if (is.null(colnames(Rmat))) {
    stop("Matrix R does not have column names (sample IDs).")
  }
  if (is.null(colnames(Gmat))) {
    stop("Matrix G does not have column names (sample IDs).")
  }
  
  # 2) Convert each matrix to a data frame, including the row names as a 'Gene' column
  dfR <- as.data.frame(Rmat)
  dfR$genes <- rownames(Rmat)
  # Pivot to a long format (key = Sample, value = Intensity)
  dfR_long <- gather(dfR, key = "Sample", value = "Intensity", -genes)
  dfR_long$Channel <- "Red"
  
  dfG <- as.data.frame(Gmat)
  dfG$genes <- rownames(Gmat)
  # Pivot to a long format
  dfG_long <- gather(dfG, key = "Sample", value = "Intensity", -genes)
  dfG_long$Channel <- "Green"
  
  # 3) Combine both long data frames
  df_all <- rbind(dfR_long, dfG_long)
  
  # 4) Create the density plot
  #    - We facet by Sample so that each sample has its own panel
  #    - We use color/fill to distinguish Red vs Green
  p <- ggplot(df_all, aes(x = log2(Intensity), color = Channel, group = interaction(Sample, Channel))) +
    geom_density(linewidth = .3) +
    scale_color_manual(values = c("Red" = "red", "Green" = "green")) +
    labs(
      x = expression(log[2](Intensity)),
      y = "Density"
    ) +
    theme_minimal()+
    theme(
      panel.grid = element_blank() # Remove all grid lines
    )
  
  return(p)
}

5.2 Visualisation of vulcanopplot for the ratio stats

Code

## Process your data and add flags for coloring

vulcanoplot_microarray_ratio_stats<- function(comparison, adg.p.Val_i = 0.05, logFC_i = 1, interactive = F){
  comparison_title = paste0(gsub("_Rep_Bio_\\d+", "", gsub("_rep_bio_\\d+", "", str_replace_all(comparison, "_[0-9]+\\.[0-9]+", ""))), " (FDR:", adg.p.Val_i, " ; logFC:",logFC_i,")")
  
  results <- list_ratio_stats_global[[comparison]] %>%
    dplyr::rename(adj.P.Val = BH, 
                  ID = id_probe) %>% 
    mutate(adj.P.Val = ifelse(adj.P.Val == 0, 0.00000000000001,adj.P.Val)) %>% 
    mutate(
      isSignificant = (adj.P.Val < 0.05) & (abs(logFC) > logFC_i),                         # Flag for significance
      
      colorCategory = case_when(
        !isSignificant ~ "Not significant",                      # Non-significant points
        logFC < -abs(logFC_i) ~ "Negative",                                  # Negative fold change
        logFC > abs(logFC_i) ~ "Positive"                                   # Positive fold change
      )
    ) %>%
    arrange(adj.P.Val)
  
  sum(results$isSignificant)
  
  # Title for the plot
  
  # ggplot for static rendering
  p1 <- ggplot(results, aes(x = logFC, y = -log10(adj.P.Val), color = colorCategory)) +
    geom_point(size = 1, alpha = .6) +                                     # Adjust point size
    geom_text_repel(
      data = subset(results, colorCategory == "Highlighted"),    # Label only the highlighted point
      aes(label = sub("^symbols:", "", CAT)),
      color = "green",                                           # Label in green
      size = 3,
      max.overlaps = 10
    ) +
    scale_color_manual(
      values = c(
        "Not significant" = "gray50",                            # Black for non-significant points
        "Negative" = "#6697EA",                                    # Blue for negative fold change
        "Positive" = "#B02428"                                     # Red for positive fold change
      )
    ) +
    theme_minimal() +
    labs(
      title = comparison_title,
      x = expression(log[2]~Fold~Change),
      y = expression(-log[10](adj.P.Val))
    )
  
  # ggplot for interactive rendering with ggplotly
  p1_interactive <- ggplot(results, aes(x = logFC, y = -log10(adj.P.Val), color = colorCategory)) +
    geom_point(aes(text = paste("ID:", ID, "<br>CAT:", CAT, "<br>logFC:", logFC, "<br>P.Value:", adj.P.Val)), size = 1, alpha = .6) +
    scale_color_manual(
      values = c(
        "Not significant" = "gray50",                            # Black for non-significant points
        "Negative" = "#6697EA",                                    # Blue for negative fold change
        "Positive" = "#B02428"                                     # Red for positive fold change
      )
    ) +
    theme_minimal() +
    labs(
      title = comparison_title,
      x = "log2 Fold Change",  # Replace expression with plain text
      y = "-log10(p-value)"
    )
  if(interactive){
    # Convert to interactive plot
    return(ggplotly(p1_interactive, tooltip = "text"))
  }else{
    return(p1_interactive)
  }
}
# names(list_ratio_stats)

5.3 Visualisation of vulcanopplot when i select different condition

Code

vulcanoplot_microarray_selection<- function(diff_by = "condition", levels_comparison_1, levels_comparison_2, logFC_i = 0){
  load(file =here::here("data/microarray/output/normalized_data_microarray_leaf_PeaSulf.RData"))
  
  df_info_sample_clean_compile <- 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) %>% 
    mutate(
      splitted = str_split_fixed(sample_id, "_", 6)
    ) %>%
    mutate(
      genotype_simplify = splitted[,1]
    ) %>%
    dplyr::select(-splitted) %>% 
    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")
    )
  if (!diff_by %in% c("sulfur_condition", "condition", "genotype", "genotype_simplify")) {
    stop("Invalid value for 'diff'. It must be one of 'sulfur_condition', 'condition', 'genotype' or 'genotype_simplify'.")
  }
  
  if(diff_by =="sulfur_condition"){
    df_to_select_s <- df_info_sample_clean_compile %>% 
      filter(sulfur_condition %in% c(levels_comparison_1, levels_comparison_2))
  }else if(diff_by == "condition"){
    df_to_select_s <- df_info_sample_clean_compile %>% 
      filter(condition %in% c(levels_comparison_1, levels_comparison_2))
  }else if(diff_by == "genotype"){
    df_to_select_s <- df_info_sample_clean_compile %>% 
      filter(genotype %in% c(levels_comparison_1, levels_comparison_2))
  }else if(diff_by == "genotype_simplify"){
    df_to_select_s <- df_info_sample_clean_compile %>% 
      filter(genotype_simplify %in% c(levels_comparison_1, levels_comparison_2)) %>% 
      filter(!sample_id %in% c("WT_SD_1.1_Rep1_Red_17", "WT_SS_17.1_Rep1_Green_17", "WT_SD_4.4_Rep2_Red_18", "WT_SS_19.3_Rep2_Green_18", "WT_SD_9.5_Rep3_Red_19", "WT_SS_25.5_Rep3_Green_19", "WT_SD_12.8_Rep4_Red_20", "WT_SS_28.8_Rep4_Green_20"))
  }
  
  levels_comparison = c(levels_comparison_1, levels_comparison_2)
  
  RG.pAq = RG.MA(MA.pAq)
  colnames(RG.pAq$G)<- colnames(RG$G)
  
  RG_all <- cbind(RG.pAq$R, RG.pAq$G)
  
  RG.pAq_sub = RG.pAq
  
  if(diff_by =="sulfur_condition"){
    RG.pAq_sub$R <- RG_all[, df_to_select_s %>% filter(sulfur_condition == levels_comparison[1]) %>% pull(sample_id), drop = FALSE]
    RG.pAq_sub$G <- RG_all[, df_to_select_s %>% filter(sulfur_condition == levels_comparison[2]) %>% pull(sample_id), drop = FALSE]
  }else if(diff_by == "condition"){
    RG.pAq_sub$R <- RG_all[, df_to_select_s %>% filter(condition == levels_comparison[1]) %>% pull(sample_id), drop = FALSE]
    RG.pAq_sub$G <- RG_all[, df_to_select_s %>% filter(condition == levels_comparison[2]) %>% pull(sample_id), drop = FALSE]
  }else if(diff_by == "genotype"){
    RG.pAq_sub$R <- RG_all[, df_to_select_s %>% filter(genotype == levels_comparison[1]) %>% pull(sample_id), drop = FALSE]
    RG.pAq_sub$G <- RG_all[, df_to_select_s %>% filter(genotype == levels_comparison[2]) %>% pull(sample_id), drop = FALSE]
  }else if(diff_by == "genotype_simplify"){
    RG.pAq_sub$R <- RG_all[, df_to_select_s %>% filter(genotype_simplify == levels_comparison[1]) %>% pull(sample_id), drop = FALSE]
    RG.pAq_sub$G <- RG_all[, df_to_select_s %>% filter(genotype_simplify == levels_comparison[2]) %>% pull(sample_id), drop = FALSE]
  }
  
  
  MA.pAq_sub <- MA.RG(RG.pAq_sub)
  
  design <- matrix(1, nrow = length(colnames(MA.pAq_sub$M)), ncol = 1)
  
  colnames(design) <- "Intercept"
  
  # Fit the linear model
  fit <- lmFit(MA.pAq_sub, design)
  
  # Apply empirical Bayes moderation
  fit.eb <- eBayes(fit)
  
  results <- topTable(fit.eb, coef = "Intercept", number = Inf, adjust.method = "BH") %>%
    mutate(
      isSignificant = (adj.P.Val < 0.05) & (abs(logFC) > logFC_i),                         # Flag for significance
      
      colorCategory = case_when(
        ID %in% c("PsCam042688", "PsCam042745") ~ "Highlighted",                     # Highlighted point
        !isSignificant ~ "Not significant",                      # Non-significant points
        logFC < -abs(logFC_i) ~ "Negative",                                  # Negative fold change
        logFC > abs(logFC_i) ~ "Positive"                                   # Positive fold change
      )
    ) %>%
    arrange(adj.P.Val) %>%
    arrange(adj.P.Val) %>%
    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(ID = id_probe),
      by = "ID"
    )
  
  n_gene_significant <- sum(results$isSignificant)
  n_gene_significant_up <-sum(results$isSignificant & results$colorCategory =="Positive")
  n_gene_significant_down <-sum(results$isSignificant & results$colorCategory =="Negative")
  n_gene_significant
  n_gene_significant_up
  n_gene_significant_down
  
  # Title for the plot
  title_i <- paste0(levels_comparison[1], " vs ", levels_comparison[2]," (Normalized)", "(FDR:0.05;logFC:", logFC_i, ")\n", n_gene_significant," genes significant (",n_gene_significant_up," Up;",n_gene_significant_down," Down)" )
  
  # ggplot for static rendering
  p1 <- ggplot(results, aes(x = logFC, y = -log10(P.Value), color = colorCategory)) +
    geom_point(size = 1, alpha = .6) +                                     # Adjust point size
    geom_text_repel(
      data = subset(results, colorCategory == "Highlighted"),    # Label only the highlighted point
      aes(label = sub("^symbols:", "", CAT)),
      color = sulfate_pallet[2],                                           # Label in green
      size = 3,
      max.overlaps = 10
    ) +
    scale_color_manual(
      values = c(
        "Not significant" = "gray50",                            # Black for non-significant points
        "Negative" = "#6697EA",                                    # Blue for negative fold change
        "Positive" = "#B02428",                                     # Red for positive fold change
        "Highlighted" = as.character(sulfate_pallet[2])                                 # Green for highlighted point
      )
    ) +
    theme_minimal() +
    labs(
      title = title_i,
      x = expression(log[2]~Fold~Change),
      y = expression(-log[10](p-value))
    )
  
  # ggplot for interactive rendering with ggplotly
  p1_interactive <- ggplot(results, aes(x = logFC, y = -log10(P.Value), color = colorCategory)) +
    geom_point(aes(text = paste("ID:", ID, "<br>CAT:", CAT, "<br>logFC:", logFC, "<br>P.Value:", P.Value)), size = 1, alpha = .6) +
    scale_color_manual(
      values = c(
        "Not significant" = "gray50",                            # Black for non-significant points
        "Negative" = "#6697EA",                                    # Blue for negative fold change
        "Positive" = "#B02428",                                     # Red for positive fold change
        "Highlighted" = as.character(sulfate_pallet[2])                                 # Green for highlighted point
      )
    ) +
    theme_minimal() +
    labs(
      title = title_i,
      x = "log2 Fold Change",  # Replace expression with plain text
      y = "-log10(p-value)"
    )
  return(p1_interactive)
}