01_normDiffana.Rmd

---
params:
  projectName: "GSE107401"
  designPath: "./project_GSE107401/deseq2_GSE107401-deseq2Design.txt"
  comparisonPath: "./project_GSE107401/deseq2_GSE107401-comparisonFile.txt"
  diffanaTest: TRUE
  expHeader: TRUE
  deseqModel: "~Condition"
  sizeFactorType: "ratio"
  fitType: "parametric"
  statisticTest: "Wald"
  weightContrast: FALSE
  prefix: "./project_GSE107401/deseq2_"
  plotInteractive: TRUE
  logoUrl: "https://www.outils.genomique.biologie.ens.fr/eoulsan/images/eoulsan.png"
  authorName: "Eoulsan"
  authorEmail: "eoulsan@bio.ens.psl.eu"
  leaveOnError: TRUE
title: "<center>Differential expression analysis</center>"
subtitle: "<center>`r paste0('Project: ', params$projectName)`</center>"
date: "`r invisible(Sys.setlocale('LC_TIME', 'C')) ; paste0('Analysis run on: ', format(Sys.time(), '%d %B %Y'))`
<br><span class='glyphicon glyphicon-pushpin'></span>`r paste0(' Project: ', params$projectName)`
<br><span class='glyphicon glyphicon-menu-right'></span> Version 2.1
<br><span class='glyphicon glyphicon-link'></span> Generated by <a href='https://www.outils.genomique.biologie.ens.fr/eoulsan/'>Eoulsan</a>"
# above is a small hack to add new lines in the postample without editing the template
author:
  - name: "`r params$authorName`"
    email: "`r params$authorEmail`"
output:
  rmdformats::readthedown:
    code_folding: show
    code_download: true
    number_sections: true
    toc_depth: 3
---

<!-- CSS to customize the colors -->
```{css style_section1, echo = FALSE}
body {
text-align: justify
}

/* Width of the main div */
#content {
max-width: 80%;
}
```

```{r style_section2, echo = FALSE, results = "asis"}
# NOTE: A CSS chunk cannot interpret the R input parameter. In order to parse the 'params$logoUrl' parameter, we require an R chunk. Therefore, the CSS content is displayed through an R chunk.

# Header in the sidebar
cat(paste0("
<style>
#sidebar h2 {
background-color:#A3ADB8;
background-size: 32%;
background-repeat: no-repeat;
background-position: center bottom;
background-image:url('", params$logoUrl, "');
background-origin: content-box;
height: 13%;
}
</style>
"))
```

```{css style_section3, echo = FALSE}
/* Sidebar body */
#sidebar {
background:#535F6A;
}

/* Small space about the authoring section */
#postamble {
border-top:solid 10px #A3ADB8;
}

/* Titles */
h1, h2, h3, h4, h5, h6, legend {
color: #000038;
}

/* Two-column mode */
.cols {
display: flex;
align-items: center;
}

.col-center {
justify-content: center;
}

.col-right {
margin-left: auto;
}

/* Button */
.btn-custom {
font-family: inherit;
font-size: inherit;
color: inherit;
background-color: #E3E3E3;
border-color: #E3E3E3;
padding: 1px 2px;
border-radius: 5px;
border-style: solid;
color: #000038;
}

.btn-custom:hover {
background-color: #535F6A;
border-color: #535F6A;
color: #FCFCFC;
}
```

<!-- Automatically computes and prints in the output the running time for any code chunk -->
```{r, echo = FALSE}
hooks = knitr::knit_hooks$get()
hook_foldable = function(type) {
  force(type)
  function(x, options) {
    res = hooks[[type]](x, options)
    
    if (isFALSE(options[[paste0("fold_", type)]])) return(res)
    
    paste0(
      "<details><summary>", "\U25BC show", "</summary>\n\n",
      res,
      "\n\n</details>"
    )
  }
}
knitr::knit_hooks$set(
  output = hook_foldable("output"),
  plot = hook_foldable("plot")
)
```

<!-- knit setup -->
```{r knit_setup, echo = FALSE}
fig_save_folder = paste0(params$prefix, "figures/")

knitr::opts_chunk$set(
  echo = TRUE,
  eval = TRUE,
  cache.lazy = FALSE,
  message = FALSE,
  warning = FALSE,
  fold_output = FALSE,
  fold_plot = FALSE,
  fig.align = "center",
  class.source = "fold-hide",
  # R base plots on transparent background
  dev.args = list(bg = 'transparent'),
  # Export figures
  dev = c('png', 'pdf'),
  fig.path = fig_save_folder,
  pdf.options(encoding = "ISOLatin9.enc"))
```

