---
title: "scRNA Analysis - 04 - Cell Annotation"
author: "Sydney Informatics Hub"
format: html
---

This notebook is designed to annotate the cells and clusters within an integrated Seurat dataset. This includes cell cycle annotation, automatic cell type annotation, and optional manual cell type annotation with a manually-curated marker gene list.

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

The key steps in this notebook are:

1. Load in the single integrated dataset from the previous notebook
2. Annotate the cells with a predicted cell cycle state
3. Annotate the cells with a predicted cell type using the HPCA database and the `SingleR` package
4. Optionally annotate the cells with a predicted cell type using custom-defined gene marker sets

:::

## 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(scuttle)
library(ggplot2)
library(ggrepel)
library(SingleR)
library(DT)
source(here("R/shinyapps.R"))
read_tmp_rds_files <- TRUE
```

## Read in the integrated dataset

The following block assumes you have run all the notebooks up to and including `03.integration.qmd`. That notebook produces a `.Rds` file called `outputs/integrated_dataset.filtered_clustered.doublets_detected.integrated.Rds` that will be read in now:

```{r inputs}
dataset <- here("outputs/integrated_dataset.filtered_clustered.doublets_detected.integrated.Rds")

# Load data
integrated <- readRDS(dataset)
```

## Annotation 1: Cell Cycle

We will first annotate the cells with their predicted phase in the cell cycle. This is a common annotation step, particularly for cancer data sets.

```{r get_cell_cycle_genes, warning = FALSE}
# Check if the Seurat object has Ensembl IDs or gene symbols
using_ens_ids <- all(startsWith(rownames(integrated@assays$RNA), "ENS"))
if(using_ens_ids) {
  # Map Ensembl IDs to gene IDs
  s.genes <- AnnotationDbi::mapIds(
    EnsDb.Hsapiens.v86::EnsDb.Hsapiens.v86,
    keys = cc.genes$s.genes,
    column = "GENEID",
    keytype = "SYMBOL"
  )
  g2m.genes <- AnnotationDbi::mapIds(
    EnsDb.Hsapiens.v86::EnsDb.Hsapiens.v86,
    keys = cc.genes$g2m.genes,
    column = "GENEID",
    keytype = "SYMBOL"
  )
} else {
  s.genes <- cc.genes$s.genes
  g2m.genes <- cc.genes$g2m.genes
}
```

```{r cell_cycle_scoring, eval = FALSE}
integrated <- CellCycleScoring(
  integrated,
  s.features = s.genes,
  g2m.features = g2m.genes
)

gc()  # Clean up memory

dir.create(here("tmp_outputs", "04.annotation"), recursive = TRUE)
saveRDS(integrated, here("tmp_outputs", "04.annotation", "integrated.cc.Rds"))
read_tmp_rds_files <- FALSE
```

```{r read_cell_cycle_scoring}
if (read_tmp_rds_files) {
  integrated <- readRDS(here("tmp_outputs", "04.annotation", "integrated.cc.Rds"))
}

available_annots <- c("Phase")
```

### Plot cell cycle annotations

Let's plot the integrated dataset UMAP, colouring by cell cycle phase:

```{r plot_cell_cycle_scoring}
DimPlot(integrated, reduction = "umap", group.by = "Phase")
```

## Annotation 2: Automatic cell type assignment with HPCA

This next step will assign cell types to each cell using the `SingleR` and `celldex` packages and the HPCA database. We will add both the main, coarse-level annotations as well as the finer-scale annotations from the HPCA database.

**NOTE:** On first use, the function `celldex::HumanPrimaryCellAtlasData()` may create a cache in your home directory at `~/.cache/R/AnnotationHub`. This is a one-time setup step, but be aware that if you are on a HPC system such as gadi it may cause issues if your home directory is near quota. Please ensure that you have at least 1GB of free space available in your home directory.

---

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

The following block sets the lower threshold of the number of cells required to keep a given annotation. If a particular cell type annotation has been assigned to this many or fewer cells, that cell type will be discarded and the cell will be labelled as "Other".

By default, we have set this to `10`. If you are studying rare cell types, you may want to set variable `min_cell_annot_threshold` in the next block to a lower value, or even `0` to keep all HPCA annotations.

---

```{r set_min_cell_threshold_for_annotation}
min_cell_annot_threshold <- 10
```

```{r annotate_cell_types_hpca, eval = FALSE}
# Create cache dir if it doesn't already exist
hpca_cache_dir <- "~/.cache/R/AnnotationHub"
if(!dir.exists(hpca_cache_dir)) {
  dir.create(hpca_cache_dir, recursive = TRUE)
}

