---
title: "scRNA Analysis - 01 - QC"
author: "Sydney Informatics Hub"
format: html
---

The goal of this notebook is for you to be able to plot and inspect some basic quality control metrics for your single cell sequencing data. We start by loading your samples into R with the Seurat package and inspecting the general distribution of counts. We then identify read and gene count thresholds to use for filtering your data and use a clustering approach to identify groups of cells of potentially poor quality that we want to remove before proceeding to downstream analyses.

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

The key steps in this notebook are:

1. Load in all sample data
2. Add metadata to Seurat objects
3. Assess count distributions
4. Identify and apply hard filtering thresholds
5. Transform and normalise data
6. Run clustering
7. Identify poor quality clusters and remove
8. Repeat data transformation, normalisation, and clustering

:::

## Instructions

Quality control is inherently a manual process. This notebook is desiged to be worked through interactively. Run each block of code and inspect the output. In a few places, you will need to create a CSV (Comma-Separated Values) file for importing your data and metadata about your samples, including filtering thresholds. In other places, we have developed some interactive R Shiny apps that will help you quickly determine parameters to set for further analysis.

Once you have worked through the entire notebook, we recommend rendering it into an HTML file that will serve as a record of how you processed your samples.

## 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, output = FALSE, warning = FALSE}
# Imports
knitr::opts_chunk$set(echo = TRUE)
library(here)
library(tidyverse)
library(Seurat)
library(hdf5r)
library(glmGamPoi)
library(clustree)
library(celldex)
library(SingleR)
library(DT)
source(here("R/shinyapps.R"))

# Helper functions
plot_filtered <- function(all_metadata, qc_list, filtered = FALSE) {
  # Plots static image of QC with metadata thresholds
  if (filtered) {
    tmp_all_metadata <- all_metadata %>%
      dplyr::filter(nCount_RNA <= qc_list$ncount_upper & nCount_RNA >= qc_list$ncount_lower) %>%
      dplyr::filter(nFeature_RNA <= qc_list$nfeature_upper & nFeature_RNA >= qc_list$nfeature_lower) %>%
      dplyr::filter(percent.mt <= qc_list$mtpercent_upper & percent.mt >= qc_list$mtpercent_lower)
    
    title <- "After Filtering"
  } else {
    tmp_all_metadata <- all_metadata
    title <- "Before Filtering"
  }
  tmp_all_metadata %>%
    ggplot(aes(x = nCount_RNA, y = nFeature_RNA, colour = percent.mt)) +
    geom_point(size = 0.1) +
    scale_color_viridis_c() +
    facet_wrap(~ orig.ident) +
    scale_x_log10() +
    scale_y_log10() +
    theme_light() +
    annotation_logticks(colour = "lightgrey") +
    ggtitle(title)
}

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)
}

# Create XDG_DATA_HOME directory in case it doesn't exist
data_dir <- Sys.getenv("XDG_DATA_HOME")
if(data_dir != "" && !dir.exists(data_dir)) {
  dir.create(data_dir, recursive = TRUE)
}

read_tmp_rds_files <- TRUE
```

## Define inputs {#define_inputs}

---

**❱❱❱ INPUT DATA ❰❰❰**

This notebook has been designed to be run after initial alignment and counting with the [nf-core/scrnaseq](https://nf-co.re/scrnaseq) Nextflow pipeline. This pipeline outputs one `.Rds` file containing a Seurat data object per sample. **Importantly**, this pipeline expects the **filtered** Seurat data, which has had non-cell barcodes removed. For version `4.0.0` of the `nf-core/scrnaseq` pipeline (the latest version at the time of publishing these notebooks), these RDS files can be found in the `results/` output directory at `results/cellranger/mtx_conversions/<SAMPLE_ID>/<SAMPLE_ID>_filtered_matrix.seurat.rds`.

---

In the following block, we load a CSV sample sheet `inputs/samplesheet.csv` which defines your samples and any metadata you would like to annotate them with. A template file already exists at this location; you can modify it to point to your single cell data, add the appropriate sample IDs and metadata columns. Then run the following block of code to read in your data and annotate it. The CSV file should have at least two columns, `sample` and `rds_path`. Additional columns will be treated as metadata variables that will be added to the Seurat objects. The value under the `sample` column will be used as the sample's name in the Seurat object. The value under the `rds_path` column should be the path to the RDS file generated by `nf-core/scrnaseq`. For example:

```
sample,rds_path,tumour
normal_sample,/path/to/normal_sample.rds,neg
tumour_sample,/path/to/tumour_sample.rds,pos
```

---

**❱❱❱ ACTION ❰❰❰**

1. Open the file `inputs/samplesheet.csv` in a text editor. To do this within RStudio:
    1. In the bottom right pane of RStudio, go to the `Files` tab.
    2. Navigate to the `inputs` folder.
    3. Click on `samplesheet.csv`. If a menu comes up, click on the `View File` option.
2. Update the contents of this file to define your samples and metadata.
    1. In the first row, ensure that the first two columns headers remain `sample` and `rds_path`. Add any other columns relevant to your data.
    2. In the remaining rows, ensure that the first two columns are the sample name (no spaces) and the path to the `.Rds` file.
    3. Ensure each sample has the same number of columns as the header row and that all metadata is filled in.
3. Save the file and close.
4. Run the next code block.

---

```{r inputs, warning = FALSE}
samples <- read_csv(here("inputs/samplesheet.csv"), show_col_types = FALSE)
cols <- colnames(samples)
cols[1:2] <- c("Sample", "RDS Path")
datatable(samples, colnames = cols)
```

### Generate input files for downstream steps

The following code chunks will generate all the input files for the downstream processing steps in this notebook and the following notebooks.

Most of these files will be templates filled with placeholder values that you will need to update based on the QC results below.

```{r define_input_files}
cluster_filter_samplesheet_path <- here("inputs/cluster_filter.csv")
cluster_res_samplesheet_path <- here("inputs/cluster_resolutions.csv")
doublet_samplesheet_path <- here("inputs/doublet_samples.csv")
```

```{r generate_input_files, eval = FALSE}
tmp_df <- samples
tmp_df <- tmp_df %>%
  mutate(
    res = 1.0,
    clusters_to_remove = "0;1;2",
    multiplet_rate = ""
  )

