Adding option to increase clusters included in damaged cell detection

R/limiric.R

@@ -10,7 +10,6 @@
 #' @param filtered_path Directory of filtered alignment output
 #' @param seurat_input 'Seurat' object to be used as input over raw files. Default NULL
 #' @param min_cells In how many cells should a gene be expressed to be kept
-#' @param resolution Numeric between 0 and 1.6 describing cluster division. Default 1
 #' @param soupx Perform ambient RNA correction, if TRUE raw_path must be given. Default is FALSE
 #' @param raw_path Directory of unfiltered alignment output
 #' @param droplet_qc Verify output with droplet_qc, if TRUE velocyto_path must be given. Default is FALSE
@@ -20,6 +19,8 @@
 #' @param filter_output Should output contain no damaged cells. Default is TRUE
 #' @param output_path Directory where 'limiric' output should be generated
 #' @param organism "Hsap" if human sample or "Mmus" if mouse sample
+#' @param resolution Numeric between 0 and 1.6 describing cluster division. Default 1
+#' @param cluster_ranks Numeric describing the number of top ranking clusters to be included as damaged cells. Default 1.
 #' @param sample_list Input multiple samples in list. Default is FALSE
 #'
 #' @return (list) Output storing the final 'Seurat' object
@@ -66,7 +67,6 @@ limiric <- function(
   filtered_path = NULL,
   seurat_input  = NULL,
   min_cells     = NULL,
-  resolution    = NULL,
   soupx         = NULL,
   raw_path      = NULL,
   droplet_qc    = NULL,
@@ -76,6 +76,8 @@ limiric <- function(
   filter_output = NULL,
   output_path   = NULL,
   organism      = NULL,
+  resolution    = NULL,
+  cluster_ranks = NULL,
   sample_list   = NULL
 
 ){
@@ -90,6 +92,7 @@ limiric <- function(
     if (is.null(seurat_input))  {seurat_input = NULL}
     if (is.null(min_cells))     {min_cells = 0}
     if (is.null(resolution))    {resolution = 1}
+    if (is.null(cluster_ranks)) {cluster_ranks = 1}
     if (is.null(soupx))         {soupx = FALSE}
     if (is.null(droplet_qc))    {droplet_qc = FALSE}
     if (is.null(filter_rbc))    {filter_rbc = TRUE}
@@ -105,6 +108,7 @@ limiric <- function(
       seurat_input  = seurat_input,
       min_cells     = min_cells,
       resolution    = resolution,
+      cluster_ranks = cluster_ranks,
       soupx         = soupx,
       raw_path      = raw_path,
       droplet_qc    = droplet_qc,
@@ -163,6 +167,7 @@ limiric <- function(
       seurat_input  <- sample$seurat_input
       min_cells     <- sample$min_cells
       resolution    <- sample$resolution
+      cluster_ranks <- sample$cluster_ranks
       soupx         <- sample$soupx
       raw_path      <- sample$raw_path
       droplet_qc    <- sample$droplet_qc
@@ -177,6 +182,7 @@ limiric <- function(
       if (is.null(seurat_input))  {seurat_input = NULL}
       if (is.null(min_cells))     {min_cells = 0}
       if (is.null(resolution))    {resolution = 1}
+      if (is.null(cluster_ranks)) {cluster_ranks = 1}
       if (is.null(soupx))         {soupx = FALSE}
       if (is.null(droplet_qc))    {droplet_qc = FALSE}
       if (is.null(filter_rbc))    {filter_rbc = TRUE}
@@ -196,6 +202,7 @@ limiric <- function(
           seurat_input  = seurat_input,
           min_cells     = min_cells,
           resolution    = resolution,
+          cluster_ranks = cluster_ranks,
           soupx         = soupx,
           raw_path      = raw_path,
           droplet_qc    = droplet_qc,

R/limiric_calculation.R

@@ -9,6 +9,7 @@
 #' @param organism A string representing the organism ("Hsap", "Mmus").
 #' @param Seurat A 'Seurat' object containing the single-cell RNA-seq data.
 #' @param resolution Numeric between 0 and 1.6 describing cluster division. Default is 1.
+#' @param cluster_ranks Numeric describing the number of top ranking clusters to be included as damaged cells. Default 1.
 #' @param annotations A data frame containing gene annotations.
 #' @param initial_cells An integer representing the initial number of cells.
 #' @param project_name A string representing the name of the project, used for plot titles.
@@ -54,6 +55,7 @@ utils::globalVariables(c("mt_plot", "complexity_plot", "rb_plot",
 limiric_calculation <- function(organism,
                                 Seurat,
                                 resolution,
+                                cluster_ranks,
                                 annotations,
                                 initial_cells,
                                 project_name,
@@ -168,12 +170,12 @@ limiric_calculation <- function(organism,
       combined_rank = rank_mt + rank_rb
     ) %>%
     arrange(combined_rank) %>%
-    slice(1) %>%
+    slice(1:cluster_ranks) %>%
     pull(seurat_clusters) %>%
     as.character()
 
   # Label all cells belonging to this cluster as "damaged"
-  limiric$seurat_clusters <- ifelse(limiric$seurat_clusters == best_cluster, 'damaged', limiric$seurat_clusters)
+  limiric$seurat_clusters <- ifelse(limiric$seurat_clusters %in% best_cluster, 'damaged', limiric$seurat_clusters)
 
   # Add Cell QC meta data to object
   limiric$MtRb <- ifelse(limiric$seurat_clusters == "damaged", "damaged", "cell")
@@ -189,30 +191,30 @@ limiric_calculation <- function(organism,
 
   # Visualise  clusters coloured by QC metrics ------------------------------------
 
-  mt_plot <- FeaturePlot(limiric, features = c("mt.percent"), cols = c("#E1E1E1", "#5372B4"), pt.size = 1) +
+  mt_plot <- FeaturePlot(limiric, features = c("mt.percent"), cols = c("#E1E1E1", "#6765ED"), pt.size = 1) +
     NoAxes() + labs(caption = "Mitochondrial gene expression") +
     theme(
       plot.title   = element_blank(),
       panel.border = element_rect(colour = "black", fill=NA, linewidth =1),
       plot.caption = element_text(hjust = 0.5, size = 16))
 
-  rb_plot <- FeaturePlot(limiric, features = c("rb.percent"), cols = c("#E1E1E1", "#5372B4"), pt.size = 1) +
+  rb_plot <- FeaturePlot(limiric, features = c("rb.percent"), cols = c("#E1E1E1", "#6765ED"), pt.size = 1) +
     NoAxes() + labs(caption = "Ribosomal gene expression") +
     theme(
       plot.title = element_blank(),
       panel.border = element_rect(colour = "black", fill=NA, linewidth =1),
       plot.caption = element_text(hjust = 0.5, size = 16))
 
   cluster_plot <- DimPlot(
-    limiric, pt.size = 1, group.by = "MtRb", cols = c("cell" = "grey", "damaged" = "#5372B4")) +
+    limiric, pt.size = 1, group.by = "MtRb", cols = c("cell" = "grey", "damaged" = "#6765ED")) +
     labs(caption = expression("Damaged cells identified by " * italic("limiric"))) + NoAxes() +
     theme(
       plot.title = element_blank(),
       plot.caption = element_text(hjust = 0.5, size = 16),
       panel.border = element_rect(colour = "black", fill=NA, linewidth =1)
     )
 
-  complexity_plot <- FeaturePlot(limiric, features = c("complexity"), cols = c("#E1E1E1", "#5372B4"), pt.size = 1) +
+  complexity_plot <- FeaturePlot(limiric, features = c("complexity"), cols = c("#E1E1E1", "#6765ED"), pt.size = 1) +
     NoAxes() + labs(caption = "Complexity score") +
     theme(
       plot.title = element_blank(),

R/limiric_core.R

@@ -8,7 +8,6 @@
 #' @param filtered_path Directory of filtered alignment output
 #' @param seurat_input Seurat object to be used as input over raw files. Default NULL
 #' @param min_cells In how many cells should a gene be expressed to be kept
-#' @param resolution Numeric between 0 and 1.6 describing cluster division. Default 1
 #' @param soupx Perform ambient RNA correction, if TRUE raw_path must be given. Default is FALSE
 #' @param raw_path Directory of unfiltered alignment output
 #' @param droplet_qc Verify output with droplet_qc, if TRUE velocyto_path must be given. Default is FALSE
@@ -17,6 +16,8 @@
 #' @param isolate_cd45 Discard non-immune cells. Default is FALSE
 #' @param filter_output Should output contain no damaged cells. Default is TRUE
 #' @param output_path Directory where limiric output cen be generated
+#' @param resolution Numeric between 0 and 1.6 describing cluster division. Default 1
+#' @param cluster_ranks Numeric describing the number of top ranking clusters to be included as damaged cells. Default 1.
 #' @param organism "Hsap" if human sample or "Mmus" if mouse sample
 #'
 #' @return (list) Output list storing the final Seurat object
@@ -63,7 +64,6 @@ limiric_core <- function(
     filtered_path,
     seurat_input  = NULL,
     min_cells     = 0,
-    resolution    = 1,
     soupx         = FALSE,
     raw_path      = NULL,
     droplet_qc    = FALSE,
@@ -72,6 +72,8 @@ limiric_core <- function(
     isolate_cd45  = FALSE,
     filter_output = TRUE,
     output_path   = "./",
+    resolution    = 1,
+    cluster_ranks = 1,
     organism      = "Hsap"
 ){
   # Receive & prepare input ------------------------------------
@@ -339,6 +341,7 @@ limiric_core <- function(
   limiric_output <- limiric_calculation(organism = organism,
                                   Seurat = Seurat,
                                   resolution = resolution,
+                                  cluster_ranks = cluster_ranks,
                                   annotations = annotations,
                                   initial_cells = initial_cells,
                                   project_name = project_name,