if(using_ens_ids) {
  ref <- celldex::HumanPrimaryCellAtlasData(ensembl = TRUE)
} else {
  ref <- celldex::HumanPrimaryCellAtlasData(ensembl = FALSE)
}

sce <- as.SingleCellExperiment(integrated, assay = "RNA")
sceM <- logNormCounts(sce)

# Add the main-level annotations
predicted <- SingleR(test = sceM, ref = ref, labels = ref$label.main)
keep <- table(predicted$labels) > min_cell_annot_threshold
integrated$SingleR.hpca_main <- ifelse(keep[predicted$labels], predicted$labels, "Other")

gc()  # Clean up memory

# Also add the fine-level annotations
predicted <- SingleR(test = sceM, ref = ref, labels = ref$label.fine)
keep <- table(predicted$labels) > min_cell_annot_threshold
integrated$SingleR.hpca_fine <- ifelse(keep[predicted$labels], predicted$labels, "Other")

rm(sce, sceM, keep, predicted, ref)
gc()  # Clean up memory

saveRDS(integrated, here("tmp_outputs", "04.annotation", "integrated.singler.Rds"))
```

```{r read_annotated_cell_types_hpca}
if (read_tmp_rds_files) {
  integrated <- readRDS(here("tmp_outputs", "04.annotation", "integrated.singler.Rds"))
}

available_annots <- c(available_annots, "SingleR.hpca_main", "SingleR.hpca_fine")
```

### Plot cell type assignment

Again, we can plot the UMAP of our integrated dataset, grouped by both our clusters as well as our HPCA cell type annotations:

```{r plot_cell_types_hpca, message = FALSE, warning = FALSE, fig.width = 16}
cluster_res_file <- here("inputs/integrated_cluster_resolution.txt")
cluster_res <- scan(cluster_res_file, numeric())
cluster_name <- paste0("SCT_snn_res.", cluster_res)

DimPlot(
  integrated,
  reduction = "umap",
  group.by = c(cluster_name, "SingleR.hpca_main"),
  label = TRUE,
  repel = TRUE,
  label.box = TRUE
)
```

### Visualise the automatic cluster assignments

We can also plot and tabulate the number of cells assigned to each cell type for each cluster:

```{r print_cluster_assignments, message = FALSE, warning = FALSE, fig.height = 6}
integrated@meta.data %>%
  dplyr::select(any_of(cluster_name), SingleR.hpca_main) %>%
  group_by(.data[[cluster_name]], SingleR.hpca_main) %>%
  summarise(n_cells = n()) %>%
  group_by(.data[[cluster_name]]) %>%
  mutate(total_cells = sum(n_cells)) %>%
  purrr::set_names(c("cluster", "cell_type", "n_cells", "total_cells")) %>%
  mutate(prop_cells = n_cells / total_cells) %>%
  ggplot(aes(x = cell_type, y = cluster, colour = prop_cells, size = prop_cells)) +
    geom_point() +
    theme(
      axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
      panel.background = element_rect(fill = "white", colour = NA),
      panel.border = element_rect(
        fill = NA,
        colour = "grey70",
        linewidth = rel(1)
    ),
    panel.grid = element_line(colour = "grey87"),
    panel.grid.major = element_line(linewidth = rel(0.5)),
    panel.grid.minor = element_line(linewidth = rel(0.25)),
    axis.ticks = element_line(colour = "grey70", linewidth = rel(0.5)),
    strip.background = element_rect(
      fill = "grey70",
      colour = NA
    ),
    complete = TRUE
  ) +
  ggtitle("Cluster cell type composition (HPCA main)") +
  labs(x = "Cell Type", y = "Cluster", colour = "Prop. Cells", size = "Prop. Cells") +
  scale_colour_gradientn(breaks = c(0, 0.5, 1), colours = c("lightblue","beige","red"))