------------------------------------------------------------------------

The aim of this document is to present the data and differential expression analysis between samples from the project **`r params$projectName`**.

The main steps are:

- Data preparation
- Quality control
- Normalisation
- Differential expression analysis

```{r echo = FALSE, results = "asis"}
# NOTE: This section is very ugly. Due to the conditional section, we need to evaluate in advance whether the section will be executed or not, in order to create a short link to the conditionally-executed section. This requires the design file to be loaded. This section is hidden in the final rendered notebook. We only make calculations here to print a message. But, calculations are performed again in a visible way in the next sections. Not elegant and (low) resource-consuming.

text_to_print = "For a quick access to the figures, go to"

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

if (!("Reference" %in% colnames(design_tmp))) {
  # Complex design, no reference column
  design_tmp = c(0,1)
} else {
  # Simple design, check if there are diverse Reference values above 0
  design_tmp = design_tmp |>
    dplyr::filter(Reference >= 0) |>
    dplyr::pull(Reference) |>
    unique()
}

if (params$diffanaTest & length(design_tmp) > 1) {
  text_to_print = paste0(
    text_to_print, ":<br>",
    "- [Normalisation > Visualisation](#normalisation_visu) for normalised counts<br>",
    "- [Results subsections](#results_subsection) for the pairwise comparisons")
} else {
  text_to_print = paste0(
    text_to_print, " the section [Normalisation > Visualisation](#normalisation_visu).")
}

cat(paste0(text_to_print, "\n"))
```

All figures are saved in PNG and PDF format and accessible in the *`r fig_save_folder`* folder. You can download the code behind this analysis through the "Code" button in the top right corner of the document.

# Settings

The aim of this section is to:

- load the tools of interest
- define custom functions
- load the input parameters

## Environment

We load the tools to handle data and figures.

```{r library, fold_output=TRUE}
library(dplyr)
library(ggplot2)
```

```{r check_deseq2_version, echo = FALSE}
version_tested = c(1, 48, 0) # version 1.48.0

current_version = utils::packageVersion("DESeq2") %>% # A.B.C
  as.character() %>%
  base::strsplit(., split = "[^0-9]+") %>%
  unlist() %>%
  as.numeric()

if (current_version[[1]] < version_tested[[1]] |
    (current_version[[1]] <= version_tested[[1]] &
     current_version[[2]] < version_tested[[2]])) {
  warning("The notebook has been tested using version 1.48.0 of DESeq2",
          " You are using an older version (",
          utils::packageVersion("DESeq2"),
          ") so it cannot be guaranteed that the script will work...")
}
```

## Custom functions

We further define custom functions to avoid repetitions in the code.

```{r custom_functions, class.source = "fold-hide"}
# -----------------------------------------------------------------------------
# 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)
}
```

These functions are listed below.

```{r print_functions, echo = FALSE, results = 'asis'}
cat(paste0("- ", setdiff(ls(), c("hooks", "params", "hook_foldable", "design_tmp",
                                 "fig_save_folder", "text_to_print"))),
    sep = "\n")
```

## Parameters

The table below summarizes all the input parameters.

```{r input_parameter, class.source = "fold-hide"}
projectName     = params$projectName
designPath      = params$designPath
comparisonPath  = params$comparisonPath
diffanaTest     = as.logical(toupper(params$diffanaTest))
expHeader       = as.logical(toupper(params$expHeader))
deseqModel      = params$deseqModel
sizeFactorType  = params$sizeFactorType
fitType         = params$fitType
statisticTest   = params$statisticTest
weightContrast  = params$weightContrast
prefix          = params$prefix
plotInteractive = params$plotInteractive
leaveOnError    = params$leaveOnError

# Display them as a table
c("Project name (projectName)" = projectName,
  "Path to the design file (designPath)" = designPath,
  "Path to the comparison file (comparisonPath)" = comparisonPath,
  "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,
  "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,
  "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) %>%
  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)
```

# Data preparation

The purpose of this section is to:

* load the design file, containing metadata about each sample ;
* load the individual sample count data, assemble them as a gene-by-sample count matrix, and save it as a TSV file ; and,
* build a dataset in a format compatible with the differential expression analysis (DESeq2-related format).

## Load metadata

We load the design file and display the 6 first rows. 

