library(dplyr)
library(ggplot2)

# -----------------------------------------------------------------------------
# buildCountMatrix
# Create a matrix of reads count
#
# Input:
#   files : a vector of files names
#   sampleLabel : a vector of sample names
#   projectPath: path to the project directory
#
# Output:
#   countMatrix : a reads count matrix
#
# Original author : Vivien DESHAIES
# -----------------------------------------------------------------------------
buildCountMatrix = function(files, sampleLabel, expHeader) {
  
  # read first file and create countMatrix
  countMatrix = read.table(files[1],
                           header = expHeader,
                           stringsAsFactors = F,
                           quote = "")
  colnames(countMatrix) = c("id", "count1")
  
  # read and merge all remaining files with the first
  for (i in 2:length(files)) {
    
    # read files
    exp = read.table(files[i],
                     header = expHeader,
                     stringsAsFactors = F,
                     quote = "")
    
    # lowercase exp columns names
    colnames(exp) = c("id", paste0("count", i))
    
    # merge file data to count matrix by id
    countMatrix = merge(countMatrix, exp, by = "id", suffixes = "_")
  }
  
  # name rows
  rownames(countMatrix) = countMatrix[,1]
  
  # delete first row containing row names
  countMatrix = countMatrix[,-1]
  
  # name columns
  colnames(countMatrix) = sampleLabel
  return(countMatrix)
}

###############################################################################
# -----------------------------------------------------------------------------
# saveRawCountMatrix
#
#   A specific treatment is needed for the raw count matrix because the object
#   DESeq2::counts(dds) does not include the names of the samples as column names
#
#   input: dds -> DESeq object
#          fileName -> character (name of the file to output)
#   output: rawCountMatrix -> file
#
# -----------------------------------------------------------------------------

saveRawCountMatrix = function(dds,fileName) {
  
  countMatrix = DESeq2::counts(dds)
  countMatrix = cbind(countMatrix, row.names(countMatrix))
  
  # Add the samples names as column names + add column Id
  colData = SummarizedExperiment::colData(dds)
  colnames(countMatrix) = c(colData$Name, "Id")
  
  # Put Ids on the first column
  countMatrix = countMatrix[, c("Id", colData$Name)]
  
  # Write the count matrix in a file
  write.table(countMatrix,
              file = paste0(fileName),
              sep = "\t",
              row.names = FALSE,
              quote = FALSE)
}

###############################################################################
# -----------------------------------------------------------------------------
# saveTable
#
#   input: dataframe -> data frame (count matrix)
#          fileName -> character (name of the file to output)
#   output: countMatrix -> file
# -----------------------------------------------------------------------------

saveTable = function(dataframe, fileName) {
  # Add column Id
  column_names = colnames(dataframe)
  dataframe = cbind(row.names(dataframe), dataframe)
  colnames(dataframe) = c("Id", column_names)
  
  # Write the table in a file
  write.table(dataframe,
              file = paste0(fileName),
              sep = "\t",
              row.names = FALSE,
              quote = FALSE)
}