integrated@meta.data %>%
  dplyr::select(any_of(cluster_name), SingleR.hpca_main) %>%
  group_by(.data[[cluster_name]], SingleR.hpca_main) %>%
  summarise(n_cells = n()) %>%
  pivot_wider(names_from = cluster_name, values_from = n_cells) %>%
  datatable()


integrated@meta.data %>%
  dplyr::select(any_of(cluster_name), SingleR.hpca_fine) %>%
  group_by(.data[[cluster_name]], SingleR.hpca_fine) %>%
  summarise(n_cells = n()) %>%
  pivot_wider(names_from = cluster_name, values_from = n_cells) %>%
  datatable()
```

## Annotation 3: Manual cell type annotation with curated gene lists

If you have a list of cell types with curated cell markers available, you can use the template CSV file `inputs/custom_markers.csv` to define these gene programs. The CSV file contains two columns with the headers `cell_type` and `gene_ids`. The `gene_ids` column should contain all of the marker genes associated with the cell type, separated by semicolons. For example:

```
cell_type,gene_ids
T_NK,CD3D;CD3E;CD3G;CD247
B_cells,CD79A;CD79B;IGHM;IGHD
Myeloid,LYZ;S100A8;S100A9;CD14
Fibroblasts,COL1A1;COL1A2;COL3A1;ACTA2
Endothelial,PECAM1;CD34;CDH5;VWF
```

**NOTE:** Manual cell type annotation requires **at the very least** two distinct cell types to be defined, as it works by measuring each cell's alignment with each custom gene program and picking the greatest score. We **highly recommend** that you define as many cell types as reasonably possible for your samples.

---

**❱❱❱ ACTION ❰❰❰**

- If you want to run custom cell type assignment:
    1. Open the file `inputs/custom_markers.csv`
    2. Leave the first line unchanged: `cell_type,gene_ids`
    3. For every cell type you want to detect, add the cell type name (no spaces) in the first column (`cell_type`) and a semicolon-delimited list of marker genes in the second column (`gene_ids`).
        - **NOTE:**  We cannot assign a cell type unless you define at least 2 cell types, but it is highly recommended that you define as many cell types as reasonably possible. You should include all cell types you expect to observe in your dataset.
    4. Continue through this section of code
- If you don't want to run custom cell type assignment, delete or rename the template file (`inputs/custom_markers.csv`), then skip ahead to the section "[Assign cell type identities to clusters](#assign_cluster_types)".

---

```{r define_custom_programs}
# Determine whether or not to run custom cell type assignemnt depending on whether the custom_markers.csv file exists
custom_marker_genes_file <- here("inputs/custom_markers.csv")
run_custom_annotation <- file.exists(custom_marker_genes_file)
```

Read in the custom marker genes.

```{r read_custom_programs, echo = run_custom_annotation, eval = run_custom_annotation}
custom_marker_genes <- read_csv(custom_marker_genes_file, show_col_types = FALSE)

custom_programs <- custom_marker_genes$cell_type
custom_marker_genes <- strsplit(custom_marker_genes$gene_ids, ";")
names(custom_marker_genes) <- custom_programs

if(!(length(custom_marker_genes) > 1)) {
  warning("You have only defined one custom gene program. Custom cell type annotation requires at least two gene programs. Please define additional cell types, or skip ahead to the section 'Assign cell type identities to clusters'.")
  run_custom_annotation <- FALSE
}

print("Custom marker genes:")
print(custom_marker_genes)

