From e40a10484e2400ad497b115d9083637ddc7250f7 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Tue, 12 Aug 2025 20:09:38 +1000
Subject: [PATCH 01/46] Update emd vertical when it cannot be calculated

---
 src/metrics/emd/helper.py | 73 +++++++++++++++++++++++----------------
 src/metrics/emd/script.py | 21 +++++++----
 2 files changed, 58 insertions(+), 36 deletions(-)

diff --git a/src/metrics/emd/helper.py b/src/metrics/emd/helper.py
index 94881d94..dd4a481d 100644
--- a/src/metrics/emd/helper.py
+++ b/src/metrics/emd/helper.py
@@ -50,35 +50,48 @@ def calculate_vertical_emd(
         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()
-    )
+    # safeguard
+    mean_emd_global = np.nan
+    max_emd_global = np.nan
+    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 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()
-    )
+        # 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,
@@ -122,10 +135,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()
 
diff --git a/src/metrics/emd/script.py b/src/metrics/emd/script.py
index e953b757..48a388f1 100644
--- a/src/metrics/emd/script.py
+++ b/src/metrics/emd/script.py
@@ -9,9 +9,12 @@
     "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/cytonorm_data_full/cycombine_mid_out.h5ad",
+    # "input_integrated_split2": "resources_test/task_cyto_batch_integration/cytonorm_data_full/cycombine_mid_out.h5ad",
+    # "input_unintegrated": "resources_test/task_cyto_batch_integration/cytonorm_data_full/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"])
@@ -19,6 +22,8 @@
 import helper as emd_helper
 import helper_functions as global_helper
 
+# import src.metrics.emd.helper as emd_helper
+
 print("Reading input files", flush=True)
 
 input_integrated_split1 = ad.read_h5ad(par["input_integrated_split1"])
@@ -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)
 
@@ -63,9 +72,9 @@
 # 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_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()

From 53a8ae54eed851a62247fe5e008c2b4574690183 Mon Sep 17 00:00:00 2001
From: Luqui12 <luca.leomazzi@gmail.com>
Date: Tue, 12 Aug 2025 22:25:17 +0200
Subject: [PATCH 02/46] Added accessory function for subsetting (R)

---
 src/utils/helper_functions.R | 39 ++++++++++++++++++++++++++++++++++++
 1 file changed, 39 insertions(+)

diff --git a/src/utils/helper_functions.R b/src/utils/helper_functions.R
index f5730f20..443255d4 100644
--- a/src/utils/helper_functions.R
+++ b/src/utils/helper_functions.R
@@ -105,3 +105,42 @@ remove_unlabelled <- function(adata) {
         c("unlabelled", "unlabeled")
     adata[!is_unlabelled, ]
 }
+
+#' Subsets the anndata object in a stratified manner by cell type.
+#'
+#' @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) %>%
+    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, ]
+  }
+}
+

From 1d08ddf9a5775e3a7f2fd3d1e6b22f1e39bd3fb7 Mon Sep 17 00:00:00 2001
From: Luqui12 <luca.leomazzi@gmail.com>
Date: Tue, 12 Aug 2025 22:28:10 +0200
Subject: [PATCH 03/46] CMS computed on subsets

---
 src/metrics/cms/script.R | 131 +++++++++++++++++++++++----------------
 1 file changed, 77 insertions(+), 54 deletions(-)

diff --git a/src/metrics/cms/script.R b/src/metrics/cms/script.R
index b05dc8cc..699f3339 100644
--- a/src/metrics/cms/script.R
+++ b/src/metrics/cms/script.R
@@ -23,11 +23,15 @@ meta <- list(
     cpus = NULL
 )
 ## VIASH END
+t0 <- Sys.time()
 
 cores_to_use <- meta$cpus
 if (is.null(cores_to_use)) {
     cores_to_use <- min(5, parallel::detectCores() - 2)
 }
+bpparam <- BiocParallel::MulticoreParam(
+    workers = cores_to_use
+)
 
 source(paste0(meta$resources_dir, "/helper_functions.R"))
 
@@ -37,7 +41,7 @@ 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")
+cat("Fetching metadata from unintegrated\n")
 integrated_split1 <- get_obs_var_for_integrated(
     i_adata = integrated_split1,
     u_adata = unintegrated
@@ -47,78 +51,97 @@ integrated_split2 <- get_obs_var_for_integrated(
     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_distr_split1 <- list()
+medcouples_split1 <- list()
+for (i in 1:5) {
+    cat(paste("Iteration", i, "of 5\n"))
+    integrated_subset <- subset_by_celltype(
+        integrated_split1,
+        frac = 0.3,
+        seed = i
+    )
+    #Transform to SingleCellExperiment and subset markers
+    cat("Transforming to SingleCellExperiment and subsetting markers\n")
+    integrated_subset_sce <- integrated_subset$as_SingleCellExperiment()
+    integrated_subset_sce <- integrated_subset_sce[markers_to_correct, ]
+    cat("Computing Cell Mixing Scores\n")
+    integrated_subset_sce <- CellMixS::cms(
+        integrated_subset_sce,
+        group = "batch",
+        assay_name = "integrated",
+        k = par[["n_neighbors"]],
+        n_dim = par[["n_dim"]],
+        BPPARAM = bpparam
+    )
+    distr <- SingleCellExperiment::colData(integrated_subset_sce)[, "cms"]
+    cms_distr_split1[[paste0("split1_iter_", i)]] <- distr
+    medcouples_split1[[paste0("split1_iter_", i)]] <- robustbase::mc(distr)
+}
 
-cms_mc_split1 <- robustbase::mc(cms_distr_split1)
-cms_mc_split2 <- robustbase::mc(cms_distr_split2)
+cat(paste("Compute Cell Mixing Score using", cores_to_use, "cores for split 2\n"), flush = TRUE)
 
-cms_mc_mean <- mean(c(cms_mc_split1, cms_mc_split2))
+cms_distr_split2 <- list()
+medcouples_split2 <- list()
+for (i in 1:5) {
+    cat(paste("Iteration", i, "of 5\n"))
+    integrated_subset <- subset_by_celltype(
+        integrated_split2,
+        frac = 0.2,
+        seed = i
+    )
+    cat("Transforming to SingleCellExperiment and subsetting markers\n")
+    integrated_subset_sce <- integrated_subset$as_SingleCellExperiment()
+    integrated_subset_sce <- integrated_subset_sce[markers_to_correct, ]
+    cat("Computing Cell Mixing Scores\n")
+    integrated_subset_sce <- CellMixS::cms(
+        integrated_subset_sce,
+        group = "batch",
+        assay_name = "integrated",
+        k = par[["n_neighbors"]],
+        n_dim = par[["n_dim"]],
+        BPPARAM = bpparam
+    )
+    distr <- SingleCellExperiment::colData(integrated_subset_sce)[, "cms"]
+    cms_distr_split2[[paste0("split2_iter_", i)]] <- distr
+    medcouples_split2[[paste0("split2_iter_", i)]] <- robustbase::mc(distr)
+}
 
-cat("Write output AnnData to file\n")
+cat("Aggregate scores\n", flush = TRUE)
+#concat named lists
+cms_distr_list <- c(cms_distr_split1, cms_distr_split2)
+medcouples_list <- c(medcouples_split1, medcouples_split2)
+# Compute mean medcouple
+mean_medcouple_cms <- mean(unlist(medcouples_list))
+
+print("cms_list")
+print(cms_distr_list)
+print("medcouples_list")
+print(medcouples_list)
+print("mean_medcouple_cms")
+print(mean_medcouple_cms)
+
+cat("Write output AnnData to file\n", flush = TRUE)
 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,
+        metric_values = mean_medcouple_cms,
         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
-        )
+        list_medcouples = medcouples_list,
+        cms_distributions = cms_distr_list
     )
 )
 
 output$write_h5ad(par[["output"]], compression = "gzip", mode = "w")
+
+cat(sprintf("Elapsed: %.3f s\n", as.numeric(difftime(Sys.time(), t0, units = "secs"))))
\ No newline at end of file

From d9476500792c9b5301fd79f7fc98609d98846165 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Wed, 13 Aug 2025 09:58:52 +1000
Subject: [PATCH 04/46] disabling seurat rpca for now

---
 src/methods/rpca_to_goal/config.vsh.yaml | 18 +-----------------
 src/methods/rpca_to_mid/config.vsh.yaml  | 16 +---------------
 2 files changed, 2 insertions(+), 32 deletions(-)

diff --git a/src/methods/rpca_to_goal/config.vsh.yaml b/src/methods/rpca_to_goal/config.vsh.yaml
index 88153b9f..66486a29 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:
diff --git a/src/methods/rpca_to_mid/config.vsh.yaml b/src/methods/rpca_to_mid/config.vsh.yaml
index 8f4f26dd..8a76b403 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). 

From b80affe6e7014529d636be20bb67d2448dc70cdb Mon Sep 17 00:00:00 2001
From: Luqui12 <luca.leomazzi@gmail.com>
Date: Wed, 13 Aug 2025 15:40:03 +0200
Subject: [PATCH 05/46] subsetting now is stratified -per cell -per sample.

---
 src/utils/helper_functions.R | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/src/utils/helper_functions.R b/src/utils/helper_functions.R
index 443255d4..dfe5ffa0 100644
--- a/src/utils/helper_functions.R
+++ b/src/utils/helper_functions.R
@@ -106,7 +106,8 @@ remove_unlabelled <- function(adata) {
     adata[!is_unlabelled, ]
 }
 
