update transform function and broken azimuth in the pipeline

.Rbuildignore

@@ -14,4 +14,6 @@ pipeline_stores
 ^tests$
 ^dev$
 _targets
-target_framework
+target_framework
+^\.positai$
+^\.claude$

.github/workflows/rworkflows.yml

@@ -40,6 +40,12 @@ jobs:
           cont: ~
           rspm: ~
     steps:
+    - name: Prefer source installs on macOS as fallback (Bioc 3.23 mac binary gap)
+      if: runner.os == 'macOS'
+      run: |
+        mkdir -p ~/.R
+        echo 'options(install.packages.check.source = "no", pkgType = "source")' >> ~/.Rprofile
+      shell: bash
     - uses: neurogenomics/rworkflows@master
       with:
         run_bioccheck: ${{ false }}

.gitignore

@@ -38,3 +38,4 @@ fibrosis_data
 !README.md
 _target*
 meta
+.positai

DESCRIPTION

@@ -1,7 +1,9 @@
 Package: HPCell
 Title: Massively-Parallel R Native Pipeline for Single-Cell Analysis
-Version: 0.5.0
+Version: 0.6.1
 Authors@R: c(person("Stefano", "Mangiola", email = "mangiolastefano@gmail.com",
+                  role = c("aut")),
+                  person("Mengyuan", "Shen", email = "shen.m@wehi.edu.au",
                   role = c("aut", "cre")),
                   person("Jiayi", "Si", email = "si.j@wehi.edu.au",
                   role = c("aut"))
@@ -40,7 +42,6 @@ Imports:
   EnsDb.Hsapiens.v86,
   scater,
   SingleR,
-  celldex,
   scuttle,
   scDblFinder,
   magrittr,
@@ -74,6 +75,7 @@ Imports:
   rhdf5
 Suggests: 
     testthat (>= 3.0.0),
+    celldex,
     qs,
     Azimuth,
     CellChat,

NAMESPACE

@@ -156,8 +156,6 @@ importFrom(SummarizedExperiment,rowData)
 importFrom(biomaRt,getBM)
 importFrom(biomaRt,useMart)
 importFrom(callr,r)
-importFrom(celldex,BlueprintEncodeData)
-importFrom(celldex,MonacoImmuneData)
 importFrom(crew,crew_controller_local)
 importFrom(data.table,":=")
 importFrom(digest,digest)

NEWS.md

@@ -0,0 +1,58 @@
+# HPCell NEWS
+
+## HPCell 0.6.1
+
+### Bug fixes
+
+* Fixed `transform_utility()` to no longer assume the first assay of the input
+  object is named `"X"`. The canonical output assay name `"X"` is now set
+  explicitly via a single `assay_name` variable rather than being inferred from
+  the input, so inputs whose first assay carries any other name are handled
+  correctly.
+
+### Warnings
+
+* `transform_utility()` now emits a warning when the input SCE's assay is not
+  named `"X"`, reporting the original name before renaming it to `"X"` for
+  downstream consistency.
+
+## HPCell 0.6.0
+
+### Bug fixes
+
+* Fixed Azimuth-based cell-type label transfer (`annotation_label_transfer()`) failing
+  with SeuratObject >= 5.0.0. Azimuth 0.5.0 still calls the defunct
+  `GetAssayData(slot = ...)` interface; HPCell now temporarily patches the
+  Azimuth import environment so `slot` is mapped to `layer` during
+  `RunAzimuth()`.
+* Updated `GetAssayData()` calls from deprecated `slot = "counts"` to
+  `layer = "counts"` in `empty_droplet_id()`, `empty_droplet_threshold()`, and
+  `alive_identification()` for Seurat 5 compatibility.
+* Fixed `transform_utility()` so the `identity` transform also applies
+  `limit_max_to_scale()` before transformation. Previously only `expm1` was
+  pre-scaled, which could leave high-count samples unscaled and cause downstream
+  failures in the assay transformation pipeline.
+
+### Improvements
+
+* `non_batch_variation_removal()` now catches `SCTransform()` failures for edge
+  cases with very few overdispersion genes, emits a warning, and returns
+  `NULL` instead of stopping the pipeline.
+* `initialise_hpc()` gains a new `default_controller` argument to set the
+  default `crew` controller for targets that do not specify their own
+  controller via `tar_resources()`.
+
+### Documentation
+
+* Expanded and corrected roxygen documentation for `initialise_hpc()`, including
+  the new `default_controller` argument.
+
+### Development
+
+* Added CellNexus 2025 pipeline scripts under `dev/cellnexus-2025-scripts/`
+  (steps 1–10 for census download, metadata preparation, HPCell execution,
+  local-cache splitting, pseudobulk preparation, and metadata unification).
+* Updated CellNexus 2024 pipeline scripts, including renamed
+  `step9_unify_and_update_sce_metadata.R` and new steps for pseudobulk
+  preparation and CellNexus querying.
+* Added `.positai` and `.claude` to `.Rbuildignore` and `.gitignore`.

R/functions.R

@@ -48,7 +48,7 @@ empty_droplet_id <- function(input_read_RNA_assay,
 
   # Get counts
   if (inherits(input_read_RNA_assay, "Seurat")) {
-    counts <- GetAssayData(input_read_RNA_assay, assay, slot = "counts")
+    counts <- GetAssayData(input_read_RNA_assay, assay, layer = "counts")
   } else if (inherits(input_read_RNA_assay, "SingleCellExperiment")) {
     counts <- assay(input_read_RNA_assay, assay)
   }
@@ -269,7 +269,7 @@ empty_droplet_threshold<- function(input_read_RNA_assay,
 
   # Get counts
   if (inherits(input_read_RNA_assay, "Seurat")) {
-    counts <- GetAssayData(input_read_RNA_assay, assay, slot = "counts")
+    counts <- GetAssayData(input_read_RNA_assay, assay, layer = "counts")
   } else if (inherits(input_read_RNA_assay, "SingleCellExperiment")) {
     counts <- assay(input_read_RNA_assay, assay)
   }
@@ -312,9 +312,6 @@ empty_droplet_threshold<- function(input_read_RNA_assay,
 #'
 #' @return A tibble with cell-type annotation data.
 #'
-#' @importFrom celldex BlueprintEncodeData
-#' @importFrom celldex MonacoImmuneData
-#' 
 #' @importFrom Seurat CreateAssayObject
 #' @importFrom Seurat SCTransform
 #' @importFrom Seurat CreateSeuratObject
@@ -390,6 +387,12 @@ annotation_label_transfer <- function(input_read_RNA_assay,
     colnames(input_read_RNA_assay)[2]= "dummy___"
   }
 
+  if (!requireNamespace("celldex", quietly = TRUE))
+    stop(
+      "Package 'celldex' is required for SingleR-based cell-type annotation. ",
+      "Install it with: BiocManager::install('celldex')"
+    )
+
   #snapshotDate(): 2025-10-29
   blueprint <- celldex::BlueprintEncodeData(
     ensembl = feature_nomenclature == "ensembl"
@@ -479,7 +482,7 @@ annotation_label_transfer <- function(input_read_RNA_assay,
   } 
 
   if(nrow(data_annotated) <= 30 | is.null(reference_azimuth)){
-    
+
     # If too little immune cells
     return(data_annotated)
     #saveRDS(output_path)
@@ -515,11 +518,40 @@ annotation_label_transfer <- function(input_read_RNA_assay,
     azimuth_annotation = 
       tryCatch({
 
+        # This is necessary because Azimuth relies on Seurat: https://github.com/satijalab/azimuth/issues/195
+        if (!"Seurat" %in% .packages()) {
+          library(Seurat)
+        }
+
         if(ncol(input_read_RNA_assay)<200) k.weight = 25
         else k.weight = 50
 
-        input_read_RNA_assay |> RenameAssays(assay.name = assay, new.assay.name = "RNA") |> 
-          Azimuth::RunAzimuth(reference = reference_azimuth, assay = "RNA", umap.name = "refUMAP") |> 
+        # input_read_RNA_assay |> RenameAssays(assay.name = assay, new.assay.name = "RNA") |> 
+        #   Azimuth::RunAzimuth(reference = reference_azimuth, assay = "RNA", umap.name = "refUMAP")
+        {
+          tmp <- input_read_RNA_assay |> RenameAssays(assay.name = assay, new.assay.name = "RNA")
+
+          # Below is solved in Azimuth 0.5.1: https://github.com/satijalab/azimuth/issues/294
+          # # Azimuth 0.5.0 calls GetAssayData(slot = ...) in ConvertGeneNames(),
+          # # which is defunct in SeuratObject >= 5.0.0. The binding is in
+          # # `imports:Azimuth` (parent.env of Azimuth's namespace), confirmed by
+          # # inspection. Patch it there so all internal Azimuth calls are
+          # # intercepted, then restore on exit.
+          # .az_imp <- parent.env(asNamespace("Azimuth"))
+          # .orig_gad <- get("GetAssayData", envir = .az_imp, inherits = FALSE)
+          # .shim_gad <- local({
+          #   orig <- .orig_gad
+          #   function(object, slot = NULL, layer = NULL, ...) {
+          #     if (!is.null(slot) && is.null(layer)) layer <- slot
+          #     orig(object, layer = layer, ...)
+          #   }
+          # })
+          # try(unlockBinding("GetAssayData", .az_imp), silent = TRUE)
+          # assign("GetAssayData", .shim_gad, envir = .az_imp)
+          # on.exit(assign("GetAssayData", .orig_gad, envir = .az_imp), add = TRUE)
+
+          Azimuth::RunAzimuth(tmp, reference = reference_azimuth, assay = "RNA", umap.name = "refUMAP")
+        } %>% 
           as_tibble() |>
           dplyr::select(.cell, any_of(
             c(
@@ -634,7 +666,7 @@ alive_identification <- function(input_read_RNA_assay,
 
 
   if (inherits(input_read_RNA_assay, "Seurat")) {
-    counts <- GetAssayData(input_read_RNA_assay, assay = assay, slot = "counts")
+    counts <- GetAssayData(input_read_RNA_assay, assay = assay, layer = "counts")
     if (!any(str_which(colnames(input_read_RNA_assay[[]]), nFeature_name)) ||
         !any(str_which(colnames(input_read_RNA_assay[[]]), nCount_name))) {
       input_read_RNA_assay[[nFeature_name]] <-
@@ -661,21 +693,18 @@ alive_identification <- function(input_read_RNA_assay,
 
   # Returns a named vector of IDs
   # Matches the gene id's row by row and inserts NA when it can't find gene names
-  if (feature_nomenclature == "symbol") {
-    location <- mapIds(
-      EnsDb.Hsapiens.v86,
-      keys=rownames(input_read_RNA_assay),
-      column="SEQNAME",
-      keytype="SYMBOL"
-    )
-  }
-
+  location <- mapIds(
+    EnsDb.Hsapiens.v86,
+    keys = rownames(input_read_RNA_assay),
+    column = "SEQNAME",
+    keytype = if (feature_nomenclature == "symbol") "SYMBOL" else "GENEID"
+  )
 
-  which_mito = rownames(input_read_RNA_assay) |> str_which("^MT")
+  which_mito = which(location == "MT")
 
   # mitochondrion =
   #   input_read_RNA_assay |>
-  #   GetAssayData( slot = "counts", assay=assay) |>
+  #   GetAssayData( layer = "counts", assay=assay) |>
   # 
   #   # Join mitochondrion statistics
   #   # Compute per-cell quality control metrics for a count matrix or a SingleCellExperiment
@@ -726,8 +755,21 @@ alive_identification <- function(input_read_RNA_assay,
     as_tibble(rownames = ".cell") %>%
     dplyr::select(-sum, -detected)
 
-  # I HAVE TO DROP UNIQUE, AS SOON AS THE BUG IN SEURAT IS RESOLVED. UNIQUE IS BUG PRONE HERE.
-  percentage_output = PercentageFeatureSet(input_read_RNA_assay,  pattern = "^RPS|^RPL", assay = assay)
+  if (feature_nomenclature == "symbol") {
+    percentage_output = PercentageFeatureSet(input_read_RNA_assay, pattern = "^RPS|^RPL", assay = assay)
+  } else {
+    # Ensembl IDs: resolve ribo gene IDs from biomart reference
+    data(ensembl_genes_biomart)
+    ribosome_ensembl_ids <- ensembl_genes_biomart[
+      grep("^(RPL|RPS)", ensembl_genes_biomart$external_gene_name), "ensembl_gene_id"
+    ]
+    ribosome_features <- intersect(ribosome_ensembl_ids, rownames(input_read_RNA_assay))
+    percentage_output = PercentageFeatureSet(
+      input_read_RNA_assay,
+      features = if (length(ribosome_features) > 0) ribosome_features else character(0),
+      assay = assay
+    )
+  }
   percentage_output = percentage_output[!duplicated(names(percentage_output))]
   # Compute ribosome statistics
   ribosome =
@@ -765,7 +807,7 @@ alive_identification <- function(input_read_RNA_assay,
         ribosome |>
         # Only retrieve metadata so nesting in the next step won't break
         left_join(input_read_RNA_assay |> as_tibble() |> select(.cell, all_of(cell_type_column)), by = ".cell") |> as_tibble()
-        
+
     }
 
 
@@ -808,12 +850,13 @@ alive_identification <- function(input_read_RNA_assay,
   mitochondrion |>
     left_join(ribosome) |>
     mutate(alive = !high_mitochondrion) |> # & !high_ribosome ) |>
-  # Select informative columns
+    # Select informative columns
     select(.cell, {{cell_type_column}}, contains("subsets"), contains("observation"),
            contains("high"), alive)
 }
 
 
+
 #' Doublet Identification
 #'
 #' @description
@@ -1105,20 +1148,25 @@ non_batch_variation_removal <- function(input_read_RNA_assay,
   min_cells <- if (bigm_lgl) 5L else 0L
   new_min_cells = calculate_num_genes_express_in_cells(m, min_cells) # update min_cells if needed
 
-  input_read_RNA_assay <- 
+  input_read_RNA_assay <- tryCatch(
     input_read_RNA_assay |>
-          SCTransform(
-            assay = assay,
-            return.only.var.genes = FALSE,
-            residual.features = NULL,
-            vars.to.regress = factors_to_regress,
-            vst.flavor = "v2",
-            scale_factor = 2186,
-            conserve.memory = TRUE,
-            min_cells = new_min_cells
-          )
-  # |> 
-  #   GetAssayData(assay="SCT")
+      SCTransform(
+        assay = assay,
+        return.only.var.genes = FALSE,
+        residual.features = NULL,
+        vars.to.regress = factors_to_regress,
+        vst.flavor = "v2",
+        scale_factor = 2186,
+        conserve.memory = TRUE,
+        min_cells = new_min_cells
+      ),
+    error = function(e) {
+      warning("HPCell says: post transformed distribution of the sample introduced few overdispersion genes, return NULL because these are very few edge cases.")
+      NULL
+    }
+  )
+
+  if (is.null(input_read_RNA_assay)) return(NULL)
 
   sct_mat <- input_read_RNA_assay |>  GetAssayData(assay="SCT")