# -----------------------------------------------------------------------------
# plotDendrogram
# Create a ggplot object corresponding to a dendrogram
#
# Input:
#   hc_data : object obtained with the ggdendro::dendro_data function
#   fig_title : figure title
#   col_vector: named vector of colors
#
# Output:
#   p : ggplot object
#
# Original author : Audrey ONFROY
# -----------------------------------------------------------------------------
plotDendrogram = function(hc_data,
                          fig_title = "",
                          col_vector = colors_condition) {
  
  p = ggplot2::ggplot() +
    # Draw dendrogram
    ggplot2::geom_segment(data = ggdendro::segment(hcdata), 
                          aes(x = x, y = y, xend = xend, yend = yend)) +
    # Leaf label
    ggplot2::geom_label(data = ggdendro::label(hcdata), 
                        aes(x = x, y = y, label = label),
                        linewidth = 0, fill = NA,
                        size = 3, angle = 90, hjust = 1,
                        nudge_y = -0.02*max(ggdendro::segment(hcdata)$yend)) +
    # Leaf colors
    ggplot2::geom_point(data = ggdendro::label(hcdata),
                        ggplot2::aes(x = x, y = y, col = Condition),
                        size = 3, pch = 19) +
    ggplot2::scale_color_manual(values = col_vector,
                                breaks = names(col_vector)) +
    # Style
    ggplot2::labs(y = "Distance",
                  title = fig_title) +
    ggplot2::coord_cartesian(clip = "off") +
    ggplot2::theme_minimal() +
    ggplot2::theme(axis.text.x = element_blank(),
                   axis.ticks.x = element_blank(),
                   axis.title.x = element_blank(),
                   axis.line.y = element_line(),
                   axis.ticks.y = element_line(),
                   panel.grid = element_blank(),
                   plot.title = element_text(hjust = 0.5),
                   plot.margin = ggplot2::unit(c(0.5, 0.5, max(nchar(hcdata$labels$label))/6, 0.5), "cm"),
                   # Transparent background
                   panel.background = ggplot2::element_blank(),
                   plot.background = ggplot2::element_blank(),
                   legend.background = ggplot2::element_blank())
  
  return(p)
}

# -----------------------------------------------------------------------------
# plotPCA
# Create a ggplot object corresponding to the PCA plot of individuals
#
# Input:
#   pca_df : dataframe
#   fig_title : figure title
#   col_vector: named vector of colors
#
# Output:
#   p : ggplot object
# -----------------------------------------------------------------------------
plotPCA = function(pca_df,
                   fig_title = "",
                   col_vector = colors_condition) {
  
  p = ggplot2::ggplot(pca_df, aes(x = Dim.1, y = Dim.2, col = Condition, label = sample)) +
    ggplot2::geom_hline(yintercept = 0, lty = 2) +
    ggplot2::geom_vline(xintercept = 0, lty = 2) +
    # Samples
    ggplot2::geom_point(size = 3, pch = 19) +
    ggrepel::geom_label_repel(fill = NA, label.size = 0) +
    ggplot2::scale_color_manual(values = col_vector,
                                breaks = names(col_vector)) +
    # Style
    ggplot2::labs(x = paste0("Dim1 (", round(pcaCount$eig[1, 2], 2), "%)"),
                  y = paste0("Dim2 (", round(pcaCount$eig[2, 2], 2), "%)"),
                  title = fig_title) +
    ggplot2::theme_minimal() +
    ggplot2::theme(plot.title = element_text(hjust = 0.5),
                   # Transparent background
                   panel.background = ggplot2::element_blank(),
                   plot.background = ggplot2::element_blank(),
                   legend.background = ggplot2::element_blank())
  
  return(p)
}

# -----------------------------------------------------------------------------
# plotBarplot
# Create a ggplot object corresponding to a barplot
#
# Input:
#   df : dataframe
#   y_column : column in df for the Y-axis
#   y_title : Y-axis title
#   fig_title : figure title
#   col_vector: named vector of colors
#
# Output:
#   p : ggplot object
# -----------------------------------------------------------------------------
plotBarplot = function(df,
                       y_column,
                       y_title = "",
                       fig_title = "",
                       col_vector = colors_condition) {
  df$y_value = df[, y_column]
  
  p = ggplot2::ggplot(df, aes(x = sample, y = y_value, fill = Condition)) +
    ggplot2::geom_bar(stat = "identity", col = "black") +
    ggplot2::scale_fill_manual(values = col_vector,
                               breaks = names(col_vector)) +
    # Style
    ggplot2::labs(x = "",
                  y = y_title,
                  title = fig_title) +
    ggplot2::scale_y_continuous(expand = c(0,0)) +
    ggplot2::theme_classic() +
    ggplot2::theme(plot.title = element_text(hjust = 0.5),
                   axis.text.x = element_text(angle = 90, hjust = 1),
                   # Transparent background
                   panel.background = ggplot2::element_blank(),
                   plot.background = ggplot2::element_blank(),
                   legend.background = ggplot2::element_blank())
  
  return(p)
}

