8 XRF analysis

Code

#pkg
library(tidyverse)
library(here)
library(readxl)
library(ggnewscale) # to have two scale_fill
library(ggh4x)
library(ggstats)
library(patchwork)
library("FactoMineR")
library("factoextra")
library("corrplot")

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

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

Bonus: Here you can find the preliminary tests carried out with the XRF.

8.1 Test

8.1.1 Import results

Code

name_files <- list.files(
  path = here("data/XRF/test/"),
  pattern = "\\.csv$", 
  full.names = TRUE
)

df_XRF <- map_dfr(name_files, ~ read_csv(.x, skip = 0, col_names = TRUE, show_col_types = FALSE) %>%
                      mutate(
                        source_file = basename(.x),  # Nom du fichier sans le chemin
                        date = str_extract(source_file, "(?<=Results_)\\d+") %>% as.integer()  # extrait le numéro après "leaf_"
                      )) %>% 
  dplyr::rename(id_xrf = `File #`) %>% 
   separate(`DateTime`, into = c("Date", "Time"), sep = " ",remove = F) %>% 
   mutate(
    Date = as.Date(Date, format = "%m-%d-%Y"),
    Time = format(strptime(Time, format = "%H:%M"), "%H:%M")
  ) %>% 
  select(-c("Operator", "ID", "Field1", "Field2", "Application", "Method","Name","Cal Check", "Multiplier",`Alloy 1`, `Match Qual 1`, `Alloy 2`, `Match Qual 2`, `Alloy 3`, `Match Qual 3`)) %>% 
  select(-contains("Err")) %>% 
  left_join(., read_excel(here::here("data/XRF/test/plant_info_XRF.xlsx"), col_names = T), by = "id_xrf") %>% 
  left_join(., read_excel(here::here("data/plant_info.xlsx"), col_names = T), by = "plant_num") %>% 
  relocate(id_xrf, plant_num, genotype, row, line, depodding, condition, .before = DateTime) 

# export results
write_csv(x = df_XRF, here::here("data/XRF/test/output/test_XRF.csv"))

8.1.2 Analyse

Code

# select data and convert low data into 0
df_XRF <- read_csv(here::here("data/XRF/test/output/test_XRF.csv"), show_col_types = FALSE) %>% 
   dplyr::select(
    where(~ !is.character(.x) || mean(.x == "< LOD", na.rm = TRUE) <= 0.5)
  ) %>% mutate(across(where(is.character), ~ ifelse(.x == "< LOD", "0", .x))) %>%
  mutate(across(where(~ all(grepl("^\\d*\\.?\\d*$", .x))), as.numeric), 
  plant_num = as.factor(plant_num),
         leaf_num = as.factor(leaf_num),
         line = as.factor(line),
         genotype = as.factor(genotype),
         genotype=fct_relevel(genotype, "WT1", "W78*", "WT2", "E568K"), 
         Date = as.factor(Date)
  )

##### Linea mixed model for all genotype ######
contrasts(df_XRF$genotype) <- contr.sum # to say look at the big average only for genotype (juste an other representation)

element <- c("Mg","Si","P","Cl","K","Ca","Cr","Mn","Fe","Ni","Cu", "Zn","Pb")

for (element_i in element){
  cat(element_i, "\n")
  mod1 <- lm(formula = as.formula(paste(element_i, "~ genotype + leaf_num")), data = df_XRF)
  p_x<-ggcoef_model(mod1)+labs(title = paste0(element_i))
  fig_export(here::here(paste0("report/ionomic/plot/XRF/test/ggcoef_model_",element_i)), format = "png", p_x, height_i = 3, width_i = 8, res_i = 600)
}

### show boxplot
df_XRF_v <- df_XRF %>% pivot_longer(
    cols = all_of(element),
    names_to = "element",
    values_to = "value"
  )

all_possibility <- df_XRF_v %>%  distinct(element)

plots <- list()

# Boucle sur chaque combinaison
for (i in seq_len(nrow(all_possibility))) {
  
  element_i <- all_possibility$element[i]
  
  df_select <- df_XRF_v %>% 
    filter(
      element == element_i
    ) %>% 
    drop_na(value) %>% 
    as.data.frame()
  
  # Si aucune donnée n'est disponible, on passe à la combinaison suivante
  if(nrow(df_select) == 0) {
    message("Aucune donnée pour: ", ", ", element_i)
    next
  }
  
  # Définition de l'étiquette de l'axe des ordonnées
  ylab_i <- paste0("% of ", element_i)
  
  # Essayer d'exécuter stat_analyse et capturer les erreurs éventuelles
  res <- tryCatch({
      stat_analyse(
        data = df_select,
        column_value = "value",
        category_variables = "leaf_num",
        grp_var = "genotype",
        show_plot = TRUE,
        outlier_show = FALSE, 
        label_outlier = "id_xrf",
        biologist_stats = TRUE,
        Ylab_i = ylab_i,
        control_conditions = "",
        strip_normale = FALSE
      )
    },
    error = function(e) {
      message("Erreur pour: ", element_i, " -> ", e$message)
      return(NULL)
    })
  
  # Si une erreur s'est produite, on passe à l'itération suivante
  if(is.null(res)) next
  
  # Extraction du plot et ajout des labels pour la légende
  p_plot <- res[["plot"]] + labs(color = "Leaf num", fill = "Leaf num")
  
  plot_name <- paste0(element_i)
  
  fig_export(here::here(paste0("report/ionomic/plot/XRF/test/boxplot/", plot_name)), p_plot, height_i = 4, width_i = 5, res_i = 300,format = "png")
  
  # Stockage du plot dans la liste avec un nom unique
  plots[[plot_name]] <- p_plot
}