-#' Subsets the anndata object in a stratified manner by cell type.
+#' 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
@@ -121,7 +122,7 @@ subset_by_celltype <- function(adata, frac = 0.5, seed = 1, anndatar = TRUE) {
   obs$.row <- seq_len(nrow(obs))   # original order
   
   keep_ids <- obs %>%
-    group_by(cell_type) %>%
+    group_by(cell_type, sample) %>%
     slice_sample(prop = frac) %>%
     ungroup() %>%
     arrange(.row) %>%     # restore original order

From 6d3e6692883dde68a3381519ff1aed30e47d7cc8 Mon Sep 17 00:00:00 2001
From: Luqui12 <luca.leomazzi@gmail.com>
Date: Wed, 13 Aug 2025 15:43:53 +0200
Subject: [PATCH 06/46] CMS is now computed on 60% of the data.

---
 src/metrics/cms/config.vsh.yaml |   2 +
 src/metrics/cms/script.R        | 101 ++++++++++++++------------------
 2 files changed, 46 insertions(+), 57 deletions(-)

diff --git a/src/metrics/cms/config.vsh.yaml b/src/metrics/cms/config.vsh.yaml
index 92b9336c..1cec62dc 100644
--- a/src/metrics/cms/config.vsh.yaml
+++ b/src/metrics/cms/config.vsh.yaml
@@ -43,6 +43,8 @@ info:
         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.
 
+        In this implementation we subset to 60% of the total cells in each split of the technical replicates. The subset is stratified by cell type and sample.
+        The mean medcouple between the two splits of technical replicates is used as the final score.
       references:
         doi: 
           - 10.26508/lsa.202001004
diff --git a/src/metrics/cms/script.R b/src/metrics/cms/script.R
index 699f3339..ce57b830 100644
--- a/src/metrics/cms/script.R
+++ b/src/metrics/cms/script.R
@@ -23,7 +23,6 @@ meta <- list(
     cpus = NULL
 )
 ## VIASH END
-t0 <- Sys.time()
 
 cores_to_use <- meta$cpus
 if (is.null(cores_to_use)) {
@@ -58,58 +57,55 @@ cat(paste("Compute Cell Mixing Score using", cores_to_use, "cores for split 1\n"
 
 cms_distr_split1 <- list()
 medcouples_split1 <- list()
-for (i in 1:5) {
-    cat(paste("Iteration", i, "of 5\n"))
-    integrated_subset <- subset_by_celltype(
-        integrated_split1,
-        frac = 0.3,
-        seed = i
-    )
-    #Transform to SingleCellExperiment and subset markers
-    cat("Transforming to SingleCellExperiment and subsetting markers\n")
-    integrated_subset_sce <- integrated_subset$as_SingleCellExperiment()
-    integrated_subset_sce <- integrated_subset_sce[markers_to_correct, ]
-    cat("Computing Cell Mixing Scores\n")
-    integrated_subset_sce <- CellMixS::cms(
-        integrated_subset_sce,
-        group = "batch",
-        assay_name = "integrated",
-        k = par[["n_neighbors"]],
-        n_dim = par[["n_dim"]],
-        BPPARAM = bpparam
-    )
-    distr <- SingleCellExperiment::colData(integrated_subset_sce)[, "cms"]
-    cms_distr_split1[[paste0("split1_iter_", i)]] <- distr
-    medcouples_split1[[paste0("split1_iter_", i)]] <- robustbase::mc(distr)
-}
+
+integrated_subset <- subset_by_celltype(
+    integrated_split1,
+    frac = 0.6,
+    seed = 1
+)
+
+cat("Transforming to SingleCellExperiment and subsetting markers\n")
+integrated_subset_sce <- integrated_subset$as_SingleCellExperiment()
+integrated_subset_sce <- integrated_subset_sce[markers_to_correct, ]
+cat("Computing Cell Mixing Scores\n")
+integrated_subset_sce <- CellMixS::cms(
+    integrated_subset_sce,
+    group = "batch",
+    assay_name = "integrated",
+    k = par[["n_neighbors"]],
+    n_dim = par[["n_dim"]],
+    BPPARAM = bpparam
+)
+distr <- SingleCellExperiment::colData(integrated_subset_sce)[, "cms"]
+cms_distr_split1[["split1"]] <- distr
+medcouples_split1[["split1"]] <- robustbase::mc(distr)
 
 cat(paste("Compute Cell Mixing Score using", cores_to_use, "cores for split 2\n"), flush = TRUE)
 
 cms_distr_split2 <- list()
 medcouples_split2 <- list()
-for (i in 1:5) {
-    cat(paste("Iteration", i, "of 5\n"))
-    integrated_subset <- subset_by_celltype(
-        integrated_split2,
-        frac = 0.2,
-        seed = i
-    )
-    cat("Transforming to SingleCellExperiment and subsetting markers\n")
-    integrated_subset_sce <- integrated_subset$as_SingleCellExperiment()
-    integrated_subset_sce <- integrated_subset_sce[markers_to_correct, ]
-    cat("Computing Cell Mixing Scores\n")
-    integrated_subset_sce <- CellMixS::cms(
-        integrated_subset_sce,
-        group = "batch",
-        assay_name = "integrated",
-        k = par[["n_neighbors"]],
-        n_dim = par[["n_dim"]],
-        BPPARAM = bpparam
-    )
-    distr <- SingleCellExperiment::colData(integrated_subset_sce)[, "cms"]
-    cms_distr_split2[[paste0("split2_iter_", i)]] <- distr
-    medcouples_split2[[paste0("split2_iter_", i)]] <- robustbase::mc(distr)
-}
+
+integrated_subset <- subset_by_celltype(
+    integrated_split2,
+    frac = 0.6,
+    seed = 1
+)
+cat("Transforming to SingleCellExperiment and subsetting markers\n")
+integrated_subset_sce <- integrated_subset$as_SingleCellExperiment()
+integrated_subset_sce <- integrated_subset_sce[markers_to_correct, ]
+cat("Computing Cell Mixing Scores\n")
+integrated_subset_sce <- CellMixS::cms(
+    integrated_subset_sce,
+    group = "batch",
+    assay_name = "integrated",
+    k = par[["n_neighbors"]],
+    n_dim = par[["n_dim"]],
+    BPPARAM = bpparam
+)
+distr <- SingleCellExperiment::colData(integrated_subset_sce)[, "cms"]
+cms_distr_split2[["split2"]] <- distr
+medcouples_split2[["split2"]] <- robustbase::mc(distr)
+
 
 cat("Aggregate scores\n", flush = TRUE)
 #concat named lists
@@ -118,13 +114,6 @@ medcouples_list <- c(medcouples_split1, medcouples_split2)
 # Compute mean medcouple
 mean_medcouple_cms <- mean(unlist(medcouples_list))
 
-print("cms_list")
-print(cms_distr_list)
-print("medcouples_list")
-print(medcouples_list)
-print("mean_medcouple_cms")
-print(mean_medcouple_cms)
-
 cat("Write output AnnData to file\n", flush = TRUE)
 output <- anndataR::AnnData(
     shape = c(0L, 0L),
@@ -143,5 +132,3 @@ output <- anndataR::AnnData(
 )
 
 output$write_h5ad(par[["output"]], compression = "gzip", mode = "w")
-
-cat(sprintf("Elapsed: %.3f s\n", as.numeric(difftime(Sys.time(), t0, units = "secs"))))
\ No newline at end of file

From 5c70c1fa4f62e993aabfdf5905ca7289b066cd85 Mon Sep 17 00:00:00 2001
From: Luqui12 <luca.leomazzi@gmail.com>
Date: Thu, 14 Aug 2025 16:51:16 +0200
Subject: [PATCH 07/46] modified gitignore

---
 .gitignore | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

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

From 189c4fe36321bebabd1f15958beb162b5781684d Mon Sep 17 00:00:00 2001
From: Luqui12 <luca.leomazzi@gmail.com>
Date: Thu, 14 Aug 2025 17:20:48 +0200
Subject: [PATCH 08/46] n_inconsistent_peaks adapted to new schema.

---
 src/metrics/n_inconsistent_peaks/script.py | 84 ++++++++++++----------
 1 file changed, 46 insertions(+), 38 deletions(-)

diff --git a/src/metrics/n_inconsistent_peaks/script.py b/src/metrics/n_inconsistent_peaks/script.py
index 859c104a..1ede3000 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,32 +27,37 @@
     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"])
+    for marker in s1_view.var.index:
+        mexp_integrated = np.array(s1_view[:,marker].layers["integrated"])
+        mexp_validation = np.array(s2_view[:,marker].layers["integrated"])
         density_integrated = get_kde_density(mexp_integrated)
         peaks_integrated = call_peaks(density_integrated)
         density_validation = get_kde_density(mexp_validation)
@@ -64,21 +70,21 @@
 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:
+        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 integrated or validation 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"])
+        for marker in s1_view_ct.var.index:
+            mexp_integrated = np.array(s1_view_ct[:, marker].layers["integrated"])
+            mexp_validation = np.array(s2_view_ct[:, marker].layers["integrated"])
             density_integrated = get_kde_density(mexp_integrated)
             peaks_integrated = call_peaks(density_integrated)
             density_validation = get_kde_density(mexp_validation)
@@ -95,10 +101,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

From ab1c9891f4b8681985ffc5b7ca74effe792d6613 Mon Sep 17 00:00:00 2001
From: Luqui12 <luca.leomazzi@gmail.com>
Date: Thu, 14 Aug 2025 17:30:51 +0200
Subject: [PATCH 09/46] n_inconsistent_peaks adapted to new schema. (2)

---
 src/metrics/n_inconsistent_peaks/helper.py |  6 +++---
 src/metrics/n_inconsistent_peaks/script.py | 16 ++++++++--------
 2 files changed, 11 insertions(+), 11 deletions(-)

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 1ede3000..fd4cd33b 100644
--- a/src/metrics/n_inconsistent_peaks/script.py
+++ b/src/metrics/n_inconsistent_peaks/script.py
@@ -56,11 +56,11 @@
     s2_view = integrated_s2[integrated_s2.obs['donor'] == donor]
 
     for marker in s1_view.var.index:
-        mexp_integrated = np.array(s1_view[:,marker].layers["integrated"])
-        mexp_validation = np.array(s2_view[:,marker].layers["integrated"])
-        density_integrated = get_kde_density(mexp_integrated)
+        mexp_s1 = np.array(s1_view[:,marker].layers["integrated"])
+        mexp_s2 = np.array(s2_view[:,marker].layers["integrated"])
+        density_integrated = get_kde_density(mexp_s1)
         peaks_integrated = call_peaks(density_integrated)
-        density_validation = get_kde_density(mexp_validation)
+        density_validation = get_kde_density(mexp_s2)
         peaks_validation = call_peaks(density_validation)
 
         if peaks_integrated != peaks_validation:
@@ -83,11 +83,11 @@
           continue
 
         for marker in s1_view_ct.var.index:
-            mexp_integrated = np.array(s1_view_ct[:, marker].layers["integrated"])
-            mexp_validation = np.array(s2_view_ct[:, marker].layers["integrated"])
-            density_integrated = get_kde_density(mexp_integrated)
+            mexp_s1 = np.array(s1_view_ct[:, marker].layers["integrated"])
+            mexp_s2 = np.array(s2_view_ct[:, marker].layers["integrated"])
+            density_integrated = get_kde_density(mexp_s1)
             peaks_integrated = call_peaks(density_integrated)
-            density_validation = get_kde_density(mexp_validation)
+            density_validation = get_kde_density(mexp_s2)
             peaks_validation = call_peaks(density_validation)
 
             if peaks_integrated != peaks_validation:

From 4d29a78f3c905a07b8a5a2cecde98ca67a39c8be Mon Sep 17 00:00:00 2001
From: Luqui12 <luca.leomazzi@gmail.com>
Date: Thu, 14 Aug 2025 17:36:40 +0200
Subject: [PATCH 10/46] n_inconsistent_peaks adapted to new schema. (3)

---
 src/metrics/n_inconsistent_peaks/script.py | 26 +++++++++++-----------
 1 file changed, 13 insertions(+), 13 deletions(-)

diff --git a/src/metrics/n_inconsistent_peaks/script.py b/src/metrics/n_inconsistent_peaks/script.py
index fd4cd33b..02f96f86 100644
--- a/src/metrics/n_inconsistent_peaks/script.py
+++ b/src/metrics/n_inconsistent_peaks/script.py
@@ -58,13 +58,13 @@
     for marker in s1_view.var.index:
         mexp_s1 = np.array(s1_view[:,marker].layers["integrated"])
         mexp_s2 = np.array(s2_view[:,marker].layers["integrated"])
-        density_integrated = get_kde_density(mexp_s1)
-        peaks_integrated = call_peaks(density_integrated)
-        density_validation = get_kde_density(mexp_s2)
-        peaks_validation = call_peaks(density_validation)
+        density_s1 = get_kde_density(mexp_s1)
+        peaks_s1 = call_peaks(density_s1)
+        density_s2 = get_kde_density(mexp_s2)
+        peaks_s2 = call_peaks(density_s2)
 
-        if peaks_integrated != peaks_validation:
-            n_inconsistent_peaks += abs(peaks_integrated - peaks_validation)
+        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
@@ -79,19 +79,19 @@
         s2_view_ct = s2_view[s2_view.obs['cell_type'] == celltype]
 
         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 integrated or validation dataset")
+          print(donor,celltype,"skipped cause less than 100 cells are present in either split 1 or split 2 dataset")
           continue
 
         for marker in s1_view_ct.var.index:
             mexp_s1 = np.array(s1_view_ct[:, marker].layers["integrated"])
             mexp_s2 = np.array(s2_view_ct[:, marker].layers["integrated"])
-            density_integrated = get_kde_density(mexp_s1)
-            peaks_integrated = call_peaks(density_integrated)
-            density_validation = get_kde_density(mexp_s2)
-            peaks_validation = call_peaks(density_validation)
+            density_s1 = get_kde_density(mexp_s1)
+            peaks_s1 = call_peaks(density_s1)
+            density_s2 = get_kde_density(mexp_s2)
+            peaks_s2 = call_peaks(density_s2)
 
-            if peaks_integrated != peaks_validation:
-                n_inconsistent_peaks_ct += abs(peaks_integrated - peaks_validation)
+            if peaks_s1 != peaks_s2:
+                n_inconsistent_peaks_ct += abs(peaks_s1 - peaks_s2)
 
 
 

From 3074bec5472a7ac10571689f541a9226595f3792 Mon Sep 17 00:00:00 2001
From: Luqui12 <luca.leomazzi@gmail.com>
Date: Thu, 14 Aug 2025 17:51:39 +0200
Subject: [PATCH 11/46] Added z-score scaling before peak calling to reduce the
 number of false positives. Updated description in config file. metric
 re-enabled (commented statusfield  in config file).

---
 .../n_inconsistent_peaks/config.vsh.yaml      | 13 ++++-----
 src/metrics/n_inconsistent_peaks/script.py    | 28 +++++++++++++------
 2 files changed, 26 insertions(+), 15 deletions(-)

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/script.py b/src/metrics/n_inconsistent_peaks/script.py
index 02f96f86..6e1c858c 100644
--- a/src/metrics/n_inconsistent_peaks/script.py
+++ b/src/metrics/n_inconsistent_peaks/script.py
@@ -56,11 +56,17 @@
     s2_view = integrated_s2[integrated_s2.obs['donor'] == donor]
 
     for marker in s1_view.var.index:
-        mexp_s1 = np.array(s1_view[:,marker].layers["integrated"])
-        mexp_s2 = np.array(s2_view[:,marker].layers["integrated"])
-        density_s1 = get_kde_density(mexp_s1)
+        marker_expression_s1_unscaled = np.array(s1_view[:,marker].layers["integrated"])
+        marker_expression_s2_unscaled = np.array(s2_view[:,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)
+
+        density_s1 = get_kde_density(marker_expression_s1)
         peaks_s1 = call_peaks(density_s1)
-        density_s2 = get_kde_density(mexp_s2)
+        density_s2 = get_kde_density(marker_expression_s2)
         peaks_s2 = call_peaks(density_s2)
 
         if peaks_s1 != peaks_s2:
@@ -83,11 +89,17 @@
           continue
 
         for marker in s1_view_ct.var.index:
-            mexp_s1 = np.array(s1_view_ct[:, marker].layers["integrated"])
-            mexp_s2 = np.array(s2_view_ct[:, marker].layers["integrated"])
-            density_s1 = get_kde_density(mexp_s1)
+            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)
+
+            density_s1 = get_kde_density(marker_expression_s1)
             peaks_s1 = call_peaks(density_s1)
-            density_s2 = get_kde_density(mexp_s2)
+            density_s2 = get_kde_density(marker_expression_s2)
             peaks_s2 = call_peaks(density_s2)
 
             if peaks_s1 != peaks_s2:

From 5a20459e7103f57cc5a9e5ea4ab3972f2db5eb20 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Fri, 15 Aug 2025 10:41:50 +1000
Subject: [PATCH 12/46] update emd scripts

---
 src/metrics/emd/helper.py | 43 +++++++++++++++++++++++++++++++++-----
 src/metrics/emd/script.py | 44 +++++++++++----------------------------
 2 files changed, 50 insertions(+), 37 deletions(-)

diff --git a/src/metrics/emd/helper.py b/src/metrics/emd/helper.py
index dd4a481d..417fbe65 100644
--- a/src/metrics/emd/helper.py
+++ b/src/metrics/emd/helper.py
@@ -267,8 +267,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."
             )
@@ -287,9 +287,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
 
@@ -424,3 +426,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 48a388f1..80dd447b 100644
--- a/src/metrics/emd/script.py
+++ b/src/metrics/emd/script.py
@@ -9,9 +9,9 @@
     "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/cytonorm_data_full/cycombine_mid_out.h5ad",
-    # "input_integrated_split2": "resources_test/task_cyto_batch_integration/cytonorm_data_full/cycombine_mid_out.h5ad",
-    # "input_unintegrated": "resources_test/task_cyto_batch_integration/cytonorm_data_full/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/"}
@@ -19,11 +19,11 @@
 
 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
 
-# import src.metrics.emd.helper as emd_helper
-
 print("Reading input files", flush=True)
 
 input_integrated_split1 = ad.read_h5ad(par["input_integrated_split1"])
@@ -63,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,

From 69c61c3e8661e6ab8880ce313018b42311b52f27 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Tue, 19 Aug 2025 09:51:03 +1000
Subject: [PATCH 13/46] re-enabling rpca and n peaks

---
 src/methods/rpca_to_goal/config.vsh.yaml    | 2 +-
 src/methods/rpca_to_mid/config.vsh.yaml     | 2 +-
 src/workflows/run_benchmark/config.vsh.yaml | 2 +-
 src/workflows/run_benchmark/main.nf         | 2 +-
 4 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/methods/rpca_to_goal/config.vsh.yaml b/src/methods/rpca_to_goal/config.vsh.yaml
index 66486a29..d17aed78 100644
--- a/src/methods/rpca_to_goal/config.vsh.yaml
+++ b/src/methods/rpca_to_goal/config.vsh.yaml
@@ -1,6 +1,6 @@
 __merge__: ../../api/comp_method.yaml
 name: rpca_to_goal
-status: disabled
+# status: disabled
 label: Seurat RPCA (to-goal)
 summary: "Batch integrate data to a goal batch using mutual nearest neighbors identified via Seurat reciprocal PCA."
 description: |
diff --git a/src/methods/rpca_to_mid/config.vsh.yaml b/src/methods/rpca_to_mid/config.vsh.yaml
index 8a76b403..5afe2488 100644
--- a/src/methods/rpca_to_mid/config.vsh.yaml
+++ b/src/methods/rpca_to_mid/config.vsh.yaml
@@ -1,6 +1,6 @@
 __merge__: ../../api/comp_method.yaml
 name: rpca_to_mid
-status: disabled
+# status: disabled
 label: Seurat RPCA (to-middle)
 summary: "Batch integrate data to a midpoint using mutual nearest neighbors identified via Seurat reciprocal PCA."
 description: |
diff --git a/src/workflows/run_benchmark/config.vsh.yaml b/src/workflows/run_benchmark/config.vsh.yaml
index 3968de6f..d48edc12 100644
--- a/src/workflows/run_benchmark/config.vsh.yaml
+++ b/src/workflows/run_benchmark/config.vsh.yaml
@@ -109,7 +109,7 @@ dependencies:
   - name: methods/rpca_to_goal
   - name: methods/rpca_to_mid
   - 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
diff --git a/src/workflows/run_benchmark/main.nf b/src/workflows/run_benchmark/main.nf
index f4ed834b..4c443f61 100644
--- a/src/workflows/run_benchmark/main.nf
+++ b/src/workflows/run_benchmark/main.nf
@@ -40,7 +40,7 @@ methods = [
 // construct list of metrics
 metrics = [
   emd,
-  // n_inconsistent_peaks,
+  n_inconsistent_peaks,
   average_batch_r2,
   flowsom_mapping_similarity,
   cms,

From 7e7a7d107f1f61ae4d2f621dd04664f0978cb107 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Tue, 19 Aug 2025 23:36:00 +1000
Subject: [PATCH 14/46] increase resourhces to mnn

---
 src/methods/mnn/config.vsh.yaml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/methods/mnn/config.vsh.yaml b/src/methods/mnn/config.vsh.yaml
index e601c4bb..e58a86cc 100644
--- a/src/methods/mnn/config.vsh.yaml
+++ b/src/methods/mnn/config.vsh.yaml
@@ -104,4 +104,4 @@ runners:
   # Allows turning the component into a Nextflow module / pipeline.
   - type: nextflow
     directives:
-      label: [midtime,midmem,midcpu]
+      label: [hightime,midmem,midcpu]

From e15b897b9a3c7f0018cb47ba8d6ed0b8e2121010 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Tue, 19 Aug 2025 23:36:11 +1000
Subject: [PATCH 15/46] patch batchadjust

---
 src/methods/batchadjust_all_controls/script.R | 48 +++++++++++++++++--
 src/methods/batchadjust_one_control/script.R  | 48 +++++++++++++++++--
 2 files changed, 88 insertions(+), 8 deletions(-)

diff --git a/src/methods/batchadjust_all_controls/script.R b/src/methods/batchadjust_all_controls/script.R
index 350ed2d9..72077906 100644
--- a/src/methods/batchadjust_all_controls/script.R
+++ b/src/methods/batchadjust_all_controls/script.R
@@ -33,7 +33,47 @@ 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])
+
+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])
+
+}
 
 cat("Split cells\n")
 input_controls <- input[input$obs$is_control != 0, ]
@@ -49,8 +89,8 @@ print(input_no_controls)
 #avoid NA due to invalid factor level
 input_controls$obs$sample <- as.character(input_controls$obs$sample)
 # 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 +112,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,
diff --git a/src/methods/batchadjust_one_control/script.R b/src/methods/batchadjust_one_control/script.R
index 0ce78a70..ad32d08c 100644
--- a/src/methods/batchadjust_one_control/script.R
+++ b/src/methods/batchadjust_one_control/script.R
@@ -33,7 +33,47 @@ 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])
+
+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])
+
+}
 
 cat("Split cells\n")
 input_controls <- input[input$obs$is_control == 1, ]
