notebook_cell_clustering.Rmd

---
title: "Module - Clustering"
author: "BTC scRNA Pipeline"
date: "`r Sys.time()`"
output:
  html_document:
    toc: TRUE
    toc_float: TRUE
    toc_depth: 4
    code_folding: show
    self_contained: TRUE
    theme: united
    df_print: paged
params:
  project_name: 'Test'
  project_object: './data/Test_normalize_object.RDS'
  input_features_plot: 'LYZ;CCL5;IL32;PTPRCAP;FCGR3A;PF4;PTPRC'
  input_group_plot: 'source_name;Sort'
  input_integration_dimension: 'auto'
  input_cluster_step: 'main'
  thr_quantile: 'q01'
  thr_npc: 'auto'
  thr_resolution: 0.25
  thr_proportion: 0.25
  n_threads: 8
  n_memory: 8
  workdir: !r here::here()
  timestamp: !r Sys.Date()
  auto_save: TRUE
---

------------------------------------------------------------------------

# Project Name: `r params$project_name`

```{r setup, include = FALSE, message = FALSE, warning = FALSE, echo = FALSE}

# Project parameters 
project_name <- params$project_name
project_object <- params$project_object

# Inputs and thresholds
input_integration_dimension <- params$input_integration_dimension
input_features_plot <- strsplit(
  params$input_features_plot, split = ';')[[1]]
input_group_plot <- strsplit(
  params$input_group_plot, split = ';')[[1]]
input_cluster_step <- params$input_cluster_step

thr_npc <- params$thr_npc
thr_resolution <- params$thr_resolution
thr_proportion <- params$thr_proportion
thr_quantile <- params$thr_quantile

# Dataflow/Computational parameters
n_threads <- params$n_threads
n_memory <- params$n_memory

# Output parameters
work_directory <- params$workdir
timestamp <- params$timestamp
auto_save <- params$auto_save

# Rmarkdown options
knitr::opts_chunk$set(
	echo = TRUE,
	error = TRUE,
	fig.align = "center",
	fig.height = 8,
	fig.width = 12,
	message = FALSE,
	warning = FALSE,
	dpi = 300
)
knitr::opts_knit$set(
  root.dir = work_directory
  )

```

## Loading library

```{r library}

library(readr)
library(dplyr)
library(ggplot2)
library(patchwork)
library(Seurat)
library(ROGUE)

```

## Folder structure

```{r project_directory}

if(!dir.exists(work_directory)) {
  dir.create(work_directory, recursive = T)
}

for(sub_directory in c('data', 'figures', 'figures/clustering')) {
  dir.create(paste0(work_directory, '/', sub_directory))
}

# Cirro-related edition
file.create('figures/clustering/.empty')

```

## Loading Seurat object

```{r loading_rds}

# Reading seurat object
seurat_object <- readRDS(file = project_object)

# Cell number
ncell <- dim(seurat_object)[2]

# Knowledge-based parameters
if(thr_npc == 'auto') {
  thr_npc <- ifelse(ncell < 50000, 25, 50)
}

# Setting random seed
random_seed <- 2203
set.seed(random_seed)

```

## Unifying sample layers

```{r seurat_rejoin, message = FALSE, warning = FALSE, echo = FALSE}

seurat_object <- JoinLayers(seurat_object)

```

## Input validation

```{r basic_input_validation}

# Extracting gene names and metadata columns
metadata_columns <- colnames(seurat_object@meta.data)
features_names <- row.names(seurat_object)

# Checking input variables
input_features_plot <- intersect(input_features_plot, features_names)
input_group_plot <- intersect(input_group_plot, metadata_columns)

# Setting reduction variable
integration_method <- function(
    method = c('pca', 'integrated.mnn', 'integrated.rpca', 'integrated.cca', 'harmony', 'auto')) {

  # Default reduction
  reduction_variable <- 'pca'
  
  # Programmatic choices
  if(method == 'auto') {
    if(file.exists('./batch_method.selected.txt')) {
      reduction_variable <- read_lines(file = "./batch_method.selected.txt")
    } 
  }
  
  if(method %in% names(seurat_object@reductions)) {
    reduction_variable <- method
  } 
  
  return(reduction_variable)
  
}

reduction_variable <- integration_method(method = input_integration_dimension)

```

::: {.alert .alert-warning}
The clustering step will be performed over the **`r reduction_variable`** reduction.
:::

## Finding neighbors and clusters

```{r finding_neighbors}

# Finding clusters
seurat_object <- FindNeighbors(
  object = seurat_object, 
  reduction = reduction_variable, 
  dims = 1:thr_npc,
  graph.name = paste0('snn_npc', thr_npc)
)

seurat_object <- FindClusters(
  object = seurat_object,
  graph.name = paste0('snn_npc', thr_npc),
  resolution = thr_resolution
)

# Running UMAP reduction
seurat_object <- RunUMAP(
  object = seurat_object,
  reduction = reduction_variable, 
  dims = 1:thr_npc,
  seed.use = random_seed
)

```

