---
title: "scRNA Analysis - 02 - Doublet Detection"
author: "Sydney Informatics Hub"
format: html
---

This notebook is designed to detect doublets in your single cell RNA sequencing data. Doublets occur when two cells are captured in the same droplet and are given the same barcode. The RNA of a doublet is sequenced as if it came from a single cell, so the resulting transcriptome is an artificial mixture of both cells’ profiles. This can create misleading “hybrid” cell types that are challenging to annotate. Multiplets are the broader category where more than one cell (two, three, or more) ends up in the same droplet. It is important to remove putative multiplets to ensure that cells are annotated correctly.

This notebook uses `DoubletFinder` to identify and annotate likely doublets in your Seurat data objects and gives you the option at the end to remove them or keep them. Typically, you would want to remove doublets, as they will confound your results, so by default this notebook will remove them. The notebook assumes you have run the QC notebook (`01.qc.qmd`) first and that your pre-filtered sample data is saved in a directory called `outputs` in the current working directory. It requires a samplesheet `data/doublet_samples.csv` that defines the sample names and (optionally) the expected doublet rates for each sample. The notebook will then load the RDS files saved by the QC notebook and then proceed to run the R package `DoubletFinder` to detect doublets.

::: {.callout-note title="Overview"}

The key steps in this notebook are:

1. Load in the QC'd data for each sample from the previous notebook
2. Estimate the expected doublet rate
3. Run DoubletFinder to identify doublets
4. Remove doublets from the dataset (optional, but **recommended**)

:::

## Imports

While working through this notebook manually, you will need to run the following block of code to import all the necessary R libraries and helper functions:

```{r setup, include = FALSE}
# Imports
knitr::opts_chunk$set(echo = TRUE)
library(here)
library(tidyverse)
library(Seurat)
library(DoubletFinder)
library(DT)
source(here("R/shinyapps.R"))

# Run a version check on DoubletFinder - must be at least 2.0.6
stopifnot(packageVersion("DoubletFinder") >= "2.0.6")

# Helper functions
find_min_pc <- function(stdvs) {
  # Find significant PCs
  # From https://biostatsquid.com/doubletfinder-tutorial/.
  # Use this function to determine the number of PCs to include in downstream
  # analyses (i.e. other reductions (UMAP, tSNE), clustering). Requires a
  # `SeuratObject` that has a PCA reduction from `Seurat::RunPCA()`.
  percent_stdv <- (stdvs/sum(stdvs)) * 100
  cumulative <- cumsum(percent_stdv)
  co1 <- which(cumulative > 90 & percent_stdv < 5)[1]
  co2 <- sort(which((percent_stdv[1:length(percent_stdv) - 1] -
                       percent_stdv[2:length(percent_stdv)]) > 0.1),
              decreasing = T)[1] + 1
  min_pc <- min(co1, co2)
  print(min_pc)
}

SCTransform_reduce_dims <- function(so) {
  # Normalise with SCTransform and reduce dimensions with PCA.
  sct <-
    Seurat::SCTransform(so, vars.to.regress = c("percent.mt"), verbose = F) |>
    Seurat::RunPCA()

  # Calculate the minimum number of PCs that explain the majority of the variation
  min_pc <- find_min_pc(sct@reductions$pca@stdev)

  sct <-
    Seurat::RunUMAP(sct, dims = 1:min_pc, verbose = F) |>
    Seurat::FindNeighbors(dims = 1:min_pc, verbose = F)

  return(sct)
}

read_tmp_rds_files <- TRUE
```

## Define inputs

At the start of the QC notebook, we generated a template CSV sample sheet called `inputs/doublet_samples.csv`. This file is used to define your sample names and the expected multiplet rate for each sample. The column names are `sample` for the sample name column and `multiplet_rate` for the multiplet rates. Multiplet rates should be the proportion of each sample expected to be multiplets; for example, if 5% of the cells were expected to be multiplets, the value should be `0.05`. This column can be left blank if you wish to use the automatic multiplet detection method (described below). The sample names will be used to find the filtered RDS files saved by the QC notebook, which will be named like `outputs/<SAMPLE_NAME>.filtered_clustered.Rds`. An example of the samplesheet is given below; the sample `normal_sample` has been specified with a multiplet rate of `0.076`, while the other sample `tumour_sample` has had this column left blank, and so the automatic multiplet rate estimation method will be used:

```
sample,multiplet_rate
normal_sample,0.076
tumour_sample,
```

---