# If using Ensembl IDs, convert the marker gene list to them as well
if(run_custom_annotation && using_ens_ids) {
  custom_marker_genes <- lapply(custom_marker_genes, function(g) {
    gene2ens <- AnnotationDbi::mapIds(
      EnsDb.Hsapiens.v86::EnsDb.Hsapiens.v86,
      keys = g,
      column = "GENEID",
      keytype = "SYMBOL"
    )
    keep <- !is.na(gene2ens)
    gene2ens <- gene2ens[keep]
    return(gene2ens)
  })

  print("Custom marker genes - Ensembl IDs:")
  print(custom_marker_genes)
}
```

Run Seurat's `AddModuleScore` method. This scores each cell for each gene program defined in the custom marker gene file by calculating the average expression of all the genes for each program and comparing that to the average expression of genes not in those gene sets. In this way, higher scores calculated by this method represent stronger alignment between a cell and that gene program.

```{r score_cells, echo = run_custom_annotation, eval = FALSE}
integrated <- AddModuleScore(
  integrated,
  features = custom_marker_genes,
  name = names(custom_marker_genes)
)

# Remove numeric suffix from cluster names
numeric_suffix_clusters <- paste0(names(custom_marker_genes), 1:length(custom_marker_genes))
stopifnot(all(numeric_suffix_clusters %in% colnames(integrated@meta.data)))  # Sanity check
colnames(integrated@meta.data)[colnames(integrated@meta.data) %in% numeric_suffix_clusters] <- names(custom_marker_genes)
stopifnot(all(names(custom_marker_genes) %in% colnames(integrated@meta.data)))  # Sanity check

gc()  # Clean up memory

saveRDS(integrated, here("tmp_outputs", "04.annotation", "integrated.modulescore.Rds"))
```

We can tabulate a summary of the scores for each custom gene program within each cluster of cells:

```{r summarise_scores, echo = run_custom_annotation, eval = run_custom_annotation}
if (read_tmp_rds_files) {
  integrated <- readRDS(here("tmp_outputs", "04.annotation", "integrated.modulescore.Rds"))
}

scores <- integrated@meta.data %>%
  dplyr::select(Cluster = all_of(cluster_name), all_of(custom_programs)) %>%
  pivot_longer(
    cols = all_of(custom_programs),
    names_to = "cell_type",
    values_to = "score"
  ) %>%
  group_by(Cluster, cell_type) %>%
  summarise(
    avg_score    = mean(score, na.rm = TRUE),
    median_score = median(score, na.rm = TRUE),
    sd_score     = sd(score, na.rm = TRUE),
    n_cells      = n(),
    .groups      = "drop"
  ) %>%
  mutate(across(where(is.numeric), ~ round(.x, 4)))

datatable(scores, colnames = c("Cluster", "Cell Type", "Avg. Score", "Median Score", "SD Score", "# Cells"))
```

We can also plot how strongly each cluster aligns to each gene program. The following code creates a heatmap, where each row represents one of the custom gene expression programs and each column represents a cluster. The more "red" a given tile is, the more highly expressed that set of marker genes is in that cluster compared to other genes.

```{r plot_cluster_scores, echo = run_custom_annotation, eval = run_custom_annotation}
ggplot(scores, aes(x = Cluster, y = cell_type, fill = median_score)) +
  geom_tile() +
  theme_light() +
  ggtitle("Cluster Scores") +
  scale_fill_gradientn(colours = c("lightblue","beige","red"))
