scRNA_QC.qmd

---
title: "Single-cell QC Report"
author: "Harvard Chan Bioinformatics Core"
date: "`r Sys.Date()`"
format:
  html:
    number-sections: false
    default-image-extension: svg
    lightbox: true
    callout-icon: false
    format-links: true
    toc: true
    theme: sandstone
    echo: true
    eval: true
    message: false
    warning: false
    code-copy: true
    code-overflow: wrap
    code-fold: true
    code-line-numbers: true
    embed-resources: true
    standalone: true
    html-math-method: katex
    grid:
      sidebar-width: 250px
      body-width: 900px
      margin-width: 300px
    comments:
      hypothesis: true
knitr:
  opts_chunk:
    fig-width:  7
    fig-height: 4
    fig-align: center
execute:
  freeze: auto
  keep-md: false
params:
  ribosomal: FALSE
  params_file: "https://github.com/bcbio/bcbioR-test-data/raw/refs/heads/main/singlecell/parameters.R"
  project_file: ../information.R
  seurat_obj: "https://github.com/bcbio/bcbioR-test-data/raw/refs/heads/main/singlecell/tiny.rds"
  seurat_obj_filt_fn: "seurat_post-QC.rds"
---

```{r, cache = FALSE, message = FALSE, warning=FALSE}
stopifnot(R.version$major >= 4) # requires R4
if (compareVersion(R.version$minor, "3.1") < 0) warning("We recommend >= R4.3.1")
stopifnot(compareVersion(as.character(BiocManager::version()), "3.18") >= 0)
stopifnot(compareVersion(as.character(packageVersion("Seurat")), "5.0.0") >= 0)
invisible(list2env(params, environment()))
source(params_file)
source(project_file)
```

```{r sanitize-datatable}
sanitize_datatable <- function(df, ...) {
  DT::datatable(df, ...,
    rownames = gsub("-", "_", rownames(df)),
    colnames = gsub("-", "_", colnames(df))
  )
}
```