@@ -49,8 +89,8 @@ print(input_no_controls)
 #avoid NA due to invalid factor level
 input_controls$obs$sample <- as.character(input_controls$obs$sample)
 # 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 +112,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,

From a46b43dc1b2e77a8a0e1c24c95e0c2c777de62e0 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Tue, 19 Aug 2025 23:42:36 +1000
Subject: [PATCH 16/46] increase global size limit

---
 src/methods/rpca_to_goal/script.R | 2 ++
 src/methods/rpca_to_mid/script.R  | 3 +++
 2 files changed, 5 insertions(+)

diff --git a/src/methods/rpca_to_goal/script.R b/src/methods/rpca_to_goal/script.R
index 8a9fd908..f985b9a0 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 = 5 * 1024^3)  # 5 GiB
+
 cat("Reading input files\n")
 input_adata <- anndata::read_h5ad(par[["input"]])
 
diff --git a/src/methods/rpca_to_mid/script.R b/src/methods/rpca_to_mid/script.R
index d84d26e0..2fcac506 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 = 5 * 1024^3)  # 5 GiB
+
+
 cat("Reading input files\n")
 input_adata <- anndata::read_h5ad(par[["input"]])
 

From 25dd67589c1a53750005862527ccc6128fd162e4 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Wed, 20 Aug 2025 02:17:05 +1000
Subject: [PATCH 17/46] update batchadjust again

