diff --git a/CHANGELOG.md b/CHANGELOG.md
index 06459f8d..e2413da0 100644
--- a/CHANGELOG.md
+++ b/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).
diff --git a/src/methods/combat/config.vsh.yaml b/src/methods/combat/config.vsh.yaml
index cc76bf50..4b622f22 100644
--- a/src/methods/combat/config.vsh.yaml
+++ b/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
diff --git a/src/methods/combat/script.R b/src/methods/combat/script.R
new file mode 100644
index 00000000..b3ab50fa
--- /dev/null
+++ b/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")
diff --git a/src/methods/combat/script.py b/src/methods/combat/script.py
deleted file mode 100644
index 3ee98400..00000000
--- a/src/methods/combat/script.py
+++ /dev/null
@@ -1,66 +0,0 @@
-import sys
-
-import anndata as ad
-import numpy as np
-import scanpy as sc
-
-## VIASH START
-# Note: this section is auto-generated by viash at runtime. To edit it, make changes
-# in config.vsh.yaml and then run `viash config inject config.vsh.yaml`.
-par = {
-    "input": "resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/unintegrated_censored.h5ad",
-    "output": "output.h5ad",
-}
-meta = {"name": "combat"}
-## VIASH END
-
-sys.path.append(meta["resources_dir"])
-from helper_functions import subset_nocontrols
-
-print("Reading and preparing input files", flush=True)
-adata = ad.read_h5ad(par["input"])
-
-# Remove reference
-adata = subset_nocontrols(adata)
-
-print("Preprocess data", flush=True)
-markers_to_correct = adata.var[adata.var["to_correct"]].index.to_numpy()
-markers_not_correct = [
-    x for x in adata.var["to_correct"].index if x not in markers_to_correct
-]
-temp_var_index = [
-    *markers_to_correct,
-    *markers_not_correct,
-]  # temporary indexing : not corrected markers last
-
-adata_to_correct = adata[:, markers_to_correct].copy()
-adata_to_correct.X = adata_to_correct.layers[
-    "preprocessed"
-]  # sc.pp.combat operates only on .X
-
-print("Train model", flush=True)
-corrected_matrix = sc.pp.combat(adata_to_correct, key="batch", inplace=False)
-
-print("Aggregate data", flush=True)
-uncorrected_matrix = adata[:, markers_not_correct].layers["preprocessed"]
-integrated_matrix = np.concatenate([corrected_matrix, uncorrected_matrix], axis=1)
-
-print("Write output AnnData to file", flush=True)
-out_data = ad.AnnData(
-    obs=adata.obs[[]],
-    var=adata[:, temp_var_index].var[[]],
-    layers={"integrated": integrated_matrix},
-    uns={
-        "dataset_id": adata.uns["dataset_id"],
-        "method_id": meta["name"],
-        "parameters": {},
-    },
-)
-print("Re-ordering markers", flush=True)
-out_data = out_data[:, adata.var_names]
-
-out_data.write_h5ad(par["output"], compression="gzip")
-
-print(
-    "Written anndata of shape ", out_data.shape, " to file: ", par["output"], flush=True
-)