diff --git a/.dockerfile/v2.0.5/HTRfit_2.0.5.tar.gz.REMOVED.git-id b/.dockerfile/v2.0.5/HTRfit_2.0.5.tar.gz.REMOVED.git-id
index b2e8b5810d9b9d930e872408cf025f7539d932c4..15dda702d65a072dd16d1e4ed37ea424c90f23e3 100644
--- a/.dockerfile/v2.0.5/HTRfit_2.0.5.tar.gz.REMOVED.git-id
+++ b/.dockerfile/v2.0.5/HTRfit_2.0.5.tar.gz.REMOVED.git-id
@@ -1 +1 @@
-695e8e39fdf4a2a1d592e1297ee21bdf45db28c6
\ No newline at end of file
+df4fa5ce4975ee070539a659d3f18302f9a235d6
\ No newline at end of file
diff --git a/NAMESPACE b/NAMESPACE
index 251eff10ce0a8913c2c0f7f3debfda1bd2379657..79ff882843a9b0406a119a83810ee06ef48641a7 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -34,6 +34,7 @@ export(countMatrix_2longDtf)
 export(counts_plot)
 export(custom_matrix_transform)
 export(detect_categoricals_vars)
+export(detect_row_matx_bellow_threshold)
 export(diagnostic_plot)
 export(drop_randfx)
 export(endsWithDigit)
@@ -70,6 +71,7 @@ export(getColumnWithSampleID)
 export(getCovarianceMatrix)
 export(getData2computeActualFixEffect)
 export(getDataFromMvrnorm)
+export(getDataFromRnorm)
 export(getDataFromUser)
 export(getDispersionComparison)
 export(getDispersionMatrix)
@@ -95,6 +97,7 @@ export(getSettingsTable)
 export(getStandardDeviationInCorrelation)
 export(getTidyGlmmTMB)
 export(getValidDispersion)
+export(get_effects_from_rnorm)
 export(get_eval_data)
 export(get_eval_data_from_dds)
 export(get_eval_data_from_ltmb)
@@ -102,6 +105,7 @@ export(get_inference_dds)
 export(get_label_y_position)
 export(get_mad_left_threshold)
 export(get_mad_user_message)
+export(get_messages_sequencing_depth)
 export(get_ml_metrics_obj)
 export(get_performances_metrics_obj)
 export(get_pr_curve)
@@ -109,6 +113,7 @@ export(get_pr_object)
 export(get_roc_curve)
 export(get_roc_object)
 export(get_rsquare_2plot)
+export(get_scaling_factor)
 export(glance_tmb)
 export(group_logQij_per_genes_and_labels)
 export(handleAnovaError)
@@ -166,8 +171,10 @@ export(subsetGenes)
 export(subset_glance)
 export(tidy_results)
 export(tidy_tmb)
+export(trimmedMeanMvaluesNormalization)
 export(updateParallel)
 export(wald_test)
+export(warning_too_low_mu_ij_row)
 export(wrap_dds)
 export(wrapper_var_cor)
 exportClasses(performance)