# -----------------------------------------------------------------------------
# plotBoxplot
# Create a ggplot object corresponding to a boxplot
#
# Input:
#   df : dataframe
#   y_column : column in df for the Y-axis
#   y_title : Y-axis title
#   fig_title : figure title
#   colors_condition: named vector of colors
#
# Output:
#   p : ggplot object
# -----------------------------------------------------------------------------
plotBoxplot = function(df,
                       y_column,
                       y_title = "",
                       fig_title = "",
                       col_vector = colors_condition) {
  df$y_value = df[, y_column]
  
  p = ggplot2::ggplot(df, aes(x = sample, y = y_value, fill = Condition)) +
    ggplot2::geom_boxplot(col = "black") +
    ggplot2::scale_fill_manual(values = col_vector,
                               breaks = names(col_vector)) +
    # Style
    ggplot2::labs(x = "",
                  y = y_title,
                  title = fig_title) +
    ggplot2::scale_y_continuous(expand = c(0,0)) +
    ggplot2::theme_classic() +
    ggplot2::theme(plot.title = element_text(hjust = 0.5),
                   axis.text.x = element_text(angle = 90, hjust = 1),
                   # Transparent background
                   panel.background = ggplot2::element_blank(),
                   plot.background = ggplot2::element_blank(),
                   legend.background = ggplot2::element_blank())
  
  return(p)
}

# -----------------------------------------------------------------------------
# expressionColors
# Define a color panel for the heatmap
#
# Input:
#   range_values : a numerical vector with two values (min and max)
#
# Output:
#   expression_colors : a color or a function built with circlize::colorRamp2
# -----------------------------------------------------------------------------
expressionColors = function(range_values) {
  min_value = range_values[1]
  max_value = range_values[2]
  
  # min_value < 0 < max_value
  if (min_value < 0 & 0 < max_value) {
    expression_colors = circlize::colorRamp2(
      breaks = c(min_value, 0, max_value),
      colors = c("#224897", "#F7F7F7", "#D32A24"))
  } else
    
    # min_value < max_value < 0
    if (min_value < max_value & max_value < 0) {
      expression_colors = circlize::colorRamp2(
        breaks = c(min_value, 0),
        colors = c("#224897", "#F7F7F7"))
    } else
      
      # 0 < min_value < max_value
      if (0 < min_value & min_value < max_value) {
        expression_colors = circlize::colorRamp2(
          breaks = c(0, max_value),
          colors = c("#F7F7F7", "#D32A24"))
      } else {
        
        # min_value == max_value
        expression_colors = "#F7F7F7"
      }
  
  return(expression_colors)
}

# -----------------------------------------------------------------------------
# buildContrast
# Define the contrast vector to extract DE results from the dds object
#
# Input:
#   dds       : DESeqDataSet
#   group1    : list of character vectors, for instance:
#                           list(c("sexe", "F"), c("genotype", "HETERO")
#   group2    : list of character vectors, for instance:
#                           list(c("sexe", "F"), c("genotype", "HOMO")
#   weighted  : boolean     whether to weight the contrast matrix by the
#                           number of replicates per condition or not
#
# Output: a list of 3 elements
#   contrast        : a numerical vector
#   samples_group1  : a character vector containing sample names in group 1
#   samples_group2  : a character vector containing sample names in group 2
#
# Source:
# Modified from: https://www.r-bloggers.com/2024/05/a-guide-to-designs-and-contrasts-in-deseq2/
# -----------------------------------------------------------------------------
buildContrast = function(dds,
                         group1,
                         group2,
                         weighted = FALSE){
  #-- Model matrix for each sample 
  mod_mat = model.matrix(DESeq2::design(dds),
                         SummarizedExperiment::colData(dds))
  
  #-- Samples in group1
  # A list of TRUE / FALSE
  # Each element of the list is of length the number of samples
  # The i-th element of the list corresponds to the i-th condition in group1
  grp1_rows = lapply(c(1:length(group1)), FUN = function(group1_i) {
    group1_i_name = group1[[group1_i]][1]
    group1_i_values = group1[[group1_i]][2:length(group1_i)]
    
    samples_in_group1_i = (SummarizedExperiment::colData(dds)[[group1_i_name]] %in% group1_i_values) # TRUE / FALSE
    
    return(samples_in_group1_i)
  })
  grp1_rows = Reduce(function(x, y) x & y, grp1_rows)
  
  grp2_rows = lapply(c(1:length(group2)), FUN = function(group2_i) {
    group2_i_name = group2[[group2_i]][1]
    group2_i_values = group2[[group2_i]][2:length(group2_i)]
    
    samples_in_group2_i = (SummarizedExperiment::colData(dds)[[group2_i_name]] %in% group2_i_values) # TRUE / FALSE
    
    return(samples_in_group2_i)
  })
  grp2_rows = Reduce(function(x, y) x & y, grp2_rows)
  
  #-- Model matrix with only samples of interest
  mod_mat1 = mod_mat[grp1_rows, , drop = FALSE]
  mod_mat2 = mod_mat[grp2_rows, , drop = FALSE]
  
  #-- Consider the weight or not
  if(!weighted){
    mod_mat1 = mod_mat1[!duplicated(mod_mat1), , drop = FALSE]
    mod_mat2 = mod_mat2[!duplicated(mod_mat2), , drop = FALSE]
  }
  
  #-- Output
  contrast = colMeans(mod_mat1)-colMeans(mod_mat2)
  samples_group1 = colnames(dds)[grp1_rows]
  samples_group2 = colnames(dds)[grp2_rows]
  output = list(contrast = contrast,
                samples_group1 = samples_group1,
                samples_group2 = samples_group2)
  
  return(output)
}