---
 .../batchadjust_all_controls/config.vsh.yaml  |  1 +
 src/methods/batchadjust_all_controls/script.R | 51 ++++--------------
 src/methods/batchadjust_all_controls/utils.R  | 51 ++++++++++++++++++
 .../batchadjust_one_control/config.vsh.yaml   |  1 +
 src/methods/batchadjust_one_control/script.R  | 53 ++++---------------
 src/methods/batchadjust_one_control/utils.R   | 52 ++++++++++++++++++
 6 files changed, 126 insertions(+), 83 deletions(-)
 create mode 100644 src/methods/batchadjust_all_controls/utils.R
 create mode 100644 src/methods/batchadjust_one_control/utils.R

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 72077906..a8a845a4 100644
--- a/src/methods/batchadjust_all_controls/script.R
+++ b/src/methods/batchadjust_all_controls/script.R
@@ -2,8 +2,9 @@ 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 = "output.h5ad",
+  percentile = as.integer('80')
 )
 meta <- list(
   name = "batchadjust_all_controls",
@@ -11,6 +12,7 @@ meta <- list(
 )
 ## 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"))
 
@@ -35,45 +37,7 @@ 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])
 
-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])
-
-}
+input <- add_original_id(input)
 
 cat("Split cells\n")
 input_controls <- input[input$obs$is_control != 0, ]
@@ -88,6 +52,11 @@ 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] <- "Batch1_anchor"
 input_controls$obs$sample[input_controls$obs$batch == 2] <- "Batch2_anchor"
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 ad32d08c..4189456c 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"))
 
@@ -35,45 +38,7 @@ 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])
 
-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])
-
-}
+input <- add_original_id(input)
 
 cat("Split cells\n")
 input_controls <- input[input$obs$is_control == 1, ]
@@ -88,6 +53,10 @@ 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] <- "Batch1_anchor"
 input_controls$obs$sample[input_controls$obs$batch == 2] <- "Batch2_anchor"
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)
+}

From 6d5550f51203ee20a1a97245893df1f936219d85 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Wed, 20 Aug 2025 02:18:14 +1000
Subject: [PATCH 18/46] increase ram max size for rpca

---
 src/methods/rpca_to_goal/script.R | 2 +-
 src/methods/rpca_to_mid/script.R  | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/methods/rpca_to_goal/script.R b/src/methods/rpca_to_goal/script.R
index f985b9a0..4bc75551 100644
--- a/src/methods/rpca_to_goal/script.R
+++ b/src/methods/rpca_to_goal/script.R
@@ -13,7 +13,7 @@ meta <- list(
 )
 ## VIASH END
 
-options(future.globals.maxSize = 5 * 1024^3)  # 5 GiB
+options(future.globals.maxSize = 8 * 1024^3)  # 8 GiB
 
 cat("Reading input files\n")
 input_adata <- anndata::read_h5ad(par[["input"]])
diff --git a/src/methods/rpca_to_mid/script.R b/src/methods/rpca_to_mid/script.R
index 2fcac506..366d6322 100644
--- a/src/methods/rpca_to_mid/script.R
+++ b/src/methods/rpca_to_mid/script.R
@@ -13,7 +13,7 @@ meta <- list(
 )
 ## VIASH END
 
-options(future.globals.maxSize = 5 * 1024^3)  # 5 GiB
+options(future.globals.maxSize = 8 * 1024^3)  # 8 GiB
 
 
 cat("Reading input files\n")

From 8af2119ff0c6ca4ab59dcdad0867dba058099517 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Wed, 20 Aug 2025 12:26:11 +1000
Subject: [PATCH 19/46] patch bug to remove Original_ID if it wasn't there in
 the first place

---
 src/methods/batchadjust_all_controls/script.R | 13 +++++++++++--
 src/methods/batchadjust_one_control/script.R  |  7 +++++++
 2 files changed, 18 insertions(+), 2 deletions(-)

diff --git a/src/methods/batchadjust_all_controls/script.R b/src/methods/batchadjust_all_controls/script.R
index a8a845a4..74ba4e0f 100644
--- a/src/methods/batchadjust_all_controls/script.R
+++ b/src/methods/batchadjust_all_controls/script.R
@@ -3,13 +3,15 @@ library(flowCore)
 ## VIASH START
 par <- list(
   input = "resources_test/debug/batchadjust/_viash_par/input_1/censored_split1.h5ad",
-  output = "output.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"))
@@ -37,6 +39,8 @@ 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])
 
+original_id_in_var <- "Original_ID" %in% input$var_names
+
 input <- add_original_id(input)
 
 cat("Split cells\n")
@@ -105,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_one_control/script.R b/src/methods/batchadjust_one_control/script.R
index 4189456c..ca609d37 100644
--- a/src/methods/batchadjust_one_control/script.R
+++ b/src/methods/batchadjust_one_control/script.R
@@ -38,6 +38,8 @@ 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])
 
+original_id_in_var <- "Original_ID" %in% input$var_names
+
 input <- add_original_id(input)
 
 cat("Split cells\n")
@@ -105,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)],

From b2de53e01e84fe0e396ffcd62a3179c336a7c15d Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Wed, 20 Aug 2025 13:41:30 +1000
Subject: [PATCH 20/46] increase rpca global maxsize again

---
 src/methods/rpca_to_goal/script.R | 10 +++++-----
 src/methods/rpca_to_mid/script.R  | 10 +++++-----
 2 files changed, 10 insertions(+), 10 deletions(-)

diff --git a/src/methods/rpca_to_goal/script.R b/src/methods/rpca_to_goal/script.R
index 4bc75551..35f562f5 100644
--- a/src/methods/rpca_to_goal/script.R
+++ b/src/methods/rpca_to_goal/script.R
@@ -13,7 +13,7 @@ meta <- list(
 )
 ## VIASH END
 
-options(future.globals.maxSize = 8 * 1024^3)  # 8 GiB
+options(future.globals.maxSize = 25 * 1024^3)  # 25 GiB
 
 cat("Reading input files\n")
 input_adata <- anndata::read_h5ad(par[["input"]])
@@ -61,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
@@ -74,7 +74,7 @@ seurat_objs <- lapply(batches, function(batch) {
         assay = "cyto",
         npcs = par[["npcs"]],
         approx = FALSE,
-        verbose = FALSE
+        verbose = TRUE
     )
 
     return(seurat_obj)
@@ -96,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")
 )
 
@@ -110,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/script.R b/src/methods/rpca_to_mid/script.R
index 366d6322..7451eb6c 100644
--- a/src/methods/rpca_to_mid/script.R
+++ b/src/methods/rpca_to_mid/script.R
@@ -13,7 +13,7 @@ meta <- list(
 )
 ## VIASH END
 
-options(future.globals.maxSize = 8 * 1024^3)  # 8 GiB
+options(future.globals.maxSize = 25 * 1024^3)  # 25 GiB
 
 
 cat("Reading input files\n")
@@ -62,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
@@ -75,7 +75,7 @@ seurat_objs <- lapply(batches, function(batch) {
         assay = "cyto",
         npcs = par[["npcs"]],
         approx = FALSE,
-        verbose = FALSE
+        verbose = TRUE
     )
 
     return(seurat_obj)
@@ -99,7 +99,7 @@ anchors <- Seurat::FindIntegrationAnchors(
     dims = seq(npcs_computed),
     k.anchor = par[["n_neighbours"]],
     reduction = "rpca",
-    verbose = FALSE,
+    verbose = TRUE,
     reference = NULL
 )
 
@@ -109,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"

From 5ea510a6d31b1c5277ebf876bdd713a70eaabb87 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Wed, 20 Aug 2025 20:10:21 +1000
Subject: [PATCH 21/46] update rpca to hightime

---
 src/methods/rpca_to_goal/config.vsh.yaml | 2 +-
 src/methods/rpca_to_mid/config.vsh.yaml  | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/methods/rpca_to_goal/config.vsh.yaml b/src/methods/rpca_to_goal/config.vsh.yaml
index d17aed78..8970b041 100644
--- a/src/methods/rpca_to_goal/config.vsh.yaml
+++ b/src/methods/rpca_to_goal/config.vsh.yaml
@@ -84,4 +84,4 @@ runners:
   # Allows turning the component into a Nextflow module / pipeline.
   - type: nextflow
     directives:
-      label: [midtime,midmem,midcpu]
+      label: [hightime,midmem,midcpu]
diff --git a/src/methods/rpca_to_mid/config.vsh.yaml b/src/methods/rpca_to_mid/config.vsh.yaml
index 5afe2488..c0e28c29 100644
--- a/src/methods/rpca_to_mid/config.vsh.yaml
+++ b/src/methods/rpca_to_mid/config.vsh.yaml
@@ -88,4 +88,4 @@ runners:
   # Allows turning the component into a Nextflow module / pipeline.
   - type: nextflow
     directives:
-      label: [midtime,midmem,midcpu]
+      label: [hightime,midmem,midcpu]

From a584eef7089203069e6069cc9c6740db0547cd98 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Tue, 26 Aug 2025 21:32:06 +1000
Subject: [PATCH 22/46] fix bug in perfect integration get obs

---
 .../perfect_integration/script.py             |  4 ++--
 src/metrics/cms/script.R                      |  6 ++++--
 .../flowsom_mapping_similarity/script.R       | 13 ++++++------
 src/utils/helper_functions.R                  | 20 ++++++++++++++-----
 4 files changed, 28 insertions(+), 15 deletions(-)

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/metrics/cms/script.R b/src/metrics/cms/script.R
index ce57b830..b902e341 100644
--- a/src/metrics/cms/script.R
+++ b/src/metrics/cms/script.R
@@ -43,11 +43,13 @@ integrated_split2 <- anndataR::read_h5ad(par[["input_integrated_split2"]])
 cat("Fetching metadata from unintegrated\n")
 integrated_split1 <- get_obs_var_for_integrated(
     i_adata = integrated_split1,
-    u_adata = unintegrated
+    u_adata = unintegrated,
+    split_id = 1
 )
 integrated_split2 <- get_obs_var_for_integrated(
     i_adata = integrated_split2,
-    u_adata = unintegrated
+    u_adata = unintegrated,
+    split_id = 2
 )
 
 # Get markers to correct
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/utils/helper_functions.R b/src/utils/helper_functions.R
index dfe5ffa0..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)
 

From 9fbefbbd3aeca5a77c74bcbdcf659b04424c00ab Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Fri, 29 Aug 2025 08:39:40 +1000
Subject: [PATCH 23/46] increase mnn and rpca goal to very high

---
 src/methods/mnn/config.vsh.yaml          | 2 +-
 src/methods/rpca_to_goal/config.vsh.yaml | 2 +-
 src/methods/rpca_to_mid/config.vsh.yaml  | 2 +-
 3 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/methods/mnn/config.vsh.yaml b/src/methods/mnn/config.vsh.yaml
index e58a86cc..bb4cfb61 100644
--- a/src/methods/mnn/config.vsh.yaml
+++ b/src/methods/mnn/config.vsh.yaml
@@ -104,4 +104,4 @@ runners:
   # Allows turning the component into a Nextflow module / pipeline.
   - type: nextflow
     directives:
-      label: [hightime,midmem,midcpu]
+      label: [veryhightime,midmem,midcpu]
diff --git a/src/methods/rpca_to_goal/config.vsh.yaml b/src/methods/rpca_to_goal/config.vsh.yaml
index 8970b041..200dfa68 100644
--- a/src/methods/rpca_to_goal/config.vsh.yaml
+++ b/src/methods/rpca_to_goal/config.vsh.yaml
@@ -84,4 +84,4 @@ runners:
   # Allows turning the component into a Nextflow module / pipeline.
   - type: nextflow
     directives:
