Add process/script to classify with singler

modules/cell-type-neuroblastoma-04/README.md

@@ -6,4 +6,5 @@ Links to specific original scripts used in this module:
 
 * `00_convert-nbatlas.R`: <https://github.com/AlexsLemonade/OpenScPCA-analysis/blob/ad3b3c6ac6bcb7154058e4f725250dc56523caa8/analyses/cell-type-neuroblastoma-04/scripts/00_convert-nbatlas.R>
 * `01_train-singler-model.R`: <https://github.com/AlexsLemonade/OpenScPCA-analysis/blob/7c0acea43b6cfc26da75a3a1dbd3558a0d9229ed/analyses/cell-type-neuroblastoma-04/scripts/01_train-singler-model.R>
+* `02_classify-singler.R`: <https://github.com/AlexsLemonade/OpenScPCA-analysis/blob/82547028b5a9555d8cee40f6c1883015c990cc4f/analyses/cell-type-neuroblastoma-04/scripts/02_classify-singler.R>
 * `03a_train-scanvi-model.py`: <https://github.com/AlexsLemonade/OpenScPCA-analysis/blob/82547028b5a9555d8cee40f6c1883015c990cc4f/analyses/cell-type-neuroblastoma-04/scripts/03a_train-scanvi-model.py>

modules/cell-type-neuroblastoma-04/main.nf

@@ -49,7 +49,8 @@ process train_singler_model {
     train-singler-model.R \
       --nbatlas_file ${nbatlas_sce_file} \
       --gtf_file ${gtf_file} \
-      --singler_model_file ${nbatlas_singler_model}
+      --singler_model_file ${nbatlas_singler_model} \
+      --threads ${task.cpus}
     """
   stub:
     nbatlas_singler_model = "nbatlas_singler_model.rds"
@@ -83,6 +84,44 @@ process train_scanvi_model {
     """
 }
 
+process classify_singler {
+  container params.cell_type_nb_04_container
+  tag "${sample_id}"
+  label 'mem_8'
+  input:
+    path singler_model
+    tuple val(sample_id),
+          val(project_id),
+          path(library_files)
+  output:
+    tuple val(sample_id),
+          val(project_id),
+          path(singler_files)
+  script:
+    singler_files = library_files
+      .collect{
+        it.name.replaceAll(/(?i).rds$/, "_singler.tsv")
+      }
+    """
+    for file in ${library_files}; do
+      classify-singler.R \
+        --sce_file ${file} \
+        --singler_model_file ${singler_model} \
+        --singler_output_tsv \$(basename \${file%.rds}_singler.tsv) \
+        --threads ${task.cpus}
+    done
+    """
+  stub:
+    singler_files = library_files
+      .collect{
+        it.name.replaceAll(/(?i).rds$/, "_singler.tsv")
+      }
+    """
+    for file in ${library_files}; do
+      touch \$(basename \${file%.rds}_singler.tsv)
+    done
+    """
+}
 
 
 workflow cell_type_neuroblastoma_04 {
@@ -96,6 +135,10 @@ workflow cell_type_neuroblastoma_04 {
         return [sample_id, project_id, library_files]
       }
 
+    /////////////////////////////////////////////////////
+    // Prepare references for cell type classification //
+    /////////////////////////////////////////////////////
+
     // convert NBAtlas to SCE and AnnData objects
     // emits: sce, anndata, hvg_file
     convert_nbatlas(file(params.cell_type_nb_04_nbatlas_url))
@@ -107,4 +150,12 @@ workflow cell_type_neuroblastoma_04 {
     // train scANVI model
     // outputs the scanvi model directory
     train_scanvi_model(convert_nbatlas.out.anndata)
+
+    /////////////////////////////////////////////////////
+    //        Perform  cell type classification        //
+    /////////////////////////////////////////////////////
+
+    // classify with SingleR
+    classify_singler(train_singler_model.out, libraries_ch)
+
 }

modules/cell-type-neuroblastoma-04/resources/usr/bin/classify-singler.R

@@ -0,0 +1,86 @@
+#!/usr/bin/env Rscript
+#
+# This script performs cell type annotation on an SCE using a reference SingleR
+# model and exports a TSV file with annotations
+# This script was adapted from:
+# https://github.com/AlexsLemonade/OpenScPCA-analysis/blob/82547028b5a9555d8cee40f6c1883015c990cc4f/analyses/cell-type-neuroblastoma-04/scripts/02_classify-singler.R
+
+suppressWarnings({
+  suppressPackageStartupMessages({
+    library(optparse)
+    library(SingleCellExperiment)
+  })
+})
+
+option_list <- list(
+  make_option(
+    opt_str = c("--sce_file"),
+    type = "character",
+    default = "",
+    help = "Input SCE object to run SingleR on"
+  ),
+  make_option(
+    opt_str = c("--singler_model_file"),
+    type = "character",
+    default = "",
+    help = "Path to trained SingleR model"
+  ),
+  make_option(
+    opt_str = c("--singler_output_tsv"),
+    type = "character",
+    default = "",
+    help = "Path to output TSV file to save lightweight SingleR annotations"
+  ),
+  make_option(
+    opt_str = c("--threads"),
+    type = "integer",
+    default = 4,
+    help = "Number of threads for SingleR to use"
+  ),
+  make_option(
+    opt_str = c("--seed"),
+    type = "integer",
+    default = 2025,
+    help = "Random seed"
+  )
+)
+
+# Parse options and check arguments
+opts <- parse_args(OptionParser(option_list = option_list))
+stopifnot(
+  "sce_file does not exist" = file.exists(opts$sce_file),
+  "singler_model_file does not exist" = file.exists(opts$singler_model_file)
+)
+set.seed(opts$seed)
+
+if (opts$threads == 1) {
+  bp_param <- BiocParallel::SerialParam()
+} else {
+  bp_param <- BiocParallel::MulticoreParam(opts$threads)
+}
+
+# Read trained model
+singler_model <- readRDS(opts$singler_model_file)
+
+# Read query SCE and convert ids to symbols
+# note that ids that don't map will retain Ensembl ids
+sce <- readRDS(opts$sce_file) |>
+  # use the gene_symbols column in the sce object for mapping
+  rOpenScPCA::sce_to_symbols(reference = "sce")
+
+# Perform annotation
+singler_result <- SingleR::classifySingleR(
+  sce,
+  singler_model,
+  BPPARAM = bp_param
+)
+
+# Extract TSV to save separately
+singler_df <- singler_result |>
+  # keep only labels, delta.next, and pruned.labels
+  purrr::discard_at("scores") |>
+  as.data.frame() |>
+  tibble::rownames_to_column("barcodes")
+
+# Export TSV
+readr::write_tsv(singler_df, opts$singler_output_tsv)

nextflow.config

@@ -111,7 +111,7 @@ profiles {
       release_bucket = "s3://openscpca-test-data-release-public-access" // test bucket
       results_bucket = "test/stub/results" // no output
       sim_bucket = "test/stub/simulated" // local output
-      project = "SCPCP000012" // a small project
+      project = "SCPCP000012,SCPCP000004" // a small project
 
       // Override full NBAtlas with empty file for stub testing
       cell_type_nb_04_nbatlas_url = "${projectDir}/modules/cell-type-neuroblastoma-04/resources/test/NBAtlas_test.rds"