```

### Predict the cell type using the maximum score

The following code will determine the cell type that each cell most strongly aligns with. Where more than 2 cell types have been provided, we will calculate the median absolute difference (MAD) of the scores, as well as the difference between the top two scores. Where the top two scores are within 1 MAD of each other, we will call the cell "Ambiguous"; otherwise, the cell type of the highest score will be assigned.

---

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

By default, we require the top-scoring cell type to be at least 1 MAD higher than the second-top-scoring cell type; otherwise, we will call the cell 'Ambiguous'.

You can optionally change this MAD threshold to be more liberal or strict with your cell type classifications. Higher MAD thresholds will be more strict and will result in more ambiguous cell types, which is useful if you wish to be more cautious with the cell type assignments and want to manually curate the ambiguous calls. On the other hand, lower MAD thresholds are more liberal and will result in fewer ambiguous calls, although the possibility of an incorrect cell type assignment is increased.

If you want to change the MAD threshold for calling unambiguous cells, change the following value:

---

```{r set_mad_threshold, echo = run_custom_annotation, eval = run_custom_annotation}
mad_threshold <- 1.0
```

```{r get_max_scores_per_cluster, echo = run_custom_annotation, eval = FALSE}
argmax_scores <- apply(integrated@meta.data[custom_programs], 1, which.max) %>% unlist
integrated$cell_type.max_score <- custom_programs[argmax_scores]

# Also calculate how much the top two scores differ
# Can only calculate when more than 3 cell programs are present
# Determines how many MADs the top 2 scores differ by
# If they differ by less than 1 MAD, mark as ambiguous
if (length(custom_programs) > 3) {
  s <- integrated@meta.data[custom_programs]
  integrated$cell_type.top_2_score_mad_diff <- apply(s, 1, function(x) {
    xv <- unlist(x)
    mad_x <- stats::mad(xv, constant = 1)
    top_2_x <- xv[order(xv, decreasing = TRUE)][1:2]
    return((top_2_x[1] - top_2_x[2]) / mad_x)
  }) %>% unlist
  
  # If the top two scores are too close, mark the cell as ambiguous
  integrated@meta.data <- integrated@meta.data %>% mutate(
    cell_type.mark_ambiguous = case_when(
      cell_type.top_2_score_mad_diff >= mad_threshold ~ cell_type.max_score,
      .default = "Ambiguous"
    )
  )
  
  available_annots <- c(available_annots, "cell_type.max_score", "cell_type.mark_ambiguous")
} else {
  available_annots <- c(available_annots, "cell_type.max_score")
}

saveRDS(integrated, here("tmp_outputs", "04.annotation", "integrated.maxscores.Rds"))
saveRDS(available_annots, here("tmp_outputs", "04.annotation", "available_annots.maxscores.Rds"))
```

We will now plot our UMAP, grouped by our clustering, HPCA main cell type annotations, and maximum module score annotations. If we have more than 2 custom gene programs defined, we will also plot the UMAP grouped by the maximum module score annotations that include the "Ambiguous" category.

```{r plot_custom_cell_types, echo = run_custom_annotation, eval = run_custom_annotation, message = FALSE, warning = FALSE, fig.width = 16, fig.height = 16}
if (read_tmp_rds_files) {
  integrated <- readRDS(here("tmp_outputs", "04.annotation", "integrated.maxscores.Rds"))
  available_annots <- readRDS(here("tmp_outputs", "04.annotation", "available_annots.maxscores.Rds"))
}