# -----------------------------------------------------------------------------
# missingSubstring
# Identify and color the variable of interest in the context of nested models
#
# Input:
#   A : the text behind a formula
#   B : the text behind a formula
#
# Output:
#   A_color : A, with in red, the substring from A which is absent in B
# -----------------------------------------------------------------------------
missingSubstring = function(A, B) {
  # Remove the ~
  A = base::substr(x = A, start = 2, stop = nchar(A))
  B = base::substr(x = B, start = 2, stop = nchar(B))
  
  # Split by +
  A_split = stringr::str_split(string = A,
                               pattern = "\\+")[[1]]
  
  B_split = stringr::str_split(string = B,
                               pattern = "\\+")[[1]]
  
  # Difference
  id_absent_in_B = which(!(A_split %in% B_split))
  colored_substring = paste0("<span style='color:red'>",
                             A_split[id_absent_in_B],
                             "</span>")
  
  # Rebuild A
  A_color = stringr::str_replace(string = A,
                                 pattern = A_split[id_absent_in_B],
                                 replacement = colored_substring)
  A_color = paste0("~", A_color)
  
  return(A_color)
}

projectName     = params$projectName
designPath      = params$designPath
comparisonPath  = params$comparisonPath
correspPath     = params$correspPath
diffanaTest     = as.logical(toupper(params$diffanaTest))
expHeader       = as.logical(toupper(params$expHeader))
deseqModel      = params$deseqModel
nestedModel     = params$nestedModel
sizeFactorType  = params$sizeFactorType
fitType         = params$fitType
statisticTest   = params$statisticTest
weightContrast  = params$weightContrast
padjLimit       = params$padjLimit
log2FCLimit     = params$log2FCLimit
plotInteractive = params$plotInteractive
leaveOnError    = params$leaveOnError
saveRDS         = params$saveRDS

