perform_qc.R

# This script also ensures that cell names correspond to the names generated by
# cellranger agg. The latter program ensures unique barcodes across samples by
# numbering them like BARCODE-1 to BARCODE-10, using the order of samples in
# aggregation.csv. By contrast, Seurat (if data is loaded with Read10X_h5())
# uses the numbering system BARCODE-1_1 to BARCODE-1_10; sample indices
# correspond to the order of samples in c(x, y) of merge().
#
# Generates the file rna_qcpassed.rds, which contains a list of Seurat objects
# representing the QC-filtered datasets.
#
# @DEPI data_raw/rna_seq
# @DEPO rna_qcpassed.rds

library(Seurat)
library(scuttle)
library(scds)
library(tidyverse)
library(fs)
library(ggbeeswarm)
library(grid)
library(patchwork)
source("common_functions.R")



# Parameters --------------------------------------------------------------

# folder that contains cellranger results
data_dir <- "data_raw/rna_seq"

# path where the merged dataset is saved
out_file <- "data_generated/rna_qcpassed.rds"

# data frame with columns 'bsf_order', 'sample_file', 'sample_name', and 'group'
samples <-
  read_csv("metadata/sample_groups.csv", comment = "#") %>% 
  arrange(bsf_order) %>%
  select(bsf_order, sample_file = geo_rna_filtered, sample_name = sample, group)



# Load datasets -----------------------------------------------------------

#' Perform quality control checks on scRNA-seq data.
#' 
#' Method "cutoff" marks those cells for discarding whose percentage of
#' mitochondrial genes and bcds score exceed a fixed threshold, and whose number
#' of features does not fall within a fixed range. Method "outlier" uses the
#' outlier-based method implemented in `scuttle`.
#'
#' @param data A Seurat object.
#' @param method QC filtering method `"cutoff"` (default) or `"outlier"`.
#' @param max_percent_mito Maximum percentage of mitochondrial genes.
#' @param min_n_feature Minimum number of features.
#' @param max_n_feature Maximum number of features.
#' @param max_bcds_score Maximum bcds score.
#'
#' @return The Seurat object with an additional boolean metadata column
#'   `pass_qc`.
perform_qc <- function(data,
                       method = c("cutoff", "outlier"),
                       max_percent_mito = 10,
                       min_n_feature = 200,
                       max_n_feature = 5000,
                       max_bcds_score = 0.8) {
  method <- match.arg(method) 
  
  data@meta.data$bcds_score <- 
    data %>% 
    as.SingleCellExperiment() %>%
    bcds() %>% 
    colData() %>% 
    magrittr::use_series("bcds_score")
  
  if (method == "outlier") {
    qc_results <- quickPerCellQC(
      data@meta.data,
      sum.field = "library_size",
      detected.field = "n_features",
      sub.fields = c("percent_mito", "bcds_score")
    )
    pass_qc <- !qc_results$discard
  } else {
    pass_qc <- 
      data@meta.data %>% 
      mutate(
        pass_qc =
          n_features > min_n_feature &
          n_features < max_n_feature &
          percent_mito < max_percent_mito &
          bcds_score < max_bcds_score
      ) %>% 
      pull(pass_qc)
  }
  
  data %>% 
    AddMetaData(pass_qc, "pass_qc")
}


datasets <-
  pmap(
    samples,
    function(bsf_order, sample_file, sample_name, group) {
      info("Loading ", sample_file)
      data <- 
        Read10X_h5(path_join(c(data_dir, sample_file))) %>% 
        magrittr::set_colnames(
          str_replace(colnames(.), "\\d+$", as.character(bsf_order))
        ) %>% 
        CreateSeuratObject(sample_file) %>% 
        PercentageFeatureSet(pattern = "^MT-", col.name = "percent_mito")
      
      data@meta.data <-
        data@meta.data %>%
        select(!orig.ident) %>%
        dplyr::rename(library_size = nCount_RNA, n_features = nFeature_RNA) %>%
        mutate(sample = {{sample_name}}, group = {{group}})
      
      info("Calculating QC measures")
      perform_qc(data)
    }
  )



# Plot QC metrics ---------------------------------------------------------

# use black for cells passing QC and pink for discarded cells
default_color_scheme <- scale_color_manual(
  name = "QC passed",
  values = c("#dd1c77", "black"),
  guide = guide_legend(override.aes = list(size = 2, alpha = 1))
)

# tune facet labels
default_theme <- 
  theme_bw() +
  theme(
    panel.grid = element_blank(),
    strip.background = element_blank(),
    strip.text = element_text(face = "bold")
  )