**❱❱❱ OPTIONAL ACTION ❰❰❰**

1. Open the file `inputs/doublet_samples.csv` in a text editor.
2. Review the contents of this file and update **if necessary**.
    1. Ensure that the first row remains `sample,multiplet_rate`. **Do not alter this line**.
    2. In the remaining rows, first column remains unchanged as the sample name (no spaces).
    3. If you have a specific expected multiplet rate for one or more samples, update its second column accordingly. Otherwise, leave it blank to automatically estimate the rate.
3. Save the file and close.

---

Cluster resolution for each sample should have been determined while running the QC notebook, and should already be defined in the file `inputs/cluster_resolutions.csv`. See the QC notebook for further details on this file.

### A note on automated multiplet rates estimation

The doublet detection workflow requires an estimate of the expected doublet rate. This can be either specified manually in the samplesheet as described above, or automatically estimated.

The automatic multiplet rate estimation is based on the expected multiplet rate of 0.8% per 1000 cells (or `8e-6` per cell), [as published by 10X for their Chromium NextGEM v3.1 chemistry](https://cdn.10xgenomics.com/image/upload/v1722285481/support-documents/CG000315_ChromiumNextGEMSingleCell3__GeneExpression_v3.1_DualIndex__RevF.pdf). The automatic estimation method will be used if no mutliplet rate value was specified for a sample in the samplesheet. Note that the multiplet rates published by 10X only account for up to 10,000 cells, and higher cell counts may result in higher rates of multiplets. Also note that the 10X data may not apply to all cell types.

```{r inputs}
samples <- read_csv(here("inputs/doublet_samples.csv"), show_col_types = FALSE)
cluster_resolutions <- read_csv(here("inputs/cluster_resolutions.csv"), show_col_types = FALSE)

# Create a multiplet rate column if not present
if (! "multiplet_rate" %in% colnames(samples)) {
  samples <- samples %>%
    mutate(
      multiplet_rate = NA
    )
}

# Merge with the cluster resolutions file
samples <- samples %>%
  left_join(cluster_resolutions, by = c("sample"))

# Remove any additional columns if present
samples <- samples %>%
  dplyr::select(sample, res, multiplet_rate)

# Load data
sample_names <- samples$sample
all_sct_files <- paste0(sample_names, ".filtered_clustered.Rds")
all_sct <- lapply(all_sct_files, function(f) {
  readRDS(here("outputs", f))
})
names(all_sct) <- sample_names

# Get the number of cells and multiplet rates
# If multiplet rate is unavailable, estimate at 8e-6 per cell
samples <- samples %>%
  mutate(
    n_cells = pmap_int(list(sample), function(s) { nrow(all_sct[[s]]@meta.data) }),
    multiplet_rate = case_when(
      !is.na(multiplet_rate) ~ multiplet_rate,
      .default = 8e-6 * n_cells
    )
  )

samples %>%
  mutate(across(where(is.numeric), ~ round(.x, 4))) %>%
  datatable(colnames = c("Sample", "Cluster Resolution", "Multiplet Rate", "# Cells"))
```

## Define DoubletFinder function

The following function wraps up the DoubletFinder workflow into a convenient function:

```{r doubletfinder_custom, eval = FALSE}
run_doubletfinder_custom <- function(seurat_obj, cluster_res, multiplet_rate) {
  # From https://biostatsquid.com/doubletfinder-tutorial/
  full_res <- paste0("SCT_snn_res.", cluster_res)
  min_pc <- find_min_pc(seurat_obj@reductions$pca@stdev)

  # pK identification (no ground-truth) 
  #introduces artificial doublets in varying props, merges with real data set and 
  # preprocesses the data + calculates the prop of artficial neighrest neighbours, 
  # provides a list of the proportion of artificial nearest neighbours for varying
  # combinations of the pN and pK
  sweep_list <- DoubletFinder::paramSweep(seurat_obj, PCs = 1:min_pc, sct = T)
  sweep_stats <- DoubletFinder::summarizeSweep(sweep_list)
  # computes a metric to find the optimal pK value (max mean variance normalised by modality coefficient)
  bcmvn <- DoubletFinder::find.pK(sweep_stats)

  # Optimal pK is the max of the bimodality coefficient (BCmvn) distribution
  optimal_pk <- bcmvn %>% 
    dplyr::filter(BCmetric == max(BCmetric)) %>%
    dplyr::select(pK)
  optimal_pk <- as.numeric(as.character(optimal_pk[[1]]))

  # Homotypic doublet % estimate
  annotations <- seurat_obj@meta.data[[full_res]] # use the clusters as the user-defined cell types
  homotypic_prop <- DoubletFinder::modelHomotypic(annotations) # get proportions of homotypic doublets
  nExp_poi <- round(multiplet_rate * nrow(seurat_obj@meta.data)) # multiply by number of cells to get the number of expected multiplets
  nExp_poi_adj <- round(nExp_poi * (1 - homotypic_prop)) # expected number of doublets

  # run DoubletFinder
  seurat_doublets <- DoubletFinder::doubletFinder(
    seu = seurat_obj,
    PCs = 1:min_pc, 
    pN = 0.25, # default
    pK = optimal_pk, # the neighborhood size used to compute the number of artificial nearest neighbours
    nExp = nExp_poi_adj, # number of expected real doublets
    reuse.pANN = NULL,
    sct = T
  )

  # change name of metadata column with Singlet/Doublet information
  colnames(seurat_doublets@meta.data)[grepl('DF.classifications.*', colnames(seurat_doublets@meta.data))] <- "doublet_finder"

  return(seurat_doublets)
}
```

## Run doublet detection

Run the next block of code to detect doublets in your data and annotate the cells in each Seurat object with the results.

```{r run_doubletfinder, eval = FALSE}
options(future.globals.maxSize = 2600*1024^2)

all_doublets <- lapply(all_sct, function(s) {
  gc()  # Clean up memory before each run
  s_name <- as.character(s@meta.data$orig.ident[[1]])
  cluster_res <- as.numeric(samples$res[match(s_name, samples$sample)])
  mr <- as.numeric(samples$multiplet_rate[match(s_name, samples$sample)])
  run_doubletfinder_custom(s, cluster_res, mr)
})
gc()  # Clean up memory after final run

dir.create(here("tmp_outputs", "02.doublet_detection"), recursive = TRUE)
saveRDS(all_doublets, here("tmp_outputs", "02.doublet_detection", "all_doublets.Rds"))
read_tmp_rds_files <- FALSE
```

Note: when running the above code interactively, `DoubletFinder` will generate a plot for each sample. This plot graphs the `BCMVN` value - the "bimodality coefficient" - against the `pK` value - the "neighbourhood size" parameter. `DoubletFinder` uses this data internally to pick the optimal `pK` value that separates singlets and doublets effectively. See the [DoubletFinder paper](https://www.sciencedirect.com/science/article/pii/S2405471219300730) for further reading.

## Summarise results

We can print out a summary table of the doublets detected in each sample:

```{r summarise_doublets}
if(read_tmp_rds_files) {
  all_doublets <- readRDS(here("tmp_outputs", "02.doublet_detection", "all_doublets.Rds"))
}

doublet_summaries <- lapply(all_doublets, function(s) { table(s$doublet_finder) })
doublets_summary <- tibble(
  sample = names(doublet_summaries),
  n_singlets = sapply(doublet_summaries, function(s) { s[["Singlet"]] }),
  n_doublets = sapply(doublet_summaries, function(s) { s[["Doublet"]] })
)
datatable(doublets_summary, colnames = c("Sample", "# Singlets", "# Doublets"))
```

## Plot data

### UMAP

We also display the UMAP plots with cells coloured by their singlet/doublet status:

```{r plot_doublets}
doublet_plots <- lapply(all_doublets, function(s) {
  s_name <- as.character(s@meta.data$orig.ident[[1]])
  DimPlot(s, reduction = "umap", group.by = "doublet_finder") +
    ggtitle(s_name)
})

for (p in doublet_plots) {
  print(p)
}
```

### Proportion of doublets per cluster

You may notice in the above plots that some clusters are over-populated with doublets. We can further inspect this by plotting the doublet rate per cluster.

```{r plot_doublets_per_cluster, message = FALSE, warning = FALSE}
doublet_cluster_plots <- lapply(all_doublets, function(s) {
  s_name <- as.character(s@meta.data$orig.ident[[1]])
  cluster_res <- as.numeric(samples$res[match(s_name, samples$sample)])
  cluster_res <- paste0("SCT_snn_res.", cluster_res)

  s@meta.data %>%
    dplyr::select(all_of(cluster_res), doublet_finder) %>%
    group_by(across(all_of(cluster_res)), doublet_finder) %>%
    summarise(n = n()) %>%
    ggplot(aes(x = .data[[cluster_res]], y = n, fill = doublet_finder)) +
      geom_col() +
      theme_light() +
      ggtitle(s_name)
})

for (p in doublet_cluster_plots) {
  print(p)
}
```

We can also plot the proportion of doublets per cluster. In the following plot, smaller, dark blue bars represent low doublet rates for a cluster while larger, bright yellow bars represent higher doublet rates; missing bars indicate that no doublets were detected for that cluster.

```{r plot_doublet_props_per_cluster}
doublet_prop_cluster_plots <- lapply(all_doublets, function(s) {
  s_name <- as.character(s@meta.data$orig.ident[[1]])
  cluster_res <- as.numeric(samples$res[match(s_name, samples$sample)])
  cluster_res <- paste0("SCT_snn_res.", cluster_res)

  s@meta.data %>%
    dplyr::select(all_of(cluster_res), doublet_finder) %>%
    group_by(across(all_of(cluster_res)), doublet_finder) %>%
    summarise(n = n()) %>%
    pivot_wider(names_from = doublet_finder, values_from = n, values_fill = 0) %>%
    summarise(prop_doublets = Doublet / (Singlet + Doublet)) %>%
    ggplot(aes(x = .data[[cluster_res]], y = prop_doublets, fill = prop_doublets)) +
        geom_col() +
        theme_light() +
        ylim(0, 1) +
        scale_fill_continuous(type = "viridis", limits = c(0, 1)) +
        ylab("Proportion of doublets") +
        ggtitle(s_name)
})

for (p in doublet_prop_cluster_plots) {
  print(p)
}
```

If you notice some clusters that are particularly over-populated with doublets, it is very possible that this is *why* these cells were clustering together. In this case, we recommend you remove all doublets in the next step.

## Remove doublets

At this point, you can choose to either remove the doublets you have detected (recommended and the default setting for this notebook) or leave them in your dataset and simply annotate them in your metadata. If you choose to leave these barcodes in your data, you should keep this in mind when performing downstream analyses and either remove them for those analyses, or potentially take them into account as co-variates in your models.

As our recommendation is to remove these barcodes, **the downstream notebooks will make the assumption that doublets have been removed**.

---

**❱❱❱ OPTIONAL ACTION ❰❰❰**

If you want to keep the doublets in your dataset (**not recommended**), change the following code to `remove_doublets <- FALSE`:

---

```{r remove_doublets_yn}
# Set this to FALSE if you wish to leave doublets in your dataset
remove_doublets <- TRUE
```

If you are removing doublets, you can run the following code chunk to save the doublet-annotated data before removal.

```{r save_doublets, eval = FALSE}
if(remove_doublets) {
  dir.create(here("outputs"))
  for (s in all_doublets) {
    sample_name <- as.character(s@meta.data$orig.ident[[1]])
    base_name <- paste0(sample_name, ".filtered_clustered.doublets_detected_kept.Rds")
    SaveSeuratRds(s, here("outputs", base_name))
  }
}
```

The following code chunk will remove the doublets if you have chosen to do so, and it will additionally re-run the SCTransform normalisation and clustering, as these will now be invalidated by the removal of cells from the dataset.

```{r remove_doublets, eval = FALSE}
res <- scan(here("inputs/all_cluster_resolutions.txt"), numeric())
cluster_algorithm <- scan(here("inputs/cluster_algorithm.txt"))

# Set the random seed for the cluster algorithm
random_seed <- ifelse(cluster_algorithm == 1, 0, 1)

if(remove_doublets) {
  all_doublets <- lapply(all_doublets, function(s) {
    s <- subset(s, subset = (doublet_finder == "Singlet"))
    # Re-run SCTransform and clustering
    s <- SCTransform_reduce_dims(s)
    s <- FindClusters(s, resolution = res, algorithm = cluster_algorithm, random.seed = random_seed, verbose = 0)
    gc()  # Clean up memory after each run
    s
  })
  cat("Doublets have been removed.")
} else {
  warning("Doublets have not been removed. Keep this in mind when conducting downstream analyses.")
}
```

## Save data

Finally, save the Seurat objects as new `.Rds` files.

```{r save_data, eval = FALSE}
dir.create(here("outputs"))
for (s in all_doublets) {
  sample_name <- as.character(s@meta.data$orig.ident[[1]])
  base_name <- paste0(sample_name, ".filtered_clustered.doublets_detected.Rds")
  SaveSeuratRds(s, here("outputs", base_name))
}
```

## Package version summary

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