# Display them as a table
c("Project name (projectName)" = projectName,
  "Path to the design file (designPath)" = designPath,
  "Path to the comparison file (comparisonPath)" = comparisonPath,
  "Path to the correspondence file (correspPath)" = correspPath,
  "Whether to perform the differential expression analysis or not (diffanaTest)" = diffanaTest,
  "Whether the design table has a header or not (expHeader)" = expHeader,
  "DESeq2 model (deseqModel)" = deseqModel,
  "Whether to perform a test of nested models (nestedModel)" = nestedModel,
  "Size factor type (sizeFactorType)" = sizeFactorType,
  "Fit type (fitType)" = fitType,
  "Statistic test (statisticTest)" = statisticTest,
  "Whether to weight the contrast vector by the number of samples or not (weightContrast)" = weightContrast,
  "Adjusted p-value threshold to consider a gene as differentially expressed (padjLimit)" = padjLimit,
  "Absolute log2 fold-change threshold to consider a gene as differentially expressed (log2FCLimit)" = log2FCLimit,
  "Prefix to save files (prefix)" = prefix,
  "Whether to make interactive plots or not (plotInteractive)" = plotInteractive,
  "Whether to stop the rendering in case of error (leaveOnError)" = leaveOnError,
  "Whether to save the data as RDS files (saveRDS)" = saveRDS) %>%
  data.frame(Parameter = names(.),
             Value     = .,
             row.names = NULL) %>%
  knitr::kable("html", escape = FALSE) %>%
  kableExtra::kable_styling("striped", full_width = FALSE) %>%
  kableExtra::column_spec(1, bold = TRUE)

design = read.table(designPath,
                    sep = "\t",
                    header = TRUE,
                    dec = ".",
                    stringsAsFactors = FALSE)

# About a warning in DESeq2
design$Condition = stringr::str_replace_all(
  design$Condition,
  pattern = "-",
  replacement = "_")

# Add Reference if does not exist
if (!("Reference" %in% colnames(design))) {
  design$Reference = Inf
  complexMode = TRUE
} else {
  complexMode = FALSE
}

# Reorder samples based on Condition
design = design %>%
  dplyr::arrange(Condition, RepTechGroup, Name) %>%
  dplyr::mutate(Condition = factor(Condition, levels = unique(Condition))) %>%
  dplyr::mutate(RepTechGroup = factor(RepTechGroup, levels = unique(RepTechGroup))) %>%
  dplyr::mutate(Name = factor(Name, levels = unique(Name)))

# Remove columns out of interest
design = design[, !(colnames(design) %in% c("SampleId", "Description", "Date", "FastqFormat"))]

# Display
head(design) %>%
  dplyr::mutate_all(as.character) %>%
  rbind.data.frame(., "..." = rep("...", ncol(.))) %>%
  knitr::kable("html") 

colors_condition = unique(design$Condition) %>%
  setNames(nm = .,
           grDevices::hcl.colors(length(.), "Spectral"))

data.frame(Color = unname(colors_condition),
           Condition = names(colors_condition)) %>%
  knitr::kable("html",
               col.names = c("Color (hex code)", "Condition")) %>%
  kableExtra::column_spec(1, color = colors_condition) %>%
  kableExtra::kable_styling(full_width = FALSE)

count_mat = buildCountMatrix(files = design$expressionFile,
                             sampleLabel = design$Name,
                             expHeader = expHeader)

# Reorder columns
count_mat = count_mat[, levels(design$Name)]

# Define an order to display the matrix
top5_to_print = order(count_mat[, 1], decreasing = TRUE)[c(1:5)]

# Display
count_mat[top5_to_print, c(1:5)] %>%
  rbind.data.frame(., "..." = rep("...", ncol(.))) %>%
  knitr::kable("html") %>%
  kableExtra::kable_styling(full_width = FALSE)

dds = DESeq2::DESeqDataSetFromMatrix(
  countData = count_mat,
  colData = design,
  design = as.formula(deseqModel))

cat('## Identifiers to symbols\n\n')

cat('We load the table containing the correspondence between genes identifiers and symbols.\n')

correspIdSymbol = read.delim(file = correspPath) %>%
  # Force column names
  `colnames<-`(c("ID", "symbol")) %>%
  # Unique ID
  dplyr::filter(!duplicated(ID)) %>%
  # Force all ID to be present
  dplyr::left_join(x = data.frame(ID = rownames(count_mat)),
                   y = .,
                   by = "ID") %>%
  # Change missing symbols to ID
  dplyr::mutate(symbol = ifelse(is.na(symbol),
                                yes = ID,
                                no = symbol))
rownames(correspIdSymbol) = correspIdSymbol$ID

# Display
head(correspIdSymbol) %>%
  rbind.data.frame(., "..." = rep("...", ncol(.))) %>%
  knitr::kable("html",
               row.names = FALSE) %>%
  kableExtra::kable_styling(full_width = FALSE)