-      label: [hightime,midmem,midcpu]
+      label: [veryhightime,midmem,midcpu]
diff --git a/src/methods/rpca_to_mid/config.vsh.yaml b/src/methods/rpca_to_mid/config.vsh.yaml
index c0e28c29..ca01b9cb 100644
--- a/src/methods/rpca_to_mid/config.vsh.yaml
+++ b/src/methods/rpca_to_mid/config.vsh.yaml
@@ -88,4 +88,4 @@ runners:
   # Allows turning the component into a Nextflow module / pipeline.
   - type: nextflow
     directives:
-      label: [hightime,midmem,midcpu]
+      label: [veryhightime,midmem,midcpu]

From 18a7440ffb76b02ebd3f03f76a7ee6fb7140fb2a Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Sat, 30 Aug 2025 09:56:05 +1000
Subject: [PATCH 24/46] update cycombine

---
 src/methods/cycombine_all_controls_to_mid/config.vsh.yaml | 1 +
 src/methods/cycombine_no_controls_to_goal/config.vsh.yaml | 1 +
 src/methods/cycombine_no_controls_to_mid/config.vsh.yaml  | 1 +
 src/methods/cycombine_one_control_to_goal/config.vsh.yaml | 1 +
 src/methods/cycombine_one_control_to_mid/config.vsh.yaml  | 1 +
 5 files changed, 5 insertions(+)

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:

From d715b27ed45530123a59a68e5ba4f1bc766e22ef Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Sat, 30 Aug 2025 16:09:30 +1000
Subject: [PATCH 25/46] add pbmcapply as dependency

---
 src/methods/cycombine_all_controls_to_goal/config.vsh.yaml | 1 +
 1 file changed, 1 insertion(+)

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:

From 6ed68c4c7de0a56892cbd230f8090a3b631378d8 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Fri, 12 Sep 2025 09:34:15 +1000
Subject: [PATCH 26/46] update bras metric

---
 src/metrics/bras/script.py | 58 ++++++++++++++++++++------------------
 1 file changed, 30 insertions(+), 28 deletions(-)

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")

From 18ebe83b5e652dcb17b4e48d2f0bd7cddff9d729 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Tue, 23 Sep 2025 11:51:19 +1000
Subject: [PATCH 27/46] adding cytovi

---
 src/methods/cytovi/config.vsh.yaml          | 61 ++++++++++++++++++
 src/methods/cytovi/script.py                | 70 +++++++++++++++++++++
 src/workflows/run_benchmark/config.vsh.yaml |  1 +
 src/workflows/run_benchmark/main.nf         |  3 +-
 4 files changed, 134 insertions(+), 1 deletion(-)
 create mode 100644 src/methods/cytovi/config.vsh.yaml
 create mode 100644 src/methods/cytovi/script.py

diff --git a/src/methods/cytovi/config.vsh.yaml b/src/methods/cytovi/config.vsh.yaml
new file mode 100644
index 00000000..e48f12d8
--- /dev/null
+++ b/src/methods/cytovi/config.vsh.yaml
@@ -0,0 +1,61 @@
+# 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_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
+  # Additional resources your script needs (optional)
+  # - type: file
+  #   path: weights.pt
+
+engines:
+  - type: docker
+    image: openproblems/base_python:1
+    setup:
+      - type: python
+        packages:
+          - scvi-tools
+
+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/cytovi/script.py b/src/methods/cytovi/script.py
new file mode 100644
index 00000000..31dfbc14
--- /dev/null
+++ b/src/methods/cytovi/script.py
@@ -0,0 +1,70 @@
+import anndata as ad
+import numpy as np
+from scvi.external import cytovi
+
+## 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",
+}
+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()
+
+# scale data
+cytovi.scale(
+    adata=adata_to_correct, transformed_layer_key="preprocessed", batch_key="batch_str"
+)
+
+print("Run CytoVI", flush=True)
+
+cytovi.CYTOVI.setup_anndata(adata_to_correct, layer="scaled", batch_key="batch_str")
+model = cytovi.CYTOVI(adata_to_correct)
+model.train()
+
+# get batch corrected data
+corrected_data = model.get_normalized_expression()
+
+# have to add in the uncorrected markers as well
+uncorrected_data = adata[:, markers_not_correct].layers["preprocessed"]
+
+out_matrix = np.concatenate([corrected_data.to_numpy(), 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]
+
+# run umap for quick check
+# import scanpy as sc
+# test_adata = out_adata.copy()
+# test_adata = test_adata[:, markers_to_correct]
+# test_adata.X = test_adata.layers["integrated"]
+# test_adata.obs = adata.obs
+# 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/workflows/run_benchmark/config.vsh.yaml b/src/workflows/run_benchmark/config.vsh.yaml
index d48edc12..a4d6dbae 100644
--- a/src/workflows/run_benchmark/config.vsh.yaml
+++ b/src/workflows/run_benchmark/config.vsh.yaml
@@ -108,6 +108,7 @@ dependencies:
   - 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/average_batch_r2
diff --git a/src/workflows/run_benchmark/main.nf b/src/workflows/run_benchmark/main.nf
index 4c443f61..673add7e 100644
--- a/src/workflows/run_benchmark/main.nf
+++ b/src/workflows/run_benchmark/main.nf
@@ -34,7 +34,8 @@ 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

From d05262bc990d2bcd574c5cfbd1ba79b036fc2a98 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Tue, 23 Sep 2025 14:23:31 +1000
Subject: [PATCH 28/46] use gpu image for cytovi

---
 src/methods/cytovi/config.vsh.yaml | 24 ++++++++++++++----------
 src/methods/cytovi/script.py       |  6 +++++-
 2 files changed, 19 insertions(+), 11 deletions(-)

diff --git a/src/methods/cytovi/config.vsh.yaml b/src/methods/cytovi/config.vsh.yaml
index e48f12d8..7bb07279 100644
--- a/src/methods/cytovi/config.vsh.yaml
+++ b/src/methods/cytovi/config.vsh.yaml
@@ -29,12 +29,16 @@ links:
   repository: https://github.com/YosefLab/cytovi-reference-implementation
 
 # Component-specific parameters (optional)
-# arguments:
-#   - name: "--n_neighbors"
-#     type: "integer"
-#     default: 5
-#     description: Number of neighbors to use.
-
+arguments:
+  - name: --n_hidden
+    type: integer
+    default: 128
+    description: Number of hidden units.
+  - name: --n_layers
+    type: integer
+    default: 1
+    description: Number of layers.
+  
 # Resources required to run the component
 resources:
   # The script of your component (required)
@@ -46,11 +50,11 @@ resources:
 
 engines:
   - type: docker
-    image: openproblems/base_python:1
+    image: openproblems/base_pytorch_nvidia:1
     setup:
       - type: python
-        packages:
-          - scvi-tools
+        pypi:
+          - scvi-tools>=1.4.0
 
 runners:
   # This platform allows running the component natively
@@ -58,4 +62,4 @@ runners:
   # Allows turning the component into a Nextflow module / pipeline.
   - type: nextflow
     directives:
-      label: [midtime,midmem,midcpu]
+      label: [hightime, midmem, lowcpu, gpu]
diff --git a/src/methods/cytovi/script.py b/src/methods/cytovi/script.py
index 31dfbc14..89ad9028 100644
--- a/src/methods/cytovi/script.py
+++ b/src/methods/cytovi/script.py
@@ -6,6 +6,8 @@
 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,
 }
 meta = {"name": "cytovi"}
 ## VIASH END
@@ -28,7 +30,9 @@
 print("Run CytoVI", flush=True)
 
 cytovi.CYTOVI.setup_anndata(adata_to_correct, layer="scaled", batch_key="batch_str")
-model = cytovi.CYTOVI(adata_to_correct)
+model = cytovi.CYTOVI(
+    adata=adata_to_correct, n_hidden=par["n_hidden"], n_layers=par["n_layers"]
+)
 model.train()
 
 # get batch corrected data

From 1bac1fdd6333b163e65fb32b2f4209a120f80a43 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Tue, 23 Sep 2025 14:42:03 +1000
Subject: [PATCH 29/46] cannot use gpu

---
 src/methods/cytovi/config.vsh.yaml | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/methods/cytovi/config.vsh.yaml b/src/methods/cytovi/config.vsh.yaml
index 7bb07279..0e34c2a2 100644
--- a/src/methods/cytovi/config.vsh.yaml
+++ b/src/methods/cytovi/config.vsh.yaml
@@ -50,7 +50,7 @@ resources:
 
 engines:
   - type: docker
-    image: openproblems/base_pytorch_nvidia:1
+    image: openproblems/base_python:1
     setup:
       - type: python
         pypi:
@@ -62,4 +62,4 @@ runners:
   # Allows turning the component into a Nextflow module / pipeline.
   - type: nextflow
     directives:
-      label: [hightime, midmem, lowcpu, gpu]
+      label: [hightime, midmem, midcpu]

From be5f1c616ba842125f5ca25f14a10204c28fecd5 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Tue, 23 Sep 2025 23:50:37 +1000
Subject: [PATCH 30/46] used gpu for cytovi

---
 src/methods/cytovi/config.vsh.yaml | 17 ++++++++++++++---
 1 file changed, 14 insertions(+), 3 deletions(-)

diff --git a/src/methods/cytovi/config.vsh.yaml b/src/methods/cytovi/config.vsh.yaml
index 0e34c2a2..ec6fa86c 100644
--- a/src/methods/cytovi/config.vsh.yaml
+++ b/src/methods/cytovi/config.vsh.yaml
@@ -50,11 +50,22 @@ resources:
 
 engines:
   - type: docker
-    image: openproblems/base_python:1
+    image: nvcr.io/nvidia/pytorch:25.03-py3
     setup:
+      - type: apt
+        packages:
+          - procps
+          - git
       - type: python
-        pypi:
+        packages:
+          - anndata~=0.11.0
+          - scanpy~=1.11.0
+          - pyyaml
+          - requests
+          - jsonschema
           - scvi-tools>=1.4.0
+        github:
+          - openproblems-bio/core#subdirectory=packages/python/openproblems
 
 runners:
   # This platform allows running the component natively
@@ -62,4 +73,4 @@ runners:
   # Allows turning the component into a Nextflow module / pipeline.
   - type: nextflow
     directives:
-      label: [hightime, midmem, midcpu]
+      label: [hightime, highmem, midcpu, gpu]

From def6ad68b19832afe9ed1c42e4a0642ce23595c6 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Wed, 24 Sep 2025 00:04:35 +1000
Subject: [PATCH 31/46] add mnnpy

---
 src/methods/mnnpy/config.vsh.yaml           | 79 +++++++++++++++++++++
 src/methods/mnnpy/script.py                 | 71 ++++++++++++++++++
 src/workflows/run_benchmark/config.vsh.yaml |  1 +
 src/workflows/run_benchmark/main.nf         |  3 +-
 4 files changed, 153 insertions(+), 1 deletion(-)
 create mode 100644 src/methods/mnnpy/config.vsh.yaml
 create mode 100644 src/methods/mnnpy/script.py