```{r umap_main_plot, fig.height = 12}

# Plotting UMAP based on Seurat clusters
umap_main <- DimPlot(
  seurat_object, 
  reduction = 'umap',
  raster = FALSE,
  label = TRUE) +
  theme(
    legend.position = 'top',
    legend.justification = 'center'
  )

if(auto_save) {
  
  ggsave(plot = umap_main, filename = paste0('figures/clustering/UMAP_MAIN_', timestamp, '.pdf'))
  
}

umap_main

```

```{r clusters_composition}

seurat_cluster_composition <- seurat_object@meta.data %>%
  group_by(seurat_clusters, patient_id) %>%
  summarise(
    n_cell = n()
  )

seurat_cluster_composition <- seurat_cluster_composition %>%
  group_by(seurat_clusters) %>%
  mutate(
    total_cell = sum(n_cell),
    percentage = round((n_cell / total_cell) * 100, digits = 2)
  )

#

ggplot(seurat_cluster_composition, aes(x = seurat_clusters, y = percentage, fill = patient_id)) +
  geom_col() +
  labs(x = "Seurat Clusters", y = "% of Cells") +
  theme_classic() +
  theme(
    legend.text = element_text(size = 14),
    legend.position = 'top',
    axis.text = element_text(size = 18)
  )

```

```{r printing_cluster_table, message = FALSE, warning = FALSE, echo = FALSE}

seurat_cluster_composition %>%
  distinct(seurat_clusters, total_cell) %>%
  tidyr::pivot_wider(names_from = seurat_clusters, values_from = total_cell)

```

## Clustering performance {.tabset}

```{r clustering_performance}

# Sampling cells per cluster
sample_clusters_members <- seurat_object@meta.data %>%
  mutate(barcode = row.names(seurat_object@meta.data)) %>%
  group_by(seurat_clusters) %>%
  slice_sample(prop = thr_proportion)

sample_clusters_members <- sample_clusters_members %>%
  rename(
    Patient = patient_id,
    clusters = seurat_clusters
  )

# Extracting expression counts
expression_counts <- GetAssayData(seurat_object, layer = "counts")

# Sampling counts
expression_counts <- expression_counts[, sample_clusters_members$barcode] %>%
  as.matrix()
expression_counts <- matr.filter(expression_counts, min.cells = 10, min.genes = 10)

# Computing entropy background
entropy_result <- SE_fun(expression_counts)

# Avg. ROGUE score across clusters
rogue_value <- CalculateRogue(entropy_result, platform = "UMI")

# Calculating ROGUE score
rogue_result <- rogue(expression_counts, 
                      labels = sample_clusters_members$clusters,
                      samples = sample_clusters_members$Patient,
                      platform = "UMI",
                      span = 0.6)

rogue_result

```

::: {.alert .alert-warning}
The ROGUE score for the clustering composition (resolution = **`r thr_resolution`**) is equal to **`r rogue_value`**.
:::

### SEplot

```{r rogue_seplot}

# Highly variable genes
SEplot(entropy_result)

```

### Boxplot

```{r rogue_boxplot}

rogue.boxplot(rogue_result) +
  scale_x_discrete(limits = sort(unique(seurat_object@meta.data$seurat_clusters))) +
  theme(
    axis.text = element_text(size = 16)
  )

```

## Data visualization {.tabset}

### Sample-based UMAP

```{r umap_sample_plot, fig.height = 12, fig.width = 14}

umap_sample <- DimPlot(seurat_object, reduction = "umap", group.by = "sample_id", raster = FALSE)

if(auto_save) {

  ggsave(plot = umap_sample, filename = paste0('figures/clustering/UMAP_SAMPLE.pdf'))

}

umap_sample +
  guides(color = guide_legend(ncol = 2, override.aes = list(size = 5))) +
  theme(
    legend.position = 'bottom',
    legend.justification = 'center'
  )

```

### Batch-based UMAP

```{r umap_batch_plot}

umap_group <- DimPlot(seurat_object, reduction = "umap", group.by = "batch", raster = FALSE)

if(auto_save) {

  ggsave(plot = umap_group, filename = paste0('figures/clustering/UMAP_BATCH.pdf'))

}

umap_group +
  theme(
    legend.position = 'bottom',
    legend.justification = 'center'
  )

```

### Group-based UMAP

```{r umap_group_plot, fig.height = 18, fig.width = 14}

if(length(input_group_plot) > 0) {
  
  umap_group <- DimPlot(seurat_object, reduction = "umap", group.by = input_group_plot, raster = FALSE, combine = FALSE)

  if(auto_save) {
  
    for(idx in 1:length(umap_group)) {
      
      umap_temp <-  umap_group[[idx]] + 
        theme(
          legend.position = 'top',
          legend.justification = 'center'
          )

          ggsave(plot = umap_temp, filename = paste0('figures/clustering/UMAP_GROUPED_', idx, '.pdf'))
  
    }
  }
  
  patchwork::wrap_plots(umap_group, ncol = 1)
  
}


```

## Saving Seurat object

```{r object_dump}

if(auto_save) {
  
  saveRDS(seurat_object, file = paste0('./data/', project_name, '_', input_cluster_step, '_cluster_object.RDS'))
  
}

```

------------------------------------------------------------------------

## Parameters log

```{r params_log, message = FALSE, warning = FALSE, echo = FALSE}

print(params)

```

## Session info

```{r session_info, message = FALSE, warning = FALSE, echo = FALSE}

sessionInfo()

```