01-metric_summary.Rmd

---
title: "01- QC metrics summary"
author: "Marc Elosua Bayes"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output:
  BiocStyle::html_document:
    toc: yes
    toc_float: yes
    number_sections: yes
    df_print: paged
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, out.width = "100%", fig.align='center', 
                      message = FALSE, warning = FALSE)
options(width = 1200)
```

## Introduction

The objective of this notebook is to perform a basic quality control (QC) analysis of the mapping performed with [spaceranger 1.1.0](https://support.10xgenomics.com/spatial-gene-expression/software/pipelines/latest/using/count).

## Libraries
```{r warning = FALSE, message = FALSE}
library(tidyverse)
library(gt)
library(ggpubr)
library(ggrepel)
```

## Parameters
```{r}
here::dr_here(show_reason = TRUE)

source(here::here("misc/paths.R"))

"{qc}/{plt_dir}" %>%
  glue::glue() %>%
  here::here() %>%
  dir.create(path = .,
           showWarnings = FALSE,
           recursive = TRUE)

"{qc}/{robj_dir}" %>%
  glue::glue() %>%
  here::here() %>%
  dir.create(path = .,
           showWarnings = FALSE,
           recursive = TRUE)
```

## Load data
The data used in this Rmarkdown document comes from **SALAS_25** dataset.
We have 4 different datasets, 1 for each capture area from a visium slides.
```{r}
qcmetric_df <- lapply(id_sp_df$gem_id, function(id) {
  tmp_df <- readr::read_csv(file = here::here(glue::glue("{spaceranger}/{id}/outs/metrics_summary.csv")))
  tmp_df <- tmp_df %>%
    dplyr::mutate(sample_id = id_sp_df[id_sp_df$gem_id == id, ] %>% dplyr::pull(donor_id)) %>%
    dplyr::rename(gem_id = `Sample ID`)
  
  return(tmp_df)
}) %>% dplyr::bind_rows()

colnames(qcmetric_df) <- stringr::str_replace_all(string = colnames(qcmetric_df),
                                                  pattern = " ",
                                                  replacement = "_")

DT::datatable(qcmetric_df)
```

# Gene Expression QC
We will start by showing the three most relevant metrics (number of reads, estimated number of spots under tissue, fraction of reads in spots under tissue, mean reads per spot, fraction of reads mapped to exonic reads, and median genes per spot) obtained by spaceranger for each of the working libraries. This information will give us an idea of the quality of the experiment as well as the sequencing and the mapping steps.
```{r fig.width=14}
top_metrics_gex <- qcmetric_df[, c("gem_id",
                                   "sample_id",
                                   "Number_of_Reads", 
                                   "Number_of_Spots_Under_Tissue", 
                                   "Fraction_Reads_in_Spots_Under_Tissue",
                                   "Mean_Reads_per_Spot",
                                   "Reads_Mapped_Confidently_to_Exonic_Regions",
                                   "Median_Genes_per_Spot")]
top_metrics_gex %>%
  gt::gt() %>%
  gt::fmt_percent(columns = c("Fraction_Reads_in_Spots_Under_Tissue", "Reads_Mapped_Confidently_to_Exonic_Regions"), 
              scale_values = FALSE, decimals = 1) %>%
  gt::fmt_number(columns = "Number_of_Reads", scale_by = 1 / 1E6, pattern = "{x}M") %>% 
  gt::tab_header(
    title = gt::md("**GEX QC metrics**"),
    subtitle = ("spaceranger v1.1.0")
  ) %>%
  gt::cols_label(
    gem_id = gt::md("**GEM ID**"),
    sample_id = gt::md("**SAMPLE ID**"),
    Number_of_Reads = gt::md("**Number of Reads**"),
    Number_of_Spots_Under_Tissue = gt::md("**Number of_Spots Under Tissue**"),
    Fraction_Reads_in_Spots_Under_Tissue = gt::md("**Fraction Reads in Spots Under Tissue**"),
    Mean_Reads_per_Spot = gt::md("**Mean Reads per Spot**"),
    Reads_Mapped_Confidently_to_Exonic_Regions = gt::md("**Fraction of Reads Mapped to Exonic Reads**"),
    Median_Genes_per_Spot = gt::md("**Median Genes per Spot**")
  )