diff --git a/src/methods/mnnpy/config.vsh.yaml b/src/methods/mnnpy/config.vsh.yaml
new file mode 100644
index 00000000..acbfe6b9
--- /dev/null
+++ b/src/methods/mnnpy/config.vsh.yaml
@@ -0,0 +1,79 @@
+# 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: mnnpy
+# A relatively short label, used when rendering visualisations (required)
+label: mnnpy
+# A one sentence summary of how this method works (required). Used when 
+# rendering summary tables.
+summary: "Batch effect correction by matching mutual nearest neighbors, Python implementation."
+# A multi-line description of how this component works (required). Used
+# when rendering reference documentation.
+description: |
+  An implementation of MNN correct in python featuring low memory usage, full multicore support and compatibility with the scanpy framework.
+
+  Batch effect correction by matching mutual nearest neighbors (Haghverdi et al, 2018) has been implemented as a function 'mnnCorrect' in the R package scran. Sadly it's extremely slow for big datasets and doesn't make full use of the parallel architecture of modern CPUs.
+
+  This project is a python implementation of the MNN correct algorithm which takes advantage of python's extendability and hackability. It seamlessly integrates with the scanpy framework and has multicore support in its bones.
+references:
+  bibtex: |
+    @misc{Kang2022,
+      author = {Kang, Chris},
+      title = {mnnpy},
+      year = {Kang2022},
+      publisher = {GitHub},
+      journal = {GitHub repository},
+      howpublished = {\url{https://github.com/chriscainx/mnnpy}},
+      commit = {2097dec30c193f036c5ed7e1c3d1e3a6270e102b}
+    }
+links:
+  repository: https://github.com/chriscainx/mnnpy
+  documentation: https://github.com/chriscainx/mnnpy#readme
+
+# 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
+
+engines:
+  - type: docker
+    image: python:3.8
+    setup:
+      - type: apt
+        packages:
+          - procps
+          - libhdf5-dev
+      - type: python
+        pypi:
+          - h5py==3.11.0
+          - anndata~=0.8.0
+          - scanpy
+          - pyyaml
+          - requests
+          - jsonschema
+      - type: python
+        pypi:
+          - git+https://github.com/openproblems-bio/core#subdirectory=packages/python/openproblems
+          - --ignore-requires-python
+      - type: python
+        github:
+          - chriscainx/mnnpy
+runners:
+  - type: executable
+  - type: nextflow
+    directives:
+      label: [hightime, midcpu, highmem]
\ No newline at end of file
diff --git a/src/methods/mnnpy/script.py b/src/methods/mnnpy/script.py
new file mode 100644
index 00000000..a059f464
--- /dev/null
+++ b/src/methods/mnnpy/script.py
@@ -0,0 +1,71 @@
+import anndata as ad
+import mnnpy
+import numpy as np
+
+## VIASH START
+par = {
+    "input": "resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/censored_split1.h5ad",
+    "output": "resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/output_mnnpy.h5ad",
+}
+meta = {"name": "mnnpy"}
+## VIASH END
+
+print("Read input", flush=True)
+adata = ad.read_h5ad(par["input"])
+
+adata.X = adata.layers["preprocessed"]
+
+# convert batch to category as otherwise mnnpy won't work..
+adata.obs["batch_cat"] = 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("Run mnn", flush=True)
+split = []
+batch_categories = adata_to_correct.obs["batch_cat"].unique().tolist()
+
+for i in batch_categories:
+    split.append(adata_to_correct[adata_to_correct.obs["batch_cat"] == i].copy())
+
+corrected, _, _ = mnnpy.mnn_correct(
+    *split, batch_key="batch", batch_categories=batch_categories, index_unique=None
+)
+
+# have to add in the uncorrected markers as well
+uncorrected_data = adata[:, markers_not_correct].layers["preprocessed"]
+
+out_matrix = np.concatenate([corrected.X, uncorrected_data], axis=1)
+out_var_idx = np.concatenate([corrected.var.index, 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]
+
+# run umap for quick check
+# import scanpy as sc
+# test_adata = out_adata.copy()
+# test_adata = test_adata[:, markers_to_correct]
+# test_adata.X = test_adata.layers["integrated"]
+# test_adata.obs = adata.obs
+# sc.pp.neighbors(test_adata, use_rep="X")
+# sc.tl.umap(test_adata)
+# sc.pl.umap(test_adata, color="batch")
+
+
+print("Store outputs", flush=True)
+out_adata.write_h5ad(par["output"], compression="gzip")
diff --git a/src/workflows/run_benchmark/config.vsh.yaml b/src/workflows/run_benchmark/config.vsh.yaml
index a4d6dbae..3a472457 100644
--- a/src/workflows/run_benchmark/config.vsh.yaml
+++ b/src/workflows/run_benchmark/config.vsh.yaml
@@ -109,6 +109,7 @@ dependencies:
   - name: methods/rpca_to_goal
   - name: methods/rpca_to_mid
   - name: methods/cytovi
+  - name: methods/mnnpy
   - name: metrics/emd
   - name: metrics/n_inconsistent_peaks
   - name: metrics/average_batch_r2
diff --git a/src/workflows/run_benchmark/main.nf b/src/workflows/run_benchmark/main.nf
index 673add7e..3985281a 100644
--- a/src/workflows/run_benchmark/main.nf
+++ b/src/workflows/run_benchmark/main.nf
@@ -35,7 +35,8 @@ methods = [
   cytonorm_one_control_to_goal,
   rpca_to_goal,
   rpca_to_mid,
-  cytovi
+  cytovi,
+  mnnpy
 ]
 
 // construct list of metrics

From 0431d13a7e0cdf6cbf4fe06982927e7fd0b05dc3 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Wed, 24 Sep 2025 00:29:45 +1000
Subject: [PATCH 32/46] try cytovi again

---
 src/methods/cytovi/config.vsh.yaml | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/methods/cytovi/config.vsh.yaml b/src/methods/cytovi/config.vsh.yaml
index ec6fa86c..1b4cccbd 100644
--- a/src/methods/cytovi/config.vsh.yaml
+++ b/src/methods/cytovi/config.vsh.yaml
@@ -58,8 +58,8 @@ engines:
           - git
       - type: python
         packages:
-          - anndata~=0.11.0
-          - scanpy~=1.11.0
+          - anndata~=0.12.0
+          - scanpy~=1.11.4
           - pyyaml
           - requests
           - jsonschema

From 97ebe0cf643415b42746f765e397a6f32d556712 Mon Sep 17 00:00:00 2001
From: Luqui12 <luca.leomazzi@gmail.com>
Date: Tue, 23 Sep 2025 17:53:45 +0200
Subject: [PATCH 33/46] Implemented lisi

---
 src/metrics/lisi/config.vsh.yaml | 130 +++++++++++++++++++++++++++++++
 src/metrics/lisi/script.py       |  78 +++++++++++++++++++
 2 files changed, 208 insertions(+)
 create mode 100644 src/metrics/lisi/config.vsh.yaml
 create mode 100644 src/metrics/lisi/script.py

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')

From fdc33df9eb5ef372b9612199daf177d621f8ba82 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Wed, 24 Sep 2025 10:20:34 +1000
Subject: [PATCH 34/46] test cytovi setup again

---
 src/methods/cytovi/config.vsh.yaml | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/methods/cytovi/config.vsh.yaml b/src/methods/cytovi/config.vsh.yaml
index 1b4cccbd..fcc00f9b 100644
--- a/src/methods/cytovi/config.vsh.yaml
+++ b/src/methods/cytovi/config.vsh.yaml
@@ -58,12 +58,12 @@ engines:
           - git
       - type: python
         packages:
-          - anndata~=0.12.0
-          - scanpy~=1.11.4
+          - anndata~=0.11.0
+          - scanpy~=1.11.0
+          - scvi-tools~=1.4.0
           - pyyaml
           - requests
           - jsonschema
-          - scvi-tools>=1.4.0
         github:
           - openproblems-bio/core#subdirectory=packages/python/openproblems
 
@@ -73,4 +73,4 @@ runners:
   # Allows turning the component into a Nextflow module / pipeline.
   - type: nextflow
     directives:
-      label: [hightime, highmem, midcpu, gpu]
+      label: [hightime, highmem, lowcpu, gpu]

From 7c2b1eaae738ffdd60058cc27181190475dd0a39 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Wed, 24 Sep 2025 10:42:04 +1000
Subject: [PATCH 35/46] testing different dependencies

---
 src/methods/cytovi/config.vsh.yaml | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/methods/cytovi/config.vsh.yaml b/src/methods/cytovi/config.vsh.yaml
index fcc00f9b..fe176ec9 100644
--- a/src/methods/cytovi/config.vsh.yaml
+++ b/src/methods/cytovi/config.vsh.yaml
@@ -58,9 +58,9 @@ engines:
           - git
       - type: python
         packages:
-          - anndata~=0.11.0
-          - scanpy~=1.11.0
-          - scvi-tools~=1.4.0
+          - anndata>=0.11.0
+          - scanpy[skmisc]>=1.10
+          - scvi-tools==1.4.0
           - pyyaml
           - requests
           - jsonschema

From 408f24515b373ff237989e3f683f8541ea5ede92 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Wed, 24 Sep 2025 11:03:31 +1000
Subject: [PATCH 36/46] remove package versioining for now

---
 src/methods/cytovi/config.vsh.yaml | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/methods/cytovi/config.vsh.yaml b/src/methods/cytovi/config.vsh.yaml
index fe176ec9..47e4b4bb 100644
--- a/src/methods/cytovi/config.vsh.yaml
+++ b/src/methods/cytovi/config.vsh.yaml
@@ -58,9 +58,9 @@ engines:
           - git
       - type: python
         packages:
-          - anndata>=0.11.0
-          - scanpy[skmisc]>=1.10
-          - scvi-tools==1.4.0
+          - anndata
+          - scanpy
+          - scvi-tools
           - pyyaml
           - requests
           - jsonschema

From 4547726de218f7329d1297a40bab3c5069ef3dfe Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Wed, 24 Sep 2025 11:26:03 +1000
Subject: [PATCH 37/46] revert to cpu fior cytovi and update mnnpy settings

---
 src/methods/cytovi/config.vsh.yaml | 15 ++-------------
 src/methods/cytovi/script.py       |  2 ++
 src/methods/mnnpy/config.vsh.yaml  |  2 +-
 3 files changed, 5 insertions(+), 14 deletions(-)

diff --git a/src/methods/cytovi/config.vsh.yaml b/src/methods/cytovi/config.vsh.yaml
index 47e4b4bb..43289e2f 100644
--- a/src/methods/cytovi/config.vsh.yaml
+++ b/src/methods/cytovi/config.vsh.yaml
@@ -50,22 +50,11 @@ resources:
 
 engines:
   - type: docker
-    image: nvcr.io/nvidia/pytorch:25.03-py3
+    image: openproblems/base_python:1
     setup:
-      - type: apt
-        packages:
-          - procps
-          - git
       - type: python
         packages:
-          - anndata
-          - scanpy
           - scvi-tools
-          - pyyaml
-          - requests
-          - jsonschema
-        github:
-          - openproblems-bio/core#subdirectory=packages/python/openproblems
 
 runners:
   # This platform allows running the component natively
@@ -73,4 +62,4 @@ runners:
   # Allows turning the component into a Nextflow module / pipeline.
   - type: nextflow
     directives:
-      label: [hightime, highmem, lowcpu, gpu]
+      label: [veryhightime, midmem, midcpu]
diff --git a/src/methods/cytovi/script.py b/src/methods/cytovi/script.py
index 89ad9028..6290f4e0 100644
--- a/src/methods/cytovi/script.py
+++ b/src/methods/cytovi/script.py
@@ -22,6 +22,8 @@
 
 adata_to_correct = adata[:, markers_to_correct].copy()
 
+print("Scaling data", flush=True)
+
 # scale data
 cytovi.scale(
     adata=adata_to_correct, transformed_layer_key="preprocessed", batch_key="batch_str"
diff --git a/src/methods/mnnpy/config.vsh.yaml b/src/methods/mnnpy/config.vsh.yaml
index acbfe6b9..646cae4c 100644
--- a/src/methods/mnnpy/config.vsh.yaml
+++ b/src/methods/mnnpy/config.vsh.yaml
@@ -76,4 +76,4 @@ runners:
   - type: executable
   - type: nextflow
     directives:
-      label: [hightime, midcpu, highmem]
\ No newline at end of file
+      label: [veryhightime, midcpu, midmem]
\ No newline at end of file

From e18afe3c26ba25dd6e80dff13d21f1e3b5abe216 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Wed, 24 Sep 2025 12:03:58 +1000
Subject: [PATCH 38/46] testing gpu again

---
 src/methods/cytovi/config.vsh.yaml | 25 ++++++++++++++++++++++---
 1 file changed, 22 insertions(+), 3 deletions(-)

diff --git a/src/methods/cytovi/config.vsh.yaml b/src/methods/cytovi/config.vsh.yaml
index 43289e2f..723bbad6 100644
--- a/src/methods/cytovi/config.vsh.yaml
+++ b/src/methods/cytovi/config.vsh.yaml
@@ -50,11 +50,30 @@ resources:
 
 engines:
   - type: docker
-    image: openproblems/base_python:1
+    image: nvcr.io/nvidia/pytorch:25.08-py3
     setup:
+      - type: apt
+        packages:
+          - procps
+          - git
       - type: python
         packages:
-          - scvi-tools
+          - anndata>=0.11.0
+          - scanpy[skmisc]>=1.10
+          - scvi-tools==1.4.0
+          - pyyaml
+          - requests
+          - jsonschema
+        github:
+          - openproblems-bio/core#subdirectory=packages/python/openproblems
+
+# engines:
+#   - type: docker
+#     image: openproblems/base_python:1
+#     setup:
+#       - type: python
+#         packages:
+#           - scvi-tools
 
 runners:
   # This platform allows running the component natively
@@ -62,4 +81,4 @@ runners:
   # Allows turning the component into a Nextflow module / pipeline.
   - type: nextflow
     directives:
-      label: [veryhightime, midmem, midcpu]
+      label: [veryhightime, midmem, midcpu, gpu]

From d5fcc03739bf79dc9878a0c3519f296726fe5162 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Wed, 24 Sep 2025 12:44:47 +1000
Subject: [PATCH 39/46] increase memory allocation for mnnpy

---
 src/methods/mnnpy/config.vsh.yaml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/methods/mnnpy/config.vsh.yaml b/src/methods/mnnpy/config.vsh.yaml
index 646cae4c..cb36e57b 100644
--- a/src/methods/mnnpy/config.vsh.yaml
+++ b/src/methods/mnnpy/config.vsh.yaml
@@ -76,4 +76,4 @@ runners:
   - type: executable
   - type: nextflow
     directives:
-      label: [veryhightime, midcpu, midmem]
\ No newline at end of file
+      label: [veryhightime, midcpu, highmem]
\ No newline at end of file

From 7c94058b90cd9a8d4e636fc686a1fc3a43a63895 Mon Sep 17 00:00:00 2001
From: Luqui12 <luca.leomazzi@gmail.com>
Date: Wed, 24 Sep 2025 11:40:38 +0200
Subject: [PATCH 40/46] Removed CMS, changed cms with lisi in
 `workflows/run_benchmark/` files

---
 src/metrics/cms/config.vsh.yaml             | 120 -----------------
 src/metrics/cms/script.R                    | 136 --------------------
 src/workflows/run_benchmark/config.vsh.yaml |   2 +-
 src/workflows/run_benchmark/main.nf         |   2 +-
 4 files changed, 2 insertions(+), 258 deletions(-)
 delete mode 100644 src/metrics/cms/config.vsh.yaml
 delete mode 100644 src/metrics/cms/script.R

diff --git a/src/metrics/cms/config.vsh.yaml b/src/metrics/cms/config.vsh.yaml
deleted file mode 100644
index 1cec62dc..00000000
--- a/src/metrics/cms/config.vsh.yaml
+++ /dev/null
@@ -1,120 +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.
-
-        In this implementation we subset to 60% of the total cells in each split of the technical replicates. The subset is stratified by cell type and sample.
-        The mean medcouple between the two splits of technical replicates is used as the final 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 b902e341..00000000
--- a/src/metrics/cms/script.R
+++ /dev/null
@@ -1,136 +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)
-}
-bpparam <- BiocParallel::MulticoreParam(
-    workers = cores_to_use
-)
-
-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 metadata from unintegrated\n")
-integrated_split1 <- get_obs_var_for_integrated(
-    i_adata = integrated_split1,
-    u_adata = unintegrated,
-    split_id = 1
-)
-integrated_split2 <- get_obs_var_for_integrated(
-    i_adata = integrated_split2,
-    u_adata = unintegrated,
-    split_id = 2
-)
-
-# Get markers to correct
-markers_to_correct <- unintegrated$var_names[unintegrated$var$to_correct]
-
-cat(paste("Compute Cell Mixing Score using", cores_to_use, "cores for split 1\n"))
-
-cms_distr_split1 <- list()
-medcouples_split1 <- list()
-
-integrated_subset <- subset_by_celltype(
-    integrated_split1,
-    frac = 0.6,
-    seed = 1
-)
-
-cat("Transforming to SingleCellExperiment and subsetting markers\n")
-integrated_subset_sce <- integrated_subset$as_SingleCellExperiment()
-integrated_subset_sce <- integrated_subset_sce[markers_to_correct, ]
-cat("Computing Cell Mixing Scores\n")
-integrated_subset_sce <- CellMixS::cms(
-    integrated_subset_sce,
-    group = "batch",
-    assay_name = "integrated",
-    k = par[["n_neighbors"]],
-    n_dim = par[["n_dim"]],
-    BPPARAM = bpparam
-)
-distr <- SingleCellExperiment::colData(integrated_subset_sce)[, "cms"]
-cms_distr_split1[["split1"]] <- distr
-medcouples_split1[["split1"]] <- robustbase::mc(distr)
-
-cat(paste("Compute Cell Mixing Score using", cores_to_use, "cores for split 2\n"), flush = TRUE)
-
-cms_distr_split2 <- list()
-medcouples_split2 <- list()
-
-integrated_subset <- subset_by_celltype(
-    integrated_split2,
-    frac = 0.6,
-    seed = 1
-)
-cat("Transforming to SingleCellExperiment and subsetting markers\n")
-integrated_subset_sce <- integrated_subset$as_SingleCellExperiment()
-integrated_subset_sce <- integrated_subset_sce[markers_to_correct, ]
-cat("Computing Cell Mixing Scores\n")
-integrated_subset_sce <- CellMixS::cms(
-    integrated_subset_sce,
-    group = "batch",
-    assay_name = "integrated",
-    k = par[["n_neighbors"]],
-    n_dim = par[["n_dim"]],
-    BPPARAM = bpparam
-)
-distr <- SingleCellExperiment::colData(integrated_subset_sce)[, "cms"]
-cms_distr_split2[["split2"]] <- distr
-medcouples_split2[["split2"]] <- robustbase::mc(distr)
-
-
-cat("Aggregate scores\n", flush = TRUE)
-#concat named lists
-cms_distr_list <- c(cms_distr_split1, cms_distr_split2)
-medcouples_list <- c(medcouples_split1, medcouples_split2)
-# Compute mean medcouple
-mean_medcouple_cms <- mean(unlist(medcouples_list))
-
-cat("Write output AnnData to file\n", flush = TRUE)
-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 = mean_medcouple_cms,
-        cms_parameters = list(
-            n_neighbors = par[["n_neighbors"]],
-            n_dim = par[["n_dim"]]
-        ),
-        list_medcouples = medcouples_list,
-        cms_distributions = cms_distr_list
-    )
-)
-
-output$write_h5ad(par[["output"]], compression = "gzip", mode = "w")
diff --git a/src/workflows/run_benchmark/config.vsh.yaml b/src/workflows/run_benchmark/config.vsh.yaml
index 3a472457..6cc67865 100644
--- a/src/workflows/run_benchmark/config.vsh.yaml
+++ b/src/workflows/run_benchmark/config.vsh.yaml
@@ -114,7 +114,7 @@ dependencies:
   - 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 3985281a..3305e346 100644