tmp_df %>%
  dplyr::select(sample, res, clusters_to_remove) %>%
  write_csv(cluster_filter_samplesheet_path)

tmp_df %>%
  dplyr::select(sample, res) %>%
  write_csv(cluster_res_samplesheet_path)

tmp_df %>%
  dplyr::select(sample, multiplet_rate) %>%
  write_csv(doublet_samplesheet_path)

rm(tmp_df)
```

### Read in RDS data

We now read in all of your samples and store them as a list of Seurat objects. **Please note** that we require all of the Seurat objects to have consistent gene IDs, i.e. they must all be using either gene symbols or Ensembl IDs.

If you have some Seurat objects that use gene symbols and others that use Ensembl IDs, this loading step will throw an error. To make all of your samples consistent, one option is to run all of the FASTQs through the `nf-core/scrnaseq` pipeline. This will ensure that they have all been processed in an identical way and that their gene IDs are consistent with one another.

An alternative option is to convert gene symbols to Ensembl IDs for those samples that are using gene symbols. We don't recommend converting Ensembl IDs to gene symbols, as some genes can have multiple gene symbols, depending on the naming convention used; in contrast, each gene has only one Ensembl ID associated with it.

---

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

If you would like to automatically convert any gene symbols to Ensembl IDs, set the following variable to `TRUE`.

**Please be aware** that in this process, some genes may not successfully convert to Ensembl IDs and will be lost. Additionally, some genes may map to the same Ensembl ID, resulting in duplicate IDs, which will be dropped. As such, we highly recommend that you don't enable this and instead re-run `nf-core/scrnaseq` for all samples if possible.

---

```{r convert_gene_symbols_to_ensembl}
# Set to TRUE if you want to automatically convert gene symbols to Ensembl IDs
# This will only be used if there are inconsistent uses of gene symbols and Ensembl IDs in your data
convert_gene_to_ens <- FALSE
```

```{r read_inputs}
all_seurat_objects <- lapply(transpose(samples), function(x) {
  # Read in Seurat object from RDS file
  s <- readRDS(here(x$rds_path))
  # Ensure sample names are assigned
  s@project.name <- x$sample
  s$orig.ident <- as.factor(x$sample)
  # Add additional metadata (if available)
  md <- x[! names(x) %in% c("sample", "rds_path")]
  if (!is.null(md)) {
    for (m in names(md)) {
      s[[m]] <- md[[m]]
    }
  }

  # Ensure both gene symbols (`gene_symbols`) and Ensembl IDs (`gene_versions`) are present in RNA assay metadata
  ens_rownames <- all(startsWith(rownames(s@assays$RNA), "ENS"))
  if(ens_rownames) {
    ens_ids <- rownames(s@assays$RNA)
    gene_symbols <- AnnotationDbi::mapIds(
      EnsDb.Hsapiens.v86::EnsDb.Hsapiens.v86,
      keys = ens_ids,
      column = "SYMBOL",
      keytype = "GENEID"
    )
    stopifnot(all(ens_ids == names(gene_symbols)))
  } else {
    gene_symbols <- rownames(s@assays$RNA)
    ens_ids <- AnnotationDbi::mapIds(
      EnsDb.Hsapiens.v86::EnsDb.Hsapiens.v86,
      keys = gene_symbols,
      column = "GENEID",
      keytype = "SYMBOL"
    )
    stopifnot(all(gene_symbols == names(ens_ids)))
  }
  have_ens_ids <- "gene_versions" %in% colnames(s@assays$RNA@meta.data) && all(startsWith(as.character(s@assays$RNA@meta.data$gene_versions), "ENS"))
  have_gene_symbols <- "gene_symbols" %in% colnames(s@assays$RNA@meta.data) && !all(startsWith(as.character(s@assays$RNA@meta.data$gene_symbols), "ENS"))
  if(!have_ens_ids) {
    s@assays$RNA@meta.data$gene_versions <- ens_ids
  }
  if(!have_gene_symbols) {
    s@assays$RNA@meta.data$gene_symbols <- gene_symbols
  }
  
  return(s)
})
names(all_seurat_objects) <- sapply(all_seurat_objects, function(s) { as.character(s@meta.data$orig.ident[[1]]) })