```

## Sequencing QC
First, we will assess the quality of the sequenced libraries before the mapping step. To so do, we will use the "Q30" variables in our dataset, as these report the fraction of the bases with a Q-score >30 for different sequences (barcodes, reads and UMIs).
Q-score is calculated in the following way:
$$Q = -10\log10(p)$$
where *p* is the probability of the base being wrongly called. Thus, bases with a high Q-score are reliable.
```{r fig.height = 5, fig.width = 14}
qc_seq_vars <- c("Q30_Bases_in_Barcode",
                 "Q30_Bases_in_RNA_Read",
                 "Q30_Bases_in_UMI")

gg_qc_seq <- purrr::map(qc_seq_vars, function(var) {
  ggplot2::ggplot(qcmetric_df, ggplot2::aes_string(x = "sample_id", y = var, fill = "sample_id")) +
    ggplot2::geom_col() +
    ggplot2::theme_bw() +
    ggplot2::scale_fill_brewer(palette = "Set2") +
    ggplot2::ylim(0, 1) +
    ggplot2::labs(x = "Libraries (SAMPLE IDs)",
         y = stringr::str_c(stringr::str_replace_all(var, "_", " "), " (%)")) +
    ggplot2::theme(
      axis.title = ggplot2::element_text(size = 12),
      axis.text = ggplot2::element_text(size = 8),
      axis.text.x = ggplot2::element_text(hjust = 1, angle = 45),
      strip.placement = "outside",
      strip.background = ggplot2::element_rect(colour = NA),
      legend.position = "none")
})

cowplot::plot_grid(plotlist = gg_qc_seq,
                   nrow = 1,
                   ncol = 3,
                   align = "hv",
                   axis = "trbl")
```

## Mapping QC
Next, we will check the quality of the mapping step performed by `spaceranger 1.0.0` across libraries. To do so, we will compare the percentage of reads mapped to the genome, and within these mapped reads, the amount of reads mapped to intergenic regions, intronic and exonic regions. We aim to obtain libraries with a high percentage of confidently mapped reads, and specially a high percentage of exonic reads, which correspond with gene expression or RNAs. The reads mapping to intergenic regions suggest contamination of ambient DNA, whereas reads mapping to intronic regions may come from pre-mRNAs or mature spliced isoforms that retain certain introns.
```{r fig.height = 10, fig.width = 14}
qc_map_vars <- c("Reads_Mapped_Confidently_to_Genome",
                 "Reads_Mapped_Confidently_to_Intergenic_Regions",
                 "Reads_Mapped_Confidently_to_Intronic_Regions",
                 "Reads_Mapped_Confidently_to_Exonic_Regions")

gg_qc_map <- purrr::map(qc_map_vars, function(var) {
  ggplot2::ggplot(qcmetric_df, ggplot2::aes_string(x = "sample_id", y = var, fill = "sample_id")) +
    ggplot2::geom_col() +
    ggplot2::theme_bw() +
    ggplot2::scale_fill_brewer(palette = "Set2") +
    ggplot2::ylim(0, 1) +
    ggplot2::labs(x = "Libraries (Sample IDs)",
         y = stringr::str_c(stringr::str_replace_all(var, "_", " "), " (%)")) +
    ggplot2::theme(
      axis.title = ggplot2::element_text(size = 12),
          axis.text = ggplot2::element_text(size = 10),
          axis.text.x = ggplot2::element_text(hjust = 1, angle = 45),
          strip.placement = "outside",
          strip.background = ggplot2::element_rect(colour = NA),
          legend.position = "none")
})
cowplot::plot_grid(plotlist = gg_qc_map,
                   nrow = 2,
                   ncol = 2,
                   align = "hv",
                   axis = "trbl")