This code is in this ![](https://img.shields.io/badge/status-stable-green) revision.

::: {.callout-note title="READ ME FIRST"}

This is a template for scRNA QC to present to your client. The actual QC can be done using our rshiny app:

https://github.com/hbc/scRNAseq_qc_app/archive/refs/heads/main.zip

Please download the app, and execute it to identify parameters interactively

After you have decided on your QC metrics load your raw object (i.e. right after you first read data into seurat) and put the parameters.R file you got from the shiny app in the same folder as this rmd.
:::

# Overview

-   Project: `r project`
-   PI: `r PI`
-   Analyst: `r analyst`
-   Experiment: `r experiment`

```{r setup}
#| include: false
library(Seurat)
library(tidyverse)
library(ggplot2)
library(R.utils)

if (isUrl(params_file)) {
  source(url(params_file))
} else {
  source(params_file)
}
```

```{r load}
if (isUrl(seurat_obj)) {
  seurat_raw <- readRDS(url(seurat_obj))
} else {
  seurat_raw <- readRDS(seurat_obj)
}
```

```{r find doublets}
#| eval: false
# if you suspect that your dataset contains doublets, you can use this code to
# detect them and filter them out. if your seurat object contains data from more
# than one sample, it is important to pass the metadata column containing the
# sample name to scDblFinder using the "samples" argument

sce <- as.SingleCellExperiment(seurat_raw, assay = "RNA")
sce <- scDblFinder(sce, samples = "orig.ident")

meta_scdblfinder <- sce@colData@listData %>%
  as.data.frame() %>%
  dplyr::select(starts_with("scDblFinder")) %>%
  mutate(barcode = sce@colData@rownames)

seurat_raw_meta <- seurat_raw@meta.data

seurat_raw_meta_new <- seurat_raw_meta %>% left_join(meta_scdblfinder, by = "barcode")
rownames(seurat_raw_meta_new) <- seurat_raw_meta_new$barcode
seurat_raw@meta.data <- seurat_raw_meta_new

seurat_raw <- subset(seurat_raw, scDblFinder.class == "singlet")
```

```{r filter no ribo}
#| eval: !expr "!ribosomal"
#| warning: false
#| results: asis

seurat_qc <- subset(
  x = seurat_raw,
  subset = (nCount_RNA >= nCount_RNA_cutoff) &
    (nFeature_RNA >= nFeature_RNA_cutoff) &
    (mitoRatio < mitoRatio_cutoff)
  ##  & (riboRatio < riboRatio_cutoff)
  & (Log10GenesPerUMI > Log10GenesPerUMI_cutoff)
)
```

```{r filter ribo}
#| eval: !expr ribosomal
#| warning: false
#| results: asis

seurat_qc <- subset(
  x = seurat_raw,
  subset = (nCount_RNA >= nCount_RNA_cutoff) &
    (nFeature_RNA >= nFeature_RNA_cutoff) &
    (mitoRatio < mitoRatio_cutoff) &
    (riboRatio < riboRatio_cutoff) &
    (Log10GenesPerUMI > Log10GenesPerUMI_cutoff)
)
```

```{r}
## Save QC object
saveRDS(seurat_qc, file = params$seurat_obj_filt_fn)
```

```{r prep-info}
metadata0 <- seurat_raw@meta.data
metadata0 <- metadata0 %>% dplyr::rename(
  nUMI = nCount_RNA,
  nGene = nFeature_RNA
)
metadata1 <- seurat_qc@meta.data
metadata1 <- metadata1 %>% dplyr::rename(
  nUMI = nCount_RNA,
  nGene = nFeature_RNA
)
```

# QC metrics: raw data {.tabset}

In this section, we review quality control (QC) metrics for the **raw feature matrices** generated by `Cellranger`. Only a low level filter excluding cells with <100 nUMIs (= number of unique molecular identifiers, or sequenced reads per cell) was applied when uploading the data into `R`. 

## Cells per sample

```{r cells raw}
metadata0 %>%
  group_by(orig.ident) %>%
  summarize(n_cells_pre_filt = n()) %>%
  sanitize_datatable()
```

## UMIs per cell

Here, we look at the distribution of UMIs (unique molecular identifiers, or sequenced reads) per cell (droplet) in the dataset. Before QC, we expect a biomodal distribution with a first *small* peak at low numbers of UMIs (<250) corresponding to droplets that encapsulated background/dying cells, and a second higher peak centered at >1000. The line is at `r nCount_RNA_cutoff`.


```{r raw_nUMIs}
#| fig-width: 7
metadata0 %>%
  ggplot(aes(x = nUMI, color = orig.ident, fill = orig.ident)) +
  geom_density(alpha = 0.2) +
  theme_classic() +
  ylab("Cell density") +
  scale_x_log10() +
  geom_vline(xintercept = nCount_RNA_cutoff) +
  facet_wrap(. ~ orig.ident) +
  ggtitle("UMIs per cell in raw dataset")
```

```{r}
#| fig-width: 7
# Visualize the distribution of nUMIs per cell (boxplot)
metadata0 %>%
  ggplot(aes(x = orig.ident, y = log10(nUMI), fill = orig.ident)) +
  geom_violin() +
  geom_boxplot(width = 0.1, fill = alpha("white", 0.7)) +
  theme_classic() +
  geom_hline(yintercept = log10(nCount_RNA_cutoff)) +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))
```

## Genes per cell

Here, we look at the number of different genes that were detected in each cell. By "detected", we mean genes with a non-zero read count measurement. Gene detection in the range of 500 to 5000 is normal for most single-cell experiments. The line is at `r nFeature_RNA_cutoff`. 

```{r raw_nGene}
#| fig-width: 7
# Visualize the distribution of genes detected per cell (histogram)
metadata0 %>%
  ggplot(aes(x = nGene, color = orig.ident, fill = orig.ident)) +
  geom_density(alpha = 0.2) +
  theme_classic() +
  scale_x_log10() +
  geom_vline(xintercept = c(nFeature_RNA_cutoff)) +
  facet_wrap(. ~ orig.ident) +
  ggtitle("Detected genes per cell in raw dataset")
```

```{r}
#| fig-width: 7
# Visualize the distribution of nUMIs per cell (boxplot)
metadata0 %>%
  ggplot(aes(x = orig.ident, y = log10(nGene), fill = orig.ident)) +
  geom_violin() +
  geom_boxplot(width = 0.1, fill = alpha("white", 0.7)) +
  theme_classic() +
  geom_hline(yintercept = c(log10(nFeature_RNA_cutoff))) +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))
```

## Mitochondrial ratio

We evaluate overall mitochondrial gene expression as a biomarker of cellular stress during sample preparation. Typically, we expect mitochondrial genes to account for <20% of overall transcripts in each cell. The line indicates `r mitoRatio_cutoff*100` %.

```{r raw_mito}
#| fig-width: 7
#| warning: false
# Visualize the distribution of mitochondrial gene expression detected per cell
metadata0 %>%
  ggplot(aes(color = orig.ident, x = mitoRatio, fill = orig.ident)) +
  geom_density(alpha = 0.2) +
  scale_x_log10() +
  theme_classic() +
  geom_vline(xintercept = mitoRatio_cutoff) +
  facet_wrap(. ~ orig.ident) +
  ggtitle("Percentage of mitochondrial gene expression per cell in raw dataset")
```

```{r}
#| fig-width: 7
# Visualize the distribution of mitochondrial gene expression per cell (violin plot)
metadata0 %>%
  ggplot(aes(x = orig.ident, y = mitoRatio, fill = orig.ident)) +
  geom_violin() +
  geom_boxplot(width = 0.1, fill = alpha("white", 0.7)) +
  theme_classic() +
  geom_hline(yintercept = c(mitoRatio_cutoff)) +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))
```

```{r raw_ribo}
#| eval: !expr ribosomal
#| warning: false
#| results: asis

cat("## Ribosomal ratio \n")

cat("We evaluate overall ribosomal gene expression. Different cells types are expected to have different levels of ribosomal expression. Due to this, there is no suggested cutoff for ribosomal ratio. We merely expect it to be similar among samples with similar cellular composition. Note that extremely high levels can indicate low quality reads. \n")
# Visualize the distribution of mitochondrial gene expression detected per cell
metadata0 %>%
  ggplot(aes(color = orig.ident, x = riboRatio, fill = orig.ident)) +
  geom_density(alpha = 0.2) +
  scale_x_log10() +
  theme_classic() +
  facet_wrap(. ~ orig.ident) +
  ggtitle("Percentage of ribosomal gene expression per cell in raw dataset")
```

## UMIs vs. Genes

By plotting the number of UMIs per cell (x-axis) vs. the number of genes per cell (y-axis), we can visually assess whether there is a large proportion of low quality cells with low read counts and/or gene detection (bottom left quadrant of the plot). In the following representation, cells are further color-coded based on the percentage of mitochondrial genes found among total detected genes. The line for nUMI is at `r nCount_RNA_cutoff` and the line for nGene is at `r nFeature_RNA_cutoff`.

```{r raw_gene_by_umi}
#| fig-height: 6
#| fig-width: 8
#| warning: false
# Visualize the correlation between genes detected and number of UMIs and determine whether strong presence of cells with low numbers of genes/UMIs
metadata0 %>%
  ggplot(aes(x = nUMI, y = nGene, color = mitoRatio)) +
  geom_point() +
  stat_smooth(method = lm) +
  scale_x_log10() +
  scale_y_log10() +
  theme_classic() +
  geom_vline(xintercept = nCount_RNA_cutoff) +
  geom_hline(yintercept = nFeature_RNA_cutoff) +
  ggtitle("Genes vs. nUMIs in raw dataset") +
  facet_wrap(~orig.ident)
```

## Complexity

Another way to assess the quality and purity of a single-cell dataset is to look for cells that have fewer detected genes per UMI than others. Typical values for this metric are >0.8 for most cells. Cells with lower diversity in the genes they express may be low-complexity cell types such as red blood cells. With sorted populations, we expect high purity and a very similar complexity distribution across samples. The line is at `r Log10GenesPerUMI_cutoff`. 

```{r raw_novelty}
# Visualize the overall novelty of the gene expression by visualizing the genes detected per UMI
metadata0 %>%
  ggplot(aes(x = Log10GenesPerUMI, color = orig.ident, fill = orig.ident)) +
  geom_density(alpha = 0.2) +
  theme_classic() +
  geom_vline(xintercept = Log10GenesPerUMI_cutoff) +
  facet_wrap(. ~ orig.ident) +
  ggtitle("log10(Genes per UMI) in raw dataset")
```


```{r}
# Visualize the overall novelty of the gene expression by visualizing the genes detected per UMI (boxplot)
metadata0 %>%
  ggplot(aes(x = orig.ident, Log10GenesPerUMI, fill = orig.ident)) +
  geom_violin() +
  geom_boxplot(width = 0.1, fill = alpha("white", 0.7)) +
  theme_classic() +
  geom_hline(yintercept = Log10GenesPerUMI_cutoff) +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))
```


# QC metrics: Filtered data {.tabset}

Based on the above QC metrics, we filtered the dataset to isolate cells passing the following thresholds: >`r nCount_RNA_cutoff` UMIs, >`r nFeature_RNA_cutoff` genes, <`r mitoRatio_cutoff` mitochondrial gene ratio, and >`r Log10GenesPerUMI_cutoff` complexity.

In this section, we review QC metrics for our filtered dataset.

## Cells per sample

```{r cells filtered}
metadata0 %>%
  group_by(orig.ident) %>%
  summarize(n_cells_pre_filt = n()) %>%
  left_join(metadata1 %>% group_by(orig.ident) %>% summarize(n_cells_post_filt = n())) %>%
  mutate(pct_remaining = round(n_cells_post_filt / n_cells_pre_filt * 100, 3)) %>%
  sanitize_datatable()
```


## UMIs per cell

```{r qc1_nUMIs}
metadata1 %>%
  ggplot(aes(color = orig.ident, x = nUMI, fill = orig.ident)) +
  geom_density(alpha = 0.2) +
  scale_x_log10() +
  theme_classic() +
  ylab("Cell density") +
  xlab("nUMI") +
  facet_wrap(. ~ orig.ident) +
  ggtitle("UMIs per cell in filtered dataset")
```

```{r}
# Visualize the distribution of nUMIs per cell (boxplot)
metadata1 %>%
  ggplot(aes(x = orig.ident, y = log10(nUMI), fill = orig.ident)) +
  geom_violin() +
  geom_boxplot(width = 0.1, fill = alpha("white", 0.7)) +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))
```

## Genes per cell

```{r qc1_genes}
# Visualize the distribution of genes detected per cell via histogram
metadata1 %>%
  ggplot(aes(color = orig.ident, x = nGene, fill = orig.ident)) +
  geom_density(alpha = 0.2) +
  theme_classic() +
  scale_x_log10() +
  xlab("nGene") +
  facet_wrap(. ~ orig.ident) +
  ggtitle("Detected genes per cell in filtered dataset")
```

```{r}
# Visualize the distribution of genes detected per cell (boxplot)
metadata1 %>%
  ggplot(aes(x = orig.ident, y = log10(nGene), fill = orig.ident)) +
  geom_violin() +
  geom_boxplot(width = 0.1, fill = alpha("white", 0.7)) +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))
```

## Mitochondrial ratio

```{r qc1_mitoratio}
#| message: false
#| warning: false
# Visualize the distribution of mitochondrial gene expression detected per cell
metadata1 %>%
  ggplot(aes(color = orig.ident, x = mitoRatio, fill = orig.ident)) +
  geom_density(alpha = 0.2) +
  scale_x_log10() +
  theme_classic() +
  geom_vline(xintercept = 0.1)
# +
#   facet_wrap(. ~ surgery)
```

```{r}
# Visualize the distribution of mitochondrial gene expression per cell (violin plot)
metadata1 %>%
  ggplot(aes(x = orig.ident, y = mitoRatio, fill = orig.ident)) +
  geom_violin() +
  geom_boxplot(width = 0.1, fill = alpha("white", 0.7)) +
  theme_classic() +
  geom_hline(yintercept = c(mitoRatio_cutoff)) +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))
```

```{r qc1_ribo}
#| eval: !expr ribosomal
#| warning: false
#| results: asis

cat("## Ribosomal ratio \n")

# Visualize the distribution of ribosomal gene expression detected per cell
metadata1 %>%
  ggplot(aes(color = orig.ident, x = riboRatio, fill = orig.ident)) +
  geom_density(alpha = 0.2) +
  scale_x_log10() +
  theme_classic() +
  facet_wrap(. ~ orig.ident) +
  ggtitle("Percentage of ribosomal gene expression per cell in filtered dataset")
```

```{r}
#| eval: !expr ribosomal
metadata1 %>%
  ggplot(aes(x = orig.ident, y = riboRatio, fill = orig.ident)) +
  geom_violin() +
  geom_boxplot(width = 0.1, fill = alpha("white", 0.7)) +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))
```

## UMIs vs. Genes

```{r qc1_genes_per_UMI}
#| fig-height: 6
#| fig-width: 8
#| warning: false

metadata1 %>%
  ggplot(aes(x = nUMI, y = nGene, color = mitoRatio)) +
  geom_point() +
  stat_smooth(method = lm) +
  scale_x_log10() +
  scale_y_log10() +
  theme_classic() +
  ggtitle("Genes vs. nUMIs in filtered dataset") +
  xlab("nUMI") +
  ylab("nGene") +
  facet_wrap(~orig.ident)
```

## Complexity

```{r qc1_complexity}
#| fig-width: 7
# Visualize the overall novelty of the gene expression by visualizing the genes detected per UMI
metadata1 %>%
  ggplot(aes(x = Log10GenesPerUMI, color = orig.ident, fill = orig.ident)) +
  geom_density(alpha = 0.2) +
  theme_classic() +
  facet_wrap(. ~ orig.ident) +
  ggtitle("log10(Genes per UMI) in filtered dataset")
```

```{r}
#| fig-width: 7
# Visualize the distribution of nUMIs per cell (boxplot)
metadata1 %>%
  ggplot(aes(x = orig.ident, Log10GenesPerUMI, fill = orig.ident)) +
  geom_violin() +
  geom_boxplot(width = 0.1, fill = alpha("white", 0.7)) +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))
```


# R session

```{r}
sessionInfo()
```