# If there is a mixture of gene symbols and Ensembl IDs for the RNA assay rownames
# and `convert_gene_to_ens` is `TRUE`
# then convert any gene symbol rownames to Ensembl IDs
using_ens_ids <- sapply(all_seurat_objects, function(s) {
  all(startsWith(rownames(s@assays$RNA), "ENS"))
})
if(any(using_ens_ids) && !all(using_ens_ids)) {
  stopifnot(convert_gene_to_ens)
  all_seurat_objects <- lapply(all_seurat_objects, function(s) {
    if(all(startsWith(rownames(s@assays$RNA), "ENS"))) {
      return(s)
    } else {
      if(any(duplicated(s@assays$RNA@meta.data$gene_versions))) {
        warning("Duplicated Ensembl IDs have been detected when converting from gene symbols. For each duplicated ID, the first will be kept and the rest will be dropped.")
        keep <- !duplicated(s@assays$RNA@meta.data$gene_versions)
        s <- s[keep,]
      }
      if(any(is.na(s@assays$RNA@meta.data$gene_versions))) {
        warning("Missing Ensembl IDs have been detected when converting from gene symbols. These will be dropped.")
        keep <- !is.na(s@assays$RNA@meta.data$gene_versions)
        s <- s[keep,]
      }
      rownames(s@assays$RNA) <- s@assays$RNA@meta.data$gene_versions
      return(s)
    }
  })
}
```

If you are working with whole-cell data, you should expect a small percentage of your data to be reads from mitochondrial genes. In the following chunk, you need to specify the mitochondrial gene prefix so we can calculate the percentage of mitochondrial reads in your data; a high percentage is indicative of poor quality or dying cells.

Human mitochondrial genes all start with the prefix `^MT-`, while in mice this will be `^mt-`.

If you are instead working with single nuclei samples, or otherwise expect no mitochondrial genes in your dataset, set the value of `mt_pattern` to `NULL`. This will cause the `mt.percentage` annotation in each Seurat object to be set to `0`.

```{r add_mt_percentatge}
# The following assumes human data
# Change to "^mt-" for mice
# Change to NULL for nuclei data
mt_pattern <- "^MT-"

if (!is.null(mt_pattern)) {
  all_seurat_objects <- lapply(all_seurat_objects, function(s) {
    so <- s
    rn <- rownames(so@assays$RNA)
    if(all(startsWith(rn, "ENS"))) {
      rownames(so@assays$RNA) <- so@assays$RNA@meta.data$gene_symbols
      so$percent.mt <- Seurat::PercentageFeatureSet(so, pattern = mt_pattern)
      rownames(so@assays$RNA) <- rn
    } else {
      so$percent.mt <- Seurat::PercentageFeatureSet(so, pattern = mt_pattern)
    }
    return(so)
  })
} else {
  all_seurat_objects <- lapply(all_seurat_objects, function(s) {
    so <- s
    so$percent.mt <- 0  # NOTE: Setting default to 0 for now so it works with downstream code that expects percent.mt
    return(so)
  })
}
```

We will also gather all of our metadata up in a single data frame for easy plotting.

```{r get_all_metadata}
all_metadata <- lapply(all_seurat_objects, function(s) {
  s@meta.data %>%
    rownames_to_column(var = "barcode")
}) %>% bind_rows()
```

## Initial QC

First, we can check the number of cells each sample has:

```{r total_cells}
total_cells <- sapply(all_seurat_objects, function(x) length(x$orig.ident))
tibble(sample = names(total_cells), n_cells = total_cells) %>% datatable(colnames = c("Sample", "# Cells"))
```

Next, plot the distribution of nCount_RNA, nFeature_RNA, and (if available), percent.mt for each sample:

```{r metadata_distributions, warning = FALSE}
log_y <- FALSE  # Set to `TRUE` if you want log-scaled plots
jitter <- FALSE  # Set to `TRUE` if you want to plot the individual data points

md_cols <- c("nCount_RNA", "nFeature_RNA")
if ("percent.mt" %in% colnames(all_metadata)) {
  md_cols <- c(md_cols, "percent.mt")
}

p_md <- all_metadata %>%
  pivot_longer(cols = md_cols) %>%
  ggplot(aes(x = orig.ident, y = value, fill = orig.ident, colour = orig.ident)) +
  facet_wrap(~ name, scales = "free_y") +
  theme_light() +
  theme(legend.position = "none")

if (log_y) {
  p_md <- p_md + scale_y_log10()
}
if (jitter) {
  p_md <- p_md + geom_jitter(size = 0.1, alpha = 0.7, colour = "lightgrey")
}