```{r load_design}
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) %>%
  knitr::kable("html") 
```


Description of the columns:

```{r print_metadata_help, echo = FALSE, results = 'asis'}
the_message = c(Name = "- **`Name`**: Unique identifier of each technical replicate.",
                expressionFile = "- **`expressionFile`**: File path to the count expression data. The **TSV** file contains two columns: one for the gene identifiers and one for the number of reads mapping to their genomic location. For compatibility with Excel, please use the **XLSX** file instead of the TSV file ([read more](https://pmc.ncbi.nlm.nih.gov/articles/PMC8357140/)). Furthermore, the XLSX file contains additional columns, such as a description of each gene.",
                Reads = "- **`Reads`**: File path to the sequencing reads, in FASTQ format ([read more](https://en.wikipedia.org/wiki/FASTQ_format)).",
                RepTechGroup = "- **`RepTechGroup`**: Samples with the same value in the `RepTechGroup` column are **technical replicates**. This occurs when the sequencing library is split over sequencing lanes. Each lane generates a technical replicate from the same biological library.",
                Condition = "- **`Condition`**: Samples with the same value in the `Condition` column are **biological replicates**. They are produced by the experimenter.",
                Reference = "- **`Reference`**: Numerical values used to define the pairwise comparisons between conditions in the differential expression analysis.")
the_message = the_message[names(the_message) %in% colnames(design)]

cat(paste0(the_message), sep = "\n")
```

We associate a color to each condition.

```{r add_colors}
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)
```

## Build count matrix