#' Generate QC plots.
#'
#' The patchworked plot contains five plots: Violin plots of (1) mitochondrial
#' gene percentage, (2) library size, (3) number of features. and (4) bcds
#' score; as well as scatter plots of (5) mt genes vs library size and (6) mt
#' genes vs feature count.
#'
#' @param data A Seurat object after QC analysis, but before filtering.
#' @param filename Name of output file. If `NULL`, do not save the plot. If
#'   `"generate"`, automatically generate the file name from the names of the
#'   sample and the raw file.
#'
#' @return A ggplot object.
plot_qc_details <- function(data, filename = NULL) {
  vis_data <- 
    data@meta.data %>% 
    arrange(desc(pass_qc))
  
  plot_violin <- function(column, name, ylim_upper = NA) {
    ggplot(vis_data, aes(name, {{column}})) +
      geom_violin(color = "gray60", alpha = 0.2, width = 0.8) +
      geom_quasirandom(
        aes(color = pass_qc),
        width = 0.4,
        bandwidth = 1,
        alpha = 0.25,
        size = 0.1,
        show.legend = FALSE
      ) +
      ylab(NULL) +
      ylim(0, ylim_upper) +
      default_color_scheme +
      default_theme +
      theme(axis.title.x = element_blank())
  }
  
  plot_scatter <- function(x, y, x_title, y_title) {
    ggplot(vis_data, aes({{x}}, {{y}})) +
      geom_point(
        aes(color = pass_qc),
        alpha = 0.25,
        size = 0.5,
        show.legend = FALSE
      ) +
      ylim(0, 100) +
      xlab(x_title) +
      ylab(y_title) +
      default_color_scheme +
      default_theme
  }
  
  layout <- "
    AAA
    BCF
    DEG"

  p <-
    wrap_plots(
      wrap_elements(
        panel = textGrob(
          str_glue(
            "sample: {vis_data$sample[1]}\n",
            "file: {data@project.name}\n",
            "QC: passed by {sum(vis_data$pass_qc)} of {nrow(vis_data)} cells"
          ),
          x = 0,
          just = "left",
          gp = gpar(fontsize = 8)
        ),
        clip = FALSE
      ),
      plot_violin(library_size, "library size"),
      plot_violin(n_features, "number of features"),
      plot_violin(percent_mito, "% mito genes", ylim_upper = 100),
      plot_violin(bcds_score, "bcds score", ylim_upper = 1),
      plot_scatter(library_size, percent_mito, "library size", "% mito genes"),
      plot_scatter(n_features, percent_mito, "number of features", "% mito genes")
    ) +
    plot_layout(design = layout, widths = c(2, 2, 4), heights = c(1, 7, 7))
  
  if (filename == "generate") {
    file <- str_match(data@project.name , "T/(.*)_trans")[, 2]
    filename <- str_glue("qc/details/{vis_data$sample[1]}__{file}")
  }
    
  ggsave_default(filename, width = 200, height = 200)
  p
}

walk(
  datasets,
  plot_qc_details,
  filename = "generate"
)



# plot only the violin plots, facet by group
datasets %>%
  map_dfr(~.@meta.data) %>%
  arrange(desc(pass_qc)) %>% 
  pivot_longer(c(library_size, n_features, percent_mito, bcds_score)) %>% 
  mutate(
    name = recode(name,
                  library_size = "library size",
                  n_features = "number of features",
                  percent_mito = "% mito genes")
  ) %>% 
  ggplot(aes(sample, value)) +
  geom_violin(color = "gray60", alpha = 0.2, width = 0.8, scale = "width") +
  geom_quasirandom(
    aes(color = pass_qc),
    width = 0.4,
    bandwidth = 1,
    alpha = 0.25,
    size = 0.1
  ) +
  facet_grid(
    vars(name),
    vars(group),
    scales = "free",
    space = "free_x",
    switch = "y"
  ) +
  ylab(NULL) +
  default_color_scheme +
  default_theme +
  theme(strip.placement = "outside")
ggsave_default("qc/qc_overview", width = 420, height = 297)



# Discard cells and save datasets -----------------------------------------

datasets %>%
  map(
    function(dataset) {
      dataset <- subset(dataset, subset = pass_qc)
      dataset@meta.data <- 
        dataset@meta.data %>% 
        select(!c(pass_qc, nCount_RNA, nFeature_RNA))
      dataset
    }
  ) %>%
  saveRDS(out_file)