p_md + geom_violin(alpha = 0.1) 
```

### Identify filtering thresholds

The following Shiny app will help you determine thresholds to use for filtering your data. Play with the thresholds until you are happy with the data quality, then note the final thresholds for use in the following code block.

```{r app_qc_thresholds, eval = FALSE}
app_qc_thresholds(all_metadata)
```

---

**❱❱❱ ACTION ❰❰❰**

In the following chunk we have defined a list object `qc_thresholds` containing six variables. Set your desired thresholds for each of the variables. The code will also generate a plot of your filtered data for you to inspect prior to committing. The threshold variables are defined in the following table:

| Variable          | Default | Description                                                                        |
| ----------------- | ------- | ---------------------------------------------------------------------------------- |
| `ncount_lower`    | 0       | Minimum number of total transcripts required for a cell to be kept                 |
| `nfeature_lower`  | 0       | Minimum number of mapped genes required for a cell to be kept                      |
| `mtpercent_upper` | 100     | Maximum allowed percentage of mitochondrial gene counts per cell                   |
| `ncount_upper`    | Inf     | Maximum number of total transcripts allowed for a cell                             |
| `nfeature_upper`  | Inf     | Maximum number of mapped genes allowed for a cell                                  |
| `mtpercent_lower` | 0       | Minimum allowed percentage of mitochondrial gene counts per cell (**rarely used**) |

---

```{r set_thresholds}
# These defaults will result in no filtering being applied
# TODO: Set your desired thresholds
qc_thresholds <- list()
qc_thresholds$ncount_lower <- 0
qc_thresholds$nfeature_lower <- 0
qc_thresholds$mtpercent_upper <- 100

# You generally won't want to change the upper count and lower MT percentage thresholds
# unless you identify outlying cells that should be removed
qc_thresholds$ncount_upper <- Inf
qc_thresholds$nfeature_upper <- Inf
qc_thresholds$mtpercent_lower <- 0

# Plot data, before and after filtering
plot_filtered(all_metadata, qc_thresholds, filtered = FALSE)
plot_filtered(all_metadata, qc_thresholds, filtered = TRUE)
```

Apply the filtering thresholds:

```{r add_thresholds, message = FALSE}
all_seurat_objects <- lapply(all_seurat_objects, function(s) {
  s2 <- s
  s2@meta.data <- s2@meta.data %>%
  mutate(
    qc_threshold = case_when(
      percent.mt >= qc_thresholds$mtpercent_lower &
        percent.mt <= qc_thresholds$mtpercent_upper &
        nFeature_RNA >= qc_thresholds$nfeature_lower &
        nFeature_RNA <= qc_thresholds$nfeature_upper &
        nCount_RNA >= qc_thresholds$ncount_lower &
        nCount_RNA <= qc_thresholds$ncount_upper ~ "keep",
      .default = "remove"
    )
  )
  s2
})

all_filtered_seurat_objects <- lapply(all_seurat_objects, function(s) {
  subset(s, subset = (qc_threshold == "keep"))
})

n_cells_pre_filtered <- sapply(all_seurat_objects, function(s) { nrow(s@meta.data) })
n_cells_post_filtered <- sapply(all_filtered_seurat_objects, function(s) { nrow(s@meta.data) })

n_cells_pre_filtered <- tibble(sample = names(n_cells_pre_filtered), n_cells_pre_filtered = n_cells_pre_filtered)
n_cells_post_filtered <- tibble(sample = names(n_cells_post_filtered), n_cells_post_filtered = n_cells_post_filtered)

filter_summary <- n_cells_pre_filtered %>% left_join(n_cells_post_filtered) %>%
  mutate(prop_filtered = 1 - n_cells_post_filtered/n_cells_pre_filtered) %>%
  mutate(across(where(is.numeric), ~ round(.x, 4)))
datatable(filter_summary, colnames = c("Sample", "# Cells (Pre-Filter)", "# Cells (Post-Filter)", "Proportion of Cells Filtered"))

rm(all_seurat_objects)  # We no longer need the original Seurat objects
```

## SCTransform and initial clustering for further filtering

Now that we have applied some basic quality filters, we will now run the `SCTransform` method in Seurat, which is a critical pre-processing step that normalises and corrects the single cell count data.

We also want to perform some initial clustering of our single cell data to look for potentially poor-quality groups of cells that we should remove before running any analyses. To do this, we will use Seurat's `FindClusters` method. This method first constructs a nearest-neighbours graph of the data, identifying which cells are closely related to each other in terms of gene expression. Then this graph is used to detect "communities" of cells, using one of two main algorithms.

The original and default clustering algorithm that Seurat uses is the Louvain method. This has some pitfalls which has led to it being superceded by the Leiden algorithm. We recommend using this updated method, however it requires first ensuring that you have the `leidenalg` Python package installed.

---

**❱❱❱ ACTION ❰❰❰**

In the following code chunk, set the `cluster_method` variable to either `"LOUVAIN"` or `"LEIDEN"`. We will also check if you are able to use the Leiden method; if you don't have the appropriate packages installed, the following code will error when trying `cluster_method <- "LEIDEN"`.

---

```{r set_cluster_method}
cluster_method <- "LOUVAIN"  # Set this to "LEIDEN" to use the more up-to-date method

# Set the algorithm parameter value for Seurat::FindClusters
if (cluster_method == "LOUVAIN") {
  cluster_algorithm <- 1
} else {
  cluster_algorithm <- 4
}
```

```{r check_cluster_method, include = FALSE}
# Run this chunk to check whether you have configured your clustering method correctly
# and to check if the Leiden method is available to use