From the individual count data, we build a gene-by-sample count matrix. The count value for a given gene corresponds to the number (integer value) of sequencing reads that map to its genomic region. We display the top left corner of this matrix, meaning 5 genes and 5 samples. The complete matrix is saved as a TSV file ([quick access](#save_count_matrix)).

```{r count_mat}
count_mat = buildCountMatrix(files = design$expressionFile,
                             sampleLabel = design$Name,
                             expHeader = expHeader)

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

# Display
count_mat[c(1:5), c(1:5)] %>%
  knitr::kable("html") %>%
  kableExtra::kable_styling(full_width = FALSE)
```

The count matrix contains **`r nrow(count_mat)`** genes, **`r ncol(count_mat)`** samples and **`r round(100*mean(count_mat == 0),2)`**% of null values.

## Create DESeqDataSet object

We create the DESeqDataSet object that includes the count matrix and the design file.

```{r dds_object}
dds = DESeq2::DESeqDataSetFromMatrix(
  countData = count_mat,
  colData = design,
  design = as.formula(deseqModel))
```

## Save count matrix {#save_count_matrix}

We save the count matrix as a TSV file.

```{r save_count_matrix}
filename = paste0(prefix, projectName, "-normalisation_rawCountMatrix.tsv")

saveRawCountMatrix(dds,
                   file = filename)
```

Output file name:<br>*`r filename`*

# Quality control

The aim of this section is to better apprehend the variability between the technical replicates. For information on the correspondence between the identifiers of technical replicates and sample names, please refer to the design file.

## Raw counts dendrogram

We compute euclidean distances between each pair of samples, considering the raw count values.

```{r dendrogram_raw_count_matrix}
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"))
```

We visualise these distances through a dendrogram. It informs whether technical replicates or closely related conditions are closer together than with other samples.

```{r fig-normalisation_unpooled_clustering, fig.width = 10, fig.height = 6}
plotDendrogram(hc_data,
               fig_title = paste0("Dendrogram on the raw count matrix - ",
                                  projectName))
```

## Raw counts PCA

The **`r ncol(count_mat)`** samples are represented in a space of **`r nrow(count_mat)`** dimensions, corresponding to genes. We perform a Principal Component Analysis (PCA) to visualise samples on a two-dimension space such that it represents the highest variability between samples.

```{r run_pca}
pcaCount = count_mat %>%
  t() %>%
  FactoMineR::PCA(., graph = FALSE)
```

We visualise the samples on the two-dimensional projection. The coordinates of each sample depend on those of the other samples. This projection captures **`r round(pcaCount$eig[2, 3], 2)`**% of the total variability between samples. Therefore, it lacks **`r 100-round(pcaCount$eig[2, 3], 2)`**% of the initial information. In this representation, two points (samples) are close together if they have similar global transcriptomic profiles. We expect the technical and biological replicates to be closer to each other than to the other samples and conditions.

```{r fig-normalisation_unpooled_PCA, fig.width = 10, fig.height = 6}
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))
```

## Null counts barplot

This figure displays the proportion of null counts per sample, meaning, the genes that are not being expressed or not captured by the protocol, within the biological samples. We expect samples from the same condition to be similar.

```{r fig-normalisation_null_counts, fig.width = 10, fig.height = 6}
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))
```

## Raw counts barplot

This figure displays the total number of reads per sample, also known as library size. We can assess if some samples have a low library size compared to others, and whether technical replicates share a similar number of reads or not.

```{r fig-normalisation_barplot_counts, fig.width = 10, fig.height = 6}
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))
```

## Raw counts boxplot

This figure displays, as a boxplot, the number of raw counts for each gene, per sample. It informs whether the genes show a similar count distribution between samples and conditions.

```{r fig-normalisation_boxplot_count, fig.width = 10, fig.height = 6}
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))
```


<!-- Conditional execution in case technical replicates are present -->
<br>

```{r tech_rep_present, echo = FALSE}
tech_rep_present = (max(table(design$RepTechGroup, design$Condition)) > 1)

if (tech_rep_present) {
  the_message = "In this context, there are technical replicates."
} else {
  the_message = "In this context, there is no technical replicate."
  
  # Anyway, we change the sample labels to the biological meaningful ones.
  dds = DESeq2::collapseReplicates(dds,
                                   groupby = dds$RepTechGroup)
}
```

`r the_message`

```{r collapse_replicates, echo = FALSE, results = 'asis', eval = tech_rep_present}
res = knitr::knit_child(
  input = '02_child_collapse_replicates.Rmd',
  quiet = TRUE)

cat(res, sep = '\n')
```

# Normalisation

The aim of this section is to normalise the count matrix, before performing the differential expression analysis. The normalisation method considers the sequencing depth and RNA composition of each sample to compute normalised count values.

For more details, please visit:

[https://hbctraining.github.io/DGE_workshop/lessons/02_DGE_count_normalization.html](https://hbctraining.github.io/DGE_workshop/lessons/02_DGE_count_normalization.html)

## Data handling

We normalise the count matrix. We display the top left corner of this matrix, meaning 5 genes and 5 samples. Unlike count data, normalised count values are not necessarily integers.

```{r estimateSizeFactors}
dds = DESeq2::estimateSizeFactors(dds,
                                  type = sizeFactorType)

DESeq2::counts(dds, normalized = TRUE)[c(1:5), c(1:5)] %>%
  knitr::kable("html") %>%
  kableExtra::kable_styling(full_width = FALSE)
```

We save the normalised count matrix as a TSV file.

```{r save_norm_count_matrix}
filename = paste0(prefix, projectName,
                  "-normalisation_normalisedCountMatrix.tsv")

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

Output file name:<br>*`r filename`*

## Visualisation {#normalisation_visu}

### Normalised counts dendrogram

We compute euclidean distances between each pair of samples, considering the normalised count values.

```{r dendrogram_pooled_normalised_count_matrix}
ddsStabilized = dds %>%
  DESeq2::varianceStabilizingTransformation() %>%
  SummarizedExperiment::assay()
dist.mat = dist(t(ddsStabilized))
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 = unique(design[, c("RepTechGroup", "Condition")]),
  by = c("label" = "RepTechGroup"))
```

We visualise these distances through a dendrogram. It informs whether the samples share a similar transcriptomic profile per condition or not.

```{r fig-normalisation_pooled_clustering, fig.width = 10, fig.height = 6}
plotDendrogram(
  hc_data,
  fig_title = paste0("Dendrogram on the normalised count matrix - ", projectName))
```

### Normalised counts PCA

We perform a PCA for all samples using the normalised count matrix.

```{r run_pca2}
pcaCount = dds %>%
  DESeq2::counts(., normalized = TRUE) %>%
  t() %>%
  FactoMineR::PCA(., graph = FALSE)
```

We visualise the samples on the two-dimensional projection. The coordinates of each sample depend on those of the other samples. This projection captures **`r round(pcaCount$eig[2, 3], 2)`**% of the total variability between samples. Therefore, it lacks **`r 100-round(pcaCount$eig[2, 3], 2)`**% of the initial information. In this representation, two points (samples) are close together if they have similar global transcriptomic profiles. We expect samples from the same biological condition to be closer to each other than to samples from other conditions.

```{r fig-normalisation_pooled_PCA, fig.width = 10, fig.height = 6}
pcaCount$ind$coord %>%
  as.data.frame() %>%
  dplyr::mutate(sample = rownames(.)) %>%
  dplyr::left_join(x = .,
                   y = unique(design[, c("RepTechGroup", "Condition")]),
                   by = c("sample" = "RepTechGroup")) %>%
  plotPCA(.,
          fig_title =  paste0("PCA of the normalised count matrix - ", projectName))
```

### Normalised counts boxplot

This figure displays, as a boxplot, the number of normalised counts for each gene, per sample.

```{r fig-normalisation_normalised_boxplot_count, fig.width = 10, fig.height = 6}
log2(DESeq2::counts(dds, normalized = TRUE)+1) %>%
  reshape2::melt() %>%
  `colnames<-`(c("gene_ID", "sample", "log2_counts")) %>%
  dplyr::mutate(sample = factor(sample, levels = levels(design$RepTechGroup))) %>%
  dplyr::left_join(x = .,
                   y = unique(design[, c("RepTechGroup", "Condition")]),
                   by = c("sample" = "RepTechGroup")) %>%
  # Plot
  plotBoxplot(
    .,
    y_column = "log2_counts",
    y_title = "log_2 (counts+1)",
    fig_title = paste0("Normalised counts distribution per sample - ", projectName))
```

### Most expressed genes plot

We identified the most expressed gene in each sample.

```{r identify_most_expressed_genes}
# preparation of 2 data frame with the same number of column than
# the dds count matrix
maxCounts = DESeq2::counts(dds)[1,]
transcriptNames = DESeq2::counts(dds)[1,]

# for each sample (column)
for (i in 1:ncol(DESeq2::counts(dds))) {
  
  # selection of the maximum number of count
  maxCounts[i] = (max(DESeq2::counts(dds, normalized = TRUE)[,i])/
                    sum(DESeq2::counts(dds, normalized = TRUE)[,i]))*100
  
  # selection of the name of the genes this the maximum of count
  transcriptNames[i] = row.names(
    subset(DESeq2::counts(dds, normalized = TRUE),
           DESeq2::counts(dds, normalized = TRUE)[,i] == 
             max(DESeq2::counts(dds, normalized = TRUE)[,i])))
}
```

We represent the information as a barplot.

```{r fig-normalisation_most_expressed_genes, fig.width = 10, fig.height = 2+ncol(count_mat)/4}
dplyr::left_join(x = data.frame(sample = names(maxCounts),
                                max_counts = maxCounts),
                 y = data.frame(sample = names(transcriptNames),
                                gene_id = transcriptNames),
                 by = "sample") %>%
  dplyr::mutate(sample = factor(sample, levels = levels(design$RepTechGroup))) %>%
  dplyr::left_join(x = .,
                   y = unique(design[, c("RepTechGroup", "Condition")]),
                   by = c("sample" = "RepTechGroup")) %>%
  # Plot
  ggplot2::ggplot(., aes(x = sample, y = max_counts,
                         fill = Condition, label = gene_id)) +
  ggplot2::geom_bar(stat = "identity", col = "black") +
  ggplot2::geom_label(mapping = aes(x = sample, y = 0),
                      linewidth = 0, fill = NA, hjust = 0) +
  ggplot2::scale_fill_manual(values = colors_condition,
                             breaks = names(colors_condition)) +
  # Style
  ggplot2::labs(y = "Proportion of reads (%)",
                title = paste0("Most expressed genes - ", projectName)) +
  ggplot2::scale_y_continuous(expand = c(0,0)) +
  ggplot2::coord_flip() +
  ggplot2::theme_classic() +
  ggplot2::theme(plot.title = element_text(hjust = 0.5),
                 axis.title.y = element_blank(),
                 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())
```


<!-- Conditional execution in case differential expression analysis should be done -->

```{r diffana, echo = FALSE, results = 'asis', eval = diffanaTest}
res = knitr::knit_child(
  input = '03_child_differential_expression.Rmd',
  quiet = TRUE)

cat(res, sep = '\n')
```


# Versions

```{r sessioninfo, echo = FALSE, fold_output = TRUE, results = "hold"}
print("Packages location:")
.libPaths()
print("")
print("Session Information:")
sessionInfo()
```

# Reference

The differential expression analysis uses the `DESeq2` R package ([vignette](https://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html)).

> Love, M.I., Huber, W. & Anders, S., Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 15, 550 (2014). DOI: [10.1186/s13059-014-0550-8](https://doi.org/10.1186/s13059-014-0550-8)