# Script For azrael for LMM for PeaSulf

library(dplyr)
library(tidyverse)
library(lme4)
library(AICcmodavg)
library(performance)
library(foreach)
library(doParallel)
library(MuMIn)
library(igraph)
require("ggrepel")
library(stringr)

#data loading and creation
nb_cores <- 120 #20
n_rows <- 100000 #500000
seuil_correlation <- 0.8

# import df_compile

load(file="R_file/data/tmp0_peasulf_for_LMM_all.RData")

# Calcul of all comparisons with MlM

calculate_critical_value <- function(p_value) {
  # Calculate the corresponding cumulative probability for a two-way p-value
  prob_cumulative <- 1 - p_value / 2
  
  # Find the corresponding z-value in the standard normal distribution
  critical_value <- qnorm(prob_cumulative)
  
  return(critical_value)
}

# Process function and comparison

cat("Creation of all posibility \n")

#vec_variable=colnames(df_compile[6:length(df_compile)])

#cat(paste0("before ", length(vec_variable), "\n"))
#vec_variable = sample(vec_variable,20)
# cat(paste0("after ", length(vec_variable), "\n"))

#comb_vec <- combn(vec_variable, 2, simplify = F)
#comb_vec_df <- do.call(rbind, comb_vec) %>% as.data.frame()

#cat(paste0("prefiltering with pairwise.complete.obs, their is before ", length(comb_vec_df$V1) , " comparison\n"))
tictoc::tic()

cor_matrix <- cor(df_compile[, 6:ncol(df_compile)], use = "pairwise.complete.obs")

tictoc::toc()
cat("time to create cor_matrix")
# Identifier les paires avec une corrélation absolue supérieure à un seuil

cor_pairs <- which(abs(cor_matrix) > seuil_correlation, arr.ind = TRUE)

# Filtrer pour ne garder que les combinaisons uniques
cor_pairs  <- cor_pairs[cor_pairs[, 1] < cor_pairs[, 2], ]
comb_vec_df<- data.frame(
  V1 = colnames(cor_matrix)[cor_pairs[, 1]],
  V2 = colnames(cor_matrix)[cor_pairs[, 2]]
)

cat(paste0("prefiltering with pairwise.complete.obs, their is after ", length(comb_vec_df$V1) , " comparison\n"))

cat(paste0("Begin paralelisation for ", length(comb_vec_df$V1),  " possibility with cut of ",n_rows," rows \n"))

#comb_vec_df <- comb_vec_df[sample(nrow(comb_vec_df), 24000), ]

#cat(paste0("Begin paralelisation on datatet for ", length(comb_vec_df$V1),  " possibility \n"))

list_dfs <- split(comb_vec_df, (seq(nrow(comb_vec_df))-1) %/% n_rows)

for (j in 1:length(list_dfs)){
  registerDoParallel(cores = nb_cores)
  tictoc::tic()
  comb_vec_df_list <- list_dfs[[j]]
  write_csv(file = paste0("test",j,".csv"),comb_vec_df_list)   
  df_result <- foreach(i = 1:length(comb_vec_df_list$V1), .combine = rbind) %dopar% {
    
    library(dplyr)
    library(tidyr)
    library(lme4)
    library(AICcmodavg)
    library(performance)
    library(MuMIn)
    
    V1_x <- comb_vec_df_list[i, 1]
    V2_x <- comb_vec_df_list[i, 2]
    
    df_compile_select <- df_compile %>% drop_na(V1_x, V2_x) %>% as.data.frame()
    
    df_compile_select[,V1_x]<-(df_compile_select[,V1_x]-mean(df_compile_select[,V1_x]))/sd(df_compile_select[,V1_x])
    df_compile_select[,V2_x]<-(df_compile_select[,V2_x]-mean(df_compile_select[,V2_x]))/sd(df_compile_select[,V2_x])
    
    cat(paste0(round(j/length(list_dfs),4)*100,"% - ",round(i/length(comb_vec_df_list$V1),4)*100,"%") ,i, "_", V1_x, "_", V2_x, "\n")
    m1 <- lmer(formula(paste(V1_x, '~', V2_x ," + (1|genotype) + (1|sulfur_condition)")), data = df_compile_select, REML=T,control = lmerControl(check.conv.singular = "ignore"))
    
    AICc_val <- AICc(m1)
    
    df_coeff <- lmerTest:::get_coefmat(m1) %>% as.data.frame()
    intercept <- df_coeff[1,1]
    slope <- df_coeff[2,1]
    pval <- df_coeff[2,5]
    
    result <- summary(m1)
    residuals <- result$coefficients[4]
    
    CI_sup05 <- slope + residuals * calculate_critical_value(0.05)
    CI_inf05 <- slope - residuals * calculate_critical_value(0.05)
    
    CI_sup01 <- slope + residuals * calculate_critical_value(0.01)
    CI_inf01 <- slope - residuals * calculate_critical_value(0.01)
    
    CI_sup001 <- slope + residuals * calculate_critical_value(0.001)
    CI_inf001 <- slope - residuals * calculate_critical_value(0.001)
    
    r2 <- r.squaredGLMM(m1)
    r2m <- r2[1]
    r2c <- r2[2]
    
    cor=cor(df_compile_select[,V1_x], df_compile_select[,V2_x], method = 'pearson')
    
    return(data.frame(V1 = V1_x, V2 = V2_x, 
                      AICc = AICc_val,
                      intercept = intercept, 
                      slope = slope,
                      pval = pval, 
                      residuals = residuals, 
                      CI_sup05 = CI_sup05,
                      CI_inf05 = CI_inf05,
                      CI_sup01 = CI_sup01,
                      CI_inf01 = CI_inf01,
                      CI_sup001 = CI_sup001,
                      CI_inf001 = CI_inf001,
                      r2m = r2m, 
                      r2c = r2c,
		      cor= cor) %>% 
             mutate(between_IC05 = ifelse(CI_sup05 > CI_inf05, 
                                          ifelse(0 >= CI_inf05 & 0 <= CI_sup05, "yes", "no"), 
                                          "no")) %>% 
             mutate(between_IC01 = ifelse(CI_sup01 > CI_inf01, 
                                          ifelse(0 >= CI_inf01 & 0 <= CI_sup01, "yes", "no"), 
                                          "no")) %>% 
             mutate(between_IC001 = ifelse(CI_sup001 > CI_inf001, 
                                           ifelse(0 >= CI_inf001 & 0 <= CI_sup001, "yes", "no"), 
                                           "no"))
           
    )
  }
  cat("stopImplicitCluster")
  tictoc::toc()
  stopImplicitCluster()
  write.csv(df_result, file = paste0("R_file/data/cut_df/tmp1_test1_",j,"_",length(list_dfs),".csv"))
  # write.csv(df_result, file = here::here(paste0("data/multi_omics/output/cut_df/tmp1_test1_",j,"_",length(list_dfs),".csv")))
  
}

cat("before saving \n")
#save(df_result,file = "R_file/data/tmp1_test1.RData")
cat(paste0("End of the code and save time:",tictoc::toc(), "\n"))