stopifnot(cluster_method %in% c("LOUVAIN", "LEIDEN"))

can_use_leiden <- system("python3 -c 'import leidenalg'", ignore.stderr = TRUE) == 0

if (cluster_method == "LOUVAIN" && can_use_leiden) {
  warning("You have chosen to use the older Louvain clustering algorithm, but you can use the more up-to-date Leiden method, which we recommend. Please consider updating the `cluster_method` variable to 'LEIDEN'.")
} else if (cluster_method == "LEIDEN" && !can_use_leiden) {
  stop("You have chosen to use the up-to-date Leiden clustering algorithm, but you do not have the required `leidenalg` Python package installed. If you wish to use the Leiden method, please ensure you have Python installed as well as the `leidenalg` package. You can install this package with: `pip install leidenalg`. Otherwise, please change the `cluster_method` variable to 'LOUVAIN'.")
}
```

We will now save your cluster algorithm choice to `inputs/cluster_algorithm.txt` for use in the later notebooks.

```{r save_cluster_method, eval = FALSE}
sink(here("inputs/cluster_algorithm.txt"))
cat(cluster_algorithm)
sink()
```

For clustering with either algorithm, you will need to specify a resolution parameter. This controls how fine or coarse the resulting clusters are; smaller values result in fewer, larger clusters, while higher values result in more smaller clusters. Typically, values between 0.4 and 1.2 will return good results [as noted by the authors of Seurat](https://satijalab.org/seurat/articles/pbmc3k_tutorial.html#cluster-the-cells), but it is a good idea to explore a wider range. Seurat handily supports clustering at multiple resolutions and storing the results of each clustering pass as separate metadata columns. In the following code block we set the range of resolution parameter values to between `0.2` and `2.0`, increasing in steps of `0.2` which should cover a large enough range of values for most datasets. Feel free to change this to suit your data.

---

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

If you want to change the range of clustering resolutions to try, do so in the next block.

**Note** that these resolutions will be used for the remainder of the analyses in these notebooks.

---

```{r sct_resolutions}
# Set the range of clustering resolutions to try
min_res <- 0.6
max_res <- 2.4
step_size <- 0.2

res <- seq(min_res, max_res, step_size)
```

We will write this range of resolutions to `inputs/all_cluster_resolutions.txt` so we can use them in later notebooks.

```{r write_sct_resolutions, eval = FALSE}
# Write the cluster resolutions to file for use later on
all_cluster_res_file <- here("inputs/all_cluster_resolutions.txt")
sink(all_cluster_res_file)
cat(paste(res, collapse = "\n"))
sink()
```

Now run the following code block to run both the `SCTransform` normalisation and the clustering of your data.

```{r init_sct_reduce_dims, eval = FALSE, output = FALSE, warning = FALSE}
options(future.globals.maxSize = 1000*1024^2)

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

all_init_sct <- lapply(all_filtered_seurat_objects, function(s) {
  # Run SCTransform and find clusters
  s_sct <- SCTransform_reduce_dims(s)
  s_sct <- FindClusters(s_sct, resolution = res, algorithm = cluster_algorithm, random.seed = random_seed, verbose = 0)
  gc()  # Clean up memory after each run
  s_sct
})
```

In order to determine the best cluster resolution to use, we can use the `clustree` package to create a clustering tree graph. This summarises how the cells move between clusters when we change the clustering resolution. We are looking for a range of clustering resolutions where the clusters remain relatively stable; in contrast, we want to avoid choosing unstable resolutions that cause large changes in the cluster assignments.

The following block will generate these clustering trees.

```{r init_create_clustertrees, eval = FALSE}
all_init_clustrees <- lapply(all_init_sct, function(s) {
  s_name <- as.character(s@meta.data$orig.ident[[1]])
  # Plot cluster trees
  clustree::clustree(s, prefix = "SCT_snn_res.") +
    ggtitle(s_name) +
    theme(
      legend.position = "right",
      legend.box = "vertical",
      legend.justification.right = "top"
    )
})

dir.create(here("tmp_outputs", "01.qc"), recursive = TRUE)
saveRDS(all_init_clustrees, here("tmp_outputs", "01.qc", "all_init_clustrees.Rds"))
read_tmp_rds_files <- FALSE
```

Now we display the cluster trees to aid in picking a cluster resolution:

```{r init_plot_clustertrees, fig.width = 10, fig.height = 10}
if(read_tmp_rds_files) {
  all_init_clustrees <- readRDS(here("tmp_outputs", "01.qc", "all_init_clustrees.Rds"))
}

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

::: {.callout-tip title="Interpreting clustree plots"}

There is a lot going on in these clustree graphs.

First, each row represents a distinct clustering resolution, with the larger, coarse clusters at the top and the smaller, fine clusters at the bottom. Each row is colour-coded, with the corresponding clustering resolution shown in the legend under `SCT_snn_res.`.

Each circle represents a cluster, with the circle's size corresponding to the number of cells in that cluster. In the legend, the `size` section shows the relationship between cluster size and the size of the circles.