--- a/src/workflows/run_benchmark/main.nf
+++ b/src/workflows/run_benchmark/main.nf
@@ -45,7 +45,7 @@ metrics = [
   n_inconsistent_peaks,
   average_batch_r2,
   flowsom_mapping_similarity,
-  cms,
+  lisi,
   bras
 ]
 

From 1387ec858867f478c5c5afcb7552ed79fb26cc2b Mon Sep 17 00:00:00 2001
From: Luqui12 <luca.leomazzi@gmail.com>
Date: Thu, 25 Sep 2025 17:13:51 +0200
Subject: [PATCH 41/46] Removed mnn.

---
 src/methods/mnn/config.vsh.yaml             | 107 --------------------
 src/methods/mnn/script.R                    |  57 -----------
 src/workflows/run_benchmark/config.vsh.yaml |   1 -
 src/workflows/run_benchmark/main.nf         |   1 -
 4 files changed, 166 deletions(-)
 delete mode 100644 src/methods/mnn/config.vsh.yaml
 delete mode 100644 src/methods/mnn/script.R

diff --git a/src/methods/mnn/config.vsh.yaml b/src/methods/mnn/config.vsh.yaml
deleted file mode 100644
index bb4cfb61..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: [veryhightime,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/workflows/run_benchmark/config.vsh.yaml b/src/workflows/run_benchmark/config.vsh.yaml
index 6cc67865..fc5f863b 100644
--- a/src/workflows/run_benchmark/config.vsh.yaml
+++ b/src/workflows/run_benchmark/config.vsh.yaml
@@ -103,7 +103,6 @@ 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
diff --git a/src/workflows/run_benchmark/main.nf b/src/workflows/run_benchmark/main.nf
index 3305e346..4c2e5234 100644
--- a/src/workflows/run_benchmark/main.nf
+++ b/src/workflows/run_benchmark/main.nf
@@ -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,

From 6810ed0898246438cd4e29c31688df948362e251 Mon Sep 17 00:00:00 2001
From: Luqui12 <luca.leomazzi@gmail.com>
Date: Mon, 29 Sep 2025 15:59:06 +0200
Subject: [PATCH 42/46] Re-implemented control methods
 `shuffle_integration_by_batch` and `shuffle_integration_by_cell_type`

---
 .../config.vsh.yaml                           |  2 +-
 .../shuffle_integration_by_batch/script.py    | 39 ++++++++++++++----
 .../script.py                                 | 41 ++++++++++++++-----
 src/workflows/run_benchmark/config.vsh.yaml   |  4 +-
 src/workflows/run_benchmark/main.nf           |  4 +-
 5 files changed, 65 insertions(+), 25 deletions(-)

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/workflows/run_benchmark/config.vsh.yaml b/src/workflows/run_benchmark/config.vsh.yaml
index fc5f863b..6ed316b0 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
diff --git a/src/workflows/run_benchmark/main.nf b/src/workflows/run_benchmark/main.nf
index 4c2e5234..b14cf770 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,

From 381a140091acb1e1285cf34383d6321d78a8cea9 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Tue, 30 Sep 2025 22:04:27 +1000
Subject: [PATCH 43/46] test cytovi

---
 src/methods/cytovi/config.vsh.yaml | 18 +++++-----
 src/methods/cytovi/script.py       | 56 ++++++++++++++++++++++++------
 2 files changed, 55 insertions(+), 19 deletions(-)

diff --git a/src/methods/cytovi/config.vsh.yaml b/src/methods/cytovi/config.vsh.yaml
index 723bbad6..91095326 100644
--- a/src/methods/cytovi/config.vsh.yaml
+++ b/src/methods/cytovi/config.vsh.yaml
@@ -38,6 +38,14 @@ arguments:
     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:
@@ -67,18 +75,10 @@ engines:
         github:
           - openproblems-bio/core#subdirectory=packages/python/openproblems
 
-# engines:
-#   - type: docker
-#     image: openproblems/base_python:1
-#     setup:
-#       - type: python
-#         packages:
-#           - scvi-tools
-
 runners:
   # This platform allows running the component natively
   - type: executable
   # Allows turning the component into a Nextflow module / pipeline.
   - type: nextflow
     directives:
-      label: [veryhightime, midmem, midcpu, gpu]
+      label: [veryhightime, lowmem, lowcpu, gpu]
diff --git a/src/methods/cytovi/script.py b/src/methods/cytovi/script.py
index 6290f4e0..6a6fd73f 100644
--- a/src/methods/cytovi/script.py
+++ b/src/methods/cytovi/script.py
@@ -1,6 +1,8 @@
 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 = {
@@ -8,6 +10,8 @@
     "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
@@ -24,26 +28,54 @@
 
 print("Scaling data", flush=True)
 
-# scale data
+# 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"
+    adata=adata_to_correct,
+    transformed_layer_key="preprocessed",
+    batch_key="batch_str",
+    inplace=True,
 )
 
-print("Run CytoVI", flush=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_to_correct, layer="scaled", batch_key="batch_str")
+cytovi.CYTOVI.setup_anndata(adata_subsampled, layer="scaled", batch_key="batch_str")
 model = cytovi.CYTOVI(
-    adata=adata_to_correct, n_hidden=par["n_hidden"], n_layers=par["n_layers"]
+    adata=adata_subsampled, n_hidden=par["n_hidden"], n_layers=par["n_layers"]
 )
 model.train()
 
 # get batch corrected data
-corrected_data = model.get_normalized_expression()
+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.to_numpy(), uncorrected_data], axis=1)
+out_matrix = np.concatenate([corrected_data, uncorrected_data], axis=1)
 out_var_idx = np.concatenate([corrected_data.columns, markers_not_correct])
 
 # create new anndata
