diff --git a/.gitignore b/.gitignore
index ae5a6576..3a57f454 100644
--- a/.gitignore
+++ b/.gitignore
@@ -10,6 +10,6 @@ trace-*
 .ipynb_checkpoints
 /temp
 /.vscode
-
+*_iterative_test.*
 /resources_raw
 /*.tar
\ No newline at end of file
diff --git a/scripts/run_benchmark/run_full_seqeracloud.sh b/scripts/run_benchmark/run_full_seqeracloud.sh
index 313f7514..979440c9 100755
--- a/scripts/run_benchmark/run_full_seqeracloud.sh
+++ b/scripts/run_benchmark/run_full_seqeracloud.sh
@@ -17,16 +17,16 @@ cat > /tmp/params.yaml << HERE
 input_states: s3://openproblems-data/resources/task_cyto_batch_integration/datasets/**/state.yaml
 rename_keys: 'input_censored_split1:output_censored_split1;input_censored_split2:output_censored_split2;input_unintegrated:output_unintegrated'
 output_state: "state.yaml"
-# settings: '{"methods_include": ["combat", "gaussnorm"]}'
+settings: '{"metrics_exclude": ["cms"], "methods_include": ["mnnpy", "cytovi"]}'
 publish_dir: "$publish_dir"
 HERE
 
 tw launch https://github.com/openproblems-bio/task_cyto_batch_integration.git \
-  --revision build/main \
+  --revision build/fix_failed_stuff \
   --pull-latest \
   --main-script target/nextflow/workflows/run_benchmark/main.nf \
   --workspace 53907369739130 \
   --params-file /tmp/params.yaml \
   --entry-name auto \
   --config common/nextflow_helpers/labels_tw.config \
-  --labels task_cyto_batch_integration,full
+  --labels task_cyto_batch_integration,mnnnpy
diff --git a/src/control_methods/perfect_integration/script.py b/src/control_methods/perfect_integration/script.py
index f79e34e3..6e5869bf 100644
--- a/src/control_methods/perfect_integration/script.py
+++ b/src/control_methods/perfect_integration/script.py
@@ -4,8 +4,8 @@
 # The following code has been auto-generated by Viash.
 par = {
     "input_unintegrated": "resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/unintegrated.h5ad",
-    "output_integrated_split1": "resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/integrated_split1.h5ad",
-    "output_integrated_split2": "resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/integrated_split2.h5ad",
+    "output_integrated_split1": "resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/perfect_integrated_split1.h5ad",
+    "output_integrated_split2": "resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/perfect_integrated_split2.h5ad",
 }
 meta = {"name": "perfect_integration"}
 
diff --git a/src/control_methods/shuffle_integration_by_batch/config.vsh.yaml b/src/control_methods/shuffle_integration_by_batch/config.vsh.yaml
index 0f07e06a..dfa15d45 100644
--- a/src/control_methods/shuffle_integration_by_batch/config.vsh.yaml
+++ b/src/control_methods/shuffle_integration_by_batch/config.vsh.yaml
@@ -3,7 +3,7 @@ name: shuffle_integration_by_batch
 label: Shuffle integration by batch
 summary: Integrations are randomly permuted within each batch
 description: Integrations are randomly permuted within each batch
-status: disabled
+# status: disabled
 resources:
   - type: python_script
     path: script.py
diff --git a/src/control_methods/shuffle_integration_by_batch/script.py b/src/control_methods/shuffle_integration_by_batch/script.py
index 2a3aaee1..0f7fbe25 100644
--- a/src/control_methods/shuffle_integration_by_batch/script.py
+++ b/src/control_methods/shuffle_integration_by_batch/script.py
@@ -16,26 +16,47 @@
 
 print("Reading and preparing input files", flush=True)
 adata = ad.read_h5ad(par["input_unintegrated"])
+adata_split1 = adata[(adata.obs.is_control > 0) | (adata.obs.batch == 1)].copy()
+adata_split2 = adata[(adata.obs.is_control > 0) | (adata.obs.split == 2)].copy()
 
-adata.obs["batch_str"] = adata.obs["batch"].astype(str)
+print("Randomise features - split 1", flush=True)
+adata_split1.obs["batch_str"] = adata_split1.obs["batch"].astype(str)
+integrated = _randomize_features(
+    adata_split1.layers["preprocessed"],
+    partition=adata_split1.obs["batch"],
+)
+
+# create new anndata
+output_split1 = ad.AnnData(
+    obs=adata_split1.obs[[]],
+    var=adata_split1.var[[]],
+    layers={"integrated": integrated},
+    uns={
+        "dataset_id": adata_split1.uns["dataset_id"],
+        "method_id": meta["name"],
+        "parameters": {},
+    },
+)
 
-print("Randomise features", flush=True)
+print("Randomise features - split 2", flush=True)
+adata_split2.obs["batch_str"] = adata_split2.obs["batch"].astype(str)
 integrated = _randomize_features(
-    adata.layers["preprocessed"],
-    partition=adata.obs["batch"],
+    adata_split2.layers["preprocessed"],
+    partition=adata_split2.obs["batch"],
 )
 
 # create new anndata