The arrows between clusters of neighbouring rows represent the movement of cells between clusters when the resolution increases. The transparency of the arrows indicates the "in-proportion" - the proportion of cells moving into a cluster relative to that cluster's overall size. This is labelled as `in_prop` in the legend. The arrows are also colour-coded by the raw number of cells moving between clusters; this is labelled in the legend as `count`.

Stable, well-defined clusters tend to remain unchanged between clustering resolutions, with a single dark arrow connecting the clusters.

Unstable clusters tend to have multiple, more transparent incoming arrows, indicating that the cells are difficult to cluster and "jump around" between clusters when moving from one resolution to the next. This is an indication of over-clustering and that a smaller clustering resolution is more appropriate.

When choosing a clustering resolution, you want to select one where the clusters are relatively stable, with clusters having singular, darker incoming arrows and as few cases of multiple faint incoming arrows as possible.

:::

Running the next block will open up an interactive Shiny app that you can use to inspect each of your clustering resolutions and the effect they have on the data:

```{r init_app_explore_clusters, eval = FALSE}
app_explore_clusters(all_init_sct, all_init_clustrees, res)
```

### Define cluster resolutions and clusters to remove

Based on the cluster tree graphs and the exploratory analysis of the clusters at each resolution, you should now pick a stable cluster resolution for each sample and define any cluster IDs that should be removed. For example, clusters that have high percentage MT counts or consistently low read counts are likely artefacts that you will want to exclude from downstream analyses.

You can specify your cluster resolutions and IDs to remove in the file `inputs/cluster_filter.csv` that was generated earlier. The file contains three columns: `sample`, `res`, and `clusters_to_remove`, e.g.:

```         
sample,res,clusters_to_remove
normal_sample,1.2,15;20
tumour_sample,0.8,16
```

Multiple clusters to remove should be separated by semicolons. If you have no clusters to remove, simply leave the `clusters_to_remove` column blank.

---

**❱❱❱ ACTION ❰❰❰**

1. Open `inputs/cluster_filter.csv` in RStudio or an external text editor.
2. Update the cluster resolutions in the second column.
3. Update the clusters to remove in the third column, ensuring multiple clusters are separated by semicolons (`;`).

---

```{r read_filter_clusters}
cluster_filter_samplesheet <- read_csv(cluster_filter_samplesheet_path, show_col_types = FALSE)
datatable(cluster_filter_samplesheet, colnames = c("Sample", "Cluster Resolution", "Clusters to Remove"))
```

The following code block will generate the per-cluster count plots for each sample, using the resolutions specified in the above samplesheet.

```{r prep_plot_init_clusters, eval = FALSE}
all_init_sct_cluster_plots <- lapply(all_init_sct, function(s) {
  s_name <- as.character(s@meta.data$orig.ident[[1]])
  s_res <- cluster_filter_samplesheet$res[match(s_name, cluster_filter_samplesheet$sample)]
  s_res_named <- paste0("SCT_snn_res.", s_res)
  s@meta.data %>%
    ggplot(aes(x = nCount_RNA, y = nFeature_RNA, col = percent.mt)) +
      geom_point(size = 0.3) +
      facet_wrap(s_res_named) +
      scale_x_log10() +
      scale_y_log10() +
      theme_light() +
      viridis::scale_color_viridis() +
      annotation_logticks(side = "lb", colour = "lightgrey") +
      ggtitle(paste0(s_name, ": ", s_res_named))
})

saveRDS(all_init_sct_cluster_plots, here("tmp_outputs", "01.qc", "all_init_sct_cluster_plots.Rds"))
```

Plot the per-cluster count plots for each sample:

```{r plot_init_clusters}
if(read_tmp_rds_files) {
  all_init_sct_cluster_plots <- readRDS(here("tmp_outputs", "01.qc", "all_init_sct_cluster_plots.Rds"))
}

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

Now, before removing our suspect clusters, we want to make sure that these cells aren't clustering due to real biological variation. So, we will look at the top 20 most differentially expressed genes in each cluster. Take note of any clusters where multiple genes may be involved in the same biological process - in particular, biological processes of interest to your research. This could indicate real biological variation and you may wish to keep these cells.

```{r get_bad_cluster_degs, eval = FALSE}
bad_cluster_degs <- lapply(all_init_sct, function(s) {
  s_name <- as.character(s@meta.data$orig.ident[[1]])
  cluster_res <- as.numeric(cluster_filter_samplesheet$res[match(s_name, cluster_filter_samplesheet$sample)])
  cluster_res <- paste0("SCT_snn_res.", cluster_res)
  stopifnot(cluster_res %in% names(s@meta.data))
  clusters_to_remove <- as.character(cluster_filter_samplesheet$clusters_to_remove[match(s_name, cluster_filter_samplesheet$sample)])
  if (clusters_to_remove == "" | is.na(clusters_to_remove)) {
    return(NULL)
  }
  clusters_to_remove <- as.integer(strsplit(clusters_to_remove, ";")[[1]])
  Idents(s) <- cluster_res

  degs_list <- lapply(clusters_to_remove, function(cluster) {
    FindMarkers(s, ident.1 = cluster) %>%
      head(20) %>%
      arrange(desc(avg_log2FC)) %>%
      mutate(sample = s_name, cluster = cluster)
  })

  degs_combined <- do.call(rbind, degs_list)

  # Get `gene_symbol` and `ensembl_id` columns
  idx <- match(rownames(degs_combined), rownames(s@assays$RNA))
  if(all(startsWith(rownames(degs_combined), "ENS"))) {
    degs_combined$gene_symbol <- s@assays$RNA@meta.data$gene_symbols[idx]
    degs_combined$ensembl_id <- rownames(degs_combined)
  } else {
    degs_combined$gene_symbol <- rownames(degs_combined)
    degs_combined$ensembl_id <- s@assays$RNA@meta.data$gene_versions[idx]
  }
  return(degs_combined)
})