dim_groups <- c(cluster_name, "SingleR.hpca_main", "cell_type.max_score")
if (length(custom_programs) > 3) {
  dim_groups <- c(dim_groups, "cell_type.mark_ambiguous")
}
DimPlot(
  integrated,
  reduction = "umap",
  group.by = dim_groups,
  label = TRUE,
  repel = TRUE,
  label.box = TRUE
)
rm(dim_groups)
```

### Visualise the manual cluster assignments

We can again plot and tabulate the number of cells assigned to each cell type for each cluster:

```{r print_custom_cluster_assignments, message = FALSE, warning = FALSE, fig.height = 6}
integrated@meta.data %>%
  dplyr::select(any_of(cluster_name), cell_type.max_score) %>%
  group_by(.data[[cluster_name]], cell_type.max_score) %>%
  summarise(n_cells = n()) %>%
  group_by(.data[[cluster_name]]) %>%
  mutate(total_cells = sum(n_cells)) %>%
  purrr::set_names(c("cluster", "cell_type", "n_cells", "total_cells")) %>%
  mutate(prop_cells = n_cells / total_cells) %>%
  ggplot(aes(x = cell_type, y = cluster, colour = prop_cells, size = prop_cells)) +
    geom_point() +
    theme(
      axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
      panel.background = element_rect(fill = "white", colour = NA),
      panel.border = element_rect(
        fill = NA,
        colour = "grey70",
        linewidth = rel(1)
    ),
    panel.grid = element_line(colour = "grey87"),
    panel.grid.major = element_line(linewidth = rel(0.5)),
    panel.grid.minor = element_line(linewidth = rel(0.25)),
    axis.ticks = element_line(colour = "grey70", linewidth = rel(0.5)),
    strip.background = element_rect(
      fill = "grey70",
      colour = NA
    ),
    complete = TRUE
  ) +
  ggtitle("Cluster cell type composition (based on maximum score)") +
  labs(x = "Cell Type", y = "Cluster", colour = "Prop. Cells", size = "Prop. Cells") +
  scale_colour_gradientn(breaks = c(0, 0.5, 1), colours = c("lightblue","beige","red"))


integrated@meta.data %>%
  dplyr::select(any_of(cluster_name), cell_type.max_score) %>%
  group_by(.data[[cluster_name]], cell_type.max_score) %>%
  summarise(n_cells = n()) %>%
  pivot_wider(names_from = cluster_name, values_from = n_cells) %>%
  datatable()