diff --git a/R/filtering_fit.R b/R/filtering_fit.R
index f75b4cc137b40ef7d2f81e2d9140167d58c6b3b9..a94abf1ce445407408f6807c7eb24476e6816bb4 100644
--- a/R/filtering_fit.R
+++ b/R/filtering_fit.R
@@ -52,7 +52,8 @@ identifyTopFit <- function(list_tmb, metric = "AIC", filter_method = "mad", keep
   if (filter_method == "mad"){
     lft_threshold <- get_mad_left_threshold(glance_df[[metric]], mad_tolerance)
     get_mad_user_message(metric, glance_df[[metric]], lft_threshold, keep)
-    index2keep <- if (keep == "top") glance_df[[metric]] > lft_threshold else glance_df[[metric]] < lft_threshold
+    index2keep <- if (keep == "top") which(glance_df[[metric]] > lft_threshold) 
+                  else which(glance_df[[metric]] < lft_threshold)
     id2keep <- rownames(glance_df)[index2keep]
     glance_df <- glance_df[ id2keep , ]
   } ## -- feel free to implement other filtering methods .. 
@@ -81,8 +82,8 @@ identifyTopFit <- function(list_tmb, metric = "AIC", filter_method = "mad", keep
 #' # Calculate the left threshold for MAD-based filtering with tolerance 2
 #' get_mad_left_threshold(c(50, 75, 100), 2)
 get_mad_left_threshold <- function(x, tolerance){
-  med <- stats::median(x)
-  mad <- stats::mad(x, center = median(x), 
+  med <- stats::median(x, na.rm = T)
+  mad <- stats::mad(x, center = med, 
                     constant = 1, na.rm = TRUE, 
                     low = FALSE, high = FALSE)
   med - (tolerance * mad)
@@ -111,7 +112,7 @@ get_mad_user_message <- function(var, x , left_threshold, keep) {
   keep_str <- ifelse(keep == "top", "above", "bellow")
   message(paste("2. Values",  keep_str , "the threshold were identified, threshold:",  left_threshold))
   message("3. Summary of selection:")
-  nb_keep <- ifelse(keep == "top", length(x[x > left_threshold]), length(x[x < left_threshold]) )
+  nb_keep <- ifelse(keep == "top", length(x[which(x > left_threshold)]), length(x[which(x < left_threshold)]) )
   sub_msg2 <- paste("values", keep_str,  "the threshold")
   message(paste("-", nb_keep, "out of", length(x), "observations had", sub_msg2, "for the", var, "metric."))
 }
diff --git a/R/mock_rnaseq.R b/R/mock_rnaseq.R
index 1ea3751ae341d21e0873fa1ae19e437b8d2086dc..89ae03eeb62324a8decaf6e18454dfc60842f3ab 100644
--- a/R/mock_rnaseq.R
+++ b/R/mock_rnaseq.R
@@ -49,6 +49,48 @@ generateCountTable <- function(mu_ij_matx_rep, matx_dispersion_rep) {
 }
 
 
+#' Get messages related to sequencing depth
+#'
+#' This function generates informative messages regarding the scaling of counts by sequencing depth.
+#'
+#' @param scaling_factors Scaling factors obtained from sequencing depth
+#' @param threshold_cov_var Threshold coefficient of variation to detect heterogeneity in scaling
+#' @return NULL
+#' @export
+get_messages_sequencing_depth <- function(scaling_factors, threshold_cov_var = 1.5){
+    message("Scaling count table according to sequencing depth: Done")
+    message("INFO: Scaling counts by sequencing depth may exhibit some randomness due to certain parameter combinations, resulting in erratic behavior. This can be minimized by simulating more genes. We advise verifying the simulated sequencing depth to avoid drawing incorrect conclusions.\n")
+    
+    cov_var_sj <- sd(scaling_factors)/mean(scaling_factors)
+    if (cov_var_sj > threshold_cov_var) 
+        message("INFO: Heterogeneity in scaling by sequencing depth has been detected. Simulated effects may be distorted. This can occur when the number of genes is too low.\n")
+  return(NULL)
+}
+
+
+
+#' Emit a warning message for rows with low mu_ij values
+#'
+#' This function emits a warning message if any rows in the mu_ij matrix have all values below a specified threshold.
+#'
+#' @param mu_ij_matrix Matrix of mu_ij values
+#' @param threshold Threshold value
+#' @return NULL
+#' @export
+warning_too_low_mu_ij_row <- function(mu_ij_matrix, threshold = 1 ){
+  n_too_low_row <- length(which(detect_row_matx_bellow_threshold(mu_ij_matrix, threshold)))
+  if (n_too_low_row > 0){
+    msg <- paste("INFO:", n_too_low_row, "genes have all(mu_ij) < 1, indicating very low counts. Consider removing them for future analysis using prepareData2fit with row_cnt_threshold = 10. To detect them, try increasing sequencing depth.\n", 
+                 sep = " ")
+    message(msg)
+  }
+  return(NULL)
+}
+
+    
+ 
+  
+
 #' Perform RNA-seq simulation
 #'
 #' Simulates RNA-seq data based on the input variables.
@@ -60,6 +102,9 @@ generateCountTable <- function(mu_ij_matx_rep, matx_dispersion_rep) {
 #' @param sequencing_depth Sequencing depth
 #' @param basal_expression base expression gene
 #' @param dispersion User-provided dispersion vector (optional)
+#' @param normal_distr Specifies the distribution type for generating effects. Choose between 'univariate' or 'multivariate' (default).
+#' - 'univariate': Effects are drawn independently from univariate normal distributions. 
+#' - 'multivariate': Effects are drawn jointly from a multivariate normal distribution. 
 #' @return List containing the ground truth, counts, and metadata
 #' @export
 #' @examples
@@ -67,10 +112,11 @@ generateCountTable <- function(mu_ij_matx_rep, matx_dispersion_rep) {
 #'              n_genes = 1000, min_replicates = 2,   
 #'               max_replicates = 4)
 mock_rnaseq <- function(list_var, n_genes, min_replicates, max_replicates, sequencing_depth = NULL,  
-                        basal_expression = 0 , dispersion = stats::runif(n_genes, min = 0, max = 1000) ) {
+                        basal_expression = 0 , dispersion = stats::runif(n_genes, min = 0, max = 1000), 
+                        normal_distr = "multivariate") {
   
   ## -- get my effect
-  df_inputSimulation <- getInput2simulation(list_var, n_genes)
+  df_inputSimulation <- getInput2simulation(list_var, n_genes, normal_distr)
   ## -- add column logQij
   df_inputSimulation <- getLog_qij(df_inputSimulation)
   df_inputSimulation <- addBasalExpression(df_inputSimulation, n_genes, basal_expression)
@@ -95,23 +141,20 @@ mock_rnaseq <- function(list_var, n_genes, min_replicates, max_replicates, seque
   matx_dispersion_rep <- replicateMatrix(matx_dispersion, matx_bool_replication)
   
   if (!is.null(sequencing_depth)) {
-    message("Scaling count table according to sequencing depth.")
-    message("Scaling counts by sequencing depth may exhibit some randomness due to certain parameter combinations, resulting in erratic behavior. This can be minimized by simulating more genes. We advise verifying the simulated sequencing depth to avoid drawing incorrect conclusions.")
-    mu_ij_dtf_rep <- scaleCountsTable(mu_ij_matx_rep, sequencing_depth)
+    scaling_factors <- get_scaling_factor(mu_ij_matx_rep, sequencing_depth)
+    invisible(get_messages_sequencing_depth(scaling_factors))
+    mu_ij_dtf_rep <- scaleCountsTable(mu_ij_matx_rep, scaling_factors)
     mu_ij_matx_rep <- as.matrix(mu_ij_dtf_rep)
+    ## -- rescaling effect
+    df_inputSimulation$log_qij_scaled <- df_inputSimulation$log_qij_scaled + log(mean(scaling_factors, na.rm = T))
   }
   
+  invisible(warning_too_low_mu_ij_row(mu_ij_matx_rep))
   matx_countsTable <- generateCountTable(mu_ij_matx_rep, matx_dispersion_rep)
-
   message("Counts simulation: Done")
   
   
   dtf_countsTable <- matx_countsTable %>% as.data.frame()
-  #if (!is.null(sequencing_depth)) {
-  #  message("Scaling count table according to sequencing depth.")
-  #  dtf_countsTable <- scaleCountsTable(dtf_countsTable, sequencing_depth)
-  #}
-  
   checkFractionOfZero(dtf_countsTable)
   
   metaData <- getSampleMetadata(list_var, n_genes, matx_bool_replication)
diff --git a/R/prepare_data2fit.R b/R/prepare_data2fit.R
index 6500bef0b5159f9a79ed98b69da36bd0611ac9fd..03d1f1f27b48f0c6657cf9f5ba8d49b5cab0d41e 100644
--- a/R/prepare_data2fit.R
+++ b/R/prepare_data2fit.R
@@ -117,16 +117,18 @@ custom_matrix_transform <- function(count_matrix, custom_expression) {
 #'
 #' @param countMatrix Count matrix.
 #' @param metadata Metadata data frame.
-#' @param normalization A boolean value indicating whether to apply median ratio normalization. 
-#'                     If \code{TRUE} (default), the counts matrix will be normalized using median ratio normalization. 
-#'                     If \code{FALSE}, no normalization will be applied.
+#' @param response_name String referring to the target variable name that is being modeled and predicted (default: "kij").
+#' @param groupID String referring to the group variable name (default: "geneID").
+#' @param row_threshold Numeric threshold for removing rows with all counts below a specified value. Default 0.
+#'                This filtering is applied before transformation and normalization.
 #' @param transform A custom R expression to apply to each element of the countMatrix. 
 #'                  This expression should be provided as a character string.
 #'                  For example, to apply log transformation, use `"log(x)"`. 
 #'                  Note that `x` represents each element in the countMatrix. See examples for more details.
 #'                  The transformation is applied before normalization (if normalization = \code{TRUE}).
-#' @param response_name String referring to the target variable name that is being modeled and predicted (default: "kij").
-#' @param groupID String referring to the group variable name (default: "geneID").
+#' @param normalization a vector character specifying method to use (default: NULL, possible choices: c('MRN', 'TTM'))
+#'        - MRN: median ratio normalization
+#'        - TMM: Trimmed Mean of M-values
 #' @return Data frame suitable for fitting.
 #' @export
 #' @examples
@@ -138,16 +140,34 @@ custom_matrix_transform <- function(count_matrix, custom_expression) {
 #' # Prepare data for fitting with custom expression
 #' data2fit <- prepareData2fit(mock_data$counts, mock_data$metadata, transform = "sqrt(x + 1)")
 prepareData2fit <- function(countMatrix, metadata, response_name = "kij", 
-                            groupID = "geneID", transform = NULL , normalization = TRUE ) {
+                            groupID = "geneID", row_threshold = 0, transform = NULL , 
+                            normalization = NULL) {
+  
+  stopifnot( length(row_threshold) == 1 && is.numeric(row_threshold) && row_threshold >= 0  )
+  if (row_threshold > 0){          
+      message(paste("INFO: filtering", response_name, "<", row_threshold, sep = " " ))
+      idx <- detect_row_matx_bellow_threshold(countMatrix, threshold = row_threshold)
+      message(paste(length(idx), "genes removed from data.", sep = " "))
+      countMatrix <- countMatrix[ !idx , ]
+  }
   
   ## -- user transform
   if (!is.null(transform))
     countMatrix <- custom_matrix_transform(countMatrix, transform)
   
-  ## -- median ratio normalization
-  if ( isTRUE(normalization) ) {
-      message("INFO: Median ratio normalization.")
-      countMatrix <- medianRatioNormalization(countMatrix)
+  stopifnot(normalization %in% c(NULL, 'MRN', 'TTM'))
+  if(!is.null(normalization)){
+      ## -- median ratio normalization
+      if ( normalization == "MRN")  {
+          message("INFO: Median ratio normalization.")
+          countMatrix <- medianRatioNormalization(countMatrix)
+      }
+       ## -- Trimmed Mean of M-values normalization
+      if ( normalization == "TTM")  {
+          message("INFO: Trimmed Mean of M-values normalization.")
+          countMatrix <- trimmedMeanMvaluesNormalization(countMatrix)
+      }
+    
   }
 
   dtf_countsLong <- countMatrix_2longDtf(countMatrix, response_name)
@@ -206,3 +226,36 @@ medianRatioNormalization <- function(countsMatrix) {
 }
 
 
+
+#' Normalize count data using Trimmed Mean of M-values (TMM) method
+#'
+#' This function normalizes count data using the Trimmed Mean of M-values (TMM) method,
+#' which calculates scale factors to account for differences in library sizes
+#' and RNA composition between samples.
+#'
+#' @param counts_matrix A matrix of count data where rows represent genes and columns represent samples.
+#' @return Normalized count matrix
+#' @export
+trimmedMeanMvaluesNormalization <- function(counts_matrix) {
+  # Check if all counts are 0
+  if (all(counts_matrix == 0)) {
+    # If all counts are 0, return the input matrix unchanged
+    return(counts_matrix)
+  }
+  
+  # Calculate the total count for each sample
+  total_counts <- colSums(counts_matrix)
+  
+  # Calculate the effective library size (geometric mean of total counts)
+  library_sizes <- exp(mean(log(total_counts)))
+  
+  # Calculate the scale factors using TMM method
+  scale_factors <- library_sizes / total_counts
+  
+  # Apply the scale factors to normalize counts
+  normalized_counts <- t(t(counts_matrix) * scale_factors)
+  
+  return(normalized_counts)
+}
+
+
diff --git a/R/sequencing_depth_scaling.R b/R/sequencing_depth_scaling.R
index f1bad838438b9cce52e453f653685d87f7e22005..2800e133ba7a2ac078b50eea160ac313c4060b2e 100644
--- a/R/sequencing_depth_scaling.R
+++ b/R/sequencing_depth_scaling.R
@@ -6,27 +6,40 @@
 #' This function scales a counts table based on the expected sequencing depth per sample.
 #'
 #' @param countsTable A counts table containing raw read counts.
-#' @param seq_depth  sequencing depth vector
+#' @param scalingDepth_factor  sequencing depth factor vector
 #' @return A scaled counts table.
 #'
 #' @export
 #' @examples
 #' mock_data <- list(counts = matrix(c(10, 20, 30, 20, 30, 10, 10, 20, 20, 20, 30, 10), ncol = 4))
-#' scaled_counts <- scaleCountsTable(countsTable = mock_data$counts, 1000000)
+#' scaled_counts <- scaleCountsTable(countsTable = mock_data$counts, 2)
 #'
-scaleCountsTable <- function(countsTable, seq_depth){
-  seq_depth_simu <- colSums(countsTable)
+scaleCountsTable <- function(countsTable, scalingDepth_factor){
+  counts_scaled <- as.data.frame(sweep(as.matrix(countsTable), 2,  scalingDepth_factor, "*"))
+  return(counts_scaled)
+}
 
+
+
+#' Get scaling factor for count normalization
+#'
+#' Calculates the scaling factor for count normalization based on the sequencing depth.
+#'
+#' @param countsTable Matrix or data frame of counts data.
+#' @param seq_depth Numeric vector containing the sequencing depths.
+#' @return Numeric vector of scaling factors for count normalization.
+#' @export
+get_scaling_factor <- function(countsTable, seq_depth){
+  seq_depth_simu <- colSums(countsTable)
+  
   if (length(seq_depth) > length(seq_depth_simu))
     message("INFO: The length of the sequencing_depth vector exceeds the number of samples. Only the first N values will be utilized.")
   if (length(seq_depth) < length(seq_depth_simu))
     message("INFO: The length of the sequencing_depth vector is shorter than the number of samples. Values will be recycled.")
 
   scalingDepth_factor <- suppressWarnings(seq_depth/seq_depth_simu)
-  counts_scaled <- as.data.frame(sweep(as.matrix(countsTable), 2,  scalingDepth_factor, "*"))
-  return(counts_scaled)
+  return(scalingDepth_factor)
 }
 
 
 
-
diff --git a/R/simulation.R b/R/simulation.R
index 6e25b6f7d65a8c0e7468f1c69b2ab549632907e7..9bf7d01e1f643c9c3676df3b0f904ade90bea64a 100644
--- a/R/simulation.R
+++ b/R/simulation.R
@@ -1,11 +1,15 @@
 # WARNING - Generated by {fusen} from dev/flat_full.Rmd: do not edit by hand
 
+
 #' Get input for simulation based on coefficients
 #'
 #' This function generates input data for simulation based on the coefficients provided in the \code{list_var} argument.
 #'
 #' @param list_var A list of variables (already initialized)
 #' @param n_genes Number of genes to simulate (default: 1)
+#' @param normal_distr Specifies the distribution type for generating effects. Choose between 'univariate' or 'multivariate' (default).
+#' - 'univariate': Effects are drawn independently from univariate normal distributions. 
+#' - 'multivariate': Effects are drawn jointly from a multivariate normal distribution. 
 #' @param input2mvrnorm Input to the \code{mvrnorm} function for simulating data from multivariate normal distribution (default: NULL)
 #' @return A data frame with input coefficients for simulation
 #' @export
@@ -13,17 +17,83 @@
 #' # Example usage
 #' list_var <- init_variable()
 #' getInput2simulation(list_var, n_genes = 10)
-getInput2simulation <- function(list_var, n_genes = 1, input2mvrnorm = NULL) {
+getInput2simulation <- function(list_var, n_genes = 1, normal_distr = "multivariate", input2mvrnorm = NULL) {
   
-  # Use default input to mvrnorm if not provided by the user
-  if (is.null(input2mvrnorm)) input2mvrnorm = getInput2mvrnorm(list_var)  
+  stopifnot( normal_distr %in% c("multivariate", "univariate") )
 
-  l_dataFromMvrnorm = getDataFromMvrnorm(list_var, input2mvrnorm, n_genes)
+  if (normal_distr == "multivariate"){
+      if (is.null(input2mvrnorm)) input2mvrnorm = getInput2mvrnorm(list_var)    
+      l_dataFrom_normdistr <- getDataFromMvrnorm(list_var, input2mvrnorm, n_genes)
+    } 
+  if (normal_distr == "univariate"){
+      l_dataFrom_normdistr <- getDataFromRnorm(list_var, n_genes)
+    }
+  
   l_dataFromUser = getDataFromUser(list_var)
-  df_input2simu <- getCoefficients(list_var, l_dataFromMvrnorm, l_dataFromUser, n_genes)
+  
+  df_input2simu <- getCoefficients(list_var, l_dataFrom_normdistr, l_dataFromUser, n_genes)
+  
   return(df_input2simu)
 }
 
+
+
+#' Prepare data using effects from a normal distribution
+#'
+#' Prepares the data by generating effects from a normal distribution for each gene.
+#'
+#' @param list_var A list of variables (already initialized)
+#' @param n_genes Number of genes to generate data for.
+#' @return A dataframe containing gene metadata and effects generated from a normal distribution.
+#' @export
+getDataFromRnorm <- function(list_var, n_genes){
+    ## -- check if all data have been provided by user
+    if (is.null(getInput2mvrnorm(list_var)$covMatrix))
+        return(list())
+
+    metadata <- getGeneMetadata(list_var , n_genes)
+    df_effects <- get_effects_from_rnorm(list_var, metadata)
+    data <- cbind(metadata, df_effects)  
+    return(list(data))
+   
+}
+
+#' Generate effects from a normal distribution
+#'
+#' Generates effects from a normal distribution for each gene.
+#'
+#' @param list_var A list of variables (already initialized)
+#' @param metadata Gene metadata.
+#' @return A dataframe containing effects generated from a normal distribution.
+#' @export
+get_effects_from_rnorm <- function(list_var, metadata){
+  
+  variable_standard_dev <- getGivenAttribute(list_var, attribute = "sd") %>% unlist()
+  interaction_standard_dev <- getGivenAttribute(list_var$interactions, attribute = "sd") %>% unlist()
+  list_stdev <- c(variable_standard_dev, interaction_standard_dev)
+  
+  # -- mu
+  variable_mu <- getGivenAttribute(list_var, attribute = "mu") %>% unlist()
+  interaction_mu <- getGivenAttribute(list_var$interactions, attribute = "mu") %>% unlist()
+  list_mu <- c(variable_mu, interaction_mu)
+  
+  variable_2rnorm <- names(list_stdev)
+  l_effects <- lapply(stats::setNames(variable_2rnorm, variable_2rnorm) , function(var){
+    col_labels <- paste("label", var, sep = "_")
+    list_combinations <- paste(metadata[["geneID"]], metadata[[col_labels]])
+    list_effects <- unique(list_combinations)
+    list_beta <- rnorm(length(list_effects), mean = list_mu[var], sd = list_stdev[var])
+    names(list_beta) <- list_effects
+    unname(list_beta[list_combinations])
+  })
+  
+  df_effects <- do.call("cbind", l_effects)
+  return(df_effects)
+}
+
+
+
+
 #' getCoefficients
 #'
 #' Get the coefficients.
diff --git a/R/utils.R b/R/utils.R
index 147ee8d7b3d35a801d1ca1fd1ac9e5fae901b6cb..7e2bb8f55c728571122d0713e2a0fa1dfa2f90e9 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -43,6 +43,21 @@ first_non_null_index <- function(lst) {
   return(NULL)
 }
 
+
+
+#' Detect rows in a matrix with all values below a given threshold
+#'
+#' This function detects rows in a matrix where all values are below a specified threshold.
+#'
+#' @param matrix The input matrix
+#' @param threshold The threshold value
+#' @return A logical vector indicating rows below the threshold
+#' @export
+detect_row_matx_bellow_threshold <- function(matrix, threshold) {
+             apply(matrix, 1, function(row) all(row < threshold))
+}
+
+
 #' Clean Variable Name
 #'
 #' This function removes digits, spaces, and special characters from a variable name.
diff --git a/dev/flat_full.Rmd b/dev/flat_full.Rmd
index 7af307ee486c87e38b49208f0d7c77b28df4b86e..b15c08e7ed2e3e56b48ae332642a933413e987c6 100644
--- a/dev/flat_full.Rmd
+++ b/dev/flat_full.Rmd
@@ -68,6 +68,21 @@ first_non_null_index <- function(lst) {
   return(NULL)
 }
 
+
+
+#' Detect rows in a matrix with all values below a given threshold
+#'
+#' This function detects rows in a matrix where all values are below a specified threshold.
+#'
+#' @param matrix The input matrix
+#' @param threshold The threshold value
+#' @return A logical vector indicating rows below the threshold
+#' @export
+detect_row_matx_bellow_threshold <- function(matrix, threshold) {
+             apply(matrix, 1, function(row) all(row < threshold))
+}
+
+
 #' Clean Variable Name
 #'
 #' This function removes digits, spaces, and special characters from a variable name.
@@ -477,6 +492,19 @@ test_that("removeDuplicatedWord returns expected output", {
   expect_equal(cleaned_words, c("hello", "world", "programming", "R is great"))
 })
 
+
+# Test for detect_row_matx_bellow_threshold function
+test_that("detect_row_matx_bellow_threshold detects rows below threshold", {
+  # Create a sample matrix
+  matrix <- matrix(c(0.5, 0.7, 1.2, 0.2, 0.9, 0.9), nrow = 2)
+  # Test with threshold 1
+  expect_equal(detect_row_matx_bellow_threshold(matrix, 1), c(FALSE, TRUE))
+  # Test with threshold 0.5
+  expect_equal(detect_row_matx_bellow_threshold(matrix, 0.5), c(FALSE, FALSE))
+  expect_equal(detect_row_matx_bellow_threshold(matrix, 2), c(TRUE, TRUE))
+})
+
+
 test_that("reorderColumns returns expected output",{
     df <- data.frame(A = 1:3, B = 4:6, C = 7:9)
     # Define the desired column order
@@ -1489,12 +1517,16 @@ test_that("set_correlation sets the correlation between variables correctly", {
 ```
 
 ```{r function-simulation , filename = "simulation"}
+
 #' Get input for simulation based on coefficients
 #'
 #' This function generates input data for simulation based on the coefficients provided in the \code{list_var} argument.
 #'
 #' @param list_var A list of variables (already initialized)
 #' @param n_genes Number of genes to simulate (default: 1)
+#' @param normal_distr Specifies the distribution type for generating effects. Choose between 'univariate' or 'multivariate' (default).
+#' - 'univariate': Effects are drawn independently from univariate normal distributions. 
+#' - 'multivariate': Effects are drawn jointly from a multivariate normal distribution. 
 #' @param input2mvrnorm Input to the \code{mvrnorm} function for simulating data from multivariate normal distribution (default: NULL)
 #' @return A data frame with input coefficients for simulation
 #' @export
@@ -1502,17 +1534,83 @@ test_that("set_correlation sets the correlation between variables correctly", {
 #' # Example usage
 #' list_var <- init_variable()
 #' getInput2simulation(list_var, n_genes = 10)
-getInput2simulation <- function(list_var, n_genes = 1, input2mvrnorm = NULL) {
+getInput2simulation <- function(list_var, n_genes = 1, normal_distr = "multivariate", input2mvrnorm = NULL) {
   
-  # Use default input to mvrnorm if not provided by the user
-  if (is.null(input2mvrnorm)) input2mvrnorm = getInput2mvrnorm(list_var)  
+  stopifnot( normal_distr %in% c("multivariate", "univariate") )
 
-  l_dataFromMvrnorm = getDataFromMvrnorm(list_var, input2mvrnorm, n_genes)
+  if (normal_distr == "multivariate"){
+      if (is.null(input2mvrnorm)) input2mvrnorm = getInput2mvrnorm(list_var)    
+      l_dataFrom_normdistr <- getDataFromMvrnorm(list_var, input2mvrnorm, n_genes)
+    } 
+  if (normal_distr == "univariate"){
+      l_dataFrom_normdistr <- getDataFromRnorm(list_var, n_genes)
+    }
+  
   l_dataFromUser = getDataFromUser(list_var)
-  df_input2simu <- getCoefficients(list_var, l_dataFromMvrnorm, l_dataFromUser, n_genes)
+  
+  df_input2simu <- getCoefficients(list_var, l_dataFrom_normdistr, l_dataFromUser, n_genes)
+  
   return(df_input2simu)
 }
 
+
+
+#' Prepare data using effects from a normal distribution
+#'
+#' Prepares the data by generating effects from a normal distribution for each gene.
+#'
+#' @param list_var A list of variables (already initialized)
+#' @param n_genes Number of genes to generate data for.
+#' @return A dataframe containing gene metadata and effects generated from a normal distribution.
+#' @export
+getDataFromRnorm <- function(list_var, n_genes){
+    ## -- check if all data have been provided by user
+    if (is.null(getInput2mvrnorm(list_var)$covMatrix))
+        return(list())
+
+    metadata <- getGeneMetadata(list_var , n_genes)
+    df_effects <- get_effects_from_rnorm(list_var, metadata)
+    data <- cbind(metadata, df_effects)  
+    return(list(data))
+   
+}
+
+#' Generate effects from a normal distribution
+#'
+#' Generates effects from a normal distribution for each gene.
+#'
+#' @param list_var A list of variables (already initialized)
+#' @param metadata Gene metadata.
+#' @return A dataframe containing effects generated from a normal distribution.
+#' @export
+get_effects_from_rnorm <- function(list_var, metadata){
+  
+  variable_standard_dev <- getGivenAttribute(list_var, attribute = "sd") %>% unlist()
+  interaction_standard_dev <- getGivenAttribute(list_var$interactions, attribute = "sd") %>% unlist()
+  list_stdev <- c(variable_standard_dev, interaction_standard_dev)
+  
+  # -- mu
+  variable_mu <- getGivenAttribute(list_var, attribute = "mu") %>% unlist()
+  interaction_mu <- getGivenAttribute(list_var$interactions, attribute = "mu") %>% unlist()
+  list_mu <- c(variable_mu, interaction_mu)
+  
+  variable_2rnorm <- names(list_stdev)
+  l_effects <- lapply(stats::setNames(variable_2rnorm, variable_2rnorm) , function(var){
+    col_labels <- paste("label", var, sep = "_")
+    list_combinations <- paste(metadata[["geneID"]], metadata[[col_labels]])
+    list_effects <- unique(list_combinations)
+    list_beta <- rnorm(length(list_effects), mean = list_mu[var], sd = list_stdev[var])
+    names(list_beta) <- list_effects
+    unname(list_beta[list_combinations])
+  })
+  
+  df_effects <- do.call("cbind", l_effects)
+  return(df_effects)
+}
+
+
+
+
 #' getCoefficients
 #'
 #' Get the coefficients.
@@ -1838,6 +1936,29 @@ getSampleID <- function(list_var){
 ```
 
 ```{r test-simulation}
+
+test_that("getDataFromRnorm generates correct data frame", {
+  input_var_list <- init_variable(name = "varA", mu = 10, sd = 0.1, level = 3) %>%
+                    init_variable(name = "varB", mu = 1, sd = 2, level = 2)
+  metadata <- getGeneMetadata(input_var_list , n_genes = 5)
+  df <- getDataFromRnorm(input_var_list, n_genes = 5)
+  expect_is(df[[1]], "data.frame")
+  expect_equal(nrow(df[[1]]), 30)
+  expect_equal(colnames(df[[1]]), c("geneID", "label_varA", "label_varB", "varA", "varB"))  
+})
+
+test_that("get_effects_from_rnorm generates correct effects", {
+  input_var_list <- init_variable(name = "varA", mu = 10, sd = 0.1, level = 3) %>%
+                    init_variable(name = "varB", mu = 1, sd = 2, level = 2)
+  metadata <- getGeneMetadata(input_var_list , n_genes = 5)
+  df_effects <- get_effects_from_rnorm(input_var_list, metadata)
+  
+  expect_is(df_effects, "matrix")
+  expect_equal(nrow(df_effects), nrow(metadata))
+  expect_equal(colnames(df_effects), c("varA", "varB"))
+})
+
+
 # Test case 1: Check if the function returns a data frame
 test_that("getInput2simulation returns a data frame", {
   list_var <- init_variable()
@@ -2116,6 +2237,48 @@ generateCountTable <- function(mu_ij_matx_rep, matx_dispersion_rep) {
 }
 
 
+#' Get messages related to sequencing depth
+#'
+#' This function generates informative messages regarding the scaling of counts by sequencing depth.
+#'
+#' @param scaling_factors Scaling factors obtained from sequencing depth
+#' @param threshold_cov_var Threshold coefficient of variation to detect heterogeneity in scaling
+#' @return NULL
+#' @export
+get_messages_sequencing_depth <- function(scaling_factors, threshold_cov_var = 1.5){
+    message("Scaling count table according to sequencing depth: Done")
+    message("INFO: Scaling counts by sequencing depth may exhibit some randomness due to certain parameter combinations, resulting in erratic behavior. This can be minimized by simulating more genes. We advise verifying the simulated sequencing depth to avoid drawing incorrect conclusions.\n")
+    
+    cov_var_sj <- sd(scaling_factors)/mean(scaling_factors)
+    if (cov_var_sj > threshold_cov_var) 
+        message("INFO: Heterogeneity in scaling by sequencing depth has been detected. Simulated effects may be distorted. This can occur when the number of genes is too low.\n")
+  return(NULL)
+}
+
+
+
+#' Emit a warning message for rows with low mu_ij values
+#'
+#' This function emits a warning message if any rows in the mu_ij matrix have all values below a specified threshold.
+#'
+#' @param mu_ij_matrix Matrix of mu_ij values
+#' @param threshold Threshold value
+#' @return NULL
+#' @export
+warning_too_low_mu_ij_row <- function(mu_ij_matrix, threshold = 1 ){
+  n_too_low_row <- length(which(detect_row_matx_bellow_threshold(mu_ij_matrix, threshold)))
+  if (n_too_low_row > 0){
+    msg <- paste("INFO:", n_too_low_row, "genes have all(mu_ij) < 1, indicating very low counts. Consider removing them for future analysis using prepareData2fit with row_cnt_threshold = 10. To detect them, try increasing sequencing depth.\n", 
+                 sep = " ")
+    message(msg)
+  }
+  return(NULL)
+}
+
+    
+ 
+  
+
 #' Perform RNA-seq simulation
 #'
 #' Simulates RNA-seq data based on the input variables.
@@ -2127,6 +2290,9 @@ generateCountTable <- function(mu_ij_matx_rep, matx_dispersion_rep) {
 #' @param sequencing_depth Sequencing depth
 #' @param basal_expression base expression gene
 #' @param dispersion User-provided dispersion vector (optional)
+#' @param normal_distr Specifies the distribution type for generating effects. Choose between 'univariate' or 'multivariate' (default).
+#' - 'univariate': Effects are drawn independently from univariate normal distributions. 
+#' - 'multivariate': Effects are drawn jointly from a multivariate normal distribution. 
 #' @return List containing the ground truth, counts, and metadata
 #' @export
 #' @examples
@@ -2134,10 +2300,11 @@ generateCountTable <- function(mu_ij_matx_rep, matx_dispersion_rep) {
 #'              n_genes = 1000, min_replicates = 2,   
 #'               max_replicates = 4)
 mock_rnaseq <- function(list_var, n_genes, min_replicates, max_replicates, sequencing_depth = NULL,  
-                        basal_expression = 0 , dispersion = stats::runif(n_genes, min = 0, max = 1000) ) {
+                        basal_expression = 0 , dispersion = stats::runif(n_genes, min = 0, max = 1000), 
+                        normal_distr = "multivariate") {
   
   ## -- get my effect
-  df_inputSimulation <- getInput2simulation(list_var, n_genes)
+  df_inputSimulation <- getInput2simulation(list_var, n_genes, normal_distr)
   ## -- add column logQij
   df_inputSimulation <- getLog_qij(df_inputSimulation)
   df_inputSimulation <- addBasalExpression(df_inputSimulation, n_genes, basal_expression)
@@ -2162,23 +2329,20 @@ mock_rnaseq <- function(list_var, n_genes, min_replicates, max_replicates, seque
   matx_dispersion_rep <- replicateMatrix(matx_dispersion, matx_bool_replication)
   
   if (!is.null(sequencing_depth)) {
-    message("Scaling count table according to sequencing depth.")
-    message("Scaling counts by sequencing depth may exhibit some randomness due to certain parameter combinations, resulting in erratic behavior. This can be minimized by simulating more genes. We advise verifying the simulated sequencing depth to avoid drawing incorrect conclusions.")
-    mu_ij_dtf_rep <- scaleCountsTable(mu_ij_matx_rep, sequencing_depth)
+    scaling_factors <- get_scaling_factor(mu_ij_matx_rep, sequencing_depth)
+    invisible(get_messages_sequencing_depth(scaling_factors))
+    mu_ij_dtf_rep <- scaleCountsTable(mu_ij_matx_rep, scaling_factors)
     mu_ij_matx_rep <- as.matrix(mu_ij_dtf_rep)
+    ## -- rescaling effect
+    df_inputSimulation$log_qij_scaled <- df_inputSimulation$log_qij_scaled + log(mean(scaling_factors, na.rm = T))
   }
   
+  invisible(warning_too_low_mu_ij_row(mu_ij_matx_rep))
   matx_countsTable <- generateCountTable(mu_ij_matx_rep, matx_dispersion_rep)
-
   message("Counts simulation: Done")
   
   
   dtf_countsTable <- matx_countsTable %>% as.data.frame()
-  #if (!is.null(sequencing_depth)) {
-  #  message("Scaling count table according to sequencing depth.")
-  #  dtf_countsTable <- scaleCountsTable(dtf_countsTable, sequencing_depth)
-  #}
-  
   checkFractionOfZero(dtf_countsTable)
   
   metaData <- getSampleMetadata(list_var, n_genes, matx_bool_replication)
@@ -2364,34 +2528,13 @@ test_that("replicateMatrix replicates matrix correctly", {
   expect_equal(dim(replicated_matrix), c(3, 9))
 })
 
-
-# Test for mock_rnaseq
-#test_that("mock_rnaseq returns expected output", {
-  # Set up input variables
-#  list_var <- NULL
-#  n_genes <- 3
-#  min_replicates <- 2
-#  max_replicates <- 4
-#  df_inputSimulation <- data.frame(gene_id = 1:3, coef_value = c(0.5, 0.3, 0.2))
-#  matx_dispersion <- matrix(1:9, nrow = 3, ncol = 3)
-
-  # Run the function
-#  expect_error(mock_rnaseq(list_var, n_genes, min_replicates, max_replicates, df_inputSimulation, 
-#                           matx_dispersion))
-  
-  
-  #list_var <- init_variable(name = "my_var", mu = c(10, 20), level = 2 )
-  #n_genes <- 10
-  #min_replicates <- 2
-  #max_replicates <- 4
-  #scaling_factor <- 1
-  #df_inputSimulation <- getInput2simulation(list_var, n_genes)
-  #dispersion <- getDispersionMatrix(list_var, n_genes, c(1000, 1000, 1000, 1000, 1000, 1, 1, 1, 1, 1))
-  #mock_rnaseq(list_var, n_genes, min_replicates, 
-  #            max_replicates, 
-  #            df_inputSimulation, dispersion)
-  #ERROOR
-#})
+# Test for warning_too_low_mu_ij_row function
+test_that("warning_too_low_mu_ij_row emits warning for low mu_ij values", {
+  # Create a sample mu_ij matrix
+  mu_ij_matrix <- matrix(c(0.5, 0.7, 1.2, 0.2, 0.9, 1.1), nrow = 2)
+  # Verify if warning message is emitted
+  expect_message(warning_too_low_mu_ij_row(mu_ij_matrix, 2))
+})
 
 
 # Test for generateReplicationMatrix
@@ -2521,16 +2664,18 @@ custom_matrix_transform <- function(count_matrix, custom_expression) {
 #'
 #' @param countMatrix Count matrix.
 #' @param metadata Metadata data frame.
-#' @param normalization A boolean value indicating whether to apply median ratio normalization. 
-#'                     If \code{TRUE} (default), the counts matrix will be normalized using median ratio normalization. 
-#'                     If \code{FALSE}, no normalization will be applied.
+#' @param response_name String referring to the target variable name that is being modeled and predicted (default: "kij").
+#' @param groupID String referring to the group variable name (default: "geneID").
+#' @param row_threshold Numeric threshold for removing rows with all counts below a specified value. Default 0.
+#'                This filtering is applied before transformation and normalization.
 #' @param transform A custom R expression to apply to each element of the countMatrix. 
 #'                  This expression should be provided as a character string.
 #'                  For example, to apply log transformation, use `"log(x)"`. 
 #'                  Note that `x` represents each element in the countMatrix. See examples for more details.
 #'                  The transformation is applied before normalization (if normalization = \code{TRUE}).
-#' @param response_name String referring to the target variable name that is being modeled and predicted (default: "kij").
-#' @param groupID String referring to the group variable name (default: "geneID").
+#' @param normalization a vector character specifying method to use (default: NULL, possible choices: c('MRN', 'TTM'))
+#'        - MRN: median ratio normalization
+#'        - TMM: Trimmed Mean of M-values
 #' @return Data frame suitable for fitting.
 #' @export
 #' @examples
@@ -2542,16 +2687,34 @@ custom_matrix_transform <- function(count_matrix, custom_expression) {
 #' # Prepare data for fitting with custom expression
 #' data2fit <- prepareData2fit(mock_data$counts, mock_data$metadata, transform = "sqrt(x + 1)")
 prepareData2fit <- function(countMatrix, metadata, response_name = "kij", 
-                            groupID = "geneID", transform = NULL , normalization = TRUE ) {
+                            groupID = "geneID", row_threshold = 0, transform = NULL , 
+                            normalization = NULL) {
+  
+  stopifnot( length(row_threshold) == 1 && is.numeric(row_threshold) && row_threshold >= 0  )
+  if (row_threshold > 0){          
+      message(paste("INFO: filtering", response_name, "<", row_threshold, sep = " " ))
+      idx <- detect_row_matx_bellow_threshold(countMatrix, threshold = row_threshold)
+      message(paste(length(idx), "genes removed from data.", sep = " "))
+      countMatrix <- countMatrix[ !idx , ]
+  }
   
   ## -- user transform
   if (!is.null(transform))
     countMatrix <- custom_matrix_transform(countMatrix, transform)
   
-  ## -- median ratio normalization
-  if ( isTRUE(normalization) ) {
-      message("INFO: Median ratio normalization.")
-      countMatrix <- medianRatioNormalization(countMatrix)
+  stopifnot(normalization %in% c(NULL, 'MRN', 'TTM'))
+  if(!is.null(normalization)){
+      ## -- median ratio normalization
+      if ( normalization == "MRN")  {
+          message("INFO: Median ratio normalization.")
+          countMatrix <- medianRatioNormalization(countMatrix)
+      }
+       ## -- Trimmed Mean of M-values normalization
+      if ( normalization == "TTM")  {
+          message("INFO: Trimmed Mean of M-values normalization.")
+          countMatrix <- trimmedMeanMvaluesNormalization(countMatrix)
+      }
+    
   }
 
   dtf_countsLong <- countMatrix_2longDtf(countMatrix, response_name)
@@ -2610,6 +2773,39 @@ medianRatioNormalization <- function(countsMatrix) {
 }
 
 
+
+#' Normalize count data using Trimmed Mean of M-values (TMM) method
+#'
+#' This function normalizes count data using the Trimmed Mean of M-values (TMM) method,
+#' which calculates scale factors to account for differences in library sizes
+#' and RNA composition between samples.
+#'
+#' @param counts_matrix A matrix of count data where rows represent genes and columns represent samples.
+#' @return Normalized count matrix
+#' @export
+trimmedMeanMvaluesNormalization <- function(counts_matrix) {
+  # Check if all counts are 0
+  if (all(counts_matrix == 0)) {
+    # If all counts are 0, return the input matrix unchanged
+    return(counts_matrix)
+  }
+  
+  # Calculate the total count for each sample
+  total_counts <- colSums(counts_matrix)
+  
+  # Calculate the effective library size (geometric mean of total counts)
+  library_sizes <- exp(mean(log(total_counts)))
+  
+  # Calculate the scale factors using TMM method
+  scale_factors <- library_sizes / total_counts
+  
+  # Apply the scale factors to normalize counts
+  normalized_counts <- t(t(counts_matrix) * scale_factors)
+  
+  return(normalized_counts)
+}
+
+
 ```
 
 ```{r  test-prepareData2fit}
@@ -2654,7 +2850,7 @@ test_that("prepareData2fit prepares data for fitting", {
   mock_data <- mock_rnaseq(list_var, n_genes = 3, 2,2)
   
   # Prepare data for fitting
-  data2fit <- prepareData2fit(mock_data$counts, mock_data$metadata)
+  data2fit <- prepareData2fit(mock_data$counts, mock_data$metadata, normalization = 'MRN')
   
   expect_true(is.character(data2fit$sampleID))
   expect_true(is.character(data2fit$geneID))
@@ -2669,13 +2865,13 @@ test_that("prepareData2fit prepares data for fitting", {
   metadata <- data.frame(sampleID = colnames(countMatrix), condition = rep(c("A", "B"), each = 5))
   
   # Call the function with log transformation
-  data2fit <- prepareData2fit(countMatrix, metadata, normalization = F ,transform = "log(x)")
+  data2fit <- prepareData2fit(countMatrix, metadata, normalization = NULL,transform = "log(x)")
   
   # Test if log transformation has been applied correctly
   expect_equal(unique(data2fit$kij), log(c(1, 2, 3, 4)))
   
   # Call the function with sqrt transformation
-  data2fit <- prepareData2fit(countMatrix, metadata, normalization = F , transform = "sqrt(x)")
+  data2fit <- prepareData2fit(countMatrix, metadata, normalization = NULL , transform = "sqrt(x)")
   # Test if log transformation has been applied correctly
   expect_equal(unique(data2fit$kij), sqrt(c(1, 2, 3, 4)))
 })
@@ -2734,6 +2930,35 @@ test_that("Throws an error when all-zero genes are encountered", {
 })
 
 
+
+# Test case 1: Check if normalization is performed correctly
+test_that("Normalized counts are calculated correctly", {
+  # Create a mock count matrix
+  counts_matrix <- matrix(c(10, 20, 30, 40, 50, 60), nrow = 3, byrow = TRUE)
+  
+  # Expected normalized counts
+  expected_normalized <- matrix(c(11.54701, 34.64102, 57.73503, 17.32051, 34.64102, 51.96152), nrow = 3, byrow = F)
+  
+  # Perform normalization
+  normalized_counts <- trimmedMeanMvaluesNormalization(counts_matrix)
+  
+  # Check if normalized counts match the expected values
+  expect_equal(normalized_counts, expected_normalized, tolerance = 1e-2)
+})
+
+# Test case 2: Check if input matrix is unchanged when all counts are 0
+test_that("Input matrix remains unchanged when all counts are 0", {
+  # Create a mock count matrix with all counts as 0
+  counts_matrix <- matrix(0, nrow = 3, ncol = 2)
+  
+  # Perform normalization
+  normalized_counts <- trimmedMeanMvaluesNormalization(counts_matrix)
+  
+  # Check if input matrix and normalized matrix are identical
+  expect_identical(normalized_counts, counts_matrix)
+})
+
+
 ```
 
 ```{r function-fitmodel, filename = "fitmodel"}
@@ -3315,7 +3540,8 @@ identifyTopFit <- function(list_tmb, metric = "AIC", filter_method = "mad", keep
   if (filter_method == "mad"){
     lft_threshold <- get_mad_left_threshold(glance_df[[metric]], mad_tolerance)
     get_mad_user_message(metric, glance_df[[metric]], lft_threshold, keep)
-    index2keep <- if (keep == "top") glance_df[[metric]] > lft_threshold else glance_df[[metric]] < lft_threshold
+    index2keep <- if (keep == "top") which(glance_df[[metric]] > lft_threshold) 
+                  else which(glance_df[[metric]] < lft_threshold)
     id2keep <- rownames(glance_df)[index2keep]
     glance_df <- glance_df[ id2keep , ]
   } ## -- feel free to implement other filtering methods .. 
@@ -3344,8 +3570,8 @@ identifyTopFit <- function(list_tmb, metric = "AIC", filter_method = "mad", keep
 #' # Calculate the left threshold for MAD-based filtering with tolerance 2
 #' get_mad_left_threshold(c(50, 75, 100), 2)
 get_mad_left_threshold <- function(x, tolerance){
-  med <- stats::median(x)
-  mad <- stats::mad(x, center = median(x), 
+  med <- stats::median(x, na.rm = T)
+  mad <- stats::mad(x, center = med, 
                     constant = 1, na.rm = TRUE, 
                     low = FALSE, high = FALSE)
   med - (tolerance * mad)
@@ -3374,7 +3600,7 @@ get_mad_user_message <- function(var, x , left_threshold, keep) {
   keep_str <- ifelse(keep == "top", "above", "bellow")
   message(paste("2. Values",  keep_str , "the threshold were identified, threshold:",  left_threshold))
   message("3. Summary of selection:")
-  nb_keep <- ifelse(keep == "top", length(x[x > left_threshold]), length(x[x < left_threshold]) )
+  nb_keep <- ifelse(keep == "top", length(x[which(x > left_threshold)]), length(x[which(x < left_threshold)]) )
   sub_msg2 <- paste("values", keep_str,  "the threshold")
   message(paste("-", nb_keep, "out of", length(x), "observations had", sub_msg2, "for the", var, "metric."))
 }
@@ -3398,7 +3624,7 @@ test_that("identifyTopFit correctly identifies top-fitting objects", {
     ## -- prepare data & fit a model with mixed effect
     data2fit = prepareData2fit(countMatrix = mock_data$counts, 
                          metadata =  mock_data$metadata, 
-                         normalization = T)
+                         normalization = 'MRN')
     l_tmb <- fitModelParallel(formula = kij ~ myVariable, data = data2fit, 
                      group_by = "geneID", family = glmmTMB::nbinom2(link = "log"), 
                      n.cores = 1)
@@ -4650,7 +4876,7 @@ test_that("extract_tmbDispersion function extracts dispersion correctly", {
   input_var_list <- init_variable(name = "varA", mu = 10, sd = 0.1, level = 3)
   mock_data <- mock_rnaseq(input_var_list, N_GENES,
                          min_replicates = MIN_REPLICATES, max_replicates = MAX_REPLICATES)
-  data2fit <- prepareData2fit(countMatrix = mock_data$counts, metadata =  mock_data$metadata)
+  data2fit <- prepareData2fit(countMatrix = mock_data$counts, metadata =  mock_data$metadata, normalization = 'MRN')
   l_res <- fitModelParallel(formula = kij ~ varA,
                           data = data2fit, group_by = "geneID",
                           family = glmmTMB::nbinom2(link = "log"), n.cores = 1)
@@ -4683,7 +4909,7 @@ test_that("getDispersionComparison function works correctly", {
   input_var_list <- init_variable(name = "varA", mu = 10, sd = 0.1, level = 3)
   mock_data <- mock_rnaseq(input_var_list, N_GENES,
                          min_replicates = MIN_REPLICATES, max_replicates = MAX_REPLICATES)
-  data2fit <- prepareData2fit(countMatrix = mock_data$counts, metadata =  mock_data$metadata)
+  data2fit <- prepareData2fit(countMatrix = mock_data$counts, metadata =  mock_data$metadata, normalization = 'MRN')
   l_res <- fitModelParallel(formula = kij ~ varA,
                           data = data2fit, group_by = "geneID",
                           family = glmmTMB::nbinom2(link = "log"), n.cores = 1)
@@ -4708,30 +4934,43 @@ test_that("getDispersionComparison function works correctly", {
 #' This function scales a counts table based on the expected sequencing depth per sample.
 #'
 #' @param countsTable A counts table containing raw read counts.
-#' @param seq_depth  sequencing depth vector
+#' @param scalingDepth_factor  sequencing depth factor vector
 #' @return A scaled counts table.
 #'
 #' @export
 #' @examples
 #' mock_data <- list(counts = matrix(c(10, 20, 30, 20, 30, 10, 10, 20, 20, 20, 30, 10), ncol = 4))
-#' scaled_counts <- scaleCountsTable(countsTable = mock_data$counts, 1000000)
+#' scaled_counts <- scaleCountsTable(countsTable = mock_data$counts, 2)
 #'
-scaleCountsTable <- function(countsTable, seq_depth){
-  seq_depth_simu <- colSums(countsTable)
+scaleCountsTable <- function(countsTable, scalingDepth_factor){
+  counts_scaled <- as.data.frame(sweep(as.matrix(countsTable), 2,  scalingDepth_factor, "*"))
+  return(counts_scaled)
+}
+
 
+
+#' Get scaling factor for count normalization
+#'
+#' Calculates the scaling factor for count normalization based on the sequencing depth.
+#'
+#' @param countsTable Matrix or data frame of counts data.
+#' @param seq_depth Numeric vector containing the sequencing depths.
+#' @return Numeric vector of scaling factors for count normalization.
+#' @export
+get_scaling_factor <- function(countsTable, seq_depth){
+  seq_depth_simu <- colSums(countsTable)
+  
   if (length(seq_depth) > length(seq_depth_simu))
     message("INFO: The length of the sequencing_depth vector exceeds the number of samples. Only the first N values will be utilized.")
   if (length(seq_depth) < length(seq_depth_simu))
     message("INFO: The length of the sequencing_depth vector is shorter than the number of samples. Values will be recycled.")
 
   scalingDepth_factor <- suppressWarnings(seq_depth/seq_depth_simu)
-  counts_scaled <- as.data.frame(sweep(as.matrix(countsTable), 2,  scalingDepth_factor, "*"))
-  return(counts_scaled)
+  return(scalingDepth_factor)
 }
 
 
 
-
 ```
 
 ```{r  test-sequencing_depth_scaling}
@@ -4741,15 +4980,27 @@ test_that("Valid scaling of counts table", {
       # Test data
       mock_data <- list(counts = matrix(c(10, 20, 30, 20, 30, 10, 10, 20, 20, 20, 30, 10), ncol = 4))
       # Test function
-      scaled_counts <- scaleCountsTable(countsTable = mock_data$counts, 115000)
+      scaled_counts <- scaleCountsTable(countsTable = mock_data$counts, 2)
       
       # Expected scaled counts
-      expected_scaled_counts <- matrix(c(5000, 10000, 15000, 10000, 15000, 5000, 
-                                         5000, 10000, 10000, 10000, 15000, 5000), ncol = 4)
-      
+      expected_scaled_counts <- data.frame(matrix(c(20, 40, 60, 40, 60,20, 20,40,40, 40,60,20), ncol = 4))
+      colnames(expected_scaled_counts) <- c("V1", "V2", "V3", "V4")
       # Check if the scaled counts match the expected scaled counts
-      expect_true(all(colSums(scaled_counts) ==  115000))
+      expect_identical(scaled_counts, expected_scaled_counts)
+
+})
+
 
+
+
+
+test_that("get_scaling_factor returns correct scaling factors", {
+  # Mock data
+  counts <- matrix(c(10, 20, 30, 20, 30, 10, 10, 20, 20, 20, 30, 10), ncol = 4)  # Simulated counts data
+  seq_depth <- c(100, 200)  # Sequencing depths for each sample
+  scaling_factors <- suppressWarnings(get_scaling_factor(counts, seq_depth))
+  
+  expect_equal(scaling_factors, c(1.666667,3.333333,2.000000,3.333333), tolerance = 0.001)
 })
 
 ```
@@ -5154,7 +5405,7 @@ test_that("Test for subsetFixEffectInferred function", {
 
   mock_data <- mock_rnaseq(input_var_list, 10, 2, 2)
   getData2computeActualFixEffect(mock_data$groundTruth$effect)
-  data2fit <- prepareData2fit(countMatrix = mock_data$counts, metadata = mock_data$metadata, normalization = F)
+  data2fit <- prepareData2fit(countMatrix = mock_data$counts, metadata = mock_data$metadata, normalization = NULL)
 
   # Fit data
   resFit <- fitModelParallel(formula = kij ~ varA + varB + varA:varB,
@@ -5283,7 +5534,7 @@ test_that("subsetByTermLabel with non-existent term label", {
 test_that("getActualMainFixEff", {
   input_var_list <- init_variable() 
   mock_data <- mock_rnaseq(input_var_list, 2, 2, 2)
-  data2fit <- prepareData2fit(mock_data$counts, mock_data$metadata)
+  data2fit <- prepareData2fit(mock_data$counts, mock_data$metadata, normalization = 'MRN')
   inference <- fitModelParallel(kij ~ myVariable , 
                                   group_by = "geneID", data2fit, n.cores = 1)
   tidy_inference <- tidy_tmb(inference)
@@ -5309,7 +5560,7 @@ test_that("getData2computeActualFixEffect return correct output",{
   # Prepare the test data
   input_var_list <- init_variable() 
   mock_data <- mock_rnaseq(input_var_list, 2, 2, 2)
-  data2fit <- prepareData2fit(mock_data$counts, mock_data$metadata)
+  data2fit <- prepareData2fit(mock_data$counts, mock_data$metadata, normalization = 'MRN')
   inference <- fitModelParallel(kij ~ myVariable, group_by = "geneID", data2fit, n.cores = 1)
   tidy_inference <- tidy_tmb(inference)
   tidy_fix <- subsetFixEffectInferred(tidy_inference)
@@ -5330,7 +5581,7 @@ test_that("generateActualForMainFixEff returns correct values for main fixed eff
   # Prepare the test data
   input_var_list <- init_variable() 
   mock_data <- mock_rnaseq(input_var_list, 2, 2, 2)
-  data2fit <- prepareData2fit(mock_data$counts, mock_data$metadata)
+  data2fit <- prepareData2fit(mock_data$counts, mock_data$metadata, normalization = 'MRN')
   fixEff_dataActual <- getData2computeActualFixEffect(mock_data$groundTruth$effects)
   actual_intercept <- getActualIntercept(fixEff_dataActual)
   df_term <- generateActualForMainFixEff("myVariable2", actual_intercept, fixEff_dataActual$data, fixEff_dataActual$categorical_vars)
@@ -5783,7 +6034,7 @@ test_that("Generate actual interaction fixed effect correctly", {
   mock_data <- mock_rnaseq(init_var, N_GENES, min_replicates = MIN_REPLICATES, max_replicates = MAX_REPLICATES)
   
   # -- fit data
-  data2fit <- prepareData2fit(countMatrix = mock_data$counts, metadata = mock_data$metadata)
+  data2fit <- prepareData2fit(countMatrix = mock_data$counts, metadata = mock_data$metadata, normalization = 'MRN')
   
   dtf_countsLong <- countMatrix_2longDtf(mock_data$counts, "k_ij")
   metadata_columnForjoining <- getColumnWithSampleID(dtf_countsLong, mock_data$metadata)
@@ -5972,7 +6223,7 @@ test_that("Test getActualInteractionFixEff", {
   mock_data <- mock_rnaseq(init_var, N_GENES, min_replicates = MIN_REPLICATES, max_replicates = MAX_REPLICATES)
   
   # Données de fit
-  data2fit <- prepareData2fit(countMatrix = mock_data$counts, metadata = mock_data$metadata)
+  data2fit <- prepareData2fit(countMatrix = mock_data$counts, metadata = mock_data$metadata, normalization = 'MRN')
   results_fit <- fitModelParallel(formula = kij ~ varA + varB + varC + varA:varC,
                                   data = data2fit, group_by = "geneID",
                                   family = glmmTMB::nbinom2(link = "log"), n.cores = 1)
@@ -6027,7 +6278,7 @@ test_that("Test getActualInteractionFixEff", {
                            max_replicates = MAX_REPLICATES, dispersion = 100)
   
   # Données de fit
-  data2fit <- prepareData2fit(countMatrix = mock_data$counts, metadata = mock_data$metadata)
+  data2fit <- prepareData2fit(countMatrix = mock_data$counts, metadata = mock_data$metadata, normalization = 'MRN')
   results_fit <- fitModelParallel(formula = kij ~ varA + varB + varC + varA:varB + varB:varC + varA:varC + varA:varB:varC,
                                   data = data2fit, group_by = "geneID",
                                   family = glmmTMB::nbinom2(link = "log"), n.cores = 1)
@@ -7756,7 +8007,7 @@ test_that("evaluation_report returns correct output", {
   mock_data <- mock_rnaseq(input_var_list, N_GENES,
                            min_replicates = MIN_REPLICATES, max_replicates = MAX_REPLICATES)
   
-  data2fit <- prepareData2fit(countMatrix = mock_data$counts, metadata = mock_data$metadata)
+  data2fit <- prepareData2fit(countMatrix = mock_data$counts, metadata = mock_data$metadata, normalization = 'MRN')
   l_res <- fitModelParallel(formula = kij ~ varA,
                             data = data2fit, group_by = "geneID",
                             family = glmmTMB::nbinom2(link = "log"), n.cores = 1)
@@ -8048,7 +8299,7 @@ test_that("get_eval_data_from_ltmb returns correct output", {
   ## -- convert counts matrix and samples metadatas in a data frame for fitting
   data2fit = prepareData2fit(countMatrix = count_matrix,
                            metadata =  metaData,
-                           normalization = F)
+                           normalization = NULL)
   l_tmb <- fitModelParallel(formula = kij ~ varA  ,
                           data = data2fit,
                           group_by = "geneID",
@@ -8115,7 +8366,7 @@ test_that("rbind_evaldata_tmb_dds returns correct output", {
   ## -- convert counts matrix and samples metadatas in a data frame for fitting
   data2fit = prepareData2fit(countMatrix = count_matrix,
                            metadata =  metaData,
-                           normalization = F)
+                           normalization = NULL)
   l_tmb <- fitModelParallel(formula = kij ~ varA  ,
                           data = data2fit,
                           group_by = "geneID",
@@ -8188,7 +8439,7 @@ test_that("get_eval_data returns correct output", {
   ## -- convert counts matrix and samples metadatas in a data frame for fitting
   data2fit = prepareData2fit(countMatrix = count_matrix,
                            metadata =  metaData,
-                           normalization = F)
+                           normalization = NULL)
   l_tmb <- fitModelParallel(formula = kij ~ varA  ,
                           data = data2fit,
                           group_by = "geneID",
@@ -9021,7 +9272,7 @@ test_that("getCategoricalVar_inFixedEffect returns the correct result", {
                              max_replicates = MAX_REPLICATES,
                              basal_expression = 3, dispersion = 100)
     
-    data2fit = prepareData2fit(countMatrix = mock_data$counts, metadata =  mock_data$metadata, normalization = F)
+    data2fit = prepareData2fit(countMatrix = mock_data$counts, metadata =  mock_data$metadata, normalization = NULL)
     
     l_tmb <- fitModelParallel(formula = kij ~  environment  + (environment | genotype ),
                               data = data2fit, group_by = "geneID",
@@ -9048,7 +9299,7 @@ test_that("group_logQij_per_genes_and_labels returns the correct result", {
                              max_replicates = MAX_REPLICATES,
                              basal_expression = 3, dispersion = 100)
     
-    data2fit = prepareData2fit(countMatrix = mock_data$counts, metadata =  mock_data$metadata, normalization = F)
+    data2fit = prepareData2fit(countMatrix = mock_data$counts, metadata =  mock_data$metadata, normalization = NULL)
     
     l_tmb <- fitModelParallel(formula = kij ~  environment  + (environment | genotype ),
                               data = data2fit, group_by = "geneID",
@@ -9079,7 +9330,7 @@ test_that("getActualMixed_typeI returns the correct result", {
                              max_replicates = MAX_REPLICATES,
                              basal_expression = 3, dispersion = 100)
     
-    data2fit = prepareData2fit(countMatrix = mock_data$counts, metadata =  mock_data$metadata, normalization = F)
+    data2fit = prepareData2fit(countMatrix = mock_data$counts, metadata =  mock_data$metadata, normalization = NULL)
     
     l_tmb <- fitModelParallel(formula = kij ~  environment  + (environment | genotype ),
                               data = data2fit, group_by = "geneID",
@@ -9126,7 +9377,7 @@ test_that("inferenceToExpected_withMixedEff correctly compares inference to expe
                          max_replicates = MAX_REPLICATES,
                          basal_expression = 3, dispersion = 100)
   
-  data2fit <- prepareData2fit(countMatrix = mock_data$counts, metadata = mock_data$metadata, normalization = FALSE)
+  data2fit <- prepareData2fit(countMatrix = mock_data$counts, metadata = mock_data$metadata, normalization = NULL)
   
   l_tmb <- fitModelParallel(formula = kij ~ environment + (environment | genotype),
                           data = data2fit, group_by = "geneID",
@@ -9161,7 +9412,7 @@ test_that("calculate_actualMixed calculates actual mixed effects as expected", {
                          max_replicates = MAX_REPLICATES,
                          basal_expression = 3, dispersion = 100)
   
-  data2fit <- prepareData2fit(countMatrix = mock_data$counts, metadata = mock_data$metadata, normalization = FALSE)
+  data2fit <- prepareData2fit(countMatrix = mock_data$counts, metadata = mock_data$metadata, normalization = NULL)
   
   
   ground_truth_eff <- mock_data$groundTruth$effects
diff --git a/man/detect_row_matx_bellow_threshold.Rd b/man/detect_row_matx_bellow_threshold.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..533103015c6f3baaffbd961a3e5b16417785a1ef
--- /dev/null
+++ b/man/detect_row_matx_bellow_threshold.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utils.R
+\name{detect_row_matx_bellow_threshold}
+\alias{detect_row_matx_bellow_threshold}
+\title{Detect rows in a matrix with all values below a given threshold}
+\usage{
+detect_row_matx_bellow_threshold(matrix, threshold)
+}
+\arguments{
+\item{matrix}{The input matrix}
+
+\item{threshold}{The threshold value}
+}
+\value{
+A logical vector indicating rows below the threshold
+}
+\description{
+This function detects rows in a matrix where all values are below a specified threshold.
+}
diff --git a/man/getDataFromRnorm.Rd b/man/getDataFromRnorm.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..a9354eb34558fa52b66a835ba46c70aecd7257b5
--- /dev/null
+++ b/man/getDataFromRnorm.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/simulation.R
+\name{getDataFromRnorm}
+\alias{getDataFromRnorm}
+\title{Prepare data using effects from a normal distribution}
+\usage{
+getDataFromRnorm(list_var, n_genes)
+}
+\arguments{
+\item{list_var}{A list of variables (already initialized)}
+
+\item{n_genes}{Number of genes to generate data for.}
+}
+\value{
+A dataframe containing gene metadata and effects generated from a normal distribution.
+}
+\description{
+Prepares the data by generating effects from a normal distribution for each gene.
+}
diff --git a/man/getInput2simulation.Rd b/man/getInput2simulation.Rd
index 660db1d76bc901edbfae9d3503a4bf3a943050c3..f50ce204662c9713702ee8d742e7f47673f691c9 100644
--- a/man/getInput2simulation.Rd
+++ b/man/getInput2simulation.Rd
@@ -4,13 +4,24 @@
 \alias{getInput2simulation}
 \title{Get input for simulation based on coefficients}
 \usage{
-getInput2simulation(list_var, n_genes = 1, input2mvrnorm = NULL)
+getInput2simulation(
+  list_var,
+  n_genes = 1,
+  normal_distr = "multivariate",
+  input2mvrnorm = NULL
+)
 }
 \arguments{
 \item{list_var}{A list of variables (already initialized)}
 
 \item{n_genes}{Number of genes to simulate (default: 1)}
 
+\item{normal_distr}{Specifies the distribution type for generating effects. Choose between 'univariate' or 'multivariate' (default).
+\itemize{
+\item 'univariate': Effects are drawn independently from univariate normal distributions.
+\item 'multivariate': Effects are drawn jointly from a multivariate normal distribution.
+}}
+
 \item{input2mvrnorm}{Input to the \code{mvrnorm} function for simulating data from multivariate normal distribution (default: NULL)}
 }
 \value{
diff --git a/man/get_effects_from_rnorm.Rd b/man/get_effects_from_rnorm.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..afda8c72fb9d2de32b8520cd4f0a6b231b43fdd5
--- /dev/null
+++ b/man/get_effects_from_rnorm.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/simulation.R
+\name{get_effects_from_rnorm}
+\alias{get_effects_from_rnorm}
+\title{Generate effects from a normal distribution}
+\usage{
+get_effects_from_rnorm(list_var, metadata)
+}
+\arguments{
+\item{list_var}{A list of variables (already initialized)}
+
+\item{metadata}{Gene metadata.}
+}
+\value{
+A dataframe containing effects generated from a normal distribution.
+}
+\description{
+Generates effects from a normal distribution for each gene.
+}
diff --git a/man/get_messages_sequencing_depth.Rd b/man/get_messages_sequencing_depth.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..2edf59cfd5e330a59d94b2b2ae091e3456bb6c52
--- /dev/null
+++ b/man/get_messages_sequencing_depth.Rd
@@ -0,0 +1,16 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/mock_rnaseq.R
+\name{get_messages_sequencing_depth}
+\alias{get_messages_sequencing_depth}
+\title{Get messages related to sequencing depth}
+\usage{
+get_messages_sequencing_depth(scaling_factors, threshold_cov_var = 1.5)
+}
+\arguments{
+\item{scaling_factors}{Scaling factors obtained from sequencing depth}
+
+\item{threshold_cov_var}{Threshold coefficient of variation to detect heterogeneity in scaling}
+}
+\description{
+This function generates informative messages regarding the scaling of counts by sequencing depth.
+}
diff --git a/man/get_scaling_factor.Rd b/man/get_scaling_factor.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..d37edd79eeda5a18872ab04a3b5993580ac3af03
--- /dev/null
+++ b/man/get_scaling_factor.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/sequencing_depth_scaling.R
+\name{get_scaling_factor}
+\alias{get_scaling_factor}
+\title{Get scaling factor for count normalization}
+\usage{
+get_scaling_factor(countsTable, seq_depth)
+}
+\arguments{
+\item{countsTable}{Matrix or data frame of counts data.}
+
+\item{seq_depth}{Numeric vector containing the sequencing depths.}
+}
+\value{
+Numeric vector of scaling factors for count normalization.
+}
+\description{
+Calculates the scaling factor for count normalization based on the sequencing depth.
+}
diff --git a/man/mock_rnaseq.Rd b/man/mock_rnaseq.Rd
index e960338a6aed97cf713451def9271137f22b8f04..c68b463527ad693a563fadcf00fd6deb5b6a765d 100644
--- a/man/mock_rnaseq.Rd
+++ b/man/mock_rnaseq.Rd
@@ -11,7 +11,8 @@ mock_rnaseq(
   max_replicates,
   sequencing_depth = NULL,
   basal_expression = 0,
-  dispersion = stats::runif(n_genes, min = 0, max = 1000)
+  dispersion = stats::runif(n_genes, min = 0, max = 1000),
+  normal_distr = "multivariate"
 )
 }
 \arguments{
@@ -28,6 +29,12 @@ mock_rnaseq(
 \item{basal_expression}{base expression gene}
 
 \item{dispersion}{User-provided dispersion vector (optional)}
+
+\item{normal_distr}{Specifies the distribution type for generating effects. Choose between 'univariate' or 'multivariate' (default).
+\itemize{
+\item 'univariate': Effects are drawn independently from univariate normal distributions.
+\item 'multivariate': Effects are drawn jointly from a multivariate normal distribution.
+}}
 }
 \value{
 List containing the ground truth, counts, and metadata
diff --git a/man/prepareData2fit.Rd b/man/prepareData2fit.Rd
index 45fd195f0e2cc97ee232966cecb9988aab190902..6f4f9ef2462d628720a563b2755ca2679968e4e1 100644
--- a/man/prepareData2fit.Rd
+++ b/man/prepareData2fit.Rd
@@ -9,8 +9,9 @@ prepareData2fit(
   metadata,
   response_name = "kij",
   groupID = "geneID",
+  row_threshold = 0,
   transform = NULL,
-  normalization = TRUE
+  normalization = NULL
 )
 }
 \arguments{
@@ -22,15 +23,18 @@ prepareData2fit(
 
 \item{groupID}{String referring to the group variable name (default: "geneID").}
 
+\item{row_threshold}{Numeric threshold for removing rows with all counts below a specified value. Default 0.
+This filtering is applied before transformation and normalization.}
+
 \item{transform}{A custom R expression to apply to each element of the countMatrix.
 This expression should be provided as a character string.
 For example, to apply log transformation, use \code{"log(x)"}.
 Note that \code{x} represents each element in the countMatrix. See examples for more details.
 The transformation is applied before normalization (if normalization = \code{TRUE}).}
 
-\item{normalization}{A boolean value indicating whether to apply median ratio normalization.
-If \code{TRUE} (default), the counts matrix will be normalized using median ratio normalization.
-If \code{FALSE}, no normalization will be applied.}
+\item{normalization}{a vector character specifying method to use (default: NULL, possible choices: c('MRN', 'TTM'))
+- MRN: median ratio normalization
+- TMM: Trimmed Mean of M-values}
 }
 \value{
 Data frame suitable for fitting.
diff --git a/man/scaleCountsTable.Rd b/man/scaleCountsTable.Rd
index 1ad7280869fd056437cdfa9c1b4bff4332a30964..6bbd05f1ed0d42a5aaa94d43c54bcedf78a06c36 100644
--- a/man/scaleCountsTable.Rd
+++ b/man/scaleCountsTable.Rd
@@ -4,12 +4,12 @@
 \alias{scaleCountsTable}
 \title{Scale Counts Table}
 \usage{
-scaleCountsTable(countsTable, seq_depth)
+scaleCountsTable(countsTable, scalingDepth_factor)
 }
 \arguments{
 \item{countsTable}{A counts table containing raw read counts.}
 
-\item{seq_depth}{sequencing depth vector}
+\item{scalingDepth_factor}{sequencing depth factor vector}
 }
 \value{
 A scaled counts table.
@@ -19,6 +19,6 @@ This function scales a counts table based on the expected sequencing depth per s
 }
 \examples{
 mock_data <- list(counts = matrix(c(10, 20, 30, 20, 30, 10, 10, 20, 20, 20, 30, 10), ncol = 4))
-scaled_counts <- scaleCountsTable(countsTable = mock_data$counts, 1000000)
+scaled_counts <- scaleCountsTable(countsTable = mock_data$counts, 2)
 
 }
diff --git a/man/trimmedMeanMvaluesNormalization.Rd b/man/trimmedMeanMvaluesNormalization.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..fbf3f012f924715933cfbc9df0774a3171798e81
--- /dev/null
+++ b/man/trimmedMeanMvaluesNormalization.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/prepare_data2fit.R
+\name{trimmedMeanMvaluesNormalization}
+\alias{trimmedMeanMvaluesNormalization}
+\title{Normalize count data using Trimmed Mean of M-values (TMM) method}
+\usage{
+trimmedMeanMvaluesNormalization(counts_matrix)
+}
+\arguments{
+\item{counts_matrix}{A matrix of count data where rows represent genes and columns represent samples.}
+}
+\value{
+Normalized count matrix
+}
+\description{
+This function normalizes count data using the Trimmed Mean of M-values (TMM) method,
+which calculates scale factors to account for differences in library sizes
+and RNA composition between samples.
+}
diff --git a/man/warning_too_low_mu_ij_row.Rd b/man/warning_too_low_mu_ij_row.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..d296c90cc721d90980bd28294267f285547c6310
--- /dev/null
+++ b/man/warning_too_low_mu_ij_row.Rd
@@ -0,0 +1,16 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/mock_rnaseq.R
+\name{warning_too_low_mu_ij_row}
+\alias{warning_too_low_mu_ij_row}
+\title{Emit a warning message for rows with low mu_ij values}
+\usage{
+warning_too_low_mu_ij_row(mu_ij_matrix, threshold = 1)
+}
+\arguments{
+\item{mu_ij_matrix}{Matrix of mu_ij values}
+
+\item{threshold}{Threshold value}
+}
+\description{
+This function emits a warning message if any rows in the mu_ij matrix have all values below a specified threshold.
+}
diff --git a/tests/testthat/test-actual_interactionfixeffects.R b/tests/testthat/test-actual_interactionfixeffects.R
index fdd12c253c65fcc8962fd7fa2d94f64c72bffb34..b628fbd66701f64170078e58d7a312b82c20598d 100644
--- a/tests/testthat/test-actual_interactionfixeffects.R
+++ b/tests/testthat/test-actual_interactionfixeffects.R
@@ -100,7 +100,7 @@ test_that("Generate actual interaction fixed effect correctly", {
   mock_data <- mock_rnaseq(init_var, N_GENES, min_replicates = MIN_REPLICATES, max_replicates = MAX_REPLICATES)
   
   # -- fit data
-  data2fit <- prepareData2fit(countMatrix = mock_data$counts, metadata = mock_data$metadata)
+  data2fit <- prepareData2fit(countMatrix = mock_data$counts, metadata = mock_data$metadata, normalization = 'MRN')
   
   dtf_countsLong <- countMatrix_2longDtf(mock_data$counts, "k_ij")
   metadata_columnForjoining <- getColumnWithSampleID(dtf_countsLong, mock_data$metadata)
@@ -289,7 +289,7 @@ test_that("Test getActualInteractionFixEff", {
   mock_data <- mock_rnaseq(init_var, N_GENES, min_replicates = MIN_REPLICATES, max_replicates = MAX_REPLICATES)
   
   # Données de fit
-  data2fit <- prepareData2fit(countMatrix = mock_data$counts, metadata = mock_data$metadata)
+  data2fit <- prepareData2fit(countMatrix = mock_data$counts, metadata = mock_data$metadata, normalization = 'MRN')
   results_fit <- fitModelParallel(formula = kij ~ varA + varB + varC + varA:varC,
                                   data = data2fit, group_by = "geneID",
                                   family = glmmTMB::nbinom2(link = "log"), n.cores = 1)
@@ -344,7 +344,7 @@ test_that("Test getActualInteractionFixEff", {
                            max_replicates = MAX_REPLICATES, dispersion = 100)
   
   # Données de fit
-  data2fit <- prepareData2fit(countMatrix = mock_data$counts, metadata = mock_data$metadata)
+  data2fit <- prepareData2fit(countMatrix = mock_data$counts, metadata = mock_data$metadata, normalization = 'MRN')
   results_fit <- fitModelParallel(formula = kij ~ varA + varB + varC + varA:varB + varB:varC + varA:varC + varA:varB:varC,
                                   data = data2fit, group_by = "geneID",
                                   family = glmmTMB::nbinom2(link = "log"), n.cores = 1)
diff --git a/tests/testthat/test-actual_mainfixeffects.R b/tests/testthat/test-actual_mainfixeffects.R
index 308ad816dc00cba812d1b7e4eb78b37a250a3654..f13c8b15e0fe9771abe96b65759254e0b68d4543 100644
--- a/tests/testthat/test-actual_mainfixeffects.R
+++ b/tests/testthat/test-actual_mainfixeffects.R
@@ -9,7 +9,7 @@ test_that("Test for subsetFixEffectInferred function", {
 
   mock_data <- mock_rnaseq(input_var_list, 10, 2, 2)
   getData2computeActualFixEffect(mock_data$groundTruth$effect)
-  data2fit <- prepareData2fit(countMatrix = mock_data$counts, metadata = mock_data$metadata, normalization = F)
+  data2fit <- prepareData2fit(countMatrix = mock_data$counts, metadata = mock_data$metadata, normalization = NULL)
 
   # Fit data
   resFit <- fitModelParallel(formula = kij ~ varA + varB + varA:varB,
@@ -138,7 +138,7 @@ test_that("subsetByTermLabel with non-existent term label", {
 test_that("getActualMainFixEff", {
   input_var_list <- init_variable() 
   mock_data <- mock_rnaseq(input_var_list, 2, 2, 2)
-  data2fit <- prepareData2fit(mock_data$counts, mock_data$metadata)
+  data2fit <- prepareData2fit(mock_data$counts, mock_data$metadata, normalization = 'MRN')
   inference <- fitModelParallel(kij ~ myVariable , 
                                   group_by = "geneID", data2fit, n.cores = 1)
   tidy_inference <- tidy_tmb(inference)
@@ -164,7 +164,7 @@ test_that("getData2computeActualFixEffect return correct output",{
   # Prepare the test data
   input_var_list <- init_variable() 
   mock_data <- mock_rnaseq(input_var_list, 2, 2, 2)
-  data2fit <- prepareData2fit(mock_data$counts, mock_data$metadata)
+  data2fit <- prepareData2fit(mock_data$counts, mock_data$metadata, normalization = 'MRN')
   inference <- fitModelParallel(kij ~ myVariable, group_by = "geneID", data2fit, n.cores = 1)
   tidy_inference <- tidy_tmb(inference)
   tidy_fix <- subsetFixEffectInferred(tidy_inference)
@@ -185,7 +185,7 @@ test_that("generateActualForMainFixEff returns correct values for main fixed eff
   # Prepare the test data
   input_var_list <- init_variable() 
   mock_data <- mock_rnaseq(input_var_list, 2, 2, 2)
-  data2fit <- prepareData2fit(mock_data$counts, mock_data$metadata)
+  data2fit <- prepareData2fit(mock_data$counts, mock_data$metadata, normalization = 'MRN')
   fixEff_dataActual <- getData2computeActualFixEffect(mock_data$groundTruth$effects)
   actual_intercept <- getActualIntercept(fixEff_dataActual)
   df_term <- generateActualForMainFixEff("myVariable2", actual_intercept, fixEff_dataActual$data, fixEff_dataActual$categorical_vars)
diff --git a/tests/testthat/test-evaluate_dispersion.R b/tests/testthat/test-evaluate_dispersion.R
index cb7215836d1afb669e4a928ace67caae07f1776d..5537b9cdbf995605850b4bc47fd46f6924a9545b 100644
--- a/tests/testthat/test-evaluate_dispersion.R
+++ b/tests/testthat/test-evaluate_dispersion.R
@@ -10,7 +10,7 @@ test_that("extract_tmbDispersion function extracts dispersion correctly", {
   input_var_list <- init_variable(name = "varA", mu = 10, sd = 0.1, level = 3)
   mock_data <- mock_rnaseq(input_var_list, N_GENES,
                          min_replicates = MIN_REPLICATES, max_replicates = MAX_REPLICATES)
-  data2fit <- prepareData2fit(countMatrix = mock_data$counts, metadata =  mock_data$metadata)
+  data2fit <- prepareData2fit(countMatrix = mock_data$counts, metadata =  mock_data$metadata, normalization = 'MRN')
   l_res <- fitModelParallel(formula = kij ~ varA,
                           data = data2fit, group_by = "geneID",
                           family = glmmTMB::nbinom2(link = "log"), n.cores = 1)
@@ -43,7 +43,7 @@ test_that("getDispersionComparison function works correctly", {
   input_var_list <- init_variable(name = "varA", mu = 10, sd = 0.1, level = 3)
   mock_data <- mock_rnaseq(input_var_list, N_GENES,
                          min_replicates = MIN_REPLICATES, max_replicates = MAX_REPLICATES)
-  data2fit <- prepareData2fit(countMatrix = mock_data$counts, metadata =  mock_data$metadata)
+  data2fit <- prepareData2fit(countMatrix = mock_data$counts, metadata =  mock_data$metadata, normalization = 'MRN')
   l_res <- fitModelParallel(formula = kij ~ varA,
                           data = data2fit, group_by = "geneID",
                           family = glmmTMB::nbinom2(link = "log"), n.cores = 1)
diff --git a/tests/testthat/test-evaluation_withmixedeffect.R b/tests/testthat/test-evaluation_withmixedeffect.R
index 0d7af6f761ac101850dd14cf8e063718e2dde0a7..a61e1bd19d41683ea9065aebc305cad90dd21fa9 100644
--- a/tests/testthat/test-evaluation_withmixedeffect.R
+++ b/tests/testthat/test-evaluation_withmixedeffect.R
@@ -46,7 +46,7 @@ test_that("getCategoricalVar_inFixedEffect returns the correct result", {
                              max_replicates = MAX_REPLICATES,
                              basal_expression = 3, dispersion = 100)
     
-    data2fit = prepareData2fit(countMatrix = mock_data$counts, metadata =  mock_data$metadata, normalization = F)
+    data2fit = prepareData2fit(countMatrix = mock_data$counts, metadata =  mock_data$metadata, normalization = NULL)
     
     l_tmb <- fitModelParallel(formula = kij ~  environment  + (environment | genotype ),
                               data = data2fit, group_by = "geneID",
@@ -73,7 +73,7 @@ test_that("group_logQij_per_genes_and_labels returns the correct result", {
                              max_replicates = MAX_REPLICATES,
                              basal_expression = 3, dispersion = 100)
     
-    data2fit = prepareData2fit(countMatrix = mock_data$counts, metadata =  mock_data$metadata, normalization = F)
+    data2fit = prepareData2fit(countMatrix = mock_data$counts, metadata =  mock_data$metadata, normalization = NULL)
     
     l_tmb <- fitModelParallel(formula = kij ~  environment  + (environment | genotype ),
                               data = data2fit, group_by = "geneID",
@@ -104,7 +104,7 @@ test_that("getActualMixed_typeI returns the correct result", {
                              max_replicates = MAX_REPLICATES,
                              basal_expression = 3, dispersion = 100)
     
-    data2fit = prepareData2fit(countMatrix = mock_data$counts, metadata =  mock_data$metadata, normalization = F)
+    data2fit = prepareData2fit(countMatrix = mock_data$counts, metadata =  mock_data$metadata, normalization = NULL)
     
     l_tmb <- fitModelParallel(formula = kij ~  environment  + (environment | genotype ),
                               data = data2fit, group_by = "geneID",
@@ -151,7 +151,7 @@ test_that("inferenceToExpected_withMixedEff correctly compares inference to expe
                          max_replicates = MAX_REPLICATES,
                          basal_expression = 3, dispersion = 100)
   
-  data2fit <- prepareData2fit(countMatrix = mock_data$counts, metadata = mock_data$metadata, normalization = FALSE)
+  data2fit <- prepareData2fit(countMatrix = mock_data$counts, metadata = mock_data$metadata, normalization = NULL)
   
   l_tmb <- fitModelParallel(formula = kij ~ environment + (environment | genotype),
                           data = data2fit, group_by = "geneID",
@@ -186,7 +186,7 @@ test_that("calculate_actualMixed calculates actual mixed effects as expected", {
                          max_replicates = MAX_REPLICATES,
                          basal_expression = 3, dispersion = 100)
   
-  data2fit <- prepareData2fit(countMatrix = mock_data$counts, metadata = mock_data$metadata, normalization = FALSE)
+  data2fit <- prepareData2fit(countMatrix = mock_data$counts, metadata = mock_data$metadata, normalization = NULL)
   
   
   ground_truth_eff <- mock_data$groundTruth$effects
diff --git a/tests/testthat/test-filtering_fit.R b/tests/testthat/test-filtering_fit.R
index 021f0d85d500137406fa4b2073bc246523e1e5d5..601de79b295439103e27e8c5aebe4c5bc6424fc1 100644
--- a/tests/testthat/test-filtering_fit.R
+++ b/tests/testthat/test-filtering_fit.R
@@ -15,7 +15,7 @@ test_that("identifyTopFit correctly identifies top-fitting objects", {
     ## -- prepare data & fit a model with mixed effect
     data2fit = prepareData2fit(countMatrix = mock_data$counts, 
                          metadata =  mock_data$metadata, 
-                         normalization = T)
+                         normalization = 'MRN')
     l_tmb <- fitModelParallel(formula = kij ~ myVariable, data = data2fit, 
                      group_by = "geneID", family = glmmTMB::nbinom2(link = "log"), 
                      n.cores = 1)
diff --git a/tests/testthat/test-mock_rnaseq.R b/tests/testthat/test-mock_rnaseq.R
index 6596bc42db8ffa8c6faff3855e488ef0751b6aa8..35449e93f70c4ae2e07d524200c6f2cf62201817 100644
--- a/tests/testthat/test-mock_rnaseq.R
+++ b/tests/testthat/test-mock_rnaseq.R
@@ -44,34 +44,13 @@ test_that("replicateMatrix replicates matrix correctly", {
   expect_equal(dim(replicated_matrix), c(3, 9))
 })
 
-
-# Test for mock_rnaseq
-#test_that("mock_rnaseq returns expected output", {
-  # Set up input variables
-#  list_var <- NULL
-#  n_genes <- 3
-#  min_replicates <- 2
-#  max_replicates <- 4
-#  df_inputSimulation <- data.frame(gene_id = 1:3, coef_value = c(0.5, 0.3, 0.2))
-#  matx_dispersion <- matrix(1:9, nrow = 3, ncol = 3)
-
-  # Run the function
-#  expect_error(mock_rnaseq(list_var, n_genes, min_replicates, max_replicates, df_inputSimulation, 
-#                           matx_dispersion))
-  
-  
-  #list_var <- init_variable(name = "my_var", mu = c(10, 20), level = 2 )
-  #n_genes <- 10
-  #min_replicates <- 2
-  #max_replicates <- 4
-  #scaling_factor <- 1
-  #df_inputSimulation <- getInput2simulation(list_var, n_genes)
-  #dispersion <- getDispersionMatrix(list_var, n_genes, c(1000, 1000, 1000, 1000, 1000, 1, 1, 1, 1, 1))
-  #mock_rnaseq(list_var, n_genes, min_replicates, 
-  #            max_replicates, 
-  #            df_inputSimulation, dispersion)
-  #ERROOR
-#})
+# Test for warning_too_low_mu_ij_row function
+test_that("warning_too_low_mu_ij_row emits warning for low mu_ij values", {
+  # Create a sample mu_ij matrix
+  mu_ij_matrix <- matrix(c(0.5, 0.7, 1.2, 0.2, 0.9, 1.1), nrow = 2)
+  # Verify if warning message is emitted
+  expect_message(warning_too_low_mu_ij_row(mu_ij_matrix, 2))
+})
 
 
 # Test for generateReplicationMatrix
diff --git a/tests/testthat/test-prepare_data2fit.R b/tests/testthat/test-prepare_data2fit.R
index 451e19d282cd7a36f93a1050476dcb820bde6c0e..941f900ee7335cb5ebfacf762c186526e5468fa8 100644
--- a/tests/testthat/test-prepare_data2fit.R
+++ b/tests/testthat/test-prepare_data2fit.R
@@ -41,7 +41,7 @@ test_that("prepareData2fit prepares data for fitting", {
   mock_data <- mock_rnaseq(list_var, n_genes = 3, 2,2)
   
   # Prepare data for fitting
-  data2fit <- prepareData2fit(mock_data$counts, mock_data$metadata)
+  data2fit <- prepareData2fit(mock_data$counts, mock_data$metadata, normalization = 'MRN')
   
   expect_true(is.character(data2fit$sampleID))
   expect_true(is.character(data2fit$geneID))
@@ -56,13 +56,13 @@ test_that("prepareData2fit prepares data for fitting", {
   metadata <- data.frame(sampleID = colnames(countMatrix), condition = rep(c("A", "B"), each = 5))
   
   # Call the function with log transformation
-  data2fit <- prepareData2fit(countMatrix, metadata, normalization = F ,transform = "log(x)")
+  data2fit <- prepareData2fit(countMatrix, metadata, normalization = NULL,transform = "log(x)")
   
   # Test if log transformation has been applied correctly
   expect_equal(unique(data2fit$kij), log(c(1, 2, 3, 4)))
   
   # Call the function with sqrt transformation
-  data2fit <- prepareData2fit(countMatrix, metadata, normalization = F , transform = "sqrt(x)")
+  data2fit <- prepareData2fit(countMatrix, metadata, normalization = NULL , transform = "sqrt(x)")
   # Test if log transformation has been applied correctly
   expect_equal(unique(data2fit$kij), sqrt(c(1, 2, 3, 4)))
 })
@@ -121,3 +121,32 @@ test_that("Throws an error when all-zero genes are encountered", {
 })
 
 
+
+# Test case 1: Check if normalization is performed correctly
+test_that("Normalized counts are calculated correctly", {
+  # Create a mock count matrix
+  counts_matrix <- matrix(c(10, 20, 30, 40, 50, 60), nrow = 3, byrow = TRUE)
+  
+  # Expected normalized counts
+  expected_normalized <- matrix(c(11.54701, 34.64102, 57.73503, 17.32051, 34.64102, 51.96152), nrow = 3, byrow = F)
+  
+  # Perform normalization
+  normalized_counts <- trimmedMeanMvaluesNormalization(counts_matrix)
+  
+  # Check if normalized counts match the expected values
+  expect_equal(normalized_counts, expected_normalized, tolerance = 1e-2)
+})
+
+# Test case 2: Check if input matrix is unchanged when all counts are 0
+test_that("Input matrix remains unchanged when all counts are 0", {
+  # Create a mock count matrix with all counts as 0
+  counts_matrix <- matrix(0, nrow = 3, ncol = 2)
+  
+  # Perform normalization
+  normalized_counts <- trimmedMeanMvaluesNormalization(counts_matrix)
+  
+  # Check if input matrix and normalized matrix are identical
+  expect_identical(normalized_counts, counts_matrix)
+})
+
+
diff --git a/tests/testthat/test-sequencing_depth_scaling.R b/tests/testthat/test-sequencing_depth_scaling.R
index 294f05fd537ae1627baf386c5e5c2433e9fb39c2..e33616dac5df00ef865db928e55b2f4676b35eb1 100644
--- a/tests/testthat/test-sequencing_depth_scaling.R
+++ b/tests/testthat/test-sequencing_depth_scaling.R
@@ -6,14 +6,26 @@ test_that("Valid scaling of counts table", {
       # Test data
       mock_data <- list(counts = matrix(c(10, 20, 30, 20, 30, 10, 10, 20, 20, 20, 30, 10), ncol = 4))
       # Test function
-      scaled_counts <- scaleCountsTable(countsTable = mock_data$counts, 115000)
+      scaled_counts <- scaleCountsTable(countsTable = mock_data$counts, 2)
       
       # Expected scaled counts
-      expected_scaled_counts <- matrix(c(5000, 10000, 15000, 10000, 15000, 5000, 
-                                         5000, 10000, 10000, 10000, 15000, 5000), ncol = 4)
-      
+      expected_scaled_counts <- data.frame(matrix(c(20, 40, 60, 40, 60,20, 20,40,40, 40,60,20), ncol = 4))
+      colnames(expected_scaled_counts) <- c("V1", "V2", "V3", "V4")
       # Check if the scaled counts match the expected scaled counts
-      expect_true(all(colSums(scaled_counts) ==  115000))
+      expect_identical(scaled_counts, expected_scaled_counts)
+
+})
+
+
+
+
 
+test_that("get_scaling_factor returns correct scaling factors", {
+  # Mock data
+  counts <- matrix(c(10, 20, 30, 20, 30, 10, 10, 20, 20, 20, 30, 10), ncol = 4)  # Simulated counts data
+  seq_depth <- c(100, 200)  # Sequencing depths for each sample
+  scaling_factors <- suppressWarnings(get_scaling_factor(counts, seq_depth))
+  
+  expect_equal(scaling_factors, c(1.666667,3.333333,2.000000,3.333333), tolerance = 0.001)
 })
 
diff --git a/tests/testthat/test-simulation.R b/tests/testthat/test-simulation.R
index 0868462c39ce176d6f00730ea482f60971fcfe4d..d43b959e0e60bdbe2c465297fe7b7c5550db2f68 100644
--- a/tests/testthat/test-simulation.R
+++ b/tests/testthat/test-simulation.R
@@ -1,5 +1,28 @@
 # WARNING - Generated by {fusen} from dev/flat_full.Rmd: do not edit by hand
 
+
+test_that("getDataFromRnorm generates correct data frame", {
+  input_var_list <- init_variable(name = "varA", mu = 10, sd = 0.1, level = 3) %>%
+                    init_variable(name = "varB", mu = 1, sd = 2, level = 2)
+  metadata <- getGeneMetadata(input_var_list , n_genes = 5)
+  df <- getDataFromRnorm(input_var_list, n_genes = 5)
+  expect_is(df[[1]], "data.frame")
+  expect_equal(nrow(df[[1]]), 30)
+  expect_equal(colnames(df[[1]]), c("geneID", "label_varA", "label_varB", "varA", "varB"))  
+})
+
+test_that("get_effects_from_rnorm generates correct effects", {
+  input_var_list <- init_variable(name = "varA", mu = 10, sd = 0.1, level = 3) %>%
+                    init_variable(name = "varB", mu = 1, sd = 2, level = 2)
+  metadata <- getGeneMetadata(input_var_list , n_genes = 5)
+  df_effects <- get_effects_from_rnorm(input_var_list, metadata)
+  
+  expect_is(df_effects, "matrix")
+  expect_equal(nrow(df_effects), nrow(metadata))
+  expect_equal(colnames(df_effects), c("varA", "varB"))
+})
+
+
 # Test case 1: Check if the function returns a data frame
 test_that("getInput2simulation returns a data frame", {
   list_var <- init_variable()
diff --git a/tests/testthat/test-simulation_report.R b/tests/testthat/test-simulation_report.R
index 39e538c18ad6ea83eb2ebe57e397776ca359b142..36192272f5d8556855aa02b9183f682d7bef9623 100644
--- a/tests/testthat/test-simulation_report.R
+++ b/tests/testthat/test-simulation_report.R
@@ -72,7 +72,7 @@ test_that("evaluation_report returns correct output", {
   mock_data <- mock_rnaseq(input_var_list, N_GENES,
                            min_replicates = MIN_REPLICATES, max_replicates = MAX_REPLICATES)
   
-  data2fit <- prepareData2fit(countMatrix = mock_data$counts, metadata = mock_data$metadata)
+  data2fit <- prepareData2fit(countMatrix = mock_data$counts, metadata = mock_data$metadata, normalization = 'MRN')
   l_res <- fitModelParallel(formula = kij ~ varA,
                             data = data2fit, group_by = "geneID",
                             family = glmmTMB::nbinom2(link = "log"), n.cores = 1)
@@ -219,7 +219,7 @@ test_that("get_eval_data_from_ltmb returns correct output", {
   ## -- convert counts matrix and samples metadatas in a data frame for fitting
   data2fit = prepareData2fit(countMatrix = count_matrix,
                            metadata =  metaData,
-                           normalization = F)
+                           normalization = NULL)
   l_tmb <- fitModelParallel(formula = kij ~ varA  ,
                           data = data2fit,
                           group_by = "geneID",
@@ -286,7 +286,7 @@ test_that("rbind_evaldata_tmb_dds returns correct output", {
   ## -- convert counts matrix and samples metadatas in a data frame for fitting
   data2fit = prepareData2fit(countMatrix = count_matrix,
                            metadata =  metaData,
-                           normalization = F)
+                           normalization = NULL)
   l_tmb <- fitModelParallel(formula = kij ~ varA  ,
                           data = data2fit,
                           group_by = "geneID",
@@ -359,7 +359,7 @@ test_that("get_eval_data returns correct output", {
   ## -- convert counts matrix and samples metadatas in a data frame for fitting
   data2fit = prepareData2fit(countMatrix = count_matrix,
                            metadata =  metaData,
-                           normalization = F)
+                           normalization = NULL)
   l_tmb <- fitModelParallel(formula = kij ~ varA  ,
                           data = data2fit,
                           group_by = "geneID",
diff --git a/tests/testthat/test-utils.R b/tests/testthat/test-utils.R
index 58c27db7c0137da0da518a78db10f020eb1c8b36..4bff341b792282cc1131588f5b1bb232653db6f6 100644
--- a/tests/testthat/test-utils.R
+++ b/tests/testthat/test-utils.R
@@ -85,6 +85,19 @@ test_that("removeDuplicatedWord returns expected output", {
   expect_equal(cleaned_words, c("hello", "world", "programming", "R is great"))
 })
 
+
+# Test for detect_row_matx_bellow_threshold function
+test_that("detect_row_matx_bellow_threshold detects rows below threshold", {
+  # Create a sample matrix
+  matrix <- matrix(c(0.5, 0.7, 1.2, 0.2, 0.9, 0.9), nrow = 2)
+  # Test with threshold 1
+  expect_equal(detect_row_matx_bellow_threshold(matrix, 1), c(FALSE, TRUE))
+  # Test with threshold 0.5
+  expect_equal(detect_row_matx_bellow_threshold(matrix, 0.5), c(FALSE, FALSE))
+  expect_equal(detect_row_matx_bellow_threshold(matrix, 2), c(TRUE, TRUE))
+})
+
+
 test_that("reorderColumns returns expected output",{
     df <- data.frame(A = 1:3, B = 4:6, C = 7:9)
     # Define the desired column order
diff --git a/vignettes/02-tutorial.Rmd b/vignettes/02-tutorial.Rmd
index 16c05f1beccc7203c0b5c6c9cfea62df33600def..04ee73cfb515691d8a21f08ccf18f4736c121207 100644
--- a/vignettes/02-tutorial.Rmd
+++ b/vignettes/02-tutorial.Rmd
@@ -48,7 +48,7 @@ To estimate plausible parameter distributions, we will fit a generalized linear
 ## -- prepare data &  fit a model 
 data2fit = prepareData2fit(countMatrix = pub_data, 
                            metadata = pub_metadata, 
-                           normalization = F,   
+                           normalization = 'MRN',   
                            response_name = "kij")
 l_tmb <- fitModelParallel(
             formula = kij ~ genotype + environment + genotype:environment,
@@ -149,7 +149,7 @@ We fit a model to the simulated RNAseq data using the `fitModelParallel()` funct
 ## -- prepare data & fit a model with mixed effect
 data2fit = prepareData2fit(countMatrix = mock_data$counts, 
                            metadata =  mock_data$metadata, 
-                           normalization = T,   
+                           normalization = 'MRN',   
                            response_name = "kij")
 l_tmb <- fitModelParallel(
           formula = kij ~ genotype + environment + genotype:environment,
@@ -275,7 +275,7 @@ We fit a model to the simulated RNAseq data using the `fitModelParallel()` funct
 ## -- prepare data & fit a model with mixed effect
 data2fit = prepareData2fit(countMatrix = mock_data$counts, 
                            metadata =  mock_data$metadata, 
-                           normalization = T,   
+                           normalization = 'MRN',   
                            response_name = "kij")
 l_tmb <- fitModelParallel(
             formula = kij ~ genotype + environment + genotype:environment,
@@ -326,7 +326,7 @@ The inclusion of a large number of genotypes in the experimental design allows f
 ## -- prepare data & fit a model with mixed effect
 data2fit = prepareData2fit(countMatrix = mock_data$counts, 
                            metadata =  mock_data$metadata, 
-                           normalization = T,   
+                           normalization = 'MRN',   
                            response_name = "kij")
 l_tmb <- fitModelParallel(
             formula = kij  ~environment + ( environment | genotype ),
diff --git a/vignettes/03-rnaseq_analysis.Rmd b/vignettes/03-rnaseq_analysis.Rmd
index f9c542e2eb513dabb2a7ff1bb1cf01477cda98d2..b19ba160e644b7e7615f2c25bb214aa650f9ec0f 100644
--- a/vignettes/03-rnaseq_analysis.Rmd
+++ b/vignettes/03-rnaseq_analysis.Rmd
@@ -68,7 +68,7 @@ The `prepareData2fit()` function serves the purpose of converting the counts mat
 data2fit = prepareData2fit(
              countMatrix = count_matrix, 
              metadata =  metaData, 
-             normalization = F,
+             normalization = NULL,
              response_name = "kij")
 
 
@@ -76,14 +76,14 @@ data2fit = prepareData2fit(
 data2fit = prepareData2fit(
              countMatrix = count_matrix, 
              metadata =  metaData, 
-             normalization = T, 
+             normalization = 'MRN', 
              response_name = "kij")
 
 ## -- apply + 1 transformation on each counts
 data2fit = prepareData2fit(
              countMatrix = count_matrix, 
              metadata =  metaData, 
-             normalization = T, 
+             normalization = 'MRN', 
              transform = "x+1",
              response_name = "kij")
 ```
diff --git a/vignettes/04-htrfit_vs_deseq2.Rmd b/vignettes/04-htrfit_vs_deseq2.Rmd
index f76456f96cd12bf080da96b73a35e51057310a8b..14d820a8904f6ed72a84ef9d0a7344e4ae22b4f4 100644
--- a/vignettes/04-htrfit_vs_deseq2.Rmd
+++ b/vignettes/04-htrfit_vs_deseq2.Rmd
@@ -65,7 +65,7 @@ metaData <- mock_data$metadata
 ## -- convert counts matrix and samples metadatas in a data frame for fitting
 data2fit = prepareData2fit(countMatrix = count_matrix, 
                            metadata =  metaData, 
-                           normalization = T,
+                           normalization = 'MRN',
                            response_name = "kij")
 
 l_tmb <- fitModelParallel(