saveRDS(bad_cluster_degs, here("tmp_outputs", "01.qc", "bad_cluster_degs.Rds"))
```

Run the following block to display the table of top differentially expressed genes per cluster:

```{r print_bad_cluster_degs, warning = FALSE}
if(read_tmp_rds_files) {
  bad_cluster_degs <- readRDS(here("tmp_outputs", "01.qc", "bad_cluster_degs.Rds"))
}

combined_degs <- bind_rows(bad_cluster_degs)

if (nrow(combined_degs) > 0) {
  first_cols <- c("gene_symbol", "ensembl_id")
  remaining_cols <- colnames(combined_degs)[!(colnames(combined_degs) %in% first_cols)]
  all_cols <- c(first_cols, remaining_cols)

  combined_degs <- combined_degs %>%
    dplyr::select(all_of(all_cols)) %>%
    mutate(across(where(is.numeric), ~ round(.x, 4)))

  datatable(
    combined_degs,
    colnames = c(
      "Gene Symbol",
      "Ensembl ID",
      "p-value",
      "Avg. log2FC",
      "% Cells (Cond. 1)",
      "% Cells (Cond. 2)",
      "Adjusted p-value",
      "Sample",
      "Cluster"
    )
  )
} else {
  cat("No clusters selected for removal.")
}
```

If you are happy with the current list of clusters to remove for each sample, continue. Otherwise, you should revise the list of clusters you have defined in `inputs/cluster_filter.csv` and re-run the code blocks in this section again.

---

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

If you want to revise the list of clusters to remove:

1. Open `inputs/cluster_filter.csv` in RStudio or an external text editor.
2. Update the list of clusters to remove per sample in the third column.
3. Save and close the file.
4. Go back to the start of this section ([**Define cluster resolutions and clusters to remove**](#define-cluster-resolutions-and-clusters-to-remove)).
5. Re-run all code blocks in this section.

---

### Apply cluster filtering

Now we will apply the cluster filtering to our original Seurat objects.

```{r define_bad_clusters, eval = FALSE}
cells_to_keep <- lapply(all_init_sct, function(s) {
  s_name <- as.character(s@meta.data$orig.ident[[1]])
  cluster_res <- as.numeric(cluster_filter_samplesheet$res[match(s_name, cluster_filter_samplesheet$sample)])
  cluster_res <- paste0("SCT_snn_res.", cluster_res)
  stopifnot(cluster_res %in% names(s@meta.data))
  clusters_to_remove <- as.character(cluster_filter_samplesheet$clusters_to_remove[match(s_name, cluster_filter_samplesheet$sample)])
  if (clusters_to_remove == "" | is.na(clusters_to_remove)) {
    return(NULL)
  }
  clusters_to_remove <- as.integer(strsplit(clusters_to_remove, ";")[[1]])
  Idents(s) <- cluster_res
  
  rm_cells <- s@meta.data[[cluster_res]] %in% clusters_to_remove
  keep_cells <- !rm_cells
  keep_cells <- rownames(s@meta.data[keep_cells,])
  return(keep_cells)
})

saveRDS(cells_to_keep, here("tmp_outputs", "01.qc", "cells_to_keep.Rds"))
```

```{r remove_bad_clusters, message = FALSE}
if(read_tmp_rds_files) {
  cells_to_keep <- readRDS(here("tmp_outputs", "01.qc", "cells_to_keep.Rds"))
}

all_cluster_filtered_seurat_objects <- lapply(all_filtered_seurat_objects, function(s) {
  s_name <- as.character(s@meta.data$orig.ident[[1]])
  
  keep_cells <- cells_to_keep[[s_name]]
  
  if (is.null(keep_cells)) {
    return(s)
  }
  
  s %>% subset(cells = keep_cells)
})

n_cells_post_cluster_filtered <- sapply(all_cluster_filtered_seurat_objects, function(s) { nrow(s@meta.data) })
n_cells_post_cluster_filtered <- tibble(sample = names(n_cells_post_cluster_filtered), n_cells_post_cluster_filtered = n_cells_post_cluster_filtered)

filter_summary <- n_cells_pre_filtered %>% left_join(n_cells_post_filtered) %>% left_join(n_cells_post_cluster_filtered) %>%
  mutate(prop_filtered = 1 - n_cells_post_cluster_filtered/n_cells_post_filtered) %>%
  mutate(across(where(is.numeric), ~ round(.x, 4)))