```

## Sequencing saturation and depth
After assessing mapped reads, it is important to test which is the sequencing saturation and depth for each library. The sequencing saturation is dependent on the library complexity and sequencing depth. The library complexity is the total number of different transcripts present in the library and it varies between the cell types/tissues, whereas the sequencing depth is the number of paired reads per cell. For this reason, we will plot the number of detected genes as a function of depth (sequenced reads). As sequencing depth increases, more genes are detected, but this function reaches a plateau, whereby more sequenced reads does not result in more detected genes; therefore, at this point we assure we sequenced until saturation. More specifically, the sequencing saturation the fraction of confidently mapped, valid cell-barcode, valid UMI reads that had a non-unique (cell-barcode, UMI, gene).
```{r fig.height = 10, fig.width = 14}
gg_lib_size <- qcmetric_df %>%
  dplyr::mutate(Number_of_Reads_mil = Number_of_Reads / 1000000) %>%
  ggplot2::ggplot(ggplot2::aes(x = sample_id, y = Number_of_Reads_mil, fill = sample_id)) +
  ggplot2::geom_bar(stat = "identity") +
  ggplot2::theme_bw() +
  ggplot2::scale_fill_brewer(palette = "Set2") +
  ggplot2::labs(x = "Libraries (SAMPLE IDs)", y = "Library size (in millions)") +
  ggplot2::theme(
    axis.title = ggplot2::element_text(size = 12),
    axis.text = ggplot2::element_text(size = 8),
    axis.text.x = ggplot2::element_text(hjust = 1, angle = 45),
    strip.placement = "outside",
    strip.background = ggplot2::element_rect(colour = NA),
    legend.position = "none")

gg_qc_seq_sat <- qcmetric_df %>%
  dplyr::mutate(Sequencing_Saturation_perc = Sequencing_Saturation,
         Mean_Reads_per_Spot_tho = Mean_Reads_per_Spot / 1000) %>%
  ggplot2::ggplot(ggplot2::aes(x = Mean_Reads_per_Spot_tho,
             y = Sequencing_Saturation_perc, color = sample_id)) +
  ggplot2::geom_point() +
  ggplot2::theme_bw() +
  ggplot2::scale_color_brewer(palette = "Set2") +
  ggplot2::ylim(0, 1) +
  ggrepel::geom_text_repel(ggplot2::aes(label = sample_id), size = 4) +
  ggplot2::labs(x = "Mean Reads per Spot (in thousands)", y = "Sequencing Saturation") +
  ggplot2::theme(
    axis.title = ggplot2::element_text(size = 12),
    axis.text = ggplot2::element_text(size = 10),
    legend.position = "none")

gg_qc_seq_depth_cell <- qcmetric_df %>%
  dplyr::mutate(Mean_Reads_per_Spot_tho = Mean_Reads_per_Spot / 1000) %>%
  ggplot2::ggplot(ggplot2::aes(x = Mean_Reads_per_Spot_tho,
             y = Median_Genes_per_Spot, color = sample_id)) +
  ggplot2::geom_point() +
  ggplot2::theme_bw() +
  ggplot2::scale_color_brewer(palette = "Set2") +
  ggrepel::geom_text_repel(ggplot2::aes(label = sample_id), size = 4) +
  ggplot2::labs(x = "Mean Reads per Spot (in thousands)",
                y = "Mean Detected Genes per Spot") +
  ggplot2::theme(
    axis.title = ggplot2::element_text(size = 12),
    axis.text = ggplot2::element_text(size = 10),
    legend.position = "none")

gg_qc_seq_depth <- qcmetric_df %>%
  dplyr::mutate(Number_of_Reads_mil = Number_of_Reads / 1000000) %>%
  ggplot2::ggplot(ggplot2::aes(x = Number_of_Reads_mil,
             y = Total_Genes_Detected, color = sample_id)) +
  ggplot2::geom_point() +
  ggplot2::theme_bw() +
  ggplot2::scale_color_brewer(palette = "Set2") +
  ggrepel::geom_text_repel(ggplot2::aes(label = sample_id), size = 4) +
  ggplot2::labs(x = "Number of Reads (in millions)", y = "Total Genes Detected") +
  ggplot2::theme(
    axis.title = ggplot2::element_text(size = 12),
    axis.text = ggplot2::element_text(size = 10),
    legend.position = "none")

cowplot::plot_grid(gg_lib_size, gg_qc_seq_sat, 
                   gg_qc_seq_depth_cell, gg_qc_seq_depth, 
                   nrow = 2,
                   ncol = 2,
                   align = "hv",
                   axis = "trbl")
```

# Session Info
```{r}
sessionInfo()
```