# Assemblage de tous les plots avec patchwork
final_plot <- wrap_plots(plots, ncol = 4) +
  plot_layout(guides = "collect") +
  plot_annotation() & theme(legend.position = 'bottom')

# Affichage du plot final
print(final_plot)
fig_export(here::here("report/ionomic/plot/XRF/test/XRF_stats"), final_plot, height_i = 8, width_i = 12, res_i = 600)

8.1.3 Import results

Code

name_files <- list.files(
  path = here("data/XRF/"),
  pattern = "\\.csv$", 
  full.names = TRUE
)

# test = read_csv("C:/Users/cmaslard/OneDrive/Documents/3_work/Dijon_FILEAS/analyse/PeaSulf_2025/data/XRF/Results_20250507.csv", skip = 0, col_names = TRUE, show_col_types = FALSE) %>% mutate(
#                         source_file = basename("C:/Users/cmaslard/OneDrive/Documents/3_work/Dijon_FILEAS/analyse/PeaSulf_2025/data/XRF/Results_20250507.csv"),  # Nom du fichier sans le chemin
#                         date = str_extract(source_file, "(?<=Results_)\\d+") %>% as.integer()  # extrait le numéro après "leaf_"
#                       )
# 
# df_XRF <- map_dfr(name_files, ~ read_csv(.x, skip = 0, col_names = TRUE, show_col_types = FALSE) %>%
#                       mutate(
#                         source_file = basename(.x),  # Nom du fichier sans le chemin
#                         date = str_extract(source_file, "(?<=Results_)\\d+") %>% as.integer()  # extrait le numéro après "leaf_"
#                       )) %>% 
#   dplyr::rename(id_xrf = `File #`) %>% 
#    separate(`DateTime`, into = c("Date", "Time"), sep = " ",remove = F) %>% 
#    mutate(
#     Date = as.Date(Date, format = "%m-%d-%Y"),
#     Time = format(strptime(Time, format = "%H:%M"), "%H:%M")
#   ) %>% 
#   select(-c("Operator", "ID", "Field1", "Field2", "Application", "Method","Name","Cal Check", "Multiplier",`Alloy 1`, `Match Qual 1`, `Alloy 2`, `Match Qual 2`, `Alloy 3`, `Match Qual 3`)) %>% 
#   select(-contains("Err")) %>% 
#   left_join(., read_excel(here::here("data/XRF/plant_info_XRF.xlsx"), col_names = T), by = "id_xrf") %>% 
#   left_join(., read_excel(here::here("data/plant_info.xlsx"), col_names = T), by = "plant_num") %>% 
#   relocate(id_xrf, plant_num, genotype, row, line, depodding, condition, .before = DateTime) 
# 


numeric_cols <- c("Cl", "Cr", "Cu","Mg", "Mn","Ni","Pb", "Si","Zn", "K", "Ca", "P", "Fe")   # complète si besoin

df_XRF <- map_dfr(
  name_files,
  function(path) {

    date_tag <- str_extract(basename(path), "(?<=Results_)\\d{8}")

    ## 1. CSV ----------------------------------------------------------------
    df_csv <- read_csv(
      path,
      col_types = cols(
        .default  = col_character(),
        `File #`  = col_character(),
        DateTime  = col_character()
      ),
      show_col_types = FALSE
    ) %>% 
      mutate(
        source_file = basename(path),
        date_tag    = date_tag
      ) %>% 
      rename(id_xrf = `File #`) %>% 
      separate(DateTime, into = c("Date", "Time"), sep = " ", remove = FALSE) %>% 
      mutate(
        Date = as.Date(Date, format = "%m-%d-%Y"),
        Time = format(strptime(Time, format = "%H:%M"), "%H:%M")
      ) %>% 
      select(-c("Operator", "ID", "Field1", "Field2", "Application", "Method",
                "Name", "Cal Check", "Multiplier",
                `Alloy 1`, `Match Qual 1`, `Alloy 2`, `Match Qual 2`,
                `Alloy 3`, `Match Qual 3`)) %>% 
      select(-contains("Err"))

    ## 2. Feuille Excel ------------------------------------------------------
    df_info <- read_excel(
      here::here("data/XRF/plant_info_XRF.xlsx"),
      sheet = date_tag,
      col_types = "text"                      # tout en texte, pour gérer "< LOD"
    )

    ## 3. Jointure locale + traitement "< LOD" ------------------------------
    df_csv %>% 
      left_join(df_info, by = "id_xrf") %>% 
      mutate(
        across(all_of(numeric_cols), ~ {
          v <- .x
          lod <- v == "< LOD"
          n_lod <- sum(lod, na.rm = TRUE)
          prop_lod <- n_lod / sum(!is.na(v))

          if (n_lod > 0 && prop_lod < 0.5) {
            v[lod] <- "0"                    # condition remplie → 0
          } else {
            v[lod] <- NA_character_          # sinon → NA
          }

          as.numeric(v)                      # conversion finale en numérique
        })
      ) 
  }
) %>% 
  ## 4. Jointure des métadonnées générales ----------------------------------
  left_join(
    read_excel(here::here("data/plant_info.xlsx"), col_names = TRUE) %>% mutate(plant_num = as.character(plant_num)),
    by = "plant_num"
  ) %>% 
  relocate(id_xrf, plant_num, genotype, row, line, depodding, condition,
           .before = DateTime)

