Adding new method: scMerge2

CHANGELOG.md

@@ -29,6 +29,8 @@
 
 * Added `method/limma_removebatcheffect` component (PR #79).
 
+* Added `methods/scmerge2` component (PR #63).
+
 ## Minor changes
 
 * Un-pin the scPRINT version and update parameters (PR #51)

src/methods/scmerge2/config.vsh.yaml

@@ -0,0 +1,59 @@
+__merge__: /src/api/comp_method.yaml
+name: scmerge2
+label: scMerge2
+summary: "scMerge2 integrates single-cell RNA-seq datasets by removing unwanted variation estimated from stably expressed genes and pseudo-replicates."
+description: |
+  scMerge2 corrects batch effects while preserving biological signal. It identifies a set of stably
+  expressed genes (SEGs), which are assumed to remain consistent across datasets, and uses them as
+  negative controls in a factor analysis model that estimates and removes unwanted variation.
+  Pseudo-replicates, constructed from pseudo-bulk profiles of cells grouped within each batch, serve
+  as anchors for the alignment.
+
+  Cell type labels are optional. When they are supplied the pseudo-replicates are built per cell
+  type; when they are not, scMerge2 identifies the groupings itself with a mutual nearest cluster
+  procedure. Both modes are exposed through `--cell_type_aware`.
+references:
+  doi:
+    - 10.1073/pnas.1820006116
+links:
+  documentation: https://sydneybiox.github.io/scMerge/articles/scMerge2.html
+  repository: https://github.com/SydneyBioX/scMerge
+info:
+  method_types: [embedding]
+  preferred_normalization: log_cp10k
+  variants:
+    scmerge2_unsupervised:
+    scmerge2_semisupervised:
+      cell_type_aware: true
+arguments:
+  - name: --cell_type_aware
+    type: boolean
+    default: false
+    description: |
+      Build the pseudo-replicates per cell type, using obs['cell_type']. When false, scMerge2
+      identifies the cell groupings itself.
+  - name: --n_control_genes
+    type: integer
+    default: 1000
+    description: Number of top-ranked stably expressed genes to use as negative controls.
+  - name: --n_dim
+    type: integer
+    default: 50
+    description: Number of principal components in the output embedding.
+resources:
+  - type: r_script
+    path: script.R
+engines:
+  - type: docker
+    image: openproblems/base_r:1
+    setup:
+      - type: apt
+        packages: cmake
+      - type: r
+        bioc:
+          - scMerge
+runners:
+  - type: executable
+  - type: nextflow
+    directives:
+      label: [hightime, highmem, midcpu]

src/methods/scmerge2/script.R

@@ -0,0 +1,86 @@
+cat("Loading dependencies\n")
+requireNamespace("anndata", quietly = TRUE)
+library(Matrix, warn.conflicts = FALSE)
+requireNamespace("scMerge", quietly = TRUE)
+requireNamespace("BiocParallel", quietly = TRUE)
+requireNamespace("BiocSingular", quietly = TRUE)
+requireNamespace("methods", quietly = TRUE)
+
+## VIASH START
+par <- list(
+  input = "resources_test/task_batch_integration/cxg_immune_cell_atlas/dataset.h5ad",
+  output = "output.h5ad",
+  cell_type_aware = FALSE,
+  n_control_genes = 1000L,
+  n_dim = 50L
+)
+meta <- list(
+  name = "scmerge2",
+  cpus = 1L
+)
+## VIASH END
+
+n_cpus <- if (is.null(meta$cpus)) 1L else meta$cpus
+bpparam <- BiocParallel::MulticoreParam(workers = n_cpus)
+
+cat("Read input\n")
+adata <- anndata::read_h5ad(par$input)
+
+# both scSEGIndex and scMerge2 want genes in the rows. scMerge2 only coerces a *dense* matrix to
+# CsparseMatrix, so a row-compressed one would slip past that check and break later on
+exprs_mat <- methods::as(Matrix::t(adata$layers[["normalized"]]), "CsparseMatrix")
+rownames(exprs_mat) <- as.character(adata$var_names)
+colnames(exprs_mat) <- as.character(adata$obs_names)
+
+cat("Select stably expressed genes\n")
+seg_df <- scMerge::scSEGIndex(exprs_mat = exprs_mat, BPPARAM = bpparam)
+seg_df <- seg_df[order(seg_df$segIdx, decreasing = TRUE), , drop = FALSE]
+ctl <- rownames(seg_df)[seq_len(min(par$n_control_genes, nrow(seg_df)))]
+
+cat("Run scMerge2\n")
+out <- scMerge::scMerge2(
+  exprsMat = exprs_mat,
+  batch = as.character(adata$obs$batch),
+  cellTypes = if (par$cell_type_aware) as.character(adata$obs$cell_type) else NULL,
+  ctl = ctl,
+  use_bpparam = bpparam,
+  use_bsparam = BiocSingular::RandomParam(),
+  verbose = TRUE
+)
+
+cat("Compute embedding\n")
+newY <- out$newY
+stopifnot(ncol(newY) == adata$n_obs)
+if (!is.null(colnames(newY))) {
+  newY <- newY[, adata$obs_names, drop = FALSE]
+}
+# subtracting the estimated unwanted variation makes the corrected matrix dense, so let
+# BiocSingular stream it rather than materialising it in one go
+corrected <- t(newY)
+n_dim <- min(par$n_dim, min(dim(corrected)) - 1L)
+embedding <- BiocSingular::runPCA(
+  corrected,
+  rank = n_dim,
+  center = TRUE,
+  scale = FALSE,
+  BSPARAM = BiocSingular::RandomParam(),
+  BPPARAM = bpparam
+)$x
+rownames(embedding) <- adata$obs_names
+
+cat("Store output\n")
+output <- anndata::AnnData(
+  obs = adata$obs[, c()],
+  var = adata$var[, c()],
+  obsm = list(
+    X_emb = embedding
+  ),
+  uns = list(
+    dataset_id = adata$uns[["dataset_id"]],
+    normalization_id = adata$uns[["normalization_id"]],
+    method_id = meta$name
+  )
+)
+
+cat("Write output to file\n")
+zzz <- output$write_h5ad(par$output, compression = "gzip")

src/workflows/run_benchmark/config.vsh.yaml

@@ -106,6 +106,7 @@ dependencies:
   - name: methods/scanorama_correct
   - name: methods/scanorama_integrate
   - name: methods/scanvi
+  - name: methods/scmerge2
   - name: methods/scgpt_finetuned
   - name: methods/scgpt_zeroshot
   - name: methods/scimilarity

src/workflows/run_benchmark/main.nf

@@ -36,6 +36,7 @@ methods = [
   scanorama_correct,
   scanorama_integrate,
   scanvi,
+  scmerge2,
   scgpt_finetuned.run(
     args: [model: file("s3://openproblems-work/cache/scGPT_human.zip")]
   ),