datatable(filter_summary, colnames = c("Sample", "# Cells (Pre-Filter)", "# Cells (Post-Hard-Filter)", "# Cells (Post-Cluster-Filter)", "Proportion of Cells Filtered"))

rm(all_filtered_seurat_objects)  # We no longer need this version of the Seurat objects
```

## Repeat SCTransform and perform final clustering

Now that we have filtered our data, we want to re-run `SCTransform` and clustering, as the filtering process will invalidate the previous round of normalisation and clustering. We will use the same range of clustering resolutions that we defined above in the `res` variable.

```{r sct_reduce_dims, eval = FALSE, output = FALSE, warning = FALSE}
options(future.globals.maxSize = 1000*1024^2)

all_sct <- lapply(all_cluster_filtered_seurat_objects, function(s) {
  # Run SCTransform and find clusters
  s_sct <- SCTransform_reduce_dims(s)
  s_sct <- FindClusters(s_sct, resolution = res, algorithm = cluster_algorithm, random.seed = random_seed, verbose = 0)
  gc()  # Clean up memory after each run
  s_sct
})

rm(all_cluster_filtered_seurat_objects)  # We no longer need this version of the Seurat objects
```

Again, we generate clustering trees to help us decide on a final resolution to use:

```{r create_clustertrees, eval = FALSE}
all_clustrees <- lapply(all_sct, function(s) {
  s_name <- as.character(s@meta.data$orig.ident[[1]])
  # Plot cluster trees
  clustree::clustree(s, prefix = "SCT_snn_res.") +
    ggtitle(s_name) +
    theme(
      legend.position = "right",
      legend.box = "vertical",
      legend.justification.right = "top"
    )
})

saveRDS(all_clustrees, here("tmp_outputs", "01.qc", "all_clustrees.Rds"))
```

```{r plot_clustertrees, fig.width = 10, fig.height = 10}
if(read_tmp_rds_files) {
  all_clustrees <- readRDS(here("tmp_outputs", "01.qc", "all_clustrees.Rds"))
}

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

We can re-run the interactive Shiny app to help us inspect the data for each sample at each clustering resolution:

```{r app_explore_clusters, eval = FALSE}
app_explore_clusters(all_sct, all_clustrees, res)
```

### Define final cluster resolutions

Once again, based on the cluster tree graphs and the exploratory analysis of the clusters at each resolution, you should now pick a final stable cluster resolution for each sample.

You can specify your cluster resolutions using the CSV file `inputs/cluster_resolutions.csv` generated earlier. It has two columns: `sample` and `res`, e.g.:

```         
sample,res
normal_sample,1.2
tumour_sample,0.8
```

Note that the `res` column values will likely be different to the values used in `inputs/cluster_filter.csv`.

---

**❱❱❱ ACTION ❰❰❰**

1. Open `inputs/cluster_resolutions.csv` in RStudio or an external text editor.
2. For each sample, update the cluster resolution to use for clustering in the second column.

---

```{r filter_clusters}
cluster_res_samplesheet <- read_csv(cluster_res_samplesheet_path, show_col_types = FALSE)
datatable(cluster_res_samplesheet, colnames = c("Sample", "Cluster Resolution"))
```

Plot the per-cluster count plots for each sample at the chosen clustering resolutions:

```{r prep_plot_clusters, eval = FALSE}
all_sct_cluster_plots <- lapply(all_sct, function(s) {
  s_name <- as.character(s@meta.data$orig.ident[[1]])
  s_res <- cluster_res_samplesheet$res[match(s_name, cluster_res_samplesheet$sample)]
  s_res_named <- paste0("SCT_snn_res.", s_res)
  s@meta.data %>%
    ggplot(aes(x = nCount_RNA, y = nFeature_RNA, col = percent.mt)) +
      geom_point(size = 0.3) +
      facet_wrap(s_res_named) +
      scale_x_log10() +
      scale_y_log10() +
      theme_light() +
      viridis::scale_color_viridis() +
      annotation_logticks(side = "lb", colour = "lightgrey") +
      ggtitle(paste0(s_name, ": ", s_res_named))
})

saveRDS(all_sct_cluster_plots, here("tmp_outputs", "01.qc", "all_sct_cluster_plots.Rds"))
```

```{r plot_clusters}
if(read_tmp_rds_files) {
  all_sct_cluster_plots <- readRDS(here("tmp_outputs", "01.qc", "all_sct_cluster_plots.Rds"))
}

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

Finally, apply the chosen clustering resolution by setting each Seurat object's Identity to that resolution.

```{r apply_cluster_resolutions, eval = FALSE}
all_sct <- lapply(all_sct, function(s) {
  s_name <- as.character(s@meta.data$orig.ident[[1]])
  cluster_res <- as.numeric(cluster_res_samplesheet$res[match(s_name, cluster_res_samplesheet$sample)])
  cluster_res <- paste0("SCT_snn_res.", cluster_res)
  Idents(s) <- cluster_res
  s
})
```

Now that we have a filtered dataset, we can save the data for each sample to a separate `.Rds` file:

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

## Package version summary

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