filename = paste0(prefix, projectName, "-normalisation_rawCountMatrix.tsv")

saveRawCountMatrix(dds,
                   file = filename)

cat("We also save the object in RDS format. The file name is the same. The extension is '.rds'.\n")

saveRDS(dds,
        file = stringr::str_replace(string = filename,
                                    pattern = "tsv$",
                                    replacement = "rds"))

dist.mat = dist(t(count_mat))
hc = stats::hclust(dist.mat)
hcdata = ggdendro::dendro_data(hc, type = "rectangle")

# Add grouping for colors
hcdata$labels = dplyr::left_join(x = hcdata$labels,
                                 y = design[, c("Name", "Condition")],
                                 by = c("label" = "Name"))

plotDendrogram(hc_data,
               fig_title = paste0("Dendrogram on the raw count matrix - ",
                                  projectName))

pcaCount = count_mat %>%
  t() %>%
  FactoMineR::PCA(., graph = FALSE)

pcaCount$ind$coord %>%
  as.data.frame() %>%
  dplyr::mutate(sample = rownames(.)) %>%
  dplyr::left_join(x = .,
                   y = unique(design[, c("Name", "Condition")]),
                   by = c("sample" = "Name")) %>%
  plotPCA(.,
          fig_title =  paste0("PCA of the raw count matrix - ", projectName))

colMeans(count_mat == 0) %>%
  data.frame(sample = factor(names(.), levels = levels(design$Name)),
             prop_null_counts = .) %>%
  dplyr::mutate(prop_null_counts = 100*prop_null_counts) %>%
  dplyr::left_join(x = .,
                   y = design[, c("Name", "Condition")],
                   by = c("sample" = "Name")) %>%
  # Plot
  plotBarplot(.,
              y_column = "prop_null_counts",
              y_title = "Proportion of null counts (%)",
              fig_title = paste0("Proportion of null counts per sample - ",
                                 projectName))

colSums(DESeq2::counts(dds)) %>%
  data.frame(sample = factor(names(.), levels = levels(design$Name)),
             read_counts = .) %>%
  dplyr::left_join(x = .,
                   y = design[, c("Name", "Condition")],
                   by = c("sample" = "Name")) %>%
  # Plot
  plotBarplot(.,
              y_column = "read_counts",
              y_title = "Total read counts",
              fig_title = paste0("Number of reads per sample - ", projectName))

log2(DESeq2::counts(dds)+1) %>%
  reshape2::melt() %>%
  `colnames<-`(c("gene_ID", "sample", "log2_counts")) %>%
  dplyr::mutate(sample = factor(sample, levels = levels(design$Name))) %>%
  dplyr::left_join(x = .,
                   y = design[, c("Name", "Condition")],
                   by = c("sample" = "Name")) %>%
  # Plot
  plotBoxplot(.,
              y_column = "log2_counts",
              y_title = "log_2 (counts+1)",
              fig_title = paste0("Raw counts distribution per sample - ", projectName))

dds = DESeq2::estimateSizeFactors(dds,
                                  type = sizeFactorType)

# Display
DESeq2::counts(dds, normalized = TRUE)[top5_to_print, c(1:5)] %>%
  rbind.data.frame(., "..." = rep("...", ncol(.))) %>%
  knitr::kable("html") %>%
  kableExtra::kable_styling(full_width = FALSE)

filename = paste0(prefix, projectName,
                  "-normalisation_normalisedCountMatrix.tsv")

saveTable(DESeq2::counts(dds, normalized = TRUE),
          fileName = filename)

cat("We also save the object in RDS format. The file name is the same. The extension is '.rds'.\n")

saveRDS(dds,
        file = stringr::str_replace(string = filename,
                                    pattern = "tsv$",
                                    replacement = "rds"))

all_interactions = combn(all_terms, m = 2) %>%
  apply(., 2, FUN = paste, collapse = ":")

# Models with all interactions
nested_model = list()
i_start = 0