# export results
write_csv(x = df_XRF, here::here("data/XRF/output/df_XRF.csv"))

8.1.4 Analyse

Leaf num = N7 is node 7 of the plant (below the first reproductive node. NR1 is the most developed leaf at the top.

Code

# select data and convert low data into 0
df_XRF <- read_csv(here::here("data/XRF/output/df_XRF.csv"), show_col_types = FALSE) %>% 
   dplyr::select(
    where(~ !is.character(.x) || mean(.x == "< LOD", na.rm = TRUE) <= 0.5)
  ) %>% mutate(across(where(is.character), ~ ifelse(.x == "< LOD", "0", .x))) %>%
  mutate(across(where(~ all(grepl("^\\d*\\.?\\d*$", .x))), as.numeric), 
  plant_num = as.factor(plant_num),
         leaf_num = as.factor(leaf_num),
         line = as.factor(line),
         genotype = as.factor(genotype),
         genotype=fct_relevel(genotype, "WT1", "W78*", "WT2", "E568K"), 
         Date = as.factor(Date)
  )

##### Linea mixed model for all genotype ######
contrasts(df_XRF$genotype) <- contr.sum # to say look at the big average only for genotype (juste an other representation)

element <- c("Mg","Si","P","Cl","K","Ca","Cr","Mn","Fe","Ni","Cu", "Zn","Pb")

for (element_i in element){
  cat(element_i, "\n")
  mod1 <- lm(formula = as.formula(paste(element_i, "~ genotype + leaf_num + row + line + Time + Date")), data = df_XRF)
  p_x<-ggcoef_model(mod1)+labs(title = paste0(element_i))
  fig_export(here::here(paste0("report/ionomic/plot/XRF/ggcoef_model_",element_i)), format = "png", p_x, height_i = 5, width_i = 8, res_i = 600)
}

### show boxplot
df_XRF_v <- df_XRF %>% pivot_longer(
    cols = all_of(element),
    names_to = "element",
    values_to = "value"
  )

all_possibility <- df_XRF_v %>%  distinct(element, Date)

plots <- list()

# Boucle sur chaque combinaison
for (i in seq_len(nrow(all_possibility))) {
  
  element_i <- all_possibility$element[i]
  Date_i <- all_possibility$Date[i]
  
  df_select <- df_XRF_v %>% 
    filter(
      element == element_i, 
      Date == Date_i
    ) %>% 
    drop_na(value) %>% 
    as.data.frame()
  
  # Si aucune donnée n'est disponible, on passe à la combinaison suivante
  if(nrow(df_select) == 0) {
    message("Aucune donnée pour: ", ", ", element_i)
    next
  }
  
  # Définition de l'étiquette de l'axe des ordonnées
  ylab_i <- paste0("% of ", element_i)
  
  # Essayer d'exécuter stat_analyse et capturer les erreurs éventuelles
  res <- tryCatch({
      stat_analyse(
        data = df_select,
        column_value = "value",
        category_variables = "leaf_num",
        grp_var = "genotype",
        show_plot = TRUE,
        outlier_show = FALSE, 
        label_outlier = "id_xrf",
        biologist_stats = TRUE,
        Ylab_i = ylab_i,
        control_conditions = "",
        strip_normale = FALSE
      )
    },
    error = function(e) {
      message("Erreur pour: ", element_i, " -> ", e$message)
      return(NULL)
    })
  
  # Si une erreur s'est produite, on passe à l'itération suivante
  if(is.null(res)) next
  
  # Extraction du plot et ajout des labels pour la légende
  p_plot <- res[["plot"]] + labs(color = "Leaf num", fill = "Leaf num", subtitle = Date_i)
  
  plot_name <- paste0(element_i, "_", Date_i)
  
  fig_export(here::here(paste0("report/ionomic/plot/XRF/boxplot/", plot_name)), p_plot, height_i = 4, width_i = 5, res_i = 300,format = "png")
  
  # Stockage du plot dans la liste avec un nom unique
  plots[[plot_name]] <- p_plot
}

# Assemblage de tous les plots avec patchwork
final_plot <- wrap_plots(plots, ncol = 3) +
  plot_layout(guides = "collect") +
  plot_annotation() & theme(legend.position = 'bottom')

# Affichage du plot final
print(final_plot)
fig_export(here::here("report/ionomic/plot/XRF/XRF_stats"), final_plot, height_i = 24, width_i = 12, res_i = 600)

# PCA
df_XRF_select = df_XRF %>% 
  dplyr::select(id_xrf, Mg,Si,P,Cl,K,Ca,Cr,Mn,Fe,Ni,Cu,Zn,Pb, genotype, leaf_num, Date) %>% 
  column_to_rownames("id_xrf") %>% 
  mutate(genotype = as.factor(genotype),
         leaf_num = as.factor(leaf_num), 
         Date = as.factor(Date)
         )

# by leaf

res_pca <- PCA(df_XRF_select, quali.sup = 14:16, graph = FALSE)
  
var <- get_pca_var(res_pca)
PCA_biplot <- fviz_pca_biplot(
  res_pca,
  geom.ind = "point",  # Affiche seulement les points
  col.ind = df_XRF_select %>% pull(Date),  # Couleur des individus selon leaf_num
  shape.ind = df_XRF_select$leaf_num,          # Forme des individus selon leaf_num
  palette = c("#1E8449", "#5E5A93"),           # Palette de couleurs
  addEllipses = TRUE,                          # Ajouter des ellipses
  ellipse.level = 0.85,
  ellipse.type = "norm",
  label = "var",                               # Afficher uniquement les flèches des variables
  repel = TRUE                                 # Évite le chevauchement des labels
) + 
  ggtitle("PCA biplot by leaf type") +
  theme_minimal()
fig_export(here::here("report/ionomic/plot/XRF/PCA/biplot_leaf_date"), format = "png", PCA_biplot, height_i = 6, width_i = 10, res_i = 600)

# by genotype 
res_pca_N7<- PCA(df_XRF_select %>% filter(leaf_num=="N7"), quali.sup = 14:15, graph = F)
  
var <- get_pca_var(res_pca_N7)
PCA_biplot <- fviz_pca_biplot(
  res_pca_N7,
  geom.ind = "point",  # Affiche seulement les points
  col.ind = df_XRF_select %>% filter(leaf_num=="N7") %>% pull(genotype),  # Couleur des individus selon leaf_num
  shape.ind = df_XRF_select %>% filter(leaf_num=="N7") %>% pull(genotype),          # Forme des individus selon leaf_num
  palette = c("#003049", "#780000", "#7FACC7", "#EC323E"),           # Palette de couleurs
  addEllipses = TRUE,                          # Ajouter des ellipses
  ellipse.level = 0.85,
  ellipse.type = "norm",
  label = "var",                               # Afficher uniquement les flèches des variables
  repel = TRUE                                 # Évite le chevauchement des labels
) + 
  ggtitle("PCA biplot by genotype for node 7") +
  theme_minimal()
fig_export(here::here("report/ionomic/plot/XRF/PCA/biplot_N7"), format = "png", PCA_biplot, height_i = 6, width_i = 10, res_i = 600)

res_pca_NR1<- PCA(df_XRF_select %>% filter(leaf_num=="NR1"), quali.sup = 14:15, graph = F)
var <- get_pca_var(res_pca_NR1)
PCA_biplot <- fviz_pca_biplot(
  res_pca_NR1,
  geom.ind = "point",  # Affiche seulement les points
  col.ind = df_XRF_select %>% filter(leaf_num=="NR1") %>% pull(genotype),  # Couleur des individus selon leaf_num
  shape.ind = df_XRF_select %>% filter(leaf_num=="NR1") %>% pull(genotype),          # Forme des individus selon leaf_num
  palette = c("#003049", "#780000", "#7FACC7", "#EC323E"),           # Palette de couleurs
  addEllipses = TRUE,                          # Ajouter des ellipses
  ellipse.level = 0.85,
  ellipse.type = "norm",
  label = "var",                               # Afficher uniquement les flèches des variables
  repel = TRUE                                 # Évite le chevauchement des labels
) + 
  ggtitle("PCA biplot by genotype for the first reproductive node") +
  theme_minimal()

fig_export(here::here("report/ionomic/plot/XRF/PCA/biplot_NR1"), format = "png", PCA_biplot, height_i = 6, width_i = 10, res_i = 600)

--- output: html_document editor_options: chunk_output_type: console --- ## XRF analysis ```{r} #pkg library(tidyverse) library(here) library(readxl) library(ggnewscale) # to have two scale_fill library(ggh4x) library(ggstats) library(patchwork) library("FactoMineR") library("factoextra") library("corrplot") # 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. # 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) ) ``` <details> <summary style="color: #007bff;"> Bonus: Here you can find the preliminary tests carried out with the XRF. </summary> ## Test ### Import results ```{r, eval = F} name_files <- list.files( path = here("data/XRF/test/"), pattern = "\\.csv$", full.names = TRUE ) df_XRF <- map_dfr(name_files, ~ read_csv(.x, skip = 0, col_names = TRUE, show_col_types = FALSE) %>% mutate( source_file = basename(.x), # Nom du fichier sans le chemin date = str_extract(source_file, "(?<=Results_)\\d+") %>% as.integer() # extrait le numéro après "leaf_" )) %>% dplyr::rename(id_xrf = `File #`) %>% separate(`DateTime`, into = c("Date", "Time"), sep = " ",remove = F) %>% mutate( Date = as.Date(Date, format = "%m-%d-%Y"), Time = format(strptime(Time, format = "%H:%M"), "%H:%M") ) %>% select(-c("Operator", "ID", "Field1", "Field2", "Application", "Method","Name","Cal Check", "Multiplier",`Alloy 1`, `Match Qual 1`, `Alloy 2`, `Match Qual 2`, `Alloy 3`, `Match Qual 3`)) %>% select(-contains("Err")) %>% left_join(., read_excel(here::here("data/XRF/test/plant_info_XRF.xlsx"), col_names = T), by = "id_xrf") %>% left_join(., read_excel(here::here("data/plant_info.xlsx"), col_names = T), by = "plant_num") %>% relocate(id_xrf, plant_num, genotype, row, line, depodding, condition, .before = DateTime) # export results write_csv(x = df_XRF, here::here("data/XRF/test/output/test_XRF.csv")) ``` ### Analyse ```{r, eval = F} # select data and convert low data into 0 df_XRF <- read_csv(here::here("data/XRF/test/output/test_XRF.csv"), show_col_types = FALSE) %>% dplyr::select( where(~ !is.character(.x) || mean(.x == "< LOD", na.rm = TRUE) <= 0.5) ) %>% mutate(across(where(is.character), ~ ifelse(.x == "< LOD", "0", .x))) %>% mutate(across(where(~ all(grepl("^\\d*\\.?\\d*$", .x))), as.numeric), plant_num = as.factor(plant_num), leaf_num = as.factor(leaf_num), line = as.factor(line), genotype = as.factor(genotype), genotype=fct_relevel(genotype, "WT1", "W78*", "WT2", "E568K"), Date = as.factor(Date) ) ##### Linea mixed model for all genotype ###### contrasts(df_XRF$genotype) <- contr.sum # to say look at the big average only for genotype (juste an other representation) element <- c("Mg","Si","P","Cl","K","Ca","Cr","Mn","Fe","Ni","Cu", "Zn","Pb") for (element_i in element){ cat(element_i, "\n") mod1 <- lm(formula = as.formula(paste(element_i, "~ genotype + leaf_num")), data = df_XRF) p_x<-ggcoef_model(mod1)+labs(title = paste0(element_i)) fig_export(here::here(paste0("report/ionomic/plot/XRF/test/ggcoef_model_",element_i)), format = "png", p_x, height_i = 3, width_i = 8, res_i = 600) } ### show boxplot df_XRF_v <- df_XRF %>% pivot_longer( cols = all_of(element), names_to = "element", values_to = "value" ) all_possibility <- df_XRF_v %>% distinct(element) plots <- list() # Boucle sur chaque combinaison for (i in seq_len(nrow(all_possibility))) { element_i <- all_possibility$element[i] df_select <- df_XRF_v %>% filter( element == element_i ) %>% drop_na(value) %>% as.data.frame() # Si aucune donnée n'est disponible, on passe à la combinaison suivante if(nrow(df_select) == 0) { message("Aucune donnée pour: ", ", ", element_i) next } # Définition de l'étiquette de l'axe des ordonnées ylab_i <- paste0("% of ", element_i) # Essayer d'exécuter stat_analyse et capturer les erreurs éventuelles res <- tryCatch({ stat_analyse( data = df_select, column_value = "value", category_variables = "leaf_num", grp_var = "genotype", show_plot = TRUE, outlier_show = FALSE, label_outlier = "id_xrf", biologist_stats = TRUE, Ylab_i = ylab_i, control_conditions = "", strip_normale = FALSE ) }, error = function(e) { message("Erreur pour: ", element_i, " -> ", e$message) return(NULL) }) # Si une erreur s'est produite, on passe à l'itération suivante if(is.null(res)) next # Extraction du plot et ajout des labels pour la légende p_plot <- res[["plot"]] + labs(color = "Leaf num", fill = "Leaf num") plot_name <- paste0(element_i) fig_export(here::here(paste0("report/ionomic/plot/XRF/test/boxplot/", plot_name)), p_plot, height_i = 4, width_i = 5, res_i = 300,format = "png") # Stockage du plot dans la liste avec un nom unique plots[[plot_name]] <- p_plot } # Assemblage de tous les plots avec patchwork final_plot <- wrap_plots(plots, ncol = 4) + plot_layout(guides = "collect") + plot_annotation() & theme(legend.position = 'bottom') # Affichage du plot final print(final_plot) fig_export(here::here("report/ionomic/plot/XRF/test/XRF_stats"), final_plot, height_i = 8, width_i = 12, res_i = 600) ``` ```{r, echo=FALSE, results = "asis"} dir_plots <- "plot/XRF/test" # chemin relatif image_files <- list.files( path = dir_plots, pattern = "^ggcoef.*\\.png$", # Ne prend que les fichiers commençant par "Mean" full.names = FALSE # pour ne garder que les noms de fichiers ) tabset_code <- ":::: {.panel-tabset}\n\n" for (file_name in image_files) { file_name_noext <- tools::file_path_sans_ext(file_name) tab_title <- sub("^ggcoef_model_", "", file_name_noext) # Création du panneau tabset_code <- paste0( tabset_code, "## ", tab_title, "\n\n" ) # Construire le chemin vers le fichier PNG et l’hypothétique version HTML png_path <- file.path(dir_plots, file_name) html_name <- sub("\\.png$", ".html", file_name) html_path <- file.path(dir_plots, html_name) # Inclusion de l’image statique tabset_code <- paste0( tabset_code, "![](", png_path, ")\n\n" ) } tabset_code <- paste0(tabset_code, "::::\n") cat(tabset_code) ``` ![](plot/XRF/test/XRF_stats.png) </details> ### Import results ```{r, eval = F} name_files <- list.files( path = here("data/XRF/"), pattern = "\\.csv$", full.names = TRUE ) # test = read_csv("C:/Users/cmaslard/OneDrive/Documents/3_work/Dijon_FILEAS/analyse/PeaSulf_2025/data/XRF/Results_20250507.csv", skip = 0, col_names = TRUE, show_col_types = FALSE) %>% mutate( # source_file = basename("C:/Users/cmaslard/OneDrive/Documents/3_work/Dijon_FILEAS/analyse/PeaSulf_2025/data/XRF/Results_20250507.csv"), # Nom du fichier sans le chemin # date = str_extract(source_file, "(?<=Results_)\\d+") %>% as.integer() # extrait le numéro après "leaf_" # ) # # df_XRF <- map_dfr(name_files, ~ read_csv(.x, skip = 0, col_names = TRUE, show_col_types = FALSE) %>% # mutate( # source_file = basename(.x), # Nom du fichier sans le chemin # date = str_extract(source_file, "(?<=Results_)\\d+") %>% as.integer() # extrait le numéro après "leaf_" # )) %>% # dplyr::rename(id_xrf = `File #`) %>% # separate(`DateTime`, into = c("Date", "Time"), sep = " ",remove = F) %>% # mutate( # Date = as.Date(Date, format = "%m-%d-%Y"), # Time = format(strptime(Time, format = "%H:%M"), "%H:%M") # ) %>% # select(-c("Operator", "ID", "Field1", "Field2", "Application", "Method","Name","Cal Check", "Multiplier",`Alloy 1`, `Match Qual 1`, `Alloy 2`, `Match Qual 2`, `Alloy 3`, `Match Qual 3`)) %>% # select(-contains("Err")) %>% # left_join(., read_excel(here::here("data/XRF/plant_info_XRF.xlsx"), col_names = T), by = "id_xrf") %>% # left_join(., read_excel(here::here("data/plant_info.xlsx"), col_names = T), by = "plant_num") %>% # relocate(id_xrf, plant_num, genotype, row, line, depodding, condition, .before = DateTime) # numeric_cols <- c("Cl", "Cr", "Cu","Mg", "Mn","Ni","Pb", "Si","Zn", "K", "Ca", "P", "Fe") # complète si besoin df_XRF <- map_dfr( name_files, function(path) { date_tag <- str_extract(basename(path), "(?<=Results_)\\d{8}") ## 1. CSV ---------------------------------------------------------------- df_csv <- read_csv( path, col_types = cols( .default = col_character(), `File #` = col_character(), DateTime = col_character() ), show_col_types = FALSE ) %>% mutate( source_file = basename(path), date_tag = date_tag ) %>% rename(id_xrf = `File #`) %>% separate(DateTime, into = c("Date", "Time"), sep = " ", remove = FALSE) %>% mutate( Date = as.Date(Date, format = "%m-%d-%Y"), Time = format(strptime(Time, format = "%H:%M"), "%H:%M") ) %>% select(-c("Operator", "ID", "Field1", "Field2", "Application", "Method", "Name", "Cal Check", "Multiplier", `Alloy 1`, `Match Qual 1`, `Alloy 2`, `Match Qual 2`, `Alloy 3`, `Match Qual 3`)) %>% select(-contains("Err")) ## 2. Feuille Excel ------------------------------------------------------ df_info <- read_excel( here::here("data/XRF/plant_info_XRF.xlsx"), sheet = date_tag, col_types = "text" # tout en texte, pour gérer "< LOD" ) ## 3. Jointure locale + traitement "< LOD" ------------------------------ df_csv %>% left_join(df_info, by = "id_xrf") %>% mutate( across(all_of(numeric_cols), ~ { v <- .x lod <- v == "< LOD" n_lod <- sum(lod, na.rm = TRUE) prop_lod <- n_lod / sum(!is.na(v)) if (n_lod > 0 && prop_lod < 0.5) { v[lod] <- "0" # condition remplie → 0 } else { v[lod] <- NA_character_ # sinon → NA } as.numeric(v) # conversion finale en numérique }) ) } ) %>% ## 4. Jointure des métadonnées générales ---------------------------------- left_join( read_excel(here::here("data/plant_info.xlsx"), col_names = TRUE) %>% mutate(plant_num = as.character(plant_num)), by = "plant_num" ) %>% relocate(id_xrf, plant_num, genotype, row, line, depodding, condition, .before = DateTime) # export results write_csv(x = df_XRF, here::here("data/XRF/output/df_XRF.csv")) ``` ### Analyse Leaf num = N7 is node 7 of the plant (below the first reproductive node. NR1 is the most developed leaf at the top. ```{r, eval = F} # select data and convert low data into 0 df_XRF <- read_csv(here::here("data/XRF/output/df_XRF.csv"), show_col_types = FALSE) %>% dplyr::select( where(~ !is.character(.x) || mean(.x == "< LOD", na.rm = TRUE) <= 0.5) ) %>% mutate(across(where(is.character), ~ ifelse(.x == "< LOD", "0", .x))) %>% mutate(across(where(~ all(grepl("^\\d*\\.?\\d*$", .x))), as.numeric), plant_num = as.factor(plant_num), leaf_num = as.factor(leaf_num), line = as.factor(line), genotype = as.factor(genotype), genotype=fct_relevel(genotype, "WT1", "W78*", "WT2", "E568K"), Date = as.factor(Date) ) ##### Linea mixed model for all genotype ###### contrasts(df_XRF$genotype) <- contr.sum # to say look at the big average only for genotype (juste an other representation) element <- c("Mg","Si","P","Cl","K","Ca","Cr","Mn","Fe","Ni","Cu", "Zn","Pb") for (element_i in element){ cat(element_i, "\n") mod1 <- lm(formula = as.formula(paste(element_i, "~ genotype + leaf_num + row + line + Time + Date")), data = df_XRF) p_x<-ggcoef_model(mod1)+labs(title = paste0(element_i)) fig_export(here::here(paste0("report/ionomic/plot/XRF/ggcoef_model_",element_i)), format = "png", p_x, height_i = 5, width_i = 8, res_i = 600) } ### show boxplot df_XRF_v <- df_XRF %>% pivot_longer( cols = all_of(element), names_to = "element", values_to = "value" ) all_possibility <- df_XRF_v %>% distinct(element, Date) plots <- list() # Boucle sur chaque combinaison for (i in seq_len(nrow(all_possibility))) { element_i <- all_possibility$element[i] Date_i <- all_possibility$Date[i] df_select <- df_XRF_v %>% filter( element == element_i, Date == Date_i ) %>% drop_na(value) %>% as.data.frame() # Si aucune donnée n'est disponible, on passe à la combinaison suivante if(nrow(df_select) == 0) { message("Aucune donnée pour: ", ", ", element_i) next } # Définition de l'étiquette de l'axe des ordonnées ylab_i <- paste0("% of ", element_i) # Essayer d'exécuter stat_analyse et capturer les erreurs éventuelles res <- tryCatch({ stat_analyse( data = df_select, column_value = "value", category_variables = "leaf_num", grp_var = "genotype", show_plot = TRUE, outlier_show = FALSE, label_outlier = "id_xrf", biologist_stats = TRUE, Ylab_i = ylab_i, control_conditions = "", strip_normale = FALSE ) }, error = function(e) { message("Erreur pour: ", element_i, " -> ", e$message) return(NULL) }) # Si une erreur s'est produite, on passe à l'itération suivante if(is.null(res)) next # Extraction du plot et ajout des labels pour la légende p_plot <- res[["plot"]] + labs(color = "Leaf num", fill = "Leaf num", subtitle = Date_i) plot_name <- paste0(element_i, "_", Date_i) fig_export(here::here(paste0("report/ionomic/plot/XRF/boxplot/", plot_name)), p_plot, height_i = 4, width_i = 5, res_i = 300,format = "png") # Stockage du plot dans la liste avec un nom unique plots[[plot_name]] <- p_plot } # Assemblage de tous les plots avec patchwork final_plot <- wrap_plots(plots, ncol = 3) + plot_layout(guides = "collect") + plot_annotation() & theme(legend.position = 'bottom') # Affichage du plot final print(final_plot) fig_export(here::here("report/ionomic/plot/XRF/XRF_stats"), final_plot, height_i = 24, width_i = 12, res_i = 600) # PCA df_XRF_select = df_XRF %>% dplyr::select(id_xrf, Mg,Si,P,Cl,K,Ca,Cr,Mn,Fe,Ni,Cu,Zn,Pb, genotype, leaf_num, Date) %>% column_to_rownames("id_xrf") %>% mutate(genotype = as.factor(genotype), leaf_num = as.factor(leaf_num), Date = as.factor(Date) ) # by leaf res_pca <- PCA(df_XRF_select, quali.sup = 14:16, graph = FALSE) var <- get_pca_var(res_pca) PCA_biplot <- fviz_pca_biplot( res_pca, geom.ind = "point", # Affiche seulement les points col.ind = df_XRF_select %>% pull(Date), # Couleur des individus selon leaf_num shape.ind = df_XRF_select$leaf_num, # Forme des individus selon leaf_num palette = c("#1E8449", "#5E5A93"), # Palette de couleurs addEllipses = TRUE, # Ajouter des ellipses ellipse.level = 0.85, ellipse.type = "norm", label = "var", # Afficher uniquement les flèches des variables repel = TRUE # Évite le chevauchement des labels ) + ggtitle("PCA biplot by leaf type") + theme_minimal() fig_export(here::here("report/ionomic/plot/XRF/PCA/biplot_leaf_date"), format = "png", PCA_biplot, height_i = 6, width_i = 10, res_i = 600) # by genotype res_pca_N7<- PCA(df_XRF_select %>% filter(leaf_num=="N7"), quali.sup = 14:15, graph = F) var <- get_pca_var(res_pca_N7) PCA_biplot <- fviz_pca_biplot( res_pca_N7, geom.ind = "point", # Affiche seulement les points col.ind = df_XRF_select %>% filter(leaf_num=="N7") %>% pull(genotype), # Couleur des individus selon leaf_num shape.ind = df_XRF_select %>% filter(leaf_num=="N7") %>% pull(genotype), # Forme des individus selon leaf_num palette = c("#003049", "#780000", "#7FACC7", "#EC323E"), # Palette de couleurs addEllipses = TRUE, # Ajouter des ellipses ellipse.level = 0.85, ellipse.type = "norm", label = "var", # Afficher uniquement les flèches des variables repel = TRUE # Évite le chevauchement des labels ) + ggtitle("PCA biplot by genotype for node 7") + theme_minimal() fig_export(here::here("report/ionomic/plot/XRF/PCA/biplot_N7"), format = "png", PCA_biplot, height_i = 6, width_i = 10, res_i = 600) res_pca_NR1<- PCA(df_XRF_select %>% filter(leaf_num=="NR1"), quali.sup = 14:15, graph = F) var <- get_pca_var(res_pca_NR1) PCA_biplot <- fviz_pca_biplot( res_pca_NR1, geom.ind = "point", # Affiche seulement les points col.ind = df_XRF_select %>% filter(leaf_num=="NR1") %>% pull(genotype), # Couleur des individus selon leaf_num shape.ind = df_XRF_select %>% filter(leaf_num=="NR1") %>% pull(genotype), # Forme des individus selon leaf_num palette = c("#003049", "#780000", "#7FACC7", "#EC323E"), # Palette de couleurs addEllipses = TRUE, # Ajouter des ellipses ellipse.level = 0.85, ellipse.type = "norm", label = "var", # Afficher uniquement les flèches des variables repel = TRUE # Évite le chevauchement des labels ) + ggtitle("PCA biplot by genotype for the first reproductive node") + theme_minimal() fig_export(here::here("report/ionomic/plot/XRF/PCA/biplot_NR1"), format = "png", PCA_biplot, height_i = 6, width_i = 10, res_i = 600) ``` ```{r, echo=FALSE, results = "asis"} dir_plots <- "plot/XRF" # chemin relatif image_files <- list.files( path = dir_plots, pattern = "^ggcoef.*\\.png$", # Ne prend que les fichiers commençant par "Mean" full.names = FALSE # pour ne garder que les noms de fichiers ) tabset_code <- ":::: {.panel-tabset}\n\n" for (file_name in image_files) { file_name_noext <- tools::file_path_sans_ext(file_name) tab_title <- sub("^ggcoef_model_", "", file_name_noext) # Création du panneau tabset_code <- paste0( tabset_code, "## ", tab_title, "\n\n" ) # Construire le chemin vers le fichier PNG et l’hypothétique version HTML png_path <- file.path(dir_plots, file_name) html_name <- sub("\\.png$", ".html", file_name) html_path <- file.path(dir_plots, html_name) # Inclusion de l’image statique tabset_code <- paste0( tabset_code, "![](", png_path, ")\n\n" ) } tabset_code <- paste0(tabset_code, "::::\n") cat(tabset_code) ``` ![](plot/XRF/XRF_stats.png) ![](plot/XRF/PCA/biplot_leaf_date.png) ![](plot/XRF/PCA/biplot_leaf_type.png) ![](plot/XRF/PCA/biplot_N7.png) ![](plot/XRF/PCA/biplot_NR1.png)