Adapt combat

CHANGELOG.md

@@ -209,6 +209,8 @@
 
 * Remove `metrics/bras` , `metrics/n_inconsitent_peaks` and `control_methods/shuffle_by_cell_type` and `control_methods/batch`
 
+* Switch to `sva::` implementation for `methods/combat` (PR #149)
+
 ## BUG FIXES
 
 * Change n_inconsistent_peaks output to float and add R2 to main.nf (PR #40).

src/methods/combat/config.vsh.yaml

@@ -1,40 +1,50 @@
-# The API specifies which type of component this is.
-# It contains specifications for:
-#   - The input/output files
-#   - Common parameters
-#   - A unit test
 __merge__: ../../api/comp_method.yaml
 
-# A unique identifier for your component (required).
-# Can contain only lowercase letters or underscores.
 name: combat
-# A relatively short label, used when rendering visualisations (required)
-label: Combat
-# A one sentence summary of how this method works (required). Used when 
-# rendering summary tables.
-summary: "ComBat batch correction for single-cell data, implemented in the scanpy package"
-# A multi-line description of how this component works (required). Used
-# when rendering reference documentation.
+label: ComBat
+summary: "ComBat batch correction using an empirical Bayes framework, parameterised by correction target"
 description: |
-  Corrects for batch effects by fitting linear models, gains statistical power via an EB framework where information is borrowed across genes. 
-  This uses the implementation combat.py
+  Popular method using empirical Bayes frameworks to correct for batch effects in biological data.
+  Parametric adjustment enabled. This uses the implementation in the R package sva.
+
+  `--target` selects whether batches are corrected towards reference batch 1 (`goal`,
+  using `ref.batch = "1"`) or towards the pooled mean and variance across batches (`mid`,
+  the default `ref.batch = NULL`).
 references:
   doi: 10.1093/biostatistics/kxj037
 links:
-  repository: https://github.com/brentp/combat.py
-  documentation: https://scanpy.readthedocs.io/en/latest/api/generated/scanpy.pp.combat.html
+  documentation: https://bioconductor.org/packages/release/bioc/html/sva.html
+  repository: https://github.com/bioc/sva
 
-resources:
-  # The script of your component (required)
-  - type: python_script
-    path: script.py
-  - path: /src/utils/helper_functions.py
+info:
+  variants:
+    no_controls_to_goal:
+      target: goal
+    no_controls_to_mid:
+      target: mid
 
+argument_groups:
+  - name: Parameter set
+    arguments:
+      - type: string
+        name: --target
+        choices: [goal, mid]
+        default: mid
+        description: Whether to correct batches towards reference batch 1 (`goal`) or towards the pooled mean across batches (`mid`).
+
+resources:
+  - type: r_script
+    path: script.R
+  - type: r_script
+    path: /src/utils/helper_functions.R
 
 engines:
   - type: docker
-    image: openproblems/base_python:1
-
+    image: openproblems/base_r:1
+    setup:
+      - type: r
+        bioc: sva
+        packages: [ lifecycle, vctrs, dplyr, rlang ]
 
 runners:
   - type: executable

src/methods/combat/script.R

@@ -0,0 +1,57 @@
+library(anndata)
+requireNamespace("sva", quietly = TRUE)
+
+## VIASH START
+par <- list(
+  input = "resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/censored_split1.h5ad",
+  output = "output.h5ad",
+  target = "mid"
+)
+meta <- list(
+  name = "combat",
+  resources_dir = "src/utils"
+)
+## VIASH END
+
+source(paste0(meta$resources_dir, "/helper_functions.R"))
+
+cat("Reading input files\n")
+input <- anndata::read_h5ad(par[["input"]]) |>
+  subset_nocontrols()
+
+cat("Subset data\n")
+data_not_correct <- input[, !input$var$to_correct]
+data_to_correct <- input[, input$var$to_correct]
+
+# correct towards batch 1 (goal) or towards the pooled mean across batches (mid)
+ref_batch <- if (par[["target"]] == "goal") "1" else NULL
+
+cat("Run ComBat\n")
+output <- sva::ComBat(
+  dat = Matrix::t(data_to_correct$layers[["preprocessed"]]),
+  batch = data_to_correct$obs$batch,
+  ref.batch = ref_batch
+)
+
+cat("Preparing output Anndata\n")
+output <- Matrix::t(output)
+output <- cbind(output, data_not_correct$layers[["preprocessed"]])
+
+cat("Write output AnnData to file\n")
+output <- anndata::AnnData(
+  obs = input$obs[, integer(0)],
+  var = input$var[colnames(output), integer(0)],
+  layers = list(integrated = output),
+  uns = list(
+    dataset_id = input$uns$dataset_id,
+    method_id = meta$name,
+    parameters = list(target = par[["target"]])
+  )
+)
+
+# reorder the var in output
+output <- output[, input$var_names]
+
+output$write_h5ad(par[["output"]], compression = "gzip")
+
+cat("Written anndata of shape ", dim(output), " to file: ", par[["output"]], "\n")