if (length(custom_programs) > 3) {
  p <- integrated@meta.data %>%
    dplyr::select(any_of(cluster_name), cell_type.mark_ambiguous) %>%
    group_by(.data[[cluster_name]], cell_type.mark_ambiguous) %>%
    summarise(n_cells = n()) %>%
    group_by(.data[[cluster_name]]) %>%
    mutate(total_cells = sum(n_cells)) %>%
    purrr::set_names(c("cluster", "cell_type", "n_cells", "total_cells")) %>%
    mutate(prop_cells = n_cells / total_cells) %>%
    ggplot(aes(x = cell_type, y = cluster, colour = prop_cells, size = prop_cells)) +
      geom_point() +
      theme(
        axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
        panel.background = element_rect(fill = "white", colour = NA),
        panel.border = element_rect(
          fill = NA,
          colour = "grey70",
          linewidth = rel(1)
      ),
      panel.grid = element_line(colour = "grey87"),
      panel.grid.major = element_line(linewidth = rel(0.5)),
      panel.grid.minor = element_line(linewidth = rel(0.25)),
      axis.ticks = element_line(colour = "grey70", linewidth = rel(0.5)),
      strip.background = element_rect(
        fill = "grey70",
        colour = NA
      ),
      complete = TRUE
    ) +
    ggtitle("Cluster cell type composition (including ambiguous cell types)") +
    labs(x = "Cell Type", y = "Cluster", colour = "Prop. Cells", size = "Prop. Cells") +
    scale_colour_gradientn(breaks = c(0, 0.5, 1), colours = c("lightblue","beige","red"))
  print(p)


  integrated@meta.data %>%
    dplyr::select(any_of(cluster_name), cell_type.mark_ambiguous) %>%
    group_by(.data[[cluster_name]], cell_type.mark_ambiguous) %>%
    summarise(n_cells = n()) %>%
    pivot_wider(names_from = cluster_name, values_from = n_cells) %>%
    datatable()
}
```

Additionally, it is helpful to plot the module scores themselves for each cell type definition. These plots demonstrate how some cell clusters will align very strongly with one cell type, while others may more weakly align with two or more cell types.

```{r plot_scores_per_cell_type, echo = run_custom_annotation, eval = run_custom_annotation, message = FALSE, warning = FALSE}
for (ct in custom_programs) {
  p <- FeaturePlot(
    integrated,
    reduction = "umap",
    features = c(ct)
  )
  
  print(
    p +
    scale_colour_gradient2(
      low = 'royalblue',
      mid = 'lightgrey',
      high = 'indianred2',
      midpoint = 0
    )
  )
}
```

## Assign cell type identities to clusters {#assign_cluster_types}

So far, our cell types have been assigned to each individual cell. Going forward, we will want to annotate each cluster with a cell type. First, you should define the cell type annotation that you want to use. You currently have the following options to choose from:

```{r print_available_cell_type_annotations, echo = FALSE}
cat(paste(available_annots, collapse = ", "))
```


---

**❱❱❱ ACTION ❰❰❰**

Define your chosen annotation in the following block:

---

```{r define_annotation_for_clusters}
cluster_annotation <- "SingleR.hpca_main"
```

The next code block will automatically determined cluster cell type assignments. This process involves simply looking at the proportions of each assigned cell type in the cluster and uses the most prevalent cell type as the representative cell type for the cluster. We also define a proportion threshold (`cell_type_proportion_threshold`) at the top of the code block; if the most prevalent cell type does not exceed this proportion, we will consider the cluster "Ambiguous". We have set this value to 67% by default, so 2/3 cells must be assigned to the same cell type for the cluster to recieve the same cell type annotation. Change this value depending on how strict or lenient you wish to be, although values less than 0.67 will likely result in clusters that are too heterogeneous to give meaningful results in downstream analyses.

---

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

By default, to assign a cluster with a definitive cell type, we require 2/3 of the cells in a cluster to be of the same cell type.

You can optionally change this threshold to be more liberal or strict with the cluster cell type assignments. Higher values will result in more ambiguous clusters, but the ones that are assigned a definitive cell type will be more homogeneous. This is useful if you wish to be more cautious with the assignments and want to manually review the ambiguous clusters. On the other hand, lower thresholds will result in fewer ambiguous calls, but this may hide considerable heterogeneity within the clusters that could confound downstream analyses.

If you want to change the threshold for calling unambiguous clusters, change the following value:

---

```{r set_ambiguous_cluster_threshold}
# TODO: Define the cell type proportion threshold to use to define a cluster's overall cell type
#       Anything below this value will be classed as ambiguous
cell_type_proportion_threshold <- 0.67
```

```{r automatic_cluster_assignment, message = FALSE, warning = FALSE, fig.height = 6}
# Summarise cluster cell types
cluster_annotations <- integrated@meta.data[c(cluster_name, cluster_annotation)]
colnames(cluster_annotations) <- c("cluster", "cell_type")
cluster_annotations$cell_id <- rownames(cluster_annotations)
cluster_sizes <- cluster_annotations %>% as_tibble %>%
  group_by(cluster) %>%
  summarise(cluster_size = n())
cluster_annotations_summary <- cluster_annotations %>% as_tibble %>%
  group_by(cluster, cell_type) %>%
  summarise(n_cells = n()) %>%
  left_join(cluster_sizes, by = "cluster") %>%
  mutate(prop_cells = n_cells / cluster_size)

# Create a dot plot of the cell type annotations per cluster
cluster_annotations_summary %>%
  arrange(prop_cells) %>%
  ggplot(aes(x = cell_type, y = cluster, colour = prop_cells, size = prop_cells)) +
  geom_point() +
  theme(
    axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
    panel.background = element_rect(fill = "white", colour = NA),
    panel.border = element_rect(
      fill = NA, 
      colour = "grey70",
      linewidth = rel(1)
    ),
    panel.grid = element_line(colour = "grey87"),
    panel.grid.major = element_line(linewidth = rel(0.5)),
    panel.grid.minor = element_line(linewidth = rel(0.25)),
    axis.ticks = element_line(colour = "grey70", linewidth = rel(0.5)),
    strip.background = element_rect(
      fill = "grey70",
      colour = NA
    ),
    complete = TRUE
  ) +
  ggtitle("Cluster cell type composition") +
  labs(x = "Cell Type", y = "Cluster", colour = "Prop. Cells", size = "Prop. Cells") +
  scale_colour_gradientn(breaks = c(0, 0.5, 1), colours = c("lightblue","beige","red"))