for (i in c(1:length(all_interactions))) {
  full_formula = c(all_terms, all_interactions) %>%
    paste(., collapse = "+")
  reduced_formula = c(all_terms, all_interactions[-i]) %>%
    paste(., collapse = "+")
  
  nested_model[[i_start + i]] = c(full_formula,
                                  reduced_formula)
}

# Models with a single interaction
i_start = i_start + length(all_interactions)

for (i in c(1:length(all_interactions))) {
  full_formula = c(all_terms, all_interactions[i]) %>%
    paste(., collapse = "+")
  reduced_formula = all_terms %>%
    paste(., collapse = "+")
  
  nested_model[[i_start + i]] = c(full_formula,
                                  reduced_formula)
}

# Models with all terms
i_start = i_start + length(all_terms)

for (i in c(1:length(all_terms))) {
  full_formula = all_terms %>%
    paste(., collapse = "+")
  reduced_formula = all_terms[-i] %>%
    paste(., collapse = "+")
  
  nested_model[[i_start + i]] = c(full_formula,
                                  reduced_formula)
}

# Models with 1 term
i_start = i_start + length(all_terms)

for (i in c(1:length(all_terms))) {
  full_formula = all_terms[i]
  reduced_formula = "1"
  
  nested_model[[i_start + i]] = c(full_formula,
                                  reduced_formula)
}

# Make as a dataframe
nested_model = nested_model %>%
  do.call(rbind.data.frame, .) %>%
  unique %>%
  `colnames<-`(c("full_formula", "reduced_formula")) %>%
  dplyr::mutate(full_formula = paste0("~", full_formula)) %>%
  dplyr::mutate(reduced_formula = paste0("~", reduced_formula))

nested_model[, c("nb_genes", "decision")] = apply(nested_model, 1, FUN = function(models) {
  # Get the formula
  full_formula    = as.formula(models[[1]])
  reduced_formula = as.formula(models[[2]])
  
  # DESeq2 object with full formula
  dds = DESeq2::DESeqDataSetFromMatrix(
    countData = DESeq2::counts(dds, normalized = FALSE),
    colData = design,
    design = full_formula)
  
  # Nested model
  dds_res = DESeq2::DESeq(dds,
                          test = "LRT",
                          full = full_formula,
                          reduced = reduced_formula)
  
  # Number of DE genes between models
  res_LRT = DESeq2::results(dds_res)
  signif_or_not = (res_LRT$padj < padjLimit) # TRUE / FALSE
  nb_genes = sum(signif_or_not, na.rm = TRUE)
  prop_genes = mean(signif_or_not, na.rm = TRUE)  # between 0 and 1
  percent_genes = 100*(round(prop_genes, 2))      # between 0 and 100
  
  # Decision
  decision = "irrelevant"
  if (prop_genes > 0.05) {
    decision = "relevant"
  }
  
  # Output
  nb_prop_genes = paste0(nb_genes, " (", percent_genes, "%)")
  
  return(c(nb_prop_genes, decision))
}) %>% t() %>%
  rbind.data.frame()

# Color the variable being evaluated
nested_model$full_formula = apply(
  X = nested_model,
  MARGIN = 1,
  FUN = function(row_i) {
    missingSubstring(A = row_i[[1]],
                     B = row_i[[2]]) 
  })

# Display
nested_model %>%
  knitr::kable("html",
               col.names = c("Full formula",
                             "Reduced formula",
                             "Nb DE genes",
                             "The colored term is..."),
               escape = FALSE,
               row.names = FALSE) %>%
  kableExtra::kable_styling(full_width = FALSE)

dds = DESeq2::estimateDispersions(dds,
                                  fitType = fitType)

DESeq2::plotDispEsts(
  dds,
  main = paste0("Dispersion estimation scatter plot\n", projectName))

dds = DESeq2::DESeq(dds,
                    test = statisticTest,
                    fitType = fitType,
                    betaPrior = FALSE)

cat("## Save object\n\n")

cat("We save the dds object in RDS format.\n")

filename = paste0(prefix, projectName,
                  "-processed_dds_object.rds")

saveRDS(dds,
        file = filename)

cat(paste0("Output file name:<br>*", filename, "*\n"))