@@ -61,12 +93,16 @@
 # 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 = out_adata.copy()
+
+# test_adata = ad.AnnData(
+#     X=out_adata.layers["integrated"].toarray(),
+#     obs=adata.obs,
+#     var=adata.var,
+# )
 # test_adata = test_adata[:, markers_to_correct]
-# test_adata.X = test_adata.layers["integrated"]
-# test_adata.obs = adata.obs
 # sc.pp.neighbors(test_adata, use_rep="X")
 # sc.tl.umap(test_adata)
 # sc.pl.umap(test_adata, color="batch")

From 6b1d40a91b004c1dcd05df1b133e3c2b68aa767b Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Wed, 1 Oct 2025 00:28:50 +1000
Subject: [PATCH 44/46] removed global emd

---
 scripts/run_benchmark/run_full_seqeracloud.sh |   6 +-
 src/metrics/emd/config.vsh.yaml               | 172 ------------------
 src/metrics/emd/helper.py                     |  70 +------
 src/metrics/emd/script.py                     |   8 -
 4 files changed, 7 insertions(+), 249 deletions(-)

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/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 417fbe65..a8a3253e 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
@@ -51,8 +45,6 @@ def calculate_vertical_emd(
     )
 
     # safeguard
-    mean_emd_global = np.nan
-    max_emd_global = np.nan
     mean_emd_ct = np.nan
     max_emd_ct = np.nan
 
@@ -65,20 +57,6 @@ def calculate_vertical_emd(
     if len(emd_long) > 0:
         emd_long = pd.concat(emd_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"]
@@ -94,8 +72,6 @@ def calculate_vertical_emd(
         )
 
     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,
@@ -150,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]
@@ -239,10 +204,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
@@ -252,8 +213,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]
@@ -305,19 +264,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(
@@ -327,24 +274,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,
     }
 
 
diff --git a/src/metrics/emd/script.py b/src/metrics/emd/script.py
index 80dd447b..6586070d 100644
--- a/src/metrics/emd/script.py
+++ b/src/metrics/emd/script.py
@@ -94,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],
         ],

From 98f26a943d84df817553199478a35495ea8974b8 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Wed, 1 Oct 2025 01:27:54 +1000
Subject: [PATCH 45/46] simplify emd vertical output

---
 src/metrics/emd/helper.py | 52 ++++++++++++++++++++-------------------
 1 file changed, 27 insertions(+), 25 deletions(-)

diff --git a/src/metrics/emd/helper.py b/src/metrics/emd/helper.py
index a8a3253e..c2f78a2c 100644
--- a/src/metrics/emd/helper.py
+++ b/src/metrics/emd/helper.py
@@ -36,11 +36,11 @@ 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
     )
 
@@ -74,8 +74,8 @@ def calculate_vertical_emd(
     return {
         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,
     }
 
 
@@ -153,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
@@ -160,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]
+
+    #         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, emd_wide
+    return emd_vals
 
 
 def calculate_horizontal_emd(

From 62d88631836d0f01b05c3049d5f62eec32d7f2ed Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Wed, 1 Oct 2025 22:08:39 +1000
Subject: [PATCH 46/46] remove mnnpy

---
 src/methods/mnnpy/config.vsh.yaml           | 79 ---------------------
 src/methods/mnnpy/script.py                 | 71 ------------------
 src/workflows/run_benchmark/config.vsh.yaml |  1 -
 src/workflows/run_benchmark/main.nf         |  3 +-
 4 files changed, 1 insertion(+), 153 deletions(-)
 delete mode 100644 src/methods/mnnpy/config.vsh.yaml
 delete mode 100644 src/methods/mnnpy/script.py

diff --git a/src/methods/mnnpy/config.vsh.yaml b/src/methods/mnnpy/config.vsh.yaml
deleted file mode 100644
index cb36e57b..00000000
--- a/src/methods/mnnpy/config.vsh.yaml
+++ /dev/null
@@ -1,79 +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: mnnpy
-# A relatively short label, used when rendering visualisations (required)
-label: mnnpy
-# A one sentence summary of how this method works (required). Used when 
-# rendering summary tables.
-summary: "Batch effect correction by matching mutual nearest neighbors, Python implementation."
-# A multi-line description of how this component works (required). Used
-# when rendering reference documentation.
-description: |
-  An implementation of MNN correct in python featuring low memory usage, full multicore support and compatibility with the scanpy framework.
-
-  Batch effect correction by matching mutual nearest neighbors (Haghverdi et al, 2018) has been implemented as a function 'mnnCorrect' in the R package scran. Sadly it's extremely slow for big datasets and doesn't make full use of the parallel architecture of modern CPUs.
-
-  This project is a python implementation of the MNN correct algorithm which takes advantage of python's extendability and hackability. It seamlessly integrates with the scanpy framework and has multicore support in its bones.
-references:
-  bibtex: |
-    @misc{Kang2022,
-      author = {Kang, Chris},
-      title = {mnnpy},
-      year = {Kang2022},
-      publisher = {GitHub},
-      journal = {GitHub repository},
-      howpublished = {\url{https://github.com/chriscainx/mnnpy}},
-      commit = {2097dec30c193f036c5ed7e1c3d1e3a6270e102b}
-    }
-links:
-  repository: https://github.com/chriscainx/mnnpy
-  documentation: https://github.com/chriscainx/mnnpy#readme
-
-# 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
-
-engines:
-  - type: docker
-    image: python:3.8
-    setup:
-      - type: apt
-        packages:
-          - procps
-          - libhdf5-dev
-      - type: python
-        pypi:
-          - h5py==3.11.0
-          - anndata~=0.8.0
-          - scanpy
-          - pyyaml
-          - requests
-          - jsonschema
-      - type: python
-        pypi:
-          - git+https://github.com/openproblems-bio/core#subdirectory=packages/python/openproblems
-          - --ignore-requires-python
-      - type: python
-        github:
-          - chriscainx/mnnpy
-runners:
-  - type: executable
-  - type: nextflow
-    directives:
-      label: [veryhightime, midcpu, highmem]
\ No newline at end of file
diff --git a/src/methods/mnnpy/script.py b/src/methods/mnnpy/script.py
deleted file mode 100644
index a059f464..00000000
--- a/src/methods/mnnpy/script.py
+++ /dev/null
@@ -1,71 +0,0 @@
-import anndata as ad
-import mnnpy
-import numpy as np
-
-## VIASH START
-par = {
-    "input": "resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/censored_split1.h5ad",
-    "output": "resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/output_mnnpy.h5ad",
-}
-meta = {"name": "mnnpy"}
-## VIASH END
-
-print("Read input", flush=True)
-adata = ad.read_h5ad(par["input"])
-
-adata.X = adata.layers["preprocessed"]
-
-# convert batch to category as otherwise mnnpy won't work..
-adata.obs["batch_cat"] = 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("Run mnn", flush=True)
-split = []
-batch_categories = adata_to_correct.obs["batch_cat"].unique().tolist()
-
-for i in batch_categories:
-    split.append(adata_to_correct[adata_to_correct.obs["batch_cat"] == i].copy())
-
-corrected, _, _ = mnnpy.mnn_correct(
-    *split, batch_key="batch", batch_categories=batch_categories, index_unique=None
-)
-
-# have to add in the uncorrected markers as well
-uncorrected_data = adata[:, markers_not_correct].layers["preprocessed"]
-
-out_matrix = np.concatenate([corrected.X, uncorrected_data], axis=1)
-out_var_idx = np.concatenate([corrected.var.index, 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]
-
-# run umap for quick check
-# import scanpy as sc
-# test_adata = out_adata.copy()
-# test_adata = test_adata[:, markers_to_correct]
-# test_adata.X = test_adata.layers["integrated"]
-# test_adata.obs = adata.obs
-# sc.pp.neighbors(test_adata, use_rep="X")
-# sc.tl.umap(test_adata)
-# sc.pl.umap(test_adata, color="batch")
-
-
-print("Store outputs", flush=True)
-out_adata.write_h5ad(par["output"], compression="gzip")
diff --git a/src/workflows/run_benchmark/config.vsh.yaml b/src/workflows/run_benchmark/config.vsh.yaml
index 6ed316b0..b63c6e3e 100644
--- a/src/workflows/run_benchmark/config.vsh.yaml
+++ b/src/workflows/run_benchmark/config.vsh.yaml
@@ -108,7 +108,6 @@ dependencies:
   - name: methods/rpca_to_goal
   - name: methods/rpca_to_mid
   - name: methods/cytovi
-  - name: methods/mnnpy
   - name: metrics/emd
   - name: metrics/n_inconsistent_peaks
   - name: metrics/average_batch_r2
diff --git a/src/workflows/run_benchmark/main.nf b/src/workflows/run_benchmark/main.nf
index b14cf770..26546f04 100644
--- a/src/workflows/run_benchmark/main.nf
+++ b/src/workflows/run_benchmark/main.nf
@@ -34,8 +34,7 @@ methods = [
   cytonorm_one_control_to_goal,
   rpca_to_goal,
   rpca_to_mid,
-  cytovi,
-  mnnpy
+  cytovi
 ]
 
 // construct list of metrics