# Create a datatable of the cell type annotations per cluster
cluster_annotations_summary %>%
  mutate(across(where(is.numeric), ~ round(.x, 4))) %>%
  datatable(colnames = c("Cluster", "Cell Type", "# Cells", "Cluster Size", "Proportion of Cells"))

# Determine the consensus annotations
cluster_annotations_consensus <- cluster_annotations_summary %>%
  group_by(cluster) %>%
  dplyr::filter(prop_cells == max(prop_cells)) %>%
  mutate(
    cell_type_consensus = case_when(
      prop_cells >= cell_type_proportion_threshold ~ cell_type,
      .default = "Ambiguous"
    )
  ) %>%
  dplyr::select(cluster, cell_type_consensus) %>%
  ungroup() %>%
  dplyr::rename(cell_type = cell_type_consensus)

datatable(cluster_annotations_consensus, colnames = c("Cluster", "Cell Type"))
```

The next code block will generate an input file called `inputs/cluster_annotation.csv containing the automatically-determined cluster assignments. If you wish to manually annotate your cell clusters, you can modify this file before moving on. Refer back to the plots generated in the appropriate sections above to help guide your decision on what cell types to assign to each cluster.

```{r write_cluster_annotations_consensus, eval = FALSE}
# Write the consensus annotations to file
cluster_annotations_consensus %>% write_csv(here("inputs/cluster_annotation.csv"))
```

---

**❱❱❱ ACTION ❰❰❰**

1. Open the file `inputs/cluster_annotation.csv`
2. Review the cluster assignments
3. If you wish to change a cluster's assignment, replace the text in the second column with your desired cell type identity.
    - **CAUTION:** You should only change these assignments if there is good reason to do so, e.g. a cluster that was called "ambiguous" actually represents a cell sub-type

---

Now we will read in the new cluster annotations from `inputs/cluster_annotation.csv` and assign them to the cells. Make sure you have reviewed this file and made any manual changes you see fit.

```{r assign_cluster_annotations, eval = FALSE}
cluster_annotations_consensus <- read_csv(here("inputs/cluster_annotation.csv"), show_col_types = FALSE)

integrated@meta.data$cluster_annotation <- cluster_annotations_consensus$cell_type[match(integrated@meta.data[[cluster_name]], cluster_annotations_consensus$cluster)]

saveRDS(integrated, here("tmp_outputs", "04.annotation", "integrated.clusterannots.Rds"))
```

```{r plot_cluster_annotations}
if (read_tmp_rds_files) {
  integrated <- readRDS(here("tmp_outputs", "04.annotation", "integrated.clusterannots.Rds"))
}

DimPlot(
  integrated,
  reduction = "umap",
  group.by = "cluster_annotation",
  label = TRUE,
  repel = TRUE,
  label.box = TRUE
)

available_annots <- c(available_annots, "cluster_annotation")
```

You now have the following cell annotations to use for downstream analyses:

```{r print_cell_type_annotations, echo = FALSE}
original_samplesheet <- here("inputs/samplesheet.csv")
original_samplesheet <- read_csv(original_samplesheet, show_col_types = FALSE)
additional_annots <- colnames(original_samplesheet)
additional_annots <- additional_annots[! additional_annots %in% c("sample", "rds_path")]
available_annots <- unique(c(additional_annots, available_annots))
rm(additional_annots)

cat(paste(available_annots, collapse = ", "))
```

## Write data to file

Finally, we write our final Seurat dataset and our available annotations to file:

```{r save_data, eval = FALSE}
# Save integrated Seurat dataset
dir.create(here("outputs"))
SaveSeuratRds(integrated, here("outputs/integrated_dataset.filtered_clustered.doublets_detected.integrated.annotated.Rds"))

# Also save available annotations
sink(here("outputs/available_annotations.txt"))
cat(paste(available_annots, collapse = "\n"))
sink()
```

## Package version summary

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