Set-up

The workflow and explanations bellow are from OSCA

library(SingleCellExperiment)
library(here) #reproducible paths
library(scater) # Plot dimred
library(bluster) # Clustering 
library(clustree) # show relationship clustering
#library(Seurat) # clusternig 2
library(pals) # for palettes with large n #kelly()22, #polychrome()#36, unname(alphabet())
project<- "fire-mice"
if(!file.exists(here("processed", project, "sce_clusters_01.RDS"))){
sce <- readRDS(here("processed", project, "sce_corrected.RDS"))
} else {
  sce <- readRDS(here("processed", project, "sce_clusters_01.RDS"))
}

Motivation

Clustering is an unsupervised learning procedure that is used in scRNA-seq data analysis to empirically define groups of cells with similar expression profiles. It is worth stressing the distinction between clusters and cell types. The former is an empirical construct while the latter is a biological truth (albeit a vaguely defined one). For this reason, questions like “what is the true number of clusters?” are usually meaningless. We can define as many clusters as we like, with whatever algorithm we like - each clustering will represent its own partitioning of the high-dimensional expression space, and is as “real” as any other clustering. It is helpful to realize that clustering, like a microscope, is simply a tool to explore the data. We can zoom in and out by changing the resolution of the clustering parameters, and we can experiment with different clustering algorithms to obtain alternative perspectives of the data.

With Bioconductor

# choose the resolutions to compute
ks <- c(5, 10, 20, 30, 40, 50, 60)
names <- paste0("cluster_k", as.character(ks))

#only run if first time
if (!file.exists(here("processed", project, "clusters_SNN_k5-k60_01.RDS"))) {
  # Compute all the cluster resolutions
  #  and save in a list with an appropriate name for later use
  clusters_list <-
    mapply(function(cluster_name, k) {
      clusterRows(reducedDim(sce, "corrected"), NNGraphParam(k = k))
    },
    names,
    ks,
    USE.NAMES = TRUE,
    SIMPLIFY = FALSE)
  saveRDS(clusters_list,
          here("processed", project, "clusters_SNN_k5-k60_01.RDS"))
} else{
  clusters_list <-
    readRDS(here("processed", project, "clusters_SNN_k5-k60_01.RDS"))
}

# only save the new sce if the object is not already here
if(!file.exists((here("processed", project, "sce_clusters_01.RDS")))){
  
  # Use the list to store the info in the sce and plot result
  idx <- c(1:length(names))
  for (i in idx) {
    k <- ks[i]
    cluster_name <- names[i]
    colData(sce) <- cbind(colData(sce), clusters_list[cluster_name])
    # and plot the result
    plotDim <-
      plotReducedDim(
        sce,
        "TSNE",
        colour_by = cluster_name,
        point_size = 0.5,
        point_alpha = 0.3,
        text_by = cluster_name,
        text_size = 3
      )  + scale_color_hue() + ggtitle(paste("Clustering with k =", k))
    print(plotDim)
  }
  # removed from here, just save after all clustering methods
#  saveRDS(sce, here("processed", project, "sce_clusters_01.RDS"))
  
} else{
  # plot directly from the sce_clusters_01.RDS
  plot_list_func(sce, 
                 col_pattern = "cluster_k", 
                 plot_cols = c( unname(polychrome()), unname(watlington())),
                 label_size = 0.5,
                 reduction="TSNE",
                 
  )
  
}
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.

## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.

## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.

## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.

## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.

## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.

Cluster tree

We use this visualisation to see the relationships between the clusters. The aim is to capture the redistribution of cells from one clustering to another at progressively higher resolutions, providing a convenient depiction of how clusters merge or split apart.

clustree(sce, prefix = "cluster_k", edge_arrow = FALSE)
## Warning: The `add` argument of `group_by()` is deprecated as of dplyr 1.0.0.
## Please use the `.add` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.

With Seurat

#only run if first time
if (!file.exists( here("processed", project, "srt_clusters.RDS"))) {
srt <- as.Seurat(sce)
srt <- FindNeighbors(srt, reduction = "corrected", dims = 1:25)
srt <- FindClusters(srt, 
                    resolution = c(0.01, 0.04, 0.05, 
                                 seq(from = 0.1, to = 1, by = 0.1))
                   )
srt_cluster_names <- grep("originalexp_snn_res", names(srt@meta.data), value = TRUE)
srt_clusters_metadata <- srt[[srt_cluster_names]]

saveRDS(srt_clusters_metadata, here("processed", project, "srt_clusters.RDS"))

plot_list_func(srt, 
               col_pattern="originalexp_snn_res", 
               plot_cols = c( unname(polychrome()), unname(watlington())),
               reduction = "TSNE", 
               label_size = 4 
               )

DimPlot(srt, reduction = "TSNE", group.by = "originalexp_snn_res.1", label = TRUE)

} else{
  srt_clusters_metadata <-
    readRDS( here("processed", project, "srt_clusters.RDS"))
}
clustree(srt_clusters_metadata, prefix = "originalexp_snn_res.", edge_arrow = FALSE)

### save all clusterings in the same object

if (!file.exists( here("processed", project, "sce_clusters_01.RDS"))) {
# check the cell names are the same
identical(row.names(colData(sce)), row.names(srt_clusters_metadata))

# add them
colData(sce) <- cbind(colData(sce), srt_clusters_metadata)

## save
saveRDS(sce, here("processed", project, "sce_clusters_01.RDS"))
}
#plot seurat clustering
plot_list_func(sce, 
               col_pattern="originalexp_snn_res", 
               plot_cols = c( unname(polychrome()), unname(watlington())),
               reduction = "TSNE", 
               label_size = 0.5 
               )
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.

## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.

## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.

## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.

## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.

## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.

## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.

## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.

## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.

## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.

## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.

## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.