-output = ad.AnnData(
-    obs=adata.obs[[]],
-    var=adata.var[[]],
+output_split2 = ad.AnnData(
+    obs=adata_split2.obs[[]],
+    var=adata_split2.var[[]],
     layers={"integrated": integrated},
     uns={
-        "dataset_id": adata.uns["dataset_id"],
+        "dataset_id": adata_split2.uns["dataset_id"],
         "method_id": meta["name"],
         "parameters": {},
     },
 )
 
 print("Write output AnnData to file", flush=True)
-output.write_h5ad(par["output"], compression="gzip")
+output_split1.write_h5ad(par["output_integrated_split1"], compression="gzip")
+output_split2.write_h5ad(par["output_integrated_split2"], compression="gzip")
diff --git a/src/control_methods/shuffle_integration_by_cell_type/script.py b/src/control_methods/shuffle_integration_by_cell_type/script.py
index 3abca3c0..d4bcd218 100644
--- a/src/control_methods/shuffle_integration_by_cell_type/script.py
+++ b/src/control_methods/shuffle_integration_by_cell_type/script.py
@@ -16,28 +16,47 @@
 
 print("Reading and preparing input files", flush=True)
 adata = ad.read_h5ad(par["input_unintegrated"])
+adata_split1 = adata[(adata.obs.is_control > 0) | (adata.obs.batch == 1)].copy()
+adata_split2 = adata[(adata.obs.is_control > 0) | (adata.obs.split == 2)].copy()
 
-adata.obs["batch_str"] = adata.obs["batch"].astype(str)
-
-print("Randomise features", flush=True)
+print("Randomise features - split 1", flush=True)
+adata_split1.obs["batch_str"] = adata_split1.obs["batch"].astype(str)
 integrated = _randomize_features(
-    adata.layers["preprocessed"],
-    partition=adata.obs["cell_type"],
+    adata_split1.layers["preprocessed"],
+    partition=adata_split1.obs["cell_type"],
 )
 
 # create new anndata
-output = ad.AnnData(
-    obs=adata.obs[[]],
-    var=adata.var[[]],
+output_split1 = ad.AnnData(
+    obs=adata_split1.obs[[]],
+    var=adata_split1.var[[]],
     layers={"integrated": integrated},
     uns={
-        "dataset_id": adata.uns["dataset_id"],
+        "dataset_id": adata_split1.uns["dataset_id"],
         "method_id": meta["name"],
         "parameters": {},
     },
 )
 
-all(x==y for x,y in zip(output.var_names, adata.var_names))
+print("Randomise features - split 2", flush=True)
+adata_split2.obs["batch_str"] = adata_split2.obs["batch"].astype(str)
+integrated = _randomize_features(
+    adata_split2.layers["preprocessed"],
+    partition=adata_split2.obs["cell_type"],
+)
+
+# create new anndata
+output_split2 = ad.AnnData(
+    obs=adata_split2.obs[[]],
+    var=adata_split2.var[[]],
+    layers={"integrated": integrated},
+    uns={
+        "dataset_id": adata_split2.uns["dataset_id"],
+        "method_id": meta["name"],
+        "parameters": {},
+    },
+)
 
 print("Write output AnnData to file", flush=True)
-output.write_h5ad(par["output"], compression="gzip")
+output_split1.write_h5ad(par["output_integrated_split1"], compression="gzip")
+output_split2.write_h5ad(par["output_integrated_split2"], compression="gzip")
\ No newline at end of file
diff --git a/src/methods/batchadjust_all_controls/config.vsh.yaml b/src/methods/batchadjust_all_controls/config.vsh.yaml
index 50fa14c9..ca9bc1ad 100644
--- a/src/methods/batchadjust_all_controls/config.vsh.yaml
+++ b/src/methods/batchadjust_all_controls/config.vsh.yaml
@@ -66,6 +66,7 @@ resources:
     path: script.R
   - path: /src/utils/anndata_to_fcs.R
   - path: BatchAdjust.R
+  - path: utils.R
 
 engines:
   # Specifications for the Docker image for this component.
diff --git a/src/methods/batchadjust_all_controls/script.R b/src/methods/batchadjust_all_controls/script.R
index 350ed2d9..74ba4e0f 100644
--- a/src/methods/batchadjust_all_controls/script.R
+++ b/src/methods/batchadjust_all_controls/script.R
@@ -2,15 +2,19 @@ library(anndata)
 library(flowCore)
 ## VIASH START
 par <- list(
-  input = "resources_test/.../input.h5ad",
-  output = "output.h5ad"
+  input = "resources_test/debug/batchadjust/_viash_par/input_1/censored_split1.h5ad",
+  output = "resources_test/debug/batchadjust/output.h5ad",
+  percentile = as.integer('80')
 )
 meta <- list(
   name = "batchadjust_all_controls",
-  temp_dir = "/tmp"
+  temp_dir = "/tmp",
+  resources_dir = "src/methods/batchadjust_all_controls"
 )
+source("src/utils/anndata_to_fcs.R")
 ## VIASH END
 
+source(paste0(meta$resources_dir, "/utils.R"))
 source(paste0(meta$resources_dir, "/anndata_to_fcs.R"))
 source(paste0(meta$resources_dir, "/BatchAdjust.R"))
 
@@ -33,7 +37,11 @@ source(paste0(meta$resources_dir, "/BatchAdjust.R"))
 cat("Reading input files\n")
 input <- anndata::read_h5ad(par[["input"]])
 #use Original_ID column to restore cell order after I/O operations
-input$layers["preprocessed"][, "Original_ID"] <- seq(1, dim(input)[1])
+# input$layers["preprocessed"][, "Original_ID"] <- seq(1, dim(input)[1])
+
+original_id_in_var <- "Original_ID" %in% input$var_names
+
+input <- add_original_id(input)
 
 cat("Split cells\n")
 input_controls <- input[input$obs$is_control != 0, ]
@@ -48,9 +56,14 @@ print(input_no_controls)
 
 #avoid NA due to invalid factor level
 input_controls$obs$sample <- as.character(input_controls$obs$sample)
+
+
+# make sure there is _ after the batch1 or batch2, otherwise batchadjust won't find the fcs files.
+input_no_controls$obs$sample <- sapply(input_no_controls$obs$sample, fix_batch_underscore_anynum)
+
 # Set sample names for batch-specific control files
-input_controls$obs$sample[input_controls$obs$batch == 1] <- "0Batch1_anchor"
-input_controls$obs$sample[input_controls$obs$batch == 2] <- "0Batch2_anchor"
+input_controls$obs$sample[input_controls$obs$batch == 1] <- "Batch1_anchor"
+input_controls$obs$sample[input_controls$obs$batch == 2] <- "Batch2_anchor"
 
 cat("Writing FCS files\n")
 anndata_to_fcs(input_controls, out_dir = meta[["temp_dir"]])
@@ -72,7 +85,7 @@ BatchAdjust(
   outdir = output_dir,
   channelsFile = paste0(meta[["temp_dir"]], "/to_correct_list.txt"),
   anchorKeyword = "anchor",
-  batchKeyword = "atch", #skip 'b' to make it robust to upper/lowercase
+  batchKeyword = "Batch", #skip 'b' to make it robust to upper/lowercase
   method = perc,
   transformation = FALSE,
   addExt = NULL,
@@ -96,6 +109,11 @@ if (FALSE %in% order_check) {
   stop("Failed in restoring indexing")
 }
 
+# Remove Original_ID if it was not there in the beginning
+if (!original_id_in_var) {
+    corrected_matrix$Original_ID <- NULL
+}
+
 cat("Write output AnnData to file\n")
 output <- anndata::AnnData(
   obs = input$obs[, integer(0)],
diff --git a/src/methods/batchadjust_all_controls/utils.R b/src/methods/batchadjust_all_controls/utils.R
new file mode 100644
index 00000000..cf4f36d5
--- /dev/null
+++ b/src/methods/batchadjust_all_controls/utils.R
@@ -0,0 +1,51 @@
+fix_batch_underscore_anynum <- function(x) {
+    # Pattern: Batch followed by one or more digits, NOT followed by underscore
+    pattern <- "(Batch\\d+)(?!_)"
+    
+    # Add underscore if missing
+    x <- sub(pattern, "\\1_", x, perl = TRUE)
+    
+    return(x)
+}
+
+add_original_id <- function(input) {
+    if (!"Original_ID" %in% input$var_names) {
+        cat("Adding Original_ID to var and recreating input anndata\n")
+        old_var <- input$var
+        old_var[] <- lapply(old_var, function(x) if (is.factor(x)) as.character(x) else x)
+        old_var <- rbind(old_var,
+            data.frame(
+                numeric_id=length(input$var_names) + 1,
+                channel="Original_ID",
+                marker="",
+                marker_type="other",
+                to_correct=FALSE,
+                row.names = "Original_ID"
+            )
+        )
+        old_var[] <- lapply(old_var, function(x) if (is.character(x)) as.factor(x) else x)
+        new_mat <- Matrix::as.matrix(input$layers[["preprocessed"]])
+        new_mat <- cbind(new_mat, Original_ID = seq_len(nrow(new_mat)))
+
+        # create new anndata
+        input <- anndata::AnnData(
+            X = new_mat,
+            obs = input$obs,
+            var = old_var,
+            layers = list(preprocessed = new_mat),
+            uns = list(
+                dataset_description = input$uns$dataset_description,
+                dataset_id = input$uns$dataset_id,
+                dataset_name = input$uns$dataset_name,
+                dataset_organism = input$uns$dataset_organism,
+                dataset_reference = input$uns$dataset_reference,
+                dataset_summary = input$uns$dataset_summary,
+                dataset_url = input$uns$dataset_url
+            )
+        )
+    } else {
+        cat("Adding new Original_ID var\n")
+        input$layers["preprocessed"][, "Original_ID"] <- seq(1, dim(input)[1])
+    }
+    return(input)
+}
diff --git a/src/methods/batchadjust_one_control/config.vsh.yaml b/src/methods/batchadjust_one_control/config.vsh.yaml
index 70851434..55dbeca5 100644
--- a/src/methods/batchadjust_one_control/config.vsh.yaml
+++ b/src/methods/batchadjust_one_control/config.vsh.yaml
@@ -66,6 +66,7 @@ resources:
     path: script.R
   - path: /src/utils/anndata_to_fcs.R
   - path: BatchAdjust.R
+  - path: utils.R
 
 engines:
   # Specifications for the Docker image for this component.
diff --git a/src/methods/batchadjust_one_control/script.R b/src/methods/batchadjust_one_control/script.R
index 0ce78a70..ca609d37 100644
--- a/src/methods/batchadjust_one_control/script.R
+++ b/src/methods/batchadjust_one_control/script.R
@@ -2,15 +2,18 @@ library(anndata)
 library(flowCore)
 ## VIASH START
 par <- list(
-  input = "resources_test/.../input.h5ad",
-  output = "output.h5ad"
+  input = "resources_test/debug/batchadjust/_viash_par/input_1/censored_split1.h5ad",
+  output = "resources_test/debug/batchadjust/output.h5ad",
+  percentile = as.integer('80')
 )
 meta <- list(
   name = "batchadjust_all_controls",
-  temp_dir = "/tmp"
+  temp_dir = "resources_test/tmp",
+  resources_dir = 
 )
 ## VIASH END
 
+source(paste0(meta$resources_dir, "/utils.R"))
 source(paste0(meta$resources_dir, "/anndata_to_fcs.R"))
 source(paste0(meta$resources_dir, "/BatchAdjust.R"))
 
@@ -33,7 +36,11 @@ source(paste0(meta$resources_dir, "/BatchAdjust.R"))
 cat("Reading input files\n")
 input <- anndata::read_h5ad(par[["input"]])
 #use Original_ID column to restore cell order after I/O operations
-input$layers["preprocessed"][, "Original_ID"] <- seq(1, dim(input)[1])
+# input$layers["preprocessed"][, "Original_ID"] <- seq(1, dim(input)[1])
+
+original_id_in_var <- "Original_ID" %in% input$var_names
+
+input <- add_original_id(input)
 
 cat("Split cells\n")
 input_controls <- input[input$obs$is_control == 1, ]
@@ -48,9 +55,13 @@ print(input_no_controls)
 
 #avoid NA due to invalid factor level
 input_controls$obs$sample <- as.character(input_controls$obs$sample)
+
+# make sure there is _ after the batch1 or batch2, otherwise batchadjust won't find the fcs files.
+input_no_controls$obs$sample <- sapply(input_no_controls$obs$sample, fix_batch_underscore_anynum)
+
 # Set sample names for batch-specific control files
-input_controls$obs$sample[input_controls$obs$batch == 1] <- "0Batch1_anchor"
-input_controls$obs$sample[input_controls$obs$batch == 2] <- "0Batch2_anchor"
+input_controls$obs$sample[input_controls$obs$batch == 1] <- "Batch1_anchor"
+input_controls$obs$sample[input_controls$obs$batch == 2] <- "Batch2_anchor"
 
 cat("Writing FCS files\n")
 anndata_to_fcs(input_controls, out_dir = meta[["temp_dir"]])
@@ -72,7 +83,7 @@ BatchAdjust(
   outdir = output_dir,
   channelsFile = paste0(meta[["temp_dir"]], "/to_correct_list.txt"),
   anchorKeyword = "anchor",
-  batchKeyword = "atch", #skip 'b' to make it robust to upper/lowercase
+  batchKeyword = "Batch", #skip 'b' to make it robust to upper/lowercase
   method = perc,
   transformation = FALSE,
   addExt = NULL,
@@ -96,6 +107,11 @@ if (FALSE %in% order_check) {
   stop("Failed in restoring indexing")
 }
 
+# Remove Original_ID if it was not there in the beginning
+if (!original_id_in_var) {
+    corrected_matrix$Original_ID <- NULL
+}
+
 cat("Write output AnnData to file\n")
 output <- anndata::AnnData(
   obs = input$obs[, integer(0)],
diff --git a/src/methods/batchadjust_one_control/utils.R b/src/methods/batchadjust_one_control/utils.R
new file mode 100644
index 00000000..ac869b5c
--- /dev/null
+++ b/src/methods/batchadjust_one_control/utils.R
@@ -0,0 +1,52 @@
+library(anndata)
+fix_batch_underscore_anynum <- function(x) {
+    # Pattern: Batch followed by one or more digits, NOT followed by underscore
+    pattern <- "(Batch\\d+)(?!_)"
+    
+    # Add underscore if missing
+    x <- sub(pattern, "\\1_", x, perl = TRUE)
+    
+    return(x)
+}
+
+add_original_id <- function(input) {
+    if (!"Original_ID" %in% input$var_names) {
+        cat("Adding Original_ID to var and recreating input anndata\n")
+        old_var <- input$var
+        old_var[] <- lapply(old_var, function(x) if (is.factor(x)) as.character(x) else x)
+        old_var <- rbind(old_var,
+            data.frame(
+                numeric_id=length(input$var_names) + 1,
+                channel="Original_ID",
+                marker="",
+                marker_type="other",
+                to_correct=FALSE,
+                row.names = "Original_ID"
+            )
+        )
+        old_var[] <- lapply(old_var, function(x) if (is.character(x)) as.factor(x) else x)
+        new_mat <- Matrix::as.matrix(input$layers[["preprocessed"]])
+        new_mat <- cbind(new_mat, Original_ID = seq_len(nrow(new_mat)))
+
+        # create new anndata
+        input <- anndata::AnnData(
+            X = new_mat,
+            obs = input$obs,
+            var = old_var,
+            layers = list(preprocessed = new_mat),
+            uns = list(
+                dataset_description = input$uns$dataset_description,
+                dataset_id = input$uns$dataset_id,
+                dataset_name = input$uns$dataset_name,
+                dataset_organism = input$uns$dataset_organism,
+                dataset_reference = input$uns$dataset_reference,
+                dataset_summary = input$uns$dataset_summary,
+                dataset_url = input$uns$dataset_url
+            )
+        )
+    } else {
+        cat("Adding new Original_ID var\n")
+        input$layers["preprocessed"][, "Original_ID"] <- seq(1, dim(input)[1])
+    }
+    return(input)
+}
diff --git a/src/methods/cycombine_all_controls_to_goal/config.vsh.yaml b/src/methods/cycombine_all_controls_to_goal/config.vsh.yaml
index b567f1e2..446129cb 100644
--- a/src/methods/cycombine_all_controls_to_goal/config.vsh.yaml
+++ b/src/methods/cycombine_all_controls_to_goal/config.vsh.yaml
@@ -71,6 +71,7 @@ engines:
     setup:
       - type: r
         bioc: [sva]
+        packages: [pbmcapply]
         github: [biosurf/cyCombine]
 
 runners:
diff --git a/src/methods/cycombine_all_controls_to_mid/config.vsh.yaml b/src/methods/cycombine_all_controls_to_mid/config.vsh.yaml
index d669ef4f..c16443ba 100644
--- a/src/methods/cycombine_all_controls_to_mid/config.vsh.yaml
+++ b/src/methods/cycombine_all_controls_to_mid/config.vsh.yaml
@@ -71,6 +71,7 @@ engines:
     setup:
       - type: r
         bioc: [sva]
+        packages: [pbmcapply]
         github: [biosurf/cyCombine]
 
 runners:
diff --git a/src/methods/cycombine_no_controls_to_goal/config.vsh.yaml b/src/methods/cycombine_no_controls_to_goal/config.vsh.yaml
index 1588843f..3ed0f67e 100644
--- a/src/methods/cycombine_no_controls_to_goal/config.vsh.yaml
+++ b/src/methods/cycombine_no_controls_to_goal/config.vsh.yaml
@@ -60,6 +60,7 @@ engines:
     setup:
       - type: r
         bioc: [sva]
+        packages: [pbmcapply]
         github: [biosurf/cyCombine]
 
 runners:
diff --git a/src/methods/cycombine_no_controls_to_mid/config.vsh.yaml b/src/methods/cycombine_no_controls_to_mid/config.vsh.yaml
index c9c90f3c..6fe321e0 100644
--- a/src/methods/cycombine_no_controls_to_mid/config.vsh.yaml
+++ b/src/methods/cycombine_no_controls_to_mid/config.vsh.yaml
@@ -60,6 +60,7 @@ engines:
     setup:
       - type: r
         bioc: [sva]
+        packages: [pbmcapply]
         github: [biosurf/cyCombine]
 
 runners:
diff --git a/src/methods/cycombine_one_control_to_goal/config.vsh.yaml b/src/methods/cycombine_one_control_to_goal/config.vsh.yaml
index 93938927..79c65f8e 100644
--- a/src/methods/cycombine_one_control_to_goal/config.vsh.yaml
+++ b/src/methods/cycombine_one_control_to_goal/config.vsh.yaml
@@ -70,6 +70,7 @@ engines:
     setup:
       - type: r
         bioc: [sva]
+        packages: [pbmcapply]
         github: [biosurf/cyCombine]
 
 runners:
diff --git a/src/methods/cycombine_one_control_to_mid/config.vsh.yaml b/src/methods/cycombine_one_control_to_mid/config.vsh.yaml
index 96083a91..21eb4566 100644
--- a/src/methods/cycombine_one_control_to_mid/config.vsh.yaml
+++ b/src/methods/cycombine_one_control_to_mid/config.vsh.yaml
@@ -70,6 +70,7 @@ engines:
     setup:
       - type: r
         bioc: [sva]
+        packages: [pbmcapply]
         github: [biosurf/cyCombine]
 
 runners:
diff --git a/src/methods/cytovi/config.vsh.yaml b/src/methods/cytovi/config.vsh.yaml
new file mode 100644
index 00000000..91095326
--- /dev/null
+++ b/src/methods/cytovi/config.vsh.yaml
@@ -0,0 +1,84 @@
+# 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: cytovi
+# A relatively short label, used when rendering visualisations (required)
+label: CytoVI
+# A one sentence summary of how this method works (required). Used when 
+# rendering summary tables.
+summary: "A deep generative model for correcting batch effects"
+# A multi-line description of how this component works (required). Used
+# when rendering reference documentation.
+description: |
+  CytoVI is a deep generative model that utilizes antibody-based single-cell profiles to 
+  learn a biologically meaningful latent representation of each cell. 
+  It is part of the scvi-tools framework and is built upon the variational autoencoder (VAE) architecture.
+references:
+  doi: 
+    - 10.1101/2025.09.07.674699
+links:
+  # URL to the documentation for this method (required).
+  documentation: https://docs.scvi-tools.org/en/latest/user_guide/models/cytovi.html
+  # URL to the code repository for this method (required).
+  repository: https://github.com/YosefLab/cytovi-reference-implementation
+
+# Component-specific parameters (optional)
+arguments:
+  - name: --n_hidden
+    type: integer
+    default: 128
+    description: Number of hidden units.
+  - name: --n_layers
+    type: integer
+    default: 1
+    description: Number of layers.
+  - name: --n_clusters
+    type: integer
+    default: 20
+    description: Number of clusters to use for subsampling.
+  - name: --subsample_fraction
+    type: double
+    default: 0.5
+    description: Fraction of cells to subsample from each cluster for training.
+  
+# Resources required to run the component
+resources:
+  # The script of your component (required)
+  - type: python_script
+    path: script.py
+  # Additional resources your script needs (optional)
+  # - type: file
+  #   path: weights.pt
+
+engines:
+  - type: docker
+    image: nvcr.io/nvidia/pytorch:25.08-py3
+    setup:
+      - type: apt
+        packages:
+          - procps
+          - git
+      - type: python
+        packages:
+          - anndata>=0.11.0
+          - scanpy[skmisc]>=1.10
+          - scvi-tools==1.4.0
+          - pyyaml
+          - requests
+          - jsonschema
+        github:
+          - openproblems-bio/core#subdirectory=packages/python/openproblems
+
+runners:
+  # This platform allows running the component natively
+  - type: executable
+  # Allows turning the component into a Nextflow module / pipeline.
+  - type: nextflow
+    directives:
+      label: [veryhightime, lowmem, lowcpu, gpu]
diff --git a/src/methods/cytovi/script.py b/src/methods/cytovi/script.py
new file mode 100644
index 00000000..6a6fd73f
--- /dev/null
+++ b/src/methods/cytovi/script.py
@@ -0,0 +1,112 @@
+import anndata as ad
+import numpy as np
+import scanpy as sc
+from scvi.external import cytovi
+from sklearn.cluster import KMeans
+
+## VIASH START
+par = {
+    "input": "resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/censored_split2.h5ad",
+    "output": "resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/output_cytovi_split2.h5ad",
+    "n_hidden": 128,
+    "n_layers": 1,
+    "n_clusters": 10,
+    "subsample_fraction": 0.5,
+}
+meta = {"name": "cytovi"}
+## VIASH END
+
+print("Reading and preparing input files", flush=True)
+adata = ad.read_h5ad(par["input"])
+
+adata.obs["batch_str"] = adata.obs["batch"].astype(str)
+
+markers_to_correct = adata.var[adata.var["to_correct"]].index.to_numpy()
+markers_not_correct = adata.var[~adata.var["to_correct"]].index.to_numpy()
+
+adata_to_correct = adata[:, markers_to_correct].copy()
+
+print("Scaling data", flush=True)
+
+# scale data. this will add a layer "scaled" to the anndata
+cytovi.scale(
+    adata=adata_to_correct,
+    transformed_layer_key="preprocessed",
+    batch_key="batch_str",
+    inplace=True,
+)
+
+print("Clustering using k-means with k =", par["n_clusters"], flush=True)
+# cluster data using Kmeans
+adata_to_correct.obs["clusters"] = (
+    KMeans(n_clusters=par["n_clusters"], random_state=0)
+    .fit_predict(adata_to_correct.layers["scaled"])
+    .astype(str)
+)
+# concatenate obs so we can use it for subsampling
+adata_to_correct.obs["sample_cluster"] = (
+    adata_to_correct.obs["sample"].astype(str) + "_" + adata_to_correct.obs["clusters"]
+)
+# subsample cells without replacement
+print("Subsampling cells", flush=True)
+subsampled_cells = adata_to_correct.obs.groupby("sample_cluster")[
+    "sample_cluster"
+].apply(lambda x: x.sample(n=round(len(x) * par["subsample_fraction"]), replace=False))
+# need the cell id included in the subsample
+subsampled_cells_idx = [x[1] for x in subsampled_cells.index.to_list()]
+
+adata_subsampled = adata_to_correct[subsampled_cells_idx, :].copy()
+
+print(
+    f"Train CytoVI on subsampled data containing {adata_subsampled.shape[0]} cells",
+    flush=True,
+)
+
+cytovi.CYTOVI.setup_anndata(adata_subsampled, layer="scaled", batch_key="batch_str")
+model = cytovi.CYTOVI(
+    adata=adata_subsampled, n_hidden=par["n_hidden"], n_layers=par["n_layers"]
+)
+model.train()
+
+# get batch corrected data
+print("Correcting data", flush=True)
+corrected_data = model.get_normalized_expression(adata=adata_to_correct)
+
+# have to add in the uncorrected markers as well
+uncorrected_data = adata[:, markers_not_correct].layers["preprocessed"]
+
+out_matrix = np.concatenate([corrected_data, uncorrected_data], axis=1)
+out_var_idx = np.concatenate([corrected_data.columns, markers_not_correct])
+
+# create new anndata
+out_adata = ad.AnnData(
+    obs=adata.obs[[]],
+    var=adata.var.loc[out_var_idx][[]],
+    layers={"integrated": out_matrix},
+    uns={
+        "dataset_id": adata.uns["dataset_id"],
+        "method_id": meta["name"],
+        "parameters": {},
+    },
+)
+
+# reorder var to match input
+out_adata = out_adata[:, adata.var_names]
+
+# leave this here for debugging purposes
+# run umap for quick check
+# import scanpy as sc
+
+# test_adata = ad.AnnData(
+#     X=out_adata.layers["integrated"].toarray(),
+#     obs=adata.obs,
+#     var=adata.var,
+# )
+# test_adata = test_adata[:, markers_to_correct]
+# sc.pp.neighbors(test_adata, use_rep="X")
+# sc.tl.umap(test_adata)
+# sc.pl.umap(test_adata, color="batch")
+
+print("Write output AnnData to file", flush=True)
+
+out_adata.write_h5ad(par["output"], compression="gzip")
diff --git a/src/methods/mnn/config.vsh.yaml b/src/methods/mnn/config.vsh.yaml
deleted file mode 100644
index e601c4bb..00000000
--- a/src/methods/mnn/config.vsh.yaml
+++ /dev/null
@@ -1,107 +0,0 @@
-# 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: mnn
-# A relatively short label, used when rendering visualisations (required)
-label: mnn
-# A one sentence summary of how this method works (required). Used when 
-# rendering summary tables.
-summary: "Original implementation of the Mutual Nearest Neighbors (MNN) algorithm by Haghverdi et al."
-# A multi-line description of how this component works (required). Used
-# when rendering reference documentation.
-description: |
-  Correct for batch effect using mutual nearest neighbors (MNNs).
-  Mutual nearest neighbors are identified as pairs of cells from different batches that are
-  within each other's set of k nearest neighbors. Cell pairs that have been identified as MNNs are assumed to
-  belong to the same biological state (or cell type), so that differences in the expression profiles of MNNs can
-  be used to infer direction and degree of batch effect.
-  The method was originally developed for single-cell RNA-sequencing data, but it is often used
-  for other types of data as well. 
-  This component uses the implementation from the `batchelor` bioconductor package.
-
-references:
-  doi: 
-    - 10.1038/nbt.4091
-  bibtex:
-    - |
-      @article{haghverdi2018batch,
-        title={Batch effects in single-cell RNA-sequencing data are corrected by matching mutual nearest neighbors},
-        author={Haghverdi, Laleh and Lun, Aaron TL and Morgan, Michael D and Marioni, John C},
-        journal={Nature biotechnology},
-        volume={36},
-        number={5},
-        pages={421--427},
-        year={2018},
-        publisher={Nature Publishing Group}
-      }
-links:
-  # URL to the documentation for this method (required).
-  documentation: https://bioconductor.org/packages/devel/bioc/html/batchelor.html
-  # URL to the code repository for this method (required).
-  repository: https://github.com/LTLA/batchelor/tree/master
-
-
-# Component-specific parameters (optional)
-argument_groups:
-  - name: "Parameters"
-    arguments:
-      - name: "--num_nn"
-        type: integer
-        info:
-          optimize:
-          type: linear
-          lower: 10
-          upper: 100
-        default: 20
-        description: An integer scalar specifying the number of nearest neighbors to consider when identifying MNNs.
-      - name: "--prop_num_nn"
-        type: double
-        info:
-          optimize:
-            type: expuniform
-            lower: 0
-            upper: 1
-        default: 0
-        description: A numeric scalar in (0, 1) specifying the proportion of cells in each dataset to use for mutual nearest neighbor searching. If set, the number of nearest neighbors used for the MNN search in each batch is redefined as max(k, prop.k*N) where N is the number of cells in that batch. If 0 is set to NULL.
-      - name: "--sigma_value"
-        type: double
-        info:
-          optimize:
-            type: expuniform
-            lower: 0.01
-            upper: 100
-        default: 0.1
-        description: A numeric scalar specifying the bandwidth of the Gaussian smoothing kernel used to compute the correction vector for each cell.
-
-# Resources required to run the component
-resources:
-  # The script of your component (required)
-  - type: r_script
-    path: script.R
-  # Additional resources your script needs (optional)
-  # - type: file
-  #   path: weights.pt
-
-engines:
-  # Specifications for the Docker image for this component.
-  - type: docker
-    image: openproblems/base_r:1
-    # Add custom dependencies here (optional). For more information, see
-    # https://viash.io/reference/config/engines/docker/#setup .
-    setup:
-      - type: r
-        bioc: [batchelor]
-
-runners:
-  # This platform allows running the component natively
-  - type: executable
-  # Allows turning the component into a Nextflow module / pipeline.
-  - type: nextflow
-    directives:
-      label: [midtime,midmem,midcpu]
diff --git a/src/methods/mnn/script.R b/src/methods/mnn/script.R
deleted file mode 100644
index 0cfed88e..00000000
--- a/src/methods/mnn/script.R
+++ /dev/null
@@ -1,57 +0,0 @@
-library(anndata)
-library(batchelor)
-
-## VIASH START
-par <- list(
-  input = "resources_test/.../input.h5ad",
-  output = "output.h5ad"
-)
-meta <- list(
-  name = "mnn"
-)
-## VIASH END
-
-cat("Reading input files\n")
-input <- anndata::read_h5ad(par[["input"]])
-
-cat("Subset data\n")
-data_not_correct <- input[, !input$var$to_correct]
-data_to_correct <- input[, input$var$to_correct]
-
-cat("Run MNN\n")
-
-# If prop_num_nn is set to 0, it is set to NULL
-if (par[["prop_num_nn"]] == 0.0) {
-  print("'prop_num_nn' = 0, setting it to NULL")
-  par[["prop_num_nn"]] <- NULL
-}
-
-corrected_data <- mnnCorrect(Matrix::t(data_to_correct$layers[["preprocessed"]]),
-              batch = data_to_correct$obs$batch,
-              k= par[["num_nn"]],
-              prop.k = par[["prop_num_nn"]],
-              sigma = par[["sigma_value"]],
-              cos.norm.in = FALSE,
-              cos.norm.out = FALSE)
-
-cat("Preparing output Anndata\n")
-corrected_data <- Matrix::t(assay(corrected_data))
-corrected_data <- cbind(corrected_data, data_not_correct$layers[["preprocessed"]])
-
-cat("Write output AnnData to file\n")
-output <- anndata::AnnData(
-  obs = input$obs[, integer(0)],
-  var = input$var[colnames(corrected_data), integer(0)],
-  layers = list(integrated = corrected_data),
-  uns = list(
-    dataset_id = input$uns$dataset_id,
-    method_id = meta$name,
-    parameters = list(
-      "k" = par[["num_nn"]],
-      "prop.k" = par["prop_num_nn"],
-      "sigma" = par[["sigma_value"]]
-    )
-  )
-)
-
-output$write_h5ad(par[["output"]], compression = "gzip")
diff --git a/src/methods/rpca_to_goal/config.vsh.yaml b/src/methods/rpca_to_goal/config.vsh.yaml
index 88153b9f..200dfa68 100644
--- a/src/methods/rpca_to_goal/config.vsh.yaml
+++ b/src/methods/rpca_to_goal/config.vsh.yaml
@@ -1,22 +1,8 @@
-# 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: rpca_to_goal
-# A relatively short label, used when rendering visualisations (required)
+# status: disabled
 label: Seurat RPCA (to-goal)
-
-# A one sentence summary of how this method works (required). Used when 
-# rendering summary tables.
 summary: "Batch integrate data to a goal batch using mutual nearest neighbors identified via Seurat reciprocal PCA."
-
-# A multi-line description of how this component works (required). Used
-# when rendering reference documentation.
 description: |
   Seurat RPCA performs batch integration by projecting each query dataset into the PCA space 
   of a goal batch, and identifying anchors using reciprocal PCA (RPCA). 
@@ -48,9 +34,7 @@ references:
       }
 
 links:
-  # URL to the documentation for this method (required).
   documentation: https://satijalab.org/seurat/articles/integration_rpca.html
-  # URL to the code repository for this method (required).
   repository: https://github.com/satijalab/seurat
 
 argument_groups:
@@ -100,4 +84,4 @@ runners:
   # Allows turning the component into a Nextflow module / pipeline.
   - type: nextflow
     directives:
-      label: [midtime,midmem,midcpu]
+      label: [veryhightime,midmem,midcpu]
diff --git a/src/methods/rpca_to_goal/script.R b/src/methods/rpca_to_goal/script.R
index 8a9fd908..35f562f5 100644
--- a/src/methods/rpca_to_goal/script.R
+++ b/src/methods/rpca_to_goal/script.R
@@ -13,6 +13,8 @@ meta <- list(
 )
 ## VIASH END
 
+options(future.globals.maxSize = 25 * 1024^3)  # 25 GiB
+
 cat("Reading input files\n")
 input_adata <- anndata::read_h5ad(par[["input"]])
 
@@ -59,7 +61,7 @@ seurat_objs <- lapply(batches, function(batch) {
         object = seurat_obj,
         features = markers_to_correct,
         assay = "cyto",
-        verbose = FALSE
+        verbose = TRUE
     )
 
     # run pca. mandatory
@@ -72,7 +74,7 @@ seurat_objs <- lapply(batches, function(batch) {
         assay = "cyto",
         npcs = par[["npcs"]],
         approx = FALSE,
-        verbose = FALSE
+        verbose = TRUE
     )
 
     return(seurat_obj)
@@ -94,7 +96,7 @@ anchors <- Seurat::FindIntegrationAnchors(
     dims = seq(npcs_computed),
     k.anchor = par[["n_neighbours"]],
     reduction = "rpca",
-    verbose = FALSE,
+    verbose = TRUE,
     reference = which(names(seurat_objs) == "1")
 )
 
@@ -108,7 +110,7 @@ batch_corrected_seurat_obj <- Seurat::IntegrateData(
     features = markers_to_correct,
     features.to.integrate = markers_to_correct,
     dims = seq(npcs_computed),
-    verbose = FALSE
+    verbose = TRUE
 )
 # just to be sure!
 Seurat::DefaultAssay(batch_corrected_seurat_obj) <- "integrated"
diff --git a/src/methods/rpca_to_mid/config.vsh.yaml b/src/methods/rpca_to_mid/config.vsh.yaml
index 8f4f26dd..ca01b9cb 100644
--- a/src/methods/rpca_to_mid/config.vsh.yaml
+++ b/src/methods/rpca_to_mid/config.vsh.yaml
@@ -1,22 +1,8 @@
-# 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: rpca_to_mid
-# A relatively short label, used when rendering visualisations (required)
+# status: disabled
 label: Seurat RPCA (to-middle)
-
-# A one sentence summary of how this method works (required). Used when 
-# rendering summary tables.
 summary: "Batch integrate data to a midpoint using mutual nearest neighbors identified via Seurat reciprocal PCA."
-
-# A multi-line description of how this component works (required). Used
-# when rendering reference documentation.
 description: |
   Seurat RPCA performs batch integration by identifying mutual nearest neighbors (anchors) 
   between all batches using reciprocal PCA (RPCA). 
@@ -102,4 +88,4 @@ runners:
   # Allows turning the component into a Nextflow module / pipeline.
   - type: nextflow
     directives:
-      label: [midtime,midmem,midcpu]
+      label: [veryhightime,midmem,midcpu]
diff --git a/src/methods/rpca_to_mid/script.R b/src/methods/rpca_to_mid/script.R
index d84d26e0..7451eb6c 100644
--- a/src/methods/rpca_to_mid/script.R
+++ b/src/methods/rpca_to_mid/script.R
@@ -13,6 +13,9 @@ meta <- list(
 )
 ## VIASH END
 
+options(future.globals.maxSize = 25 * 1024^3)  # 25 GiB
+
+
 cat("Reading input files\n")
 input_adata <- anndata::read_h5ad(par[["input"]])
 
@@ -59,7 +62,7 @@ seurat_objs <- lapply(batches, function(batch) {
         object = seurat_obj,
         features = markers_to_correct,
         assay = "cyto",
-        verbose = FALSE
+        verbose = TRUE
     )
 
     # run pca. mandatory
@@ -72,7 +75,7 @@ seurat_objs <- lapply(batches, function(batch) {
         assay = "cyto",
         npcs = par[["npcs"]],
         approx = FALSE,
-        verbose = FALSE
+        verbose = TRUE
     )
 
     return(seurat_obj)
@@ -96,7 +99,7 @@ anchors <- Seurat::FindIntegrationAnchors(
     dims = seq(npcs_computed),
     k.anchor = par[["n_neighbours"]],
     reduction = "rpca",
-    verbose = FALSE,
+    verbose = TRUE,
     reference = NULL
 )
 
@@ -106,7 +109,7 @@ batch_corrected_seurat_obj <- Seurat::IntegrateData(
     anchorset = anchors,
     features.to.integrate = markers_to_correct,
     dims = seq(npcs_computed),
-    verbose = FALSE,
+    verbose = TRUE,
 )
 # just to be sure!
 Seurat::DefaultAssay(batch_corrected_seurat_obj) <- "integrated"
diff --git a/src/metrics/bras/script.py b/src/metrics/bras/script.py
index a7e88c75..98229ad0 100644
--- a/src/metrics/bras/script.py
+++ b/src/metrics/bras/script.py
@@ -1,20 +1,23 @@
+import sys
+
 import anndata as ad
 import numpy as np
-import sys
 from scib_metrics import bras
 
 ## 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`.
+# The following code has been auto-generated by Viash.
 par = {
-  'input_validation': 'resources_test/.../validation.h5ad',
-  'input_unintegrated': 'resources_test/.../unintegrated.h5ad',
-  'input_integrated': 'resources_test/.../integrated.h5ad',
-  'output': 'output.h5ad'
+    "input_unintegrated": r"resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/unintegrated.h5ad",
+    "input_integrated_split1": r"resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/integrated_split1.h5ad",
+    "input_integrated_split2": r"resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/integrated_split2.h5ad",
+    "output": r"resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/harmonypy_bras_score.h5ad",
 }
 meta = {
-  'name': 'bras'
+    "name": r"bras",
+    "functionality_name": r"bras",
+    "resources_dir": r"src/utils",
 }
+
 ## VIASH END
 
 sys.path.append(meta["resources_dir"])
@@ -45,25 +48,25 @@
 integrated_s2 = subset_nocontrols(integrated_s2)
 integrated_s2 = remove_unlabelled(integrated_s2)
 
-print('Compute metrics', flush=True)
-batch_labels_s1 = integrated_s1.obs['batch'].values
-ct_labels_s1 = integrated_s1.obs['cell_type'].values
+print("Compute metrics", flush=True)
+batch_labels_s1 = integrated_s1.obs["batch"].values
+ct_labels_s1 = integrated_s1.obs["cell_type"].values
 
 bras_s1 = bras(
     integrated_s1.layers["integrated"],
-    labels = batch_labels_s1,
-    batch = ct_labels_s1,
-    metric='cosine'
+    labels=ct_labels_s1,
+    batch=batch_labels_s1,
+    metric="euclidean",
 )
 
-batch_labels_s2 = integrated_s2.obs['batch'].values
-ct_labels_s2 = integrated_s2.obs['cell_type'].values
+batch_labels_s2 = integrated_s2.obs["batch"].values
+ct_labels_s2 = integrated_s2.obs["cell_type"].values
 
 bras_s2 = bras(
     integrated_s2.layers["integrated"],
-    labels = batch_labels_s2,
-    batch = ct_labels_s2,
-    metric='cosine'
+    labels=ct_labels_s2,
+    batch=batch_labels_s2,
+    metric="euclidean",
 )
 
 bras_score = np.mean([bras_s1, bras_s2])
@@ -71,13 +74,12 @@
 print("Write output AnnData to file", flush=True)
 output = ad.AnnData(
     uns={
-    'dataset_id': integrated_s1.uns['dataset_id'],
-    'method_id': integrated_s1.uns['method_id'],
-    'metric_ids': ['bras'],
-    'metric_values': [bras_score],
-    'bras_s1': bras_s1,
-    'bras_s2': bras_s2
-  }
+        "dataset_id": integrated_s1.uns["dataset_id"],
+        "method_id": integrated_s1.uns["method_id"],
+        "metric_ids": ["bras"],
+        "metric_values": [bras_score],
+        "bras_s1": bras_s1,
+        "bras_s2": bras_s2,
+    }
 )
-output.write_h5ad(par['output'], compression='gzip')
-
+output.write_h5ad(par["output"], compression="gzip")
diff --git a/src/metrics/cms/config.vsh.yaml b/src/metrics/cms/config.vsh.yaml
deleted file mode 100644
index 92b9336c..00000000
--- a/src/metrics/cms/config.vsh.yaml
+++ /dev/null
@@ -1,118 +0,0 @@
-__merge__: ../../api/comp_metric.yaml
-
-# A unique identifier for your component (required).
-# Can contain only lowercase letters or underscores.
-name: cms
-
-# Metadata for your component
-info:
-  metrics:
-      # A unique identifier for your metric (required).
-      # Can contain only lowercase letters or underscores.
-    - name: cms
-      # A relatively short label, used when rendering visualisarions (required)
-      label: Cell Mixing Score
-      # A one sentence summary of how this metric works (required). Used when 
-      # rendering summary tables.
-      summary: "Cellspecific Mixing Score (cms) quantifies batch effects at the cell level by computing batch-specific distance distributions towards k-nearest neighbouring cells."
-      # A multi-line description of how this component works (required). Used
-      # when rendering reference documentation.
-      description: |
-        The cellspecific mixing score cms tests for each cell the hypothesis that batch-specific distance 
-        distributions towards it's k-nearest neighbouring (knn) cells are derived from the same unspecified 
-        underlying distribution using the Anderson-Darling test. The test considers differences in the number of cells 
-        from each batch, making the cms a robust metric when evaluating batch effects in datasets with unbalanced batch sizes.
-        This implementation uses k = 50 knn cells and the first 10 principal components for distance computations.
-
-        The cms for each cell can be interpreted as a P-value, that is, the probability of observing deviations 
-        in the batch specific distance distributions by chance (assuming they are all derived from the same distribution).
-        Therefore, for a given cell:
-        - A low cms score indicates that the batch-specific distance distributions towards its knn cells are significantly different, 
-          suggesting that the cell is influenced by batch effects.
-        - A high cms score indicates that the batch-specific distance distributions towards its knn cells are similar, 
-          suggesting that the cell is not influenced by batch effects.
-        
-
-        To characterize the overall batch mixing in a dataset, we use the medcouple statistic on the distribution of cms scores.
-        The medcouple is a robust measure of skewness, which is less sensitive to outliers than the traditional skewness measure.
-        The medcouple statistic returns values between -1 and 1, where:
-        - A value close to -1 indicates a pronounced left-skewed distribution, reflecting abundance of cells with high cms scores
-        - A value close to 1 indicates a pronounced right-skewed distribution, reflecting abundance of cells with low cms scores
-        - A value around 0 indicates a symmetrical distribution of cms scores
-
-        It has been empirically observed that a uniform (thus, symmetrical) distribution of cms scores across cells in a dataset is indicative of good mixing
-        (e.g. via random shuffling of batch labels in a dataset). Therefore, a medcouple around 0 or lower is considered a good mixing score.
-
-      references:
-        doi: 
-          - 10.26508/lsa.202001004
-        bibtex:
-          - |
-            @article{lutge2021cellmixs,
-              title={CellMixS: quantifying and visualizing batch effects in single-cell RNA-seq data},
-              author={L{\"u}tge, Almut and Zyprych-Walczak, Joanna and Kunzmann, Urszula Brykczynska and Crowell, Helena L and Calini, Daniela and Malhotra, Dheeraj and Soneson, Charlotte and Robinson, Mark D},
-              journal={Life science alliance},
-              volume={4},
-              number={6},
-              year={2021},
-              publisher={Life Science Alliance}
-            }
-          - |
-            @article{brys2004robust,
-              title={A robust measure of skewness},
-              author={Brys, Guy and Hubert, Mia and Struyf, Anja},
-              journal={Journal of Computational and Graphical Statistics},
-              volume={13},
-              number={4},
-              pages={996--1017},
-              year={2004},
-              publisher={Taylor \& Francis}
-            }
-
-      links:
-        # URL to the documentation for this metric (required).
-        documentation: https://bioconductor.org/packages/release/bioc/html/CellMixS.html
-        # URL to the code repository for this metric (required).
-        repository: https://github.com/almutlue/CellMixS
-      # The minimum possible value for this metric (required)
-      min: -1
-      # The maximum possible value for this metric (required)
-      max: 1
-      # Whether a higher value represents a 'better' solution (required)
-      maximize: false
-
-# Component-specific parameters (optional)
-arguments:
-  - name: "--n_neighbors"
-    type: "integer"
-    default: 50
-    description: Number of k-nearest neighbours (knn) to use when computing cms scores.
-  - name: "--n_dim"
-    type: "integer"
-    default: 10
-    description: Number of principal components to use when computing cms scores.
-
-# Resources required to run the component
-resources:
-  # The script of your component (required)
-  - type: r_script
-    path: script.R
-  - path: /src/utils/helper_functions.R
-
-engines:
-  # Specifications for the Docker image for this component.
-  - type: docker
-    image: openproblems/base_r:1
-    setup:
-      - type: r
-        packages: [robustbase, rhdf5]
-        bioc: [CellMixS, SingleCellExperiment, Biocparallel]
-        github: [scverse/anndataR@0.99.0]
-
-runners:
-  # This platform allows running the component natively
-  - type: executable
-  # Allows turning the component into a Nextflow module / pipeline.
-  - type: nextflow
-    directives:
-      label: [midtime,midmem,midcpu]
diff --git a/src/metrics/cms/script.R b/src/metrics/cms/script.R
deleted file mode 100644
index b05dc8cc..00000000
--- a/src/metrics/cms/script.R
+++ /dev/null
@@ -1,124 +0,0 @@
-requireNamespace("anndataR", quietly = TRUE)
-requireNamespace("CellMixS", quietly = TRUE)
-requireNamespace("BiocParallel", quietly = TRUE)
-requireNamespace("robustbase", quietly = TRUE)
-requireNamespace("parallel", quietly = TRUE)
-requireNamespace("SingleCellExperiment", quietly = TRUE)
-
-## VIASH START
-# The following code has been auto-generated by Viash.
-# treat warnings as errors
-
-par <- list(
-    input_unintegrated = "resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/unintegrated.h5ad",
-    input_integrated_split1 = "resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/integrated_split1.h5ad",
-    input_integrated_split2 = "resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/integrated_split2.h5ad",
-    output = "resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/score_cms.h5ad",
-    n_neighbors = as.integer("50"),
-    n_dim = as.integer("10")
-)
-meta <- list(
-    name = "cms",
-    resources_dir = "src/utils",
-    cpus = NULL
-)
-## VIASH END
-
-cores_to_use <- meta$cpus
-if (is.null(cores_to_use)) {
-    cores_to_use <- min(5, parallel::detectCores() - 2)
-}
-
-source(paste0(meta$resources_dir, "/helper_functions.R"))
-
-
-cat("Reading input files\n")
-unintegrated <- anndataR::read_h5ad(par[["input_unintegrated"]])
-integrated_split1 <- anndataR::read_h5ad(par[["input_integrated_split1"]])
-integrated_split2 <- anndataR::read_h5ad(par[["input_integrated_split2"]])
-
-cat("Fetching some metadata from unintegrated\n")
-integrated_split1 <- get_obs_var_for_integrated(
-    i_adata = integrated_split1,
-    u_adata = unintegrated
-)
-integrated_split2 <- get_obs_var_for_integrated(
-    i_adata = integrated_split2,
-    u_adata = unintegrated
-)
-
-# Fetch batch annotations from unintegrated
-# batch_key <- input_unintegrated$obs$batch
-
-# Get markers to correct
-markers_to_correct <- unintegrated$var_names[unintegrated$var$to_correct]
-
-cat("Converting to SingleCellExperiment object\n")
-
-# Convert to SingleCellExperiment
-integrated_split1_sce <- integrated_split1$as_SingleCellExperiment()
-integrated_split1_sce <- integrated_split1_sce[markers_to_correct, ]
-
-integrated_split2_sce <- integrated_split2$as_SingleCellExperiment()
-integrated_split2_sce <- integrated_split2_sce[markers_to_correct, ]
-
-# cores_to_use <- 5
-bpparam <- BiocParallel::MulticoreParam(
-    workers = cores_to_use
-)
-cat(paste("Compute Cell Mixing Score using", cores_to_use, "cores for split 1\n"))
-
-integrated_split1_sce <- CellMixS::cms(
-    integrated_split1_sce,
-    group = "batch",
-    assay_name = "integrated",
-    k = par[["n_neighbors"]],
-    n_dim = par[["n_dim"]],
-    BPPARAM = bpparam
-)
-
-cat(paste("Compute Cell Mixing Score using", cores_to_use, "cores for split 2\n"))
-
-integrated_split2_sce <- CellMixS::cms(
-    integrated_split2_sce,
-    group = "batch",
-    assay_name = "integrated",
-    k = par[["n_neighbors"]],
-    n_dim = par[["n_dim"]],
-    BPPARAM = bpparam
-)
-
-cat("Compute Medcouple statistic\n")
-
-cms_distr_split1 <- SingleCellExperiment::colData(integrated_split1_sce)[, "cms"]
-cms_distr_split2 <- SingleCellExperiment::colData(integrated_split2_sce)[, "cms"]
-
-cms_mc_split1 <- robustbase::mc(cms_distr_split1)
-cms_mc_split2 <- robustbase::mc(cms_distr_split2)
-
-cms_mc_mean <- mean(c(cms_mc_split1, cms_mc_split2))
-
-cat("Write output AnnData to file\n")
-output <- anndataR::AnnData(
-    shape = c(0L, 0L),
-    uns = list(
-        dataset_id = integrated_split1$uns$dataset_id,
-        method_id = integrated_split1$uns$method_id,
-        metric_ids = meta$name,
-        metric_values = cms_mc_mean,
-        cms_parameters = list(
-            n_neighbors = par[["n_neighbors"]],
-            n_dim = par[["n_dim"]]
-        ),
-        cms_medcouple_score = list(
-            left = cms_mc_split1,
-            right = cms_mc_split2
-        ),
-        cms_distribution = list(
-            left = cms_distr_split1,
-            right = cms_distr_split2
-        )
-    )
-)
-
-output$write_h5ad(par[["output"]], compression = "gzip", mode = "w")
diff --git a/src/metrics/emd/config.vsh.yaml b/src/metrics/emd/config.vsh.yaml
index 8963c9b8..72cef341 100644
--- a/src/metrics/emd/config.vsh.yaml
+++ b/src/metrics/emd/config.vsh.yaml
@@ -99,178 +99,6 @@ info:
       max: .inf
       # Whether a higher value represents a 'better' solution (required)
       maximize: false
-      
-      # A unique identifier for your metric (required).
-      # Can contain only lowercase letters or underscores.
-    - name: emd_mean_global_horiz
-      # A relatively short label, used when rendering visualisarions (required)
-      label: EMD Mean Global Horizontal
-      # A one sentence summary of how this metric works (required). Used when 
-      # rendering summary tables.
-      summary: "Mean Earth Mover Distance calculated horizontally across donors for each marker."
-      # A multi-line description of how this component works (required). Used
-      # when rendering reference documentation.
-      description: |
-        Earth Mover Distance (EMD), also known as the Wasserstein metric, measures the difference 
-        between two probability distributions. 
-        
-        Here, EMD is used to compare marker expression distributions between paired samples from the same donor 
-        quantified across two different batches. 
-        For each paired sample and marker, the marker expression values are first converted into 
-        probability distributions. 
-        This is done by binning the expression values into a range from -100 to 100 with a bin width of 0.1.
-        The `wasserstein_distance` function from SciPy is then used to calculate the EMD between the two 
-        probability distributions belonging to the same cell type, marker, and a given paired samples.
-        This is then repeated for every marker and paired sample.
-        Finally, the average of all these EMD values is computed and reported as the metric score.
-
-        The key difference between this and `emd_mean_ct_horiz` is that the EMD values are
-        computed agnostic of cell types.
-
-        A high score indicates that at least one marker and cell type in a given sample pair has a 
-        large difference in distribution after batch integration.
-        A low score means that the most poorly corrected marker expression is well integrated across batches.   
-      references:
-        doi: 
-          - 10.1023/A:1026543900054
-      links:
-        # URL to the documentation for this metric (required).
-        documentation: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.wasserstein_distance.html
-        # URL to the code repository for this metric (required).
-        repository: https://github.com/scipy/scipy
-      # The minimum possible value for this metric (required)
-      min: 0
-      # The maximum possible value for this metric (required)
-      max: .inf
-      # Whether a higher value represents a 'better' solution (required)
-      maximize: false
-      
-      # A unique identifier for your metric (required).
-      # Can contain only lowercase letters or underscores.
-    - name: emd_max_global_horiz
-      # A relatively short label, used when rendering visualisarions (required)
-      label: EMD Max Global Horizontal
-      # A one sentence summary of how this metric works (required). Used when 
-      # rendering summary tables.
-      summary: "Max Earth Mover Distance calculated horizontally across donors for each marker."
-      # A multi-line description of how this component works (required). Used
-      # when rendering reference documentation.
-      description: |
-        Earth Mover Distance (EMD), also known as the Wasserstein metric, measures the difference 
-        between two probability distributions. 
-        
-        Here, EMD is used to compare marker expression distributions between paired samples from the same donor 
-        quantified across two different batches. 
-        For each paired sample and marker, the marker expression values are first converted into 
-        probability distributions. 
-        This is done by binning the expression values into a range from -100 to 100 with a bin width of 0.1.
-        The `wasserstein_distance` function from SciPy is then used to calculate the EMD between the two 
-        probability distributions belonging to the same cell type, marker, and a given paired samples.
-        This is then repeated for every cell type, marker, and paired sample.
-        Finally, the maximum of all these EMD values is computed and reported as the metric score.
-
-        The key difference between this and `emd_max_ct_horiz` is that the EMD values are
-        computed agnostic of cell types.
-
-        A high score indicates that at least one marker in a given sample pair has a large difference in 
-        distribution after batch integration.
-        A low score means that the most poorly corrected marker expression is well integrated across batches.   
-      references:
-        doi: 
-          - 10.1023/A:1026543900054
-      links:
-        # URL to the documentation for this metric (required).
-        documentation: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.wasserstein_distance.html
-        # URL to the code repository for this metric (required).
-        repository: https://github.com/scipy/scipy
-      # The minimum possible value for this metric (required)
-      min: 0
-      # The maximum possible value for this metric (required)
-      max: .inf
-      # Whether a higher value represents a 'better' solution (required)
-      maximize: false
-
-      # A unique identifier for your metric (required).
-      # Can contain only lowercase letters or underscores.
-    - name: emd_mean_global_vert
-      # A relatively short label, used when rendering visualisarions (required)
-      label: EMD Mean Global Vertical
-      # A one sentence summary of how this metric works (required). Used when 
-      # rendering summary tables.
-      summary: "Mean Earth Mover Distance across batch corrected samples and markers."
-      # A multi-line description of how this component works (required). Used
-      # when rendering reference documentation.
-      description: |
-        Earth Mover Distance (EMD), also known as the Wasserstein metric, measures the difference 
-        between two probability distributions. 
-        
-        Here, EMD is used to compare marker expression distributions between all integrated 
-        samples from the same group.
-        For each pair of samples and marker, the marker expression values are first converted into 
-        probability distributions. 
-        This is done by binning the expression values into a range from -100 to 100 with a bin width of 0.1.
-        The `wasserstein_distance` function from SciPy is then used to calculate the EMD between the two 
-        probability distributions belonging to the same cell type, marker, and a given paired samples.
-        This is then repeated for every cell type, marker, and paired sample.
-        Finally, the average of all these EMD values is computed and reported as the metric score.
-
-        A high score indicates overall, there is a large difference in distribution of marker expression after batch integration.
-        A low score means that overall, the samples are well integrated.
-      references:
-        doi: 
-          - 10.1023/A:1026543900054
-      links:
-        # URL to the documentation for this metric (required).
-        documentation: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.wasserstein_distance.html
-        # URL to the code repository for this metric (required).
-        repository: https://github.com/scipy/scipy
-      # The minimum possible value for this metric (required)
-      min: 0
-      # The maximum possible value for this metric (required)
-      max: .inf
-      # Whether a higher value represents a 'better' solution (required)
-      maximize: false
-      
-      # A unique identifier for your metric (required).
-      # Can contain only lowercase letters or underscores.
-    - name: emd_max_global_vert
-      # A relatively short label, used when rendering visualisarions (required)
-      label: EMD Max Global Vertical
-      # A one sentence summary of how this metric works (required). Used when 
-      # rendering summary tables.
-      summary: "Max Earth Mover Distance across batch corrected samples and markers."
-      # A multi-line description of how this component works (required). Used
-      # when rendering reference documentation.
-      description: |
-        Earth Mover Distance (EMD), also known as the Wasserstein metric, measures the difference 
-        between two probability distributions. 
-        
-        Here, EMD is used to compare marker expression distributions between all integrated 
-        samples from the same group.
-        For each pair of samples and marker, the marker expression values are first converted into 
-        probability distributions. 
-        This is done by binning the expression values into a range from -100 to 100 with a bin width of 0.1.
-        The `wasserstein_distance` function from SciPy is then used to calculate the EMD between the two 
-        probability distributions belonging to the same cell type, marker, and a given paired samples.
-        This is then repeated for every cell type, marker, and paired sample.
-        Finally, the maximum of all these EMD values is computed and reported as the metric score.
-
-        A high score indicates there is a pair of samples and marker which show large difference in distribution after batch integration.
-        A low score means that, the worst integrated pair of samples and marker are well integrated.
-      references:
-        doi: 
-          - 10.1023/A:1026543900054
-      links:
-        # URL to the documentation for this metric (required).
-        documentation: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.wasserstein_distance.html
-        # URL to the code repository for this metric (required).
-        repository: https://github.com/scipy/scipy
-      # The minimum possible value for this metric (required)
-      min: 0
-      # The maximum possible value for this metric (required)
-      max: .inf
-      # Whether a higher value represents a 'better' solution (required)
-      maximize: false
 
       # A unique identifier for your metric (required).
       # Can contain only lowercase letters or underscores.
diff --git a/src/metrics/emd/helper.py b/src/metrics/emd/helper.py
index 94881d94..c2f78a2c 100644
--- a/src/metrics/emd/helper.py
+++ b/src/metrics/emd/helper.py
@@ -5,8 +5,6 @@
 import pandas as pd
 from scipy.stats import wasserstein_distance
 
-KEY_MEAN_EMD_GLOBAL = "mean_emd_global"
-KEY_MAX_EMD_GLOBAL = "max_emd_global"
 KEY_MEAN_EMD_CT = "mean_emd_ct"
 KEY_MAX_EMD_CT = "max_emd_ct"
 KEY_EMD_VERT_MAT_split1 = "emd_vert_mat_split1"
@@ -27,10 +25,6 @@ def calculate_vertical_emd(
 
     Returns:
         dict: a dictionary containing the following elements.
-            "mean_emd_global": np.float32: mean emd value computed from a flattened data frame containing
-                mean emd computed for every marker across all pairing two samples from the same group.
-            "max_emd_global": np.float32: max emd value computed from a flattened data frame containing
-                max emd computed for every marker across all pairing two samples from the same group.
             "mean_emd_ct": np.float32: mean emd value computed from a flattened data frame containing
                 mean emd computed for every marker and cell type across all pairing two samples from the same group.
             "max_emd_ct": np.float32: max emd value computed from a flattened data frame containing
@@ -42,51 +36,46 @@ def calculate_vertical_emd(
                     or global, depending on what the 2d matrix represents.
     """
 
-    emd_split1_long, emd_split1_wide = get_vert_emd_for_integrated_adata(
+    emd_split1_long = get_vert_emd_for_integrated_adata(
         i_adata=i_split1_adata, markers_to_assess=markers_to_assess
     )
 
-    emd_split2_long, emd_split2_wide = get_vert_emd_for_integrated_adata(
+    emd_split2_long = get_vert_emd_for_integrated_adata(
         i_adata=i_split2_adata, markers_to_assess=markers_to_assess
     )
 
-    emd_long = pd.concat([emd_split1_long, emd_split2_long])
-
-    # mean global emd across all sample combinations, markers, and splits
-    mean_emd_global = np.nanmean(
-        emd_long[emd_long["cell_type"] == "global"]
-        .drop(columns=["cell_type", "first_sample", "second_sample"])
-        .to_numpy()
-        .flatten()
-    )
-    max_emd_global = np.nanmax(
-        emd_long[emd_long["cell_type"] == "global"]
-        .drop(columns=["cell_type", "first_sample", "second_sample"])
-        .to_numpy()
-        .flatten()
-    )
-
-    # mean cell type emd across all sample combinations, markers, and splits
-    mean_emd_ct = np.nanmean(
-        emd_long[emd_long["cell_type"] != "global"]
-        .drop(columns=["cell_type", "first_sample", "second_sample"])
-        .to_numpy()
-        .flatten()
-    )
-    max_emd_ct = np.nanmax(
-        emd_long[emd_long["cell_type"] != "global"]
-        .drop(columns=["cell_type", "first_sample", "second_sample"])
-        .to_numpy()
-        .flatten()
-    )
+    # safeguard
+    mean_emd_ct = np.nan
+    max_emd_ct = np.nan
+
+    # compute these only if we can.
+    emd_long = []
+    for df in [emd_split1_long, emd_split2_long]:
+        if isinstance(df, pd.DataFrame):
+            emd_long.append(df)
+
+    if len(emd_long) > 0:
+        emd_long = pd.concat(emd_long)
+
+        # mean cell type emd across all sample combinations, markers, and splits
+        mean_emd_ct = np.nanmean(
+            emd_long[emd_long["cell_type"] != "global"]
+            .drop(columns=["cell_type", "first_sample", "second_sample"])
+            .to_numpy()
+            .flatten()
+        )
+        max_emd_ct = np.nanmax(
+            emd_long[emd_long["cell_type"] != "global"]
+            .drop(columns=["cell_type", "first_sample", "second_sample"])
+            .to_numpy()
+            .flatten()
+        )
 
     return {
-        KEY_MEAN_EMD_GLOBAL: mean_emd_global,
-        KEY_MAX_EMD_GLOBAL: max_emd_global,
         KEY_MEAN_EMD_CT: mean_emd_ct,
         KEY_MAX_EMD_CT: max_emd_ct,
-        KEY_EMD_VERT_MAT_split1: emd_split1_wide,
-        KEY_EMD_VERT_MAT_split2: emd_split2_wide,
+        KEY_EMD_VERT_MAT_split1: emd_split1_long,
+        KEY_EMD_VERT_MAT_split2: emd_split2_long,
     }
 
 
@@ -122,10 +111,10 @@ def get_vert_emd_for_integrated_adata(i_adata: ad.AnnData, markers_to_assess: li
 
         print(
             f"{i_adata.uns['dataset_id']} from {i_adata.uns['method_id']} does not have"
-            f"at least 2 samples per group. Skipping EMD vertical calculation."
+            f" at least 2 samples per group. Skipping EMD vertical calculation."
         )
 
-        return np.nan, np.nan, np.nan
+        return np.nan, np.nan
 
     cell_types = i_adata.obs["cell_type"].unique()
 
@@ -137,17 +126,6 @@ def get_vert_emd_for_integrated_adata(i_adata: ad.AnnData, markers_to_assess: li
         first_sample_adata = i_adata[i_adata.obs["sample"] == sample_combo[0]]
         second_sample_adata = i_adata[i_adata.obs["sample"] == sample_combo[1]]
 
-        # global emd
-        emd_df = compute_emd(
-            left_sample=first_sample_adata,
-            right_sample=second_sample_adata,
-            markers_to_assess=markers_to_assess,
-        )
-        emd_df["cell_type"] = "global"
-        emd_df["first_sample"] = sample_combo[0]
-        emd_df["second_sample"] = sample_combo[1]
-        emd_vals.append(emd_df)
-
         # emd per cell type
         for cell_type in cell_types:
             # cell_type = cell_types[0]
@@ -175,6 +153,8 @@ def get_vert_emd_for_integrated_adata(i_adata: ad.AnnData, markers_to_assess: li
 
     # concatenate EMD values
     emd_vals = pd.concat(emd_vals)
+    # remove unparsable characters like "/"
+    emd_vals.columns = emd_vals.columns.str.replace("/", "_")
 
     # prepare the data to draw the heatmap in cytonorm 2 supp paper.
     # 1 row/column = 1 sample, a cell is emd for a given marker
@@ -182,31 +162,31 @@ def get_vert_emd_for_integrated_adata(i_adata: ad.AnnData, markers_to_assess: li
     # note, only run this after calculating mean, otherwise you end up having to
     # remove the sample id columns.
 
-    emd_wide = {}
+    # emd_wide = {}
 
-    emd_types = emd_vals["cell_type"].unique()
+    # emd_types = emd_vals["cell_type"].unique()
 
-    for marker in markers_to_assess:
-        # marker = markers_to_assess[0]
+    # for marker in markers_to_assess:
+    #     # marker = markers_to_assess[0]
 
-        # remove unparsable characters like "/"
-        marker_name = marker.replace("/", "_")
-        # have to initialise the dictionary..
-        emd_wide[marker_name] = {}
-
-        for emd_type in emd_types:
-            # ct = cell_types[0]
-            emd_df = emd_vals[emd_vals["cell_type"] == emd_type]
-
-            if emd_df.shape[0] > 0:
-                # safeguard. Only pivot if we computed the emd.
-                # This is a safeguard in case there is a rare cell type which we don't have
-                # any samples with at least 50 cells for.
-                emd_wide[marker_name][emd_type] = emd_df.pivot(
-                    index="second_sample", columns="first_sample", values=marker
-                )
+    #     # remove unparsable characters like "/"
+    #     marker_name = marker.replace("/", "_")
+    #     # have to initialise the dictionary..
+    #     emd_wide[marker_name] = {}
+
+    #     for emd_type in emd_types:
+    #         # ct = cell_types[0]
+    #         emd_df = emd_vals[emd_vals["cell_type"] == emd_type]
 
-    return emd_vals, emd_wide
+    #         if emd_df.shape[0] > 0:
+    #             # safeguard. Only pivot if we computed the emd.
+    #             # This is a safeguard in case there is a rare cell type which we don't have
+    #             # any samples with at least 50 cells for.
+    #             emd_wide[marker_name][emd_type] = emd_df.pivot(
+    #                 index="second_sample", columns="first_sample", values=marker
+    #             )
+
+    return emd_vals
 
 
 def calculate_horizontal_emd(
@@ -226,10 +206,6 @@ def calculate_horizontal_emd(
 
     Returns:
         dict: a dictionary containing the following elements.
-            "mean_emd_global": np.float32: mean emd value computed from a flattened data frame containing
-                mean emd computed for every marker across all pairs of samples from a given donor.
-            "max_emd_global": np.float32: max emd value computed from a flattened data frame containing
-                max emd computed for every marker across all pairs of samples from a given donor.
             "mean_emd_ct": np.float32: mean emd value computed from a flattened data frame containing
                 mean emd computed for every marker and cell type across all pairs of samples from a given donor.
             "max_emd_ct": np.float32: max emd value computed from a flattened data frame containing
@@ -239,8 +215,6 @@ def calculate_horizontal_emd(
     """
 
     emd_per_donor_per_ct = []
-    # global means agnostic of cell type labels
-    emd_per_donor_global = []
 
     for donor in donor_list:
         # donor = donor_list[0]
@@ -254,8 +228,8 @@ def calculate_horizontal_emd(
         )
         if len(cell_type_not_in_both) > 1:
             print(
-                f"In donor {donor}: some cell types are in left integrated output"
-                f" but not in right integrated output.\n"
+                f"In donor {donor}: some cell types are in split 1"
+                f" but not in split 2.\n"
                 f"Cell types missing: {''.join(cell_type_not_in_both)}]n"
                 f"Computing cell type EMD using just cell types common in both."
             )
@@ -274,9 +248,11 @@ def calculate_horizontal_emd(
             # Do not calculate if we have less than 50 cells as it does not make sense.
             if i_split1_ct.n_obs < 50 or i_split2_ct.n_obs < 50:
                 print(
-                    f"There are less than 50 cells for either left or right integrated "
-                    f"data for donor {donor} and cell type {cell_type}.\n"
-                    f"Skipping calculating EMD for this donor and cell type."
+                    f"There are less than 50 cells in either split 1 or split 2"
+                    f" for donor {donor} and cell type {cell_type}.\n"
+                    f"Split 1 {cell_type}: {i_split1_ct.n_obs} cells.\n"
+                    f"Split 2 {cell_type}: {i_split2_ct.n_obs} cells.\n"
+                    f"Skipping calculating horizontal EMD for this donor and cell type."
                 )
                 continue
 
@@ -290,19 +266,7 @@ def calculate_horizontal_emd(
 
             emd_per_donor_per_ct.append(emd_df)
 
-        # calculate EMD when combining all cell types as well.
-        emd_df = compute_emd(
-            left_sample=i_split1_donor,
-            right_sample=i_split2_donor,
-            markers_to_assess=markers_to_assess,
-        )
-        emd_df["cell_type"] = "global"
-        emd_df["donor"] = donor
-
-        emd_per_donor_global.append(emd_df)
-
     emd_per_donor_per_ct = pd.concat(emd_per_donor_per_ct)
-    emd_per_donor_global = pd.concat(emd_per_donor_global)
 
     # compute the mean and max per ct and for global.
     mean_emd_ct = np.nanmean(
@@ -312,24 +276,15 @@ def calculate_horizontal_emd(
         emd_per_donor_per_ct.drop(columns=["cell_type", "donor"]).values
     )
 
-    mean_emd_global = np.nanmean(
-        emd_per_donor_global.drop(columns=["cell_type", "donor"]).values
-    )
-    max_emd_global = np.nanmax(
-        emd_per_donor_global.drop(columns=["cell_type", "donor"]).values
-    )
-
     # concatenate the global and cell type emd
-    emd_per_donor = pd.concat([emd_per_donor_per_ct, emd_per_donor_global])
-
-    emd_per_donor.columns = emd_per_donor.columns.str.replace("/", "_", regex=False)
+    emd_per_donor_per_ct.columns = emd_per_donor_per_ct.columns.str.replace(
+        "/", "_", regex=False
+    )
 
     return {
-        KEY_MEAN_EMD_GLOBAL: mean_emd_global,
-        KEY_MAX_EMD_GLOBAL: max_emd_global,
         KEY_MEAN_EMD_CT: mean_emd_ct,
         KEY_MAX_EMD_CT: max_emd_ct,
-        KEY_EMD_HORZ_PER_DONOR: emd_per_donor,
+        KEY_EMD_HORZ_PER_DONOR: emd_per_donor_per_ct,
     }
 
 
@@ -411,3 +366,34 @@ def bin_array(values):
 
     # the 1st return value is the bin indices
     return bin_indices, bin_probabilities
+
+
+def check_donor_batches(input_integrated_split1, input_integrated_split2):
+    """
+    Ensure each donor is present in exactly one batch per split,
+    and that the batch IDs differ between splits.
+    """
+
+    donor_list = input_integrated_split1.obs["donor"].unique()
+
+    for donor in donor_list:
+        batch_split1 = input_integrated_split1.obs[
+            input_integrated_split1.obs["donor"] == donor
+        ]["batch"].unique()
+        batch_split2 = input_integrated_split2.obs[
+            input_integrated_split2.obs["donor"] == donor
+        ]["batch"].unique()
+
+        if len(batch_split1) > 1 or len(batch_split2) > 1:
+            raise ValueError(
+                f"Donor {donor} has samples in {len(batch_split1)} batches in integrated left"
+                f" and {len(batch_split2)} batches in integrated right. It should only have"
+                f" samples in exactly ONE batch in each of integrated left and integrated right."
+            )
+
+        if batch_split1[0] == batch_split2[0]:
+            raise ValueError(
+                f"Donor {donor} has samples in the same batch for both integrated left and right.\n"
+                f"Integrated left batch id: {batch_split1[0]}.\n"
+                f"Integrated right batch id: {batch_split2[0]}."
+            )
diff --git a/src/metrics/emd/script.py b/src/metrics/emd/script.py
index e953b757..6586070d 100644
--- a/src/metrics/emd/script.py
+++ b/src/metrics/emd/script.py
@@ -9,13 +9,18 @@
     "input_integrated_split1": "resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/integrated_split1.h5ad",
     "input_integrated_split2": "resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/integrated_split2.h5ad",
     "input_unintegrated": "resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/unintegrated.h5ad",
+    # "input_integrated_split1": "resources_test/task_cyto_batch_integration/human_blood_mass_cytometry_subset/integrated_split1.h5ad",
+    # "input_integrated_split2": "resources_test/task_cyto_batch_integration/human_blood_mass_cytometry_subset/integrated_split2.h5ad",
+    # "input_unintegrated": "resources_test/task_cyto_batch_integration/human_blood_mass_cytometry_subset/unintegrated.h5ad",
     "output": "resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/emd_out.h5ad",
 }
-meta = {"name": "emd", "resources_dir": "src/utils/helper_functions.py"}
+meta = {"name": "emd", "resources_dir": "src/utils/"}
 ## VIASH END
 
 sys.path.append(meta["resources_dir"])
 
+# import src.metrics.emd.helper as emd_helper
+# import src.utils.helper_functions as global_helper
 import helper as emd_helper
 import helper_functions as global_helper
 
@@ -37,11 +42,15 @@
 )
 
 # more preprocessing
-input_integrated_split1 = global_helper.subset_markers_tocorrect(input_integrated_split1)
+input_integrated_split1 = global_helper.subset_markers_tocorrect(
+    input_integrated_split1
+)
 input_integrated_split1 = global_helper.subset_nocontrols(input_integrated_split1)
 input_integrated_split1 = global_helper.remove_unlabelled(input_integrated_split1)
 
-input_integrated_split2 = global_helper.subset_markers_tocorrect(input_integrated_split2)
+input_integrated_split2 = global_helper.subset_markers_tocorrect(
+    input_integrated_split2
+)
 input_integrated_split2 = global_helper.subset_nocontrols(input_integrated_split2)
 input_integrated_split2 = global_helper.remove_unlabelled(input_integrated_split2)
 
@@ -54,36 +63,16 @@
 dataset_id = input_unintegrated.uns["dataset_id"]
 method_id = input_integrated_split1.uns["method_id"]
 
-# shouldn't need these anymore
-# del input_unintegrated
+
+# check that the data for each donor in integrated left and right are actually
+# from two different batches!
+emd_helper.check_donor_batches(
+    input_integrated_split1=input_integrated_split1,
+    input_integrated_split2=input_integrated_split2,
+)
 
 # calculate horizontal EMD for each donor across integrated left and right
 donor_list = input_integrated_split1.obs["donor"].unique()
-
-# check that the data for each donor in integrated left and right are actually from two different batches!
-for donor in donor_list:
-    # donor = donor_list[0]
-    batch_split1 = input_integrated_split1.obs[input_integrated_split1.obs["donor"] == donor][
-        "batch"
-    ].unique()
-    batch_split2 = input_integrated_split2.obs[
-        input_integrated_split2.obs["donor"] == donor
-    ]["batch"].unique()
-
-    if len(batch_split1) > 1 or len(batch_split2) > 1:
-        raise ValueError(
-            f"Donor {donor} has samples in {len(batch_split1)} batches in integrated left"
-            f" and {len(batch_split2)} batches in integrated right.It should only have"
-            f"samples in exactly ONE batch in each of integrated left and integrated right."
-        )
-
-    if batch_split1[0] == batch_split2[0]:
-        raise ValueError(
-            f"Donor {donor} has samples in the same batch for both integrated left and right.\n"
-            f"Integrated left batch id: {batch_split1[0]}.\n"
-            f"Integrated right batch id: {batch_split2[0]}."
-        )
-
 emd_horz = emd_helper.calculate_horizontal_emd(
     i_split1_adata=input_integrated_split1,
     i_split2_adata=input_integrated_split2,
@@ -105,22 +94,14 @@
         "dataset_id": dataset_id,
         "method_id": method_id,
         "metric_ids": [
-            "emd_mean_global_horiz",
-            "emd_max_global_horiz",
             "emd_mean_ct_horiz",
             "emd_max_ct_horiz",
-            "emd_mean_global_vert",
-            "emd_max_global_vert",
             "emd_mean_ct_vert",
             "emd_max_ct_vert",
         ],
         "metric_values": [
-            emd_horz[emd_helper.KEY_MEAN_EMD_GLOBAL],
-            emd_horz[emd_helper.KEY_MAX_EMD_GLOBAL],
             emd_horz[emd_helper.KEY_MEAN_EMD_CT],
             emd_horz[emd_helper.KEY_MAX_EMD_CT],
-            emd_vert[emd_helper.KEY_MEAN_EMD_GLOBAL],
-            emd_vert[emd_helper.KEY_MAX_EMD_GLOBAL],
             emd_vert[emd_helper.KEY_MEAN_EMD_CT],
             emd_vert[emd_helper.KEY_MAX_EMD_CT],
         ],
diff --git a/src/metrics/flowsom_mapping_similarity/script.R b/src/metrics/flowsom_mapping_similarity/script.R
index 832c5321..4d5fb8f5 100644
--- a/src/metrics/flowsom_mapping_similarity/script.R
+++ b/src/metrics/flowsom_mapping_similarity/script.R
@@ -6,10 +6,14 @@ par <- list(
   "input_unintegrated" = 'resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/unintegrated.h5ad',
   "input_integrated_split1" = 'resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/integrated_split1.h5ad',
   "input_integrated_split2" = 'resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/integrated_split2.h5ad',
+  # if using perfect integration
+  # input_integrated_split1 = "resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/perfect_integrated_split1.h5ad",
+  # input_integrated_split2 = "resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/perfect_integrated_split2.h5ad",
   "output" = 'resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/score.h5ad'
 )
 meta <- list(
-  "name" = 'flowsom_mapping_similarity'
+  "name" = 'flowsom_mapping_similarity',
+  "resources_dir" = "src/utils"
 )
 ## VIASH END
 
@@ -18,19 +22,16 @@ source(paste0(meta$resources_dir, "/helper_functions.R"))
 library(anndata)
 
 unintegrated <- anndata::read_h5ad(par[["input_unintegrated"]])
-integrated_s1 <- anndata::read_h5ad(par[["input_integrated_split1"]]) 
 
-print(unintegrated)
-print(integrated_s1)
 # read and filter split 1 data
 integrated_s1 <- anndata::read_h5ad(par[["input_integrated_split1"]]) |>
-  get_obs_var_for_integrated(unintegrated) |>
+  get_obs_var_for_integrated(unintegrated, split_id = 1) |>
   subset_nocontrols() |>
   remove_unlabelled()
 
 # read and filter split 2 data
 integrated_s2 <- anndata::read_h5ad(par[["input_integrated_split2"]]) |>
-  get_obs_var_for_integrated(unintegrated) |>
+  get_obs_var_for_integrated(unintegrated, split_id = 2) |>
   subset_nocontrols() |>
   remove_unlabelled()
 
diff --git a/src/metrics/lisi/config.vsh.yaml b/src/metrics/lisi/config.vsh.yaml
new file mode 100644
index 00000000..91009ecb
--- /dev/null
+++ b/src/metrics/lisi/config.vsh.yaml
@@ -0,0 +1,130 @@
+# 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_metric.yaml
+
+# A unique identifier for your component (required).
+# Can contain only lowercase letters or underscores.
+name: lisi
+
+
+
+# Metadata for your component
+info:
+  metrics:
+      # A unique identifier for your metric (required).
+      # Can contain only lowercase letters or underscores.
+    - name: iLisi
+      # A relatively short label, used when rendering visualisarions (required)
+      label: iLisi
+      summary: "Integration Lisi score."
+      description: |
+        Compute the integration local inverse simpson index (iLISI) for each cell, then return the median 
+        as the final score. 
+      references:
+        doi: 
+          - 10.1038/s41592-019-0619-0
+        bibtex:
+          - |
+            @article{korsunsky2019fast,
+            title={Fast, sensitive and accurate integration of single-cell data with Harmony},
+            author={Korsunsky, Ilya and Millard, Nghia and Fan, Jean and Slowikowski, Kamil and Zhang, Fan and Wei, Kevin and Baglaenko, Yuriy and Brenner, Michael and Loh, Po-ru and Raychaudhuri, Soumya},
+            journal={Nature methods},
+            volume={16},
+            number={12},
+            pages={1289--1296},
+            year={2019},
+            publisher={Nature Publishing Group US New York}
+            }
+      links:
+        # URL to the documentation for this metric (required).
+        documentation: https://scib-metrics.readthedocs.io/en/latest/api.html
+        # URL to the code repository for this metric (required).
+        repository: https://github.com/YosefLab/scib-metrics
+      # The minimum possible value for this metric (required)
+      min: -0.0001
+      # The maximum possible value for this metric (required)
+      max: 1.0001
+      # Whether a higher value represents a 'better' solution (required)
+      maximize: true
+
+
+    - name: cLisi
+      label: cLisi
+      summary: "Cell type Lisi score."
+      description: |
+        Compute the cell type local inverse simpson index (cLISI) for each cell, then return the median 
+        as the final score. 
+      references:
+        doi: 
+          - 10.1038/s41592-019-0619-0
+        bibtex:
+          - |
+            @article{korsunsky2019fast,
+            title={Fast, sensitive and accurate integration of single-cell data with Harmony},
+            author={Korsunsky, Ilya and Millard, Nghia and Fan, Jean and Slowikowski, Kamil and Zhang, Fan and Wei, Kevin and Baglaenko, Yuriy and Brenner, Michael and Loh, Po-ru and Raychaudhuri, Soumya},
+            journal={Nature methods},
+            volume={16},
+            number={12},
+            pages={1289--1296},
+            year={2019},
+            publisher={Nature Publishing Group US New York}
+            }
+      links:
+        # URL to the documentation for this metric (required).
+        documentation: https://scib-metrics.readthedocs.io/en/latest/api.html
+        # URL to the code repository for this metric (required).
+        repository: https://github.com/YosefLab/scib-metrics
+      # The minimum possible value for this metric (required)
+      min: -0.0001
+      # The maximum possible value for this metric (required)
+      max: 1.0001
+      # Whether a higher value represents a 'better' solution (required)
+      maximize: true
+
+
+# Component-specific parameters (optional)
+# arguments:
+#   - name: "--n_neighbors"
+#     type: "integer"
+#     default: 5
+#     description: Number of neighbors to use.
+
+# Resources required to run the component
+resources:
+  # The script of your component (required)
+  - type: python_script
+    path: script.py
+  - path: /src/utils/helper_functions.py
+
+engines:
+  # Specifications for the Docker image for this component.
+  - type: docker
+    image: python:3.11
+    setup:
+      - type: apt
+        packages:
+          - procps
+      - type: python
+        packages:
+          - jax~=0.6.2
+          - jaxlib~=0.6.2
+          - anndata~=0.11.0
+          - scanpy~=1.11.0
+          - scib-metrics~=0.5.7
+          - pyyaml
+          - requests
+          - jsonschema
+        github:
+          - "openproblems-bio/core#subdirectory=packages/python/openproblems"
+
+
+runners:
+  # This platform allows running the component natively
+  - type: executable
+  # Allows turning the component into a Nextflow module / pipeline.
+  - type: nextflow
+    directives:
+      label: [midtime,midmem,midcpu]
diff --git a/src/metrics/lisi/script.py b/src/metrics/lisi/script.py
new file mode 100644
index 00000000..f160266d
--- /dev/null
+++ b/src/metrics/lisi/script.py
@@ -0,0 +1,78 @@
+import anndata as ad
+import scib_metrics as sm
+import numpy as np
+import sys
+
+## 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_unintegrated': 'resources_test/.../unintegrated.h5ad',
+  'input_integrated_split1': 'resources_test/.../integrated_split1.h5ad',
+  'input_integrated_split2': 'resources_test/.../integrated_split2.h5ad',
+  'output': 'output.h5ad'
+}
+meta = {
+  'name': 'lisi'
+}
+## VIASH END
+
+sys.path.append(meta["resources_dir"])
+from helper_functions import (
+    get_obs_var_for_integrated,
+    subset_markers_tocorrect,
+)
+
+print('Reading input files', flush=True)
+input_unintegrated = ad.read_h5ad(par['input_unintegrated'])
+input_integrated_split1 = ad.read_h5ad(par['input_integrated_split1'])
+input_integrated_split2 = ad.read_h5ad(par['input_integrated_split2'])
+
+print("Formatting input files", flush=True)
+integrated_s1, integrated_s2 = get_obs_var_for_integrated(
+    input_integrated_split1, input_integrated_split2, input_unintegrated
+)
+integrated_s1 = subset_markers_tocorrect(integrated_s1)
+integrated_s2 = subset_markers_tocorrect(integrated_s2)
+
+print('Compute metrics', flush=True)
+n_batches = len(integrated_s1.obs.batch.unique())
+n_celltypes = len(integrated_s1.obs.cell_type.unique())
+
+print("Compute iLisi and cLisi for split 1", flush=True)
+knn = sm.nearest_neighbors.pynndescent(integrated_s1.layers['integrated'], n_neighbors=100, random_state=0)
+
+ilisi_s1_per_cell = sm.lisi_knn(knn, integrated_s1.obs.batch)
+ilisi_s1 = (np.nanmedian(ilisi_s1_per_cell) - 1) / (n_batches - 1)
+
+clisi_s1_per_cell = sm.lisi_knn(knn, integrated_s1.obs.cell_type)
+clisi_s1 = (n_celltypes - np.nanmedian(clisi_s1_per_cell)) / (n_celltypes - 1)
+
+print("Compute iLisi and cLisi for split 2", flush=True)
+knn = sm.nearest_neighbors.pynndescent(integrated_s2.layers['integrated'], n_neighbors=100, random_state=0)
+ilisi_s2_per_cell = sm.lisi_knn(knn, integrated_s2.obs.batch)
+ilisi_s2 = (np.nanmedian(ilisi_s2_per_cell) - 1) / (n_batches - 1)
+
+clisi_s2_per_cell = sm.lisi_knn(knn, integrated_s2.obs.cell_type)
+clisi_s2 = (n_celltypes - np.nanmedian(clisi_s2_per_cell)) / (n_celltypes - 1)
+
+ilisi = np.mean([ilisi_s1, ilisi_s2])
+clisi = np.mean([clisi_s1, clisi_s2])
+uns_metric_ids = [ 'ilisi', 'clisi' ]
+uns_metric_values = [ ilisi, clisi ]
+
+print("Write output AnnData to file", flush=True)
+output = ad.AnnData(
+  uns={
+        "dataset_id": integrated_s1.uns["dataset_id"],
+        "method_id": integrated_s1.uns["method_id"],
+        "metric_ids": uns_metric_ids,
+        "metric_values": uns_metric_values,
+        "ilisi_s1_index": ilisi_s1_per_cell,
+        "ilisi_s2_index": ilisi_s2_per_cell,
+        "clisi_s1_index": clisi_s1_per_cell,
+        "clisi_s2_index": clisi_s2_per_cell
+    }
+
+)
+output.write_h5ad(par['output'], compression='gzip')
diff --git a/src/metrics/n_inconsistent_peaks/config.vsh.yaml b/src/metrics/n_inconsistent_peaks/config.vsh.yaml
index 66cb1f92..569be185 100644
--- a/src/metrics/n_inconsistent_peaks/config.vsh.yaml
+++ b/src/metrics/n_inconsistent_peaks/config.vsh.yaml
@@ -2,7 +2,7 @@
 __merge__: ../../api/comp_metric.yaml
 
 name: n_inconsistent_peaks
-status: disabled
+# status: disabled
 
 info:
   metrics:
@@ -14,12 +14,11 @@ info:
       # A multi-line description of how this component works (required). Used
       # when rendering reference documentation.
       description: |
-        The metric compares the number of marker expression peaks between the validation and batch-normalized data. 
-        The number of peaks is calculated using the `scipy.signal.find_peaks` function. 
-        The metric is calculated as the absolute difference between the number of peaks in the validation and batch-normalized data.
-        The marker expression profiles are first smoothed using kernel density estimation (KDE) (`scipy.stats.gaussian_kde`),
-        and then peaks are then identified using the `scipy.signal.find_peaks` function.
-        For peak calling, the `prominence` parameter is set to 0.1 and the `height` parameter is set to 0.05*max_density.
+        The metric compares the number of marker expression peaks between batch integrated technical replicates (split 1 and split 2). 
+        The metric is calculated as the absolute difference between the number of peaks in the technical replicates.
+        The marker expression profiles are first z-score scaled (with common mu and sigma) and then smoothed using kernel density estimation (KDE) (`scipy.stats.gaussian_kde`).
+        Finally, peaks are identified using the `scipy.signal.find_peaks` function.
+        For peak calling, the `prominence` parameter is set to 0.01 and the `height` parameter is set to 0.1.
       references:
         doi: 
           - 10.1038/s41592-019-0686-2
diff --git a/src/metrics/n_inconsistent_peaks/helper.py b/src/metrics/n_inconsistent_peaks/helper.py
index 182eafec..9ffb3a57 100644
--- a/src/metrics/n_inconsistent_peaks/helper.py
+++ b/src/metrics/n_inconsistent_peaks/helper.py
@@ -16,7 +16,7 @@ def get_kde_density(expression_array):
     min_val = expression_array.min()
     max_val = expression_array.max()
     marker_values = np.reshape(expression_array, (1,-1)) # Reshape array for KDE
-    kde = gaussian_kde(marker_values, bw_method=0.3)
+    kde = gaussian_kde(marker_values, bw_method='scott')
     x_grid = np.linspace(min_val, max_val, 100)
     density = kde(x_grid)
     #Plot, for debugging
@@ -35,8 +35,8 @@ def call_peaks(density):
     peaks: array of values representing the peaks of the density
     '''
 
-    height_trsh = 0.05*density.max()
-    prom_trsh = 0.1
+    height_trsh = 0.1
+    prom_trsh = 0.01
                 
     peaks, _ = find_peaks(density, 
                         prominence=prom_trsh,
diff --git a/src/metrics/n_inconsistent_peaks/script.py b/src/metrics/n_inconsistent_peaks/script.py
index 859c104a..6e1c858c 100644
--- a/src/metrics/n_inconsistent_peaks/script.py
+++ b/src/metrics/n_inconsistent_peaks/script.py
@@ -4,17 +4,18 @@
 import numpy as np
 
 ## 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`.
+# The following code has been auto-generated by Viash.
 par = {
-  'input_validation': 'resources_test/.../validation.h5ad',
-  'input_unintegrated': 'resources_test/.../unintegrated.h5ad',
-  'input_integrated': 'resources_test/.../integrated.h5ad',
-  'output': 'output.h5ad'
+  'input_unintegrated': 'resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/unintegrated.h5ad',
+  'input_integrated_split1': 'resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/integrated_split1.h5ad',
+  'input_integrated_split2': 'resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/integrated_split2.h5ad',
+  'output': 'resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/score.h5ad'
 }
 meta = {
-  'name': 'n_inconsistent_peaks'
+  'name': 'n_inconsistent_peaks',
 }
+
+
 ## VIASH END
 
 sys.path.append(meta["resources_dir"])
@@ -26,66 +27,83 @@
     subset_nocontrols,
 )
 
-print('Reading input files', flush=True)
-input_validation = ad.read_h5ad(par['input_validation'])
-input_unintegrated = ad.read_h5ad(par['input_unintegrated'])
-input_integrated = ad.read_h5ad(par['input_integrated'])
-
-print('Formatting input files', flush=True)
-#Format data integrated data
-input_integrated = get_obs_var_for_integrated(input_integrated,input_validation,input_unintegrated)
-input_integrated = subset_markers_tocorrect(input_integrated)
-input_integrated = subset_nocontrols(input_integrated)
-input_integrated = remove_unlabelled(input_integrated)
-#Format validation data
-input_validation = subset_markers_tocorrect(input_validation)
-input_validation = remove_unlabelled(input_validation)
+print("Reading input files", flush=True)
+integrated_s1 = ad.read_h5ad(par["input_integrated_split1"])
+integrated_s2 = ad.read_h5ad(par["input_integrated_split2"])
+unintegrated = ad.read_h5ad(par["input_unintegrated"])
+
+print("Formatting input files", flush=True)
+integrated_s1, integrated_s2 = get_obs_var_for_integrated(
+    integrated_s1, integrated_s2, unintegrated
+)
+
+integrated_s1 = subset_nocontrols(integrated_s1)
+integrated_s1 = subset_markers_tocorrect(integrated_s1)
+integrated_s1 = subset_nocontrols(integrated_s1)
+integrated_s1 = remove_unlabelled(integrated_s1)
+
+integrated_s2 = subset_nocontrols(integrated_s2)
+integrated_s2 = subset_markers_tocorrect(integrated_s2)
+integrated_s2 = subset_nocontrols(integrated_s2)
+integrated_s2 = remove_unlabelled(integrated_s2)
 
 print('Compute metric (All cells)', flush=True)
-donor_list = input_integrated.obs['donor'].unique()
+donor_list = integrated_s1.obs['donor'].unique()
 n_inconsistent_peaks = 0
 
 for donor in donor_list:
-    integrated_view = input_integrated[input_integrated.obs['donor'] == donor]
-    validation_view = input_validation[input_validation.obs['donor'] == donor]
+    s1_view = integrated_s1[integrated_s1.obs['donor'] == donor]
+    s2_view = integrated_s2[integrated_s2.obs['donor'] == donor]
 
-    for marker in integrated_view.var.index:
-        mexp_integrated = np.array(integrated_view[:,marker].layers["integrated"])
-        mexp_validation = np.array(validation_view[:,marker].layers["preprocessed"])
-        density_integrated = get_kde_density(mexp_integrated)
-        peaks_integrated = call_peaks(density_integrated)
-        density_validation = get_kde_density(mexp_validation)
-        peaks_validation = call_peaks(density_validation)
+    for marker in s1_view.var.index:
+        marker_expression_s1_unscaled = np.array(s1_view[:,marker].layers["integrated"])
+        marker_expression_s2_unscaled = np.array(s2_view[:,marker].layers["integrated"])
 
-        if peaks_integrated != peaks_validation:
-            n_inconsistent_peaks += abs(peaks_integrated - peaks_validation)
+        pooled = np.concatenate([marker_expression_s1_unscaled, marker_expression_s2_unscaled])
+        mu, sd = pooled.mean(), pooled.std()
+        marker_expression_s1 = (marker_expression_s1_unscaled - mu) / (sd)
+        marker_expression_s2 = (marker_expression_s2_unscaled - mu) / (sd)
+
+        density_s1 = get_kde_density(marker_expression_s1)
+        peaks_s1 = call_peaks(density_s1)
+        density_s2 = get_kde_density(marker_expression_s2)
+        peaks_s2 = call_peaks(density_s2)
+
+        if peaks_s1 != peaks_s2:
+            n_inconsistent_peaks += abs(peaks_s1 - peaks_s2)
 
 print('Compute metric (per cell type)', flush=True)
 n_inconsistent_peaks_ct = 0
 
 for donor in donor_list:
-    integrated_view = input_integrated[input_integrated.obs['donor'] == donor]
-    validation_view = input_validation[input_validation.obs['donor'] == donor]
-    celltype_list = integrated_view.obs['cell_type'].unique()
+    s1_view = integrated_s1[integrated_s1.obs['donor'] == donor]
+    s2_view = integrated_s2[integrated_s2.obs['donor'] == donor]
+    celltype_list = s1_view.obs['cell_type'].unique()
 
     for celltype in celltype_list:
-        integrated_view_ct = integrated_view[integrated_view.obs['cell_type'] == celltype]
-        validation_view_ct = validation_view[validation_view.obs['cell_type'] == celltype]
+        s1_view_ct = s1_view[s1_view.obs['cell_type'] == celltype]
+        s2_view_ct = s2_view[s2_view.obs['cell_type'] == celltype]
 
-        if integrated_view_ct.shape[0] < 100 or validation_view_ct.shape[0] < 100:
-          print(donor,celltype,"skipped cause less than 100 cells are present in either integrated or validation dataset")
+        if s1_view_ct.shape[0] < 100 or s2_view_ct.shape[0] < 100:
+          print(donor,celltype,"skipped cause less than 100 cells are present in either split 1 or split 2 dataset")
           continue
 
-        for marker in integrated_view_ct.var.index:
-            mexp_integrated = np.array(integrated_view_ct[:, marker].layers["integrated"])
-            mexp_validation = np.array(validation_view_ct[:, marker].layers["preprocessed"])
-            density_integrated = get_kde_density(mexp_integrated)
-            peaks_integrated = call_peaks(density_integrated)
-            density_validation = get_kde_density(mexp_validation)
-            peaks_validation = call_peaks(density_validation)
+        for marker in s1_view_ct.var.index:
+            marker_expression_s1_unscaled = np.array(s1_view_ct[:, marker].layers["integrated"])
+            marker_expression_s2_unscaled = np.array(s2_view_ct[:, marker].layers["integrated"])
+
+            pooled = np.concatenate([marker_expression_s1_unscaled, marker_expression_s2_unscaled])
+            mu, sd = pooled.mean(), pooled.std()
+            marker_expression_s1 = (marker_expression_s1_unscaled - mu) / (sd)
+            marker_expression_s2 = (marker_expression_s2_unscaled - mu) / (sd)
 
-            if peaks_integrated != peaks_validation:
-                n_inconsistent_peaks_ct += abs(peaks_integrated - peaks_validation)
+            density_s1 = get_kde_density(marker_expression_s1)
+            peaks_s1 = call_peaks(density_s1)
+            density_s2 = get_kde_density(marker_expression_s2)
+            peaks_s2 = call_peaks(density_s2)
+
+            if peaks_s1 != peaks_s2:
+                n_inconsistent_peaks_ct += abs(peaks_s1 - peaks_s2)
 
 
 
@@ -95,10 +113,12 @@
 print("Write output AnnData to file", flush=True)
 output = ad.AnnData(
     uns={
-    'dataset_id': input_integrated.uns['dataset_id'],
-    'method_id': input_integrated.uns['method_id'],
+    'dataset_id': integrated_s1.uns['dataset_id'],
+    'method_id': integrated_s1.uns['method_id'],
     'metric_ids': uns_metric_ids,
     'metric_values': uns_metric_values
   }
 )
 output.write_h5ad(par['output'], compression='gzip')
+
+print(uns_metric_ids, uns_metric_values)
\ No newline at end of file
diff --git a/src/utils/helper_functions.R b/src/utils/helper_functions.R
index f5730f20..a6af5b70 100644
--- a/src/utils/helper_functions.R
+++ b/src/utils/helper_functions.R
@@ -14,9 +14,10 @@ requireNamespace("anndataR", quietly = TRUE)
 #'
 #' @param i_adata AnnData object, integrated data
 #' @param u_adata AnnData object, unintegrated dataset
+#' @param split_id numeric, split id of the integrated data
 #' @return AnnData object with .var and .obs added
 #'
-get_obs_var_for_integrated <- function(i_adata, u_adata) {
+get_obs_var_for_integrated <- function(i_adata, u_adata, split_id) {
 
     i_adata$obs <- u_adata$obs[i_adata$obs_names, ]
     i_adata$var <- u_adata$var[i_adata$var_names, ]
@@ -25,16 +26,23 @@ get_obs_var_for_integrated <- function(i_adata, u_adata) {
     # everything is from batch 1, but some samples need to be labelled to come from batch 2
     if (i_adata$uns["method_id"] == "perfect_integration") {
         cat(
-            "Control method 'perfect_integration' detected. Changing batch labels for split 2.\n"
+            "Control method 'perfect_integration' detected. Changing batch labels.\n"
         )
         cat("Computing new batch labels\n")
         # mutate is needed as donors that are used for controls, we won't have the mapping
         i_adata_new_batch_labels <- get_batch_label_perfect_integration(
             u_adata = u_adata,
             i_adata = i_adata,
-            split_id = 1
+            split_id = split_id
         )
 
+        # safeguard
+        if (! all(i_adata_new_batch_labels$donor == i_adata$obs$donor)) {
+            stop(
+                "Donor labels do not match between new batch labels and integrated data. This should not happen!"
+            )
+        }
+
         cat("Attaching new batch labels\n")
         i_adata$obs$batch <- i_adata_new_batch_labels$new_batch_label
     }
@@ -53,12 +61,14 @@ get_obs_var_for_integrated <- function(i_adata, u_adata) {
 #' @return a dataframe with donor and new batch label
 #'
 get_batch_label_perfect_integration <- function(u_adata, i_adata, split_id) {
+    # this return which batch sample we used for a donor for a given split
     actual_donor_batch_map <- unique(
         u_adata$obs[(u_adata$obs$split == split_id), c("donor", "batch")]
     )
-    # mutate is needed as donors that are used for controls, we won't have the mapping
+    # mutate is needed as donors that are used for controls, won't have batch_new as 
+    # the split id is 0
     i_adata_new_batch_labels <- i_adata$obs[, c("donor", "batch")] %>%
-        left_join(actual_donor_batch_map, by="donor", suffix = c("_old", "_new")) %>%
+        left_join(actual_donor_batch_map, by = "donor", suffix = c("_old", "_new")) %>%
         mutate(new_batch_label = ifelse(is.na(batch_new), batch_old, batch_new)) %>%
         select(donor, new_batch_label)
 
@@ -105,3 +115,43 @@ remove_unlabelled <- function(adata) {
         c("unlabelled", "unlabeled")
     adata[!is_unlabelled, ]
 }
+
+#' Subsets the anndata object in a stratified manner
+#' with 'cell type' and 'sample' as strata.
+#'
+#' @param adata AnnData object
+#' @param frac numeric, fraction of cells to keep for each cell type
+#' @param seed numeric, seed for reproducibility
+#' @param anndatar logical, whether the input is anndataR object or not
+#' @return AnnData object with only the markers to correct
+subset_by_celltype <- function(adata, frac = 0.5, seed = 1, anndatar = TRUE) {
+  set.seed(seed)
+  
+  obs <- adata$obs
+  obs$cell_id <- rownames(obs)
+  obs$.row <- seq_len(nrow(obs))   # original order
+  
+  keep_ids <- obs %>%
+    group_by(cell_type, sample) %>%
+    slice_sample(prop = frac) %>%
+    ungroup() %>%
+    arrange(.row) %>%     # restore original order
+    pull(cell_id)
+  
+  if (anndatar == TRUE){
+    keep_idx <- match(keep_ids, adata$obs_names)
+    
+    adata_sub <- anndataR::AnnData(
+      X   = NULL,
+      obs = adata$obs[keep_idx, , drop = FALSE],
+      var = adata$var,
+      uns = adata$uns,
+      layers = list(
+        "integrated" = adata$layers$integrated[keep_idx, , drop = FALSE]
+      )
+    )
+  } else{
+    adata_sub <- adata[keep_ids, ]
+  }
+}
+
diff --git a/src/workflows/run_benchmark/config.vsh.yaml b/src/workflows/run_benchmark/config.vsh.yaml
index 3968de6f..b63c6e3e 100644
--- a/src/workflows/run_benchmark/config.vsh.yaml
+++ b/src/workflows/run_benchmark/config.vsh.yaml
@@ -83,8 +83,8 @@ dependencies:
   - name: utils/extract_uns_metadata
     repository: op
   # - name: control_methods/shuffle_integration
-  # - name: control_methods/shuffle_integration_by_batch
-  # - name: control_methods/shuffle_integration_by_cell_type
+  - name: control_methods/shuffle_integration_by_batch
+  - name: control_methods/shuffle_integration_by_cell_type
   - name: control_methods/no_integration
   - name: control_methods/perfect_integration
   - name: methods/harmonypy
@@ -103,16 +103,16 @@ dependencies:
   - name: methods/cytonorm_no_controls_to_goal
   - name: methods/cytonorm_all_controls_to_goal
   - name: methods/cytonorm_one_control_to_goal
-  - name: methods/mnn
   - name: methods/batchadjust_one_control
   - name: methods/batchadjust_all_controls
   - name: methods/rpca_to_goal
   - name: methods/rpca_to_mid
+  - name: methods/cytovi
   - name: metrics/emd
-  # - name: metrics/n_inconsistent_peaks
+  - name: metrics/n_inconsistent_peaks
   - name: metrics/average_batch_r2
   - name: metrics/flowsom_mapping_similarity
-  - name: metrics/cms
+  - name: metrics/lisi
   - name: metrics/bras
 
 runners:
diff --git a/src/workflows/run_benchmark/main.nf b/src/workflows/run_benchmark/main.nf
index f4ed834b..26546f04 100644
--- a/src/workflows/run_benchmark/main.nf
+++ b/src/workflows/run_benchmark/main.nf
@@ -10,8 +10,8 @@ workflow auto {
 // construct list of methods and control methods
 methods = [
   // shuffle_integration,
-  // shuffle_integration_by_batch,
-  // shuffle_integration_by_cell_type,
+  shuffle_integration_by_batch,
+  shuffle_integration_by_cell_type,
   harmonypy,
   limma_remove_batch_effect,
   no_integration,
@@ -24,7 +24,6 @@ methods = [
   cycombine_all_controls_to_mid,
   cycombine_all_controls_to_goal,
   gaussnorm,
-  mnn,
   batchadjust_one_control,
   batchadjust_all_controls,
   cytonorm_no_controls_to_mid,
@@ -34,16 +33,17 @@ methods = [
   cytonorm_all_controls_to_goal,
   cytonorm_one_control_to_goal,
   rpca_to_goal,
-  rpca_to_mid
+  rpca_to_mid,
+  cytovi
 ]
 
 // construct list of metrics
 metrics = [
   emd,
-  // n_inconsistent_peaks,
+  n_inconsistent_peaks,
   average_batch_r2,
   flowsom_mapping_similarity,
-  cms,
+  lisi,
   bras
 ]