From 65a35a900ec83b9191b1fef03f0a04835aab147f Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Sun, 12 Oct 2025 11:05:46 +1100
Subject: [PATCH 01/31] first prototype of new implementation

---
 .../ratio_inconsistent_peaks/config.vsh.yaml  |  66 +++++
 .../ratio_inconsistent_peaks/helper.py        | 113 ++++++++
 .../ratio_inconsistent_peaks/script.py        | 263 ++++++++++++++++++
 3 files changed, 442 insertions(+)
 create mode 100644 src/metrics/ratio_inconsistent_peaks/config.vsh.yaml
 create mode 100644 src/metrics/ratio_inconsistent_peaks/helper.py
 create mode 100644 src/metrics/ratio_inconsistent_peaks/script.py

diff --git a/src/metrics/ratio_inconsistent_peaks/config.vsh.yaml b/src/metrics/ratio_inconsistent_peaks/config.vsh.yaml
new file mode 100644
index 00000000..c45d7c84
--- /dev/null
+++ b/src/metrics/ratio_inconsistent_peaks/config.vsh.yaml
@@ -0,0 +1,66 @@
+# 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: ratio_inconsistent_peaks
+
+# Metadata for your component
+info:
+  metrics:
+      # A unique identifier for your metric (required).
+      # Can contain only lowercase letters or underscores.
+    - name: ratio_inconsistent_peaks
+      label: Ratio of inconsistent peaks
+      summary: "Ratio of the number of cell‑type marker‑expression peaks between validation and batch‑normalized data."
+      description: |
+        The metric compares the number of cell type specific 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 (cell type) 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.
+      references:
+        doi: 
+          - 10.1038/s41592-019-0686-2
+      links:
+        # URL to the documentation for this metric (required).
+        documentation: https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.find_peaks.html#scipy.signal.find_peaks
+        # URL to the code repository for this metric (required).
+        repository: https://github.com/scipy/scipy/blob/v1.15.2/scipy/signal/_peak_finding.py#L0-L1
+      # 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
+
+# Resources required to run the component
+resources:
+  # The script of your component (required)
+  - type: python_script
+    path: script.py
+  - path: helper.py
+  - path: /src/utils/helper_functions.py
+
+engines:
+  # Specifications for the Docker image for this component.
+  - type: docker
+    image: openproblems/base_python:1
+    # Add custom dependencies here (optional). For more information, see
+    # https://viash.io/reference/config/engines/docker/#setup .
+    # setup:
+    #   - type: python
+    #     packages: numpy<2
+
+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/ratio_inconsistent_peaks/helper.py b/src/metrics/ratio_inconsistent_peaks/helper.py
new file mode 100644
index 00000000..427ee599
--- /dev/null
+++ b/src/metrics/ratio_inconsistent_peaks/helper.py
@@ -0,0 +1,113 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import seaborn as sns
+from ripser import ripser
+from scipy.signal import find_peaks
+from scipy.stats import gaussian_kde
+
+
+def standardise_marker_expression(dist_1, dist_2):
+    """
+    Standardises the marker expression values from two distributions.
+
+    Inputs:
+    dist_1: array of values (1D) representing the marker expression from distribution 1
+    dist_2: array of values (1D) representing the marker expression from distribution 2
+
+    Outputs:
+    std_dist_1: array of standardised values for distribution 1
+    std_dist_2: array of standardised values for distribution 2
+    """
+
+    pooled = np.concatenate([dist_1, dist_2])
+    mu, sd = pooled.mean(), pooled.std()
+    std_dist_1 = (dist_1 - mu) / (sd)
+    std_dist_2 = (dist_2 - mu) / (sd)
+
+    return std_dist_1, std_dist_2
+
+
+def get_kde_density(expression_array, return_xgrid=False, plot=False):
+    """
+    Returns the density of the array using a gaussian kernel density estimation.
+
+    Inputs:
+    expression_array: array of values (1D) representing the marker expression
+    return_xgrid: boolean, if True, also return the x_grid values used for density estimation
+    plot: boolean, if True, plot the density estimation
+
+    Outputs:
+    density: array of values representing the density of marker expression
+    x_grid (optional): array of x values where the density is evaluated
+    """
+
+    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="scott")
+    x_grid = np.linspace(min_val, max_val, 100)
+    density = kde(x_grid)
+
+    if plot:
+        fig, ax = plt.subplots()
+        sns.scatterplot(x=x_grid, y=density, ax=ax)
+        ax.set_title("KDE Density Estimation")
+        ax.set_xlabel("Marker Expression")
+        ax.set_ylabel("Density")
+        fig.tight_layout()
+        fig.show()
+
+    if return_xgrid:
+        # handy for plotting later on and maybe even save in the AnnData object
+        return density, x_grid
+    else:
+        return density
+
+
+def call_peaks(density):
+    """
+    Returns the peaks of the density using scipy.signal.find_peaks.
+
+    Inputs:
+    density: array of values representing the density of marker expression
+
+    Outputs:
+    peaks: array of values representing the peaks of the density
+    """
+
+    height_trsh = 0.1
+    prom_trsh = 0.01
+
+    peaks, _ = find_peaks(density, prominence=prom_trsh, height=height_trsh)
+    num_peaks = len(peaks)
+
+    return num_peaks
+
+
+def persistent_peak_count(ys, persistence_cutoff=0.05):
+    """
+    Counts robust peaks in a 1D dataset using persistent homology.
+
+    Args:
+        ys (np.ndarray): KDE of a marker expression (1D array)
+        persistence_cutoff (float): a threshold that decides which peaks are “significant enough” to count.
+            A large persistence peak survives over many levels of smoothing (i.e. a strong, real peak).
+            A small persistence peak quickly merges into a neighbor — likely noise.
+            0.01: very low threshold counts even weak bumps as peaks
+            0.05: moderate (default) counts clearly separated peaks
+            0.1–0.2: high threshold counts only strong, dominant peaks
+
+    Returns:
+        int: number of significant peaks
+        (diagram, persistence_values): raw persistence outputs
+    """
+
+    # Invert to turn peaks into "holes" for 0D persistence
+    Y = -ys.reshape(-1, 1)
+    diagram = ripser(Y, maxdim=0)["dgms"][0]
+    persistence = diagram[:, 1] - diagram[:, 0]
+
+    # Define significance threshold relative to data range
+    threshold = persistence_cutoff * np.ptp(ys)
+    n_peaks = np.sum(persistence > threshold)
+    return n_peaks, diagram, persistence
diff --git a/src/metrics/ratio_inconsistent_peaks/script.py b/src/metrics/ratio_inconsistent_peaks/script.py
new file mode 100644
index 00000000..5dd62766
--- /dev/null
+++ b/src/metrics/ratio_inconsistent_peaks/script.py
@@ -0,0 +1,263 @@
+import sys
+from collections import defaultdict
+
+import anndata as ad
+import numpy as np
+
+## VIASH START
+# 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",
+    "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": "ratio_inconsistent_peaks",
+}
+
+# for local testing only
+# import src.metrics.ratio_inconsistent_peaks.helper as metric_helper
+# from src.utils.helper_functions import (
+#     get_obs_var_for_integrated,
+#     remove_unlabelled,
+#     subset_markers_tocorrect,
+#     subset_nocontrols,
+# )
+
+## VIASH END
+
+sys.path.append(meta["resources_dir"])
+
+import helper as metric_helper
+
+# from helper import call_peaks, get_kde_density
+from helper_functions import (
+    get_obs_var_for_integrated,
+    remove_unlabelled,
+    subset_markers_tocorrect,
+    subset_nocontrols,
+)
+
+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)
+
+donor_list = integrated_s1.obs["donor"].unique()
+
+print("Compute metric (per cell type)", flush=True)
+
+# case 1 = consistent peaks in unintegrated and also in integrated
+# case 3 = consistent peaks in unintegrated but inconsistent in integrated
+# not recording case 2 or 4 where unintegrated is inconsistent
+n_case1 = 0
+n_case3 = 0
+
+# so we can see where each cases comes from
+case_details = defaultdict(list)
+
+
+for donor in donor_list:
+    # for testing only
+    # donor = donor_list[0]
+
+    print("Processing donor", donor, flush=True)
+
+    u_view = unintegrated[unintegrated.obs["donor"] == donor]
+
+    # process per split
+    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:
+        # for testing only
+        # celltype = celltype_list[0]
+
+        print(f"Processing celltype {celltype}", flush=True)
+
+        u_view_ct = u_view[u_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]
+
+        # safeguard checking the obsnames are exactly the same
+        print(
+            "Check that the obsnames are exactly the same between s1 and unintegrated",
+            flush=True,
+        )
+        # Check membership
+        split_names = np.concatenate([s1_view_ct.obs_names, s2_view_ct.obs_names])
+        u_names = list(u_view_ct.obs_names)
+
+        in_split_not_u = [n for n in split_names if n not in u_names]
+        in_u_not_split = [n for n in u_names if n not in split_names]
+
+        if len(in_split_not_u) > 0 or len(in_u_not_split) > 0:
+            print(
+                f"Error: the obsnames of the unintegrated and s1 and s2 do not match for {donor}, {celltype}."
+                "\nPlease check dataset as this should not happen!"
+                f"\nMissing (in s1/s2 but not in unintegrated): {in_split_not_u}."
+                f"\nExtra (in unintegrated but not in s1/s2): {in_u_not_split}.",
+                flush=True,
+            )
+            sys.exit(1)
+
+        print("Obsnames match!", flush=True)
+
+        if s1_view_ct.shape[0] < 100 or s2_view_ct.shape[0] < 100:
+            print(f"Skipping celltype {celltype} and donor {donor}.", flush=True)
+            if s1_view_ct.shape[0] < 100:
+                print(
+                    f"Because n_cells in s1 is {s1_view_ct.shape[0]}, less than 100",
+                    flush=True,
+                )
+            else:
+                print(
+                    f"Because n_cells in s2 is {s2_view_ct.shape[0]}, less than 100",
+                    flush=True,
+                )
+            # TODO uncomment me when done
+            continue
+
+        for marker in s1_view_ct.var.index:
+            # for testing only
+            # marker = u_view_ct.var.index[0]
+
+            print(f"Processing marker {marker} for celltype {celltype}", flush=True)
+
+            print("--------------------------------", flush=True)
+            print("Computing peaks for unintegrated", flush=True)
+            # unintegrated for split 1
+            u_view_ct_s1 = u_view_ct[u_view_ct.obs["split"] == 1]
+            u_view_ct_s2 = u_view_ct[u_view_ct.obs["split"] == 2]
+
+            print("Standardising marker expression", flush=True)
+            # standardise marker expression based on pooled mean and sd of
+            # unscaled marker expression for unintegrated data for split 1 and 2
+            u_s1_unscaled = np.array(u_view_ct_s1[:, marker].layers["preprocessed"])
+            u_s2_unscaled = np.array(u_view_ct_s2[:, marker].layers["preprocessed"])
+
+            u_s1_scaled, u_s2_scaled = metric_helper.standardise_marker_expression(
+                u_s1_unscaled,
+                u_s2_unscaled,
+            )
+            print("Computing KDE density", flush=True)
+            density_dist_u_s1 = metric_helper.get_kde_density(u_s1_scaled)
+            density_dist_u_s2 = metric_helper.get_kde_density(u_s2_scaled)
+
+            print("Calling peaks", flush=True)
+
+            peaks_u_s1 = metric_helper.call_peaks(density_dist_u_s1)
+            peaks_u_s2 = metric_helper.call_peaks(density_dist_u_s2)
+
+            # use persistent peak only if the peak calling method is too sensitive...
+            # persistent_peak_count_u_s1, diagram_u_s1, persistence_values_u_s1 = (
+            #     metric_helper.persistent_peak_count(density_dist_u_s1)
+            # )
+            # persistent_peak_count_u_s2, diagram_u_s2, persistence_values_u_s2 = (
+            #     metric_helper.persistent_peak_count(density_dist_u_s2)
+            # )
+
+            print("--------------------------------", flush=True)
+
+            print("\n", flush=True)
+
+            print("--------------------------------", flush=True)
+            print("Computing peaks for integrated", flush=True)
+
+            print("Standardising marker expression", flush=True)
+            # standardise marker expression based on pooled mean and sd of
+            # unscaled marker expression for unintegrated data for split 1 and 2
+            s1_unscaled = np.array(s1_view_ct[:, marker].layers["integrated"])
+            s2_unscaled = np.array(s2_view_ct[:, marker].layers["integrated"])
+
+            s1_scaled, s2_scaled = metric_helper.standardise_marker_expression(
+                s1_unscaled,
+                s2_unscaled,
+            )
+            print("Computing KDE density", flush=True)
+            density_dist_s1 = metric_helper.get_kde_density(s1_scaled)
+            density_dist_s2 = metric_helper.get_kde_density(s2_scaled)
+
+            print("Calling peaks", flush=True)
+
+            peaks_s1 = metric_helper.call_peaks(density_dist_s1)
+            peaks_s2 = metric_helper.call_peaks(density_dist_s2)
+
+            # use persistent peak only if the peak calling method is too sensitive...
+            # persistent_peak_count_s1, diagram_s1, persistence_values_s1 = (
+            #     metric_helper.persistent_peak_count(density_dist_s1)
+            # )
+            # persistent_peak_count_s2, diagram_s2, persistence_values_s2 = (
+            #     metric_helper.persistent_peak_count(density_dist_s2)
+            # )
+
+            print("--------------------------------", flush=True)
+
+            print("\n", flush=True)
+
+            print("Comparing peaks between unintegrated and integrated", flush=True)
+
+            # case 1 or 3 where we have consistent peaks in unintegrated
+            if peaks_u_s1 == peaks_u_s2:
+                if peaks_s1 != peaks_s2:
+                    n_case3 += 1
+                    case_details["case3"].append((donor, celltype, marker))
+                else:
+                    n_case1 += 1
+                    case_details["case1"].append((donor, celltype, marker))
+            else:
+                print(
+                    "WARNING! Inconsistent peaks detected in unintegrated data (case 2 or 4). Skipping calculation",
+                    flush=True,
+                )
+                case_details["case2or4"].append((donor, celltype, marker))
+
+print("Done processing all celltypes and donors", flush=True)
+print("Calculating ratio", flush=True)
+
+if n_case1 + n_case3 == 0:
+    print(
+        "Only case 2 or 4 are found!. Cannot calculate metric.",
+        flush=True,
+    )
+    metric_val = np.nan
+else:
+    metric_val = n_case3 / (n_case1 + n_case3)
+
+
+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": [meta["name"]],
+        "metric_values": [metric_val],
+        "n_cases": {
+            "case1": n_case1,
+            "case3": n_case3,
+            "case2or4": len(case_details["case2or4"]),
+        },
+        "case_details": dict(case_details),
+    }
+)
+output.write_h5ad(par["output"], compression="gzip")
+
+# print(uns_metric_ids, uns_metric_values)

From 3996bfa8389bb0e50fd468bc572d674ebd932411 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Sun, 12 Oct 2025 11:45:54 +1100
Subject: [PATCH 02/31] add scikit tda and persistent peak as alternative
 output

---
 .../ratio_inconsistent_peaks/config.vsh.yaml  |  9 +-
 .../ratio_inconsistent_peaks/helper.py        |  6 +-
 .../ratio_inconsistent_peaks/script.py        | 97 ++++++++++++-------
 3 files changed, 67 insertions(+), 45 deletions(-)

diff --git a/src/metrics/ratio_inconsistent_peaks/config.vsh.yaml b/src/metrics/ratio_inconsistent_peaks/config.vsh.yaml
index c45d7c84..bfcc00ee 100644
--- a/src/metrics/ratio_inconsistent_peaks/config.vsh.yaml
+++ b/src/metrics/ratio_inconsistent_peaks/config.vsh.yaml
@@ -51,11 +51,10 @@ engines:
   # Specifications for the Docker image for this component.
   - type: docker
     image: openproblems/base_python:1
-    # Add custom dependencies here (optional). For more information, see
-    # https://viash.io/reference/config/engines/docker/#setup .
-    # setup:
-    #   - type: python
-    #     packages: numpy<2
+    setup:
+      - type: python
+        packages:
+          - scikit-tda
 
 runners:
   # This platform allows running the component natively
diff --git a/src/metrics/ratio_inconsistent_peaks/helper.py b/src/metrics/ratio_inconsistent_peaks/helper.py
index 427ee599..008eb211 100644
--- a/src/metrics/ratio_inconsistent_peaks/helper.py
+++ b/src/metrics/ratio_inconsistent_peaks/helper.py
@@ -84,7 +84,7 @@ def call_peaks(density):
     return num_peaks
 
 
-def persistent_peak_count(ys, persistence_cutoff=0.05):
+def persistent_peak_count(ys, persistence_cutoff=0.08):
     """
     Counts robust peaks in a 1D dataset using persistent homology.
 
@@ -96,10 +96,10 @@ def persistent_peak_count(ys, persistence_cutoff=0.05):
             0.01: very low threshold counts even weak bumps as peaks
             0.05: moderate (default) counts clearly separated peaks
             0.1–0.2: high threshold counts only strong, dominant peaks
+            Default to 0.08 to biased towards strong peaks but not overly.
 
     Returns:
         int: number of significant peaks
-        (diagram, persistence_values): raw persistence outputs
     """
 
     # Invert to turn peaks into "holes" for 0D persistence
@@ -110,4 +110,4 @@ def persistent_peak_count(ys, persistence_cutoff=0.05):
     # Define significance threshold relative to data range
     threshold = persistence_cutoff * np.ptp(ys)
     n_peaks = np.sum(persistence > threshold)
-    return n_peaks, diagram, persistence
+    return n_peaks
diff --git a/src/metrics/ratio_inconsistent_peaks/script.py b/src/metrics/ratio_inconsistent_peaks/script.py
index 5dd62766..9bc26657 100644
--- a/src/metrics/ratio_inconsistent_peaks/script.py
+++ b/src/metrics/ratio_inconsistent_peaks/script.py
@@ -3,6 +3,7 @@
 
 import anndata as ad
 import numpy as np
+import pandas as pd
 
 ## VIASH START
 # The following code has been auto-generated by Viash.
@@ -72,6 +73,9 @@
 # so we can see where each cases comes from
 case_details = defaultdict(list)
 
+# for comparison only
+persistent_peaks_res = []
+
 
 for donor in donor_list:
     # for testing only
@@ -97,30 +101,6 @@
         s1_view_ct = s1_view[s1_view.obs["cell_type"] == celltype]
         s2_view_ct = s2_view[s2_view.obs["cell_type"] == celltype]
 
-        # safeguard checking the obsnames are exactly the same
-        print(
-            "Check that the obsnames are exactly the same between s1 and unintegrated",
-            flush=True,
-        )
-        # Check membership
-        split_names = np.concatenate([s1_view_ct.obs_names, s2_view_ct.obs_names])
-        u_names = list(u_view_ct.obs_names)
-
-        in_split_not_u = [n for n in split_names if n not in u_names]
-        in_u_not_split = [n for n in u_names if n not in split_names]
-
-        if len(in_split_not_u) > 0 or len(in_u_not_split) > 0:
-            print(
-                f"Error: the obsnames of the unintegrated and s1 and s2 do not match for {donor}, {celltype}."
-                "\nPlease check dataset as this should not happen!"
-                f"\nMissing (in s1/s2 but not in unintegrated): {in_split_not_u}."
-                f"\nExtra (in unintegrated but not in s1/s2): {in_u_not_split}.",
-                flush=True,
-            )
-            sys.exit(1)
-
-        print("Obsnames match!", flush=True)
-
         if s1_view_ct.shape[0] < 100 or s2_view_ct.shape[0] < 100:
             print(f"Skipping celltype {celltype} and donor {donor}.", flush=True)
             if s1_view_ct.shape[0] < 100:
@@ -168,12 +148,12 @@
             peaks_u_s2 = metric_helper.call_peaks(density_dist_u_s2)
 
             # use persistent peak only if the peak calling method is too sensitive...
-            # persistent_peak_count_u_s1, diagram_u_s1, persistence_values_u_s1 = (
-            #     metric_helper.persistent_peak_count(density_dist_u_s1)
-            # )
-            # persistent_peak_count_u_s2, diagram_u_s2, persistence_values_u_s2 = (
-            #     metric_helper.persistent_peak_count(density_dist_u_s2)
-            # )
+            persistent_peak_count_u_s1 = metric_helper.persistent_peak_count(
+                density_dist_u_s1
+            )
+            persistent_peak_count_u_s2 = metric_helper.persistent_peak_count(
+                density_dist_u_s2
+            )
 
             print("--------------------------------", flush=True)
 
@@ -202,18 +182,21 @@
             peaks_s2 = metric_helper.call_peaks(density_dist_s2)
 
             # use persistent peak only if the peak calling method is too sensitive...
-            # persistent_peak_count_s1, diagram_s1, persistence_values_s1 = (
-            #     metric_helper.persistent_peak_count(density_dist_s1)
-            # )
-            # persistent_peak_count_s2, diagram_s2, persistence_values_s2 = (
-            #     metric_helper.persistent_peak_count(density_dist_s2)
-            # )
+            persistent_peak_count_s1 = metric_helper.persistent_peak_count(
+                density_dist_s1
+            )
+            persistent_peak_count_s2 = metric_helper.persistent_peak_count(
+                density_dist_s2
+            )
 
             print("--------------------------------", flush=True)
 
             print("\n", flush=True)
 
-            print("Comparing peaks between unintegrated and integrated", flush=True)
+            print(
+                f"Comparing peaks between unintegrated and integrated for {donor}, {celltype}, {marker}",
+                flush=True,
+            )
 
             # case 1 or 3 where we have consistent peaks in unintegrated
             if peaks_u_s1 == peaks_u_s2:
@@ -228,8 +211,31 @@
                     "WARNING! Inconsistent peaks detected in unintegrated data (case 2 or 4). Skipping calculation",
                     flush=True,
                 )
+                print(
+                    f"Number of peaks in unintegrated split 1: {peaks_u_s1}, split 2: {peaks_u_s2}",
+                    flush=True,
+                )
                 case_details["case2or4"].append((donor, celltype, marker))
 
+            # for comparison only
+            persistent_peaks_res.append(
+                [
+                    donor,
+                    celltype,
+                    marker,
+                    peaks_u_s1,
+                    peaks_u_s2,
+                    persistent_peak_count_u_s1,
+                    persistent_peak_count_u_s2,
+                    peaks_s1,
+                    peaks_s2,
+                    persistent_peak_count_s1,
+                    persistent_peak_count_s2,
+                ]
+            )
+            print("Done comparing peaks.", flush=True)
+            print("\n", flush=True)
+
 print("Done processing all celltypes and donors", flush=True)
 print("Calculating ratio", flush=True)
 
@@ -242,6 +248,22 @@
 else:
     metric_val = n_case3 / (n_case1 + n_case3)
 
+persistent_peaks_res = pd.DataFrame(
+    persistent_peaks_res,
+    columns=[
+        "donor",
+        "celltype",
+        "marker",
+        "peaks_u_s1",
+        "peaks_u_s2",
+        "persistent_peaks_u_s1",
+        "persistent_peaks_u_s2",
+        "peaks_s1",
+        "peaks_s2",
+        "persistent_peaks_s1",
+        "persistent_peaks_s2",
+    ],
+)
 
 print("Write output AnnData to file", flush=True)
 output = ad.AnnData(
@@ -256,6 +278,7 @@
             "case2or4": len(case_details["case2or4"]),
         },
         "case_details": dict(case_details),
+        "peak_calling_results_comparison": persistent_peaks_res,
     }
 )
 output.write_h5ad(par["output"], compression="gzip")

From 4a793f62c9ae9eef63ebfead1843005789545900 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Sun, 12 Oct 2025 11:50:48 +1100
Subject: [PATCH 03/31] update description

---
 .../ratio_inconsistent_peaks/config.vsh.yaml      | 15 +++++++++++----
 1 file changed, 11 insertions(+), 4 deletions(-)

diff --git a/src/metrics/ratio_inconsistent_peaks/config.vsh.yaml b/src/metrics/ratio_inconsistent_peaks/config.vsh.yaml
index bfcc00ee..933c9aad 100644
--- a/src/metrics/ratio_inconsistent_peaks/config.vsh.yaml
+++ b/src/metrics/ratio_inconsistent_peaks/config.vsh.yaml
@@ -16,14 +16,21 @@ info:
       # Can contain only lowercase letters or underscores.
     - name: ratio_inconsistent_peaks
       label: Ratio of inconsistent peaks
-      summary: "Ratio of the number of cell‑type marker‑expression peaks between validation and batch‑normalized data."
+      summary: "Ratio of the number of cell‑type marker‑expression peaks between unintegrated and batch‑normalized data."
       description: |
-        The metric compares the number of cell type specific 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 metric compares the number of cell type specific marker expression peaks between unintegrated 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 unintegrated and batch-normalized data.
         The (cell type) 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.
+        Ratio of inconsistent peaks is defined as number of cases where the number of peaks differ between the two splits in the batch
+        normalized data divided by the total number of cases.
+        Cases where there are different number of peaks between the two splits in the unintegrated data are ignored from the denominator.
+        A lower score indicates better performance, means there are less cases with inconsistent peaks after batch correction.
+        An alternative peak counting method using persistent homology is also implemented for comparison because peak calling 
+        is sensitive to noise and parameter choices.
+        
       references:
         doi: 
           - 10.1038/s41592-019-0686-2

From 0bb39a034e4faf0dcea5994a69aec794e4adc3ee Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Sun, 12 Oct 2025 12:00:23 +1100
Subject: [PATCH 04/31] update workflow script

---
 src/workflows/run_benchmark/config.vsh.yaml | 1 +
 src/workflows/run_benchmark/main.nf         | 1 +
 2 files changed, 2 insertions(+)

diff --git a/src/workflows/run_benchmark/config.vsh.yaml b/src/workflows/run_benchmark/config.vsh.yaml
index c424489d..c104629c 100644
--- a/src/workflows/run_benchmark/config.vsh.yaml
+++ b/src/workflows/run_benchmark/config.vsh.yaml
@@ -110,6 +110,7 @@ dependencies:
   - name: methods/cytovi
   - name: metrics/emd
   - name: metrics/n_inconsistent_peaks
+  - name: metrics/ratio_inconsistent_peaks
   - name: metrics/average_batch_r2
   - name: metrics/flowsom_mapping_similarity
   - name: metrics/lisi
diff --git a/src/workflows/run_benchmark/main.nf b/src/workflows/run_benchmark/main.nf
index be484ce2..0138ac23 100644
--- a/src/workflows/run_benchmark/main.nf
+++ b/src/workflows/run_benchmark/main.nf
@@ -41,6 +41,7 @@ methods = [
 metrics = [
   emd,
   n_inconsistent_peaks,
+  ratio_inconsistent_peaks
   average_batch_r2,
   flowsom_mapping_similarity,
   lisi,

From d36e9c15f00c0210d6ca6249e0b7711d0c2b53e6 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Mon, 13 Oct 2025 09:48:03 +1100
Subject: [PATCH 05/31] add missing comma - facepalm

---
 src/workflows/run_benchmark/main.nf | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/workflows/run_benchmark/main.nf b/src/workflows/run_benchmark/main.nf
index 0138ac23..0cd386eb 100644
--- a/src/workflows/run_benchmark/main.nf
+++ b/src/workflows/run_benchmark/main.nf
@@ -41,7 +41,7 @@ methods = [
 metrics = [
   emd,
   n_inconsistent_peaks,
-  ratio_inconsistent_peaks
+  ratio_inconsistent_peaks,
   average_batch_r2,
   flowsom_mapping_similarity,
   lisi,

From f7be6f2d49d871412a83f1730601ef44f0d286a9 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Mon, 13 Oct 2025 23:46:45 +1100
Subject: [PATCH 06/31] add small epsilon to harmonypy to fix kmeans bug

---
 src/methods/harmonypy/script.py | 12 ++++++++----
 1 file changed, 8 insertions(+), 4 deletions(-)

diff --git a/src/methods/harmonypy/script.py b/src/methods/harmonypy/script.py
index 8fc0cccc..5ac69a5d 100644
--- a/src/methods/harmonypy/script.py
+++ b/src/methods/harmonypy/script.py
@@ -4,8 +4,8 @@
 
 ## 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_harmony_split2.h5ad",
+    "input": "/Users/putri.g/Documents/cytobenchmark/debug_general/_viash_par/input_1/censored_split1.h5ad",
+    "output": "/Users/putri.g/Documents/cytobenchmark/debug_general/_viash_par/output_1/output_harmony_split1.h5ad",
 }
 meta = {"name": "harmonypy"}
 ## VIASH END
@@ -13,6 +13,7 @@
 print("Reading and preparing input files", flush=True)
 adata = ad.read_h5ad(par["input"])
 
+# harmony can't handle integer batch labels
 adata.obs["batch_str"] = adata.obs["batch"].astype(str)
 
 markers_to_correct = adata.var[adata.var["to_correct"]].index.to_numpy()
@@ -21,10 +22,13 @@
 adata_to_correct = adata[:, markers_to_correct].copy()
 
 print("Run harmony", flush=True)
-# harmony can't handle integer batch labels
+
+# TODO numerical instability in kmeans causing problem with harmony.
+# so adding a very small value to all entries to make sure there are no zeros
+epsilon = 1e-20
 
 out = harmonypy.run_harmony(
-    data_mat=adata_to_correct.layers["preprocessed"],
+    data_mat=adata_to_correct.layers["preprocessed"] + epsilon,
     meta_data=adata_to_correct.obs,
     vars_use="batch_str",
 )

From a0740ccc0bfe897044d6e4098048acd4a1303b82 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Mon, 13 Oct 2025 23:46:54 +1100
Subject: [PATCH 07/31] increase bras chunk to 1000

---
 src/metrics/bras/script.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/metrics/bras/script.py b/src/metrics/bras/script.py
index 98229ad0..0ac51fb8 100644
--- a/src/metrics/bras/script.py
+++ b/src/metrics/bras/script.py
@@ -57,6 +57,7 @@
     labels=ct_labels_s1,
     batch=batch_labels_s1,
     metric="euclidean",
+    chunk_size=1000,
 )
 
 batch_labels_s2 = integrated_s2.obs["batch"].values

From b6beb8317c9145577d74ea0af7f06a575b57b22a Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Tue, 14 Oct 2025 13:56:57 +1100
Subject: [PATCH 08/31] testing bras with jax gpu

---
 src/metrics/bras/config.vsh.yaml | 37 +++++++++++++++++++++-----------
 src/metrics/bras/script.py       |  1 -
 2 files changed, 24 insertions(+), 14 deletions(-)

diff --git a/src/metrics/bras/config.vsh.yaml b/src/metrics/bras/config.vsh.yaml
index fed70ef3..6d7434cf 100644
--- a/src/metrics/bras/config.vsh.yaml
+++ b/src/metrics/bras/config.vsh.yaml
@@ -56,13 +56,6 @@ info:
       # 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)
@@ -73,24 +66,42 @@ resources:
 
 engines:
   # Specifications for the Docker image for this component.
+  # testing gpu jax version
   - type: docker
-    image: python:3.11
+    image: nvidia/cuda:13.0.1-cudnn-devel-ubuntu24.04
     setup:
       - type: apt
         packages:
           - procps
+          - git
       - type: python
         packages:
-          - jax~=0.6.2
-          - jaxlib~=0.6.2
           - anndata~=0.11.0
           - scanpy~=1.11.0
-          - scib-metrics~=0.5.6
           - pyyaml
           - requests
           - jsonschema
+          - scikit-learn
         github:
-          - "openproblems-bio/core#subdirectory=packages/python/openproblems"
+          - openproblems-bio/core#subdirectory=packages/python/openproblems
+  # - 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.6
+  #         - pyyaml
+  #         - requests
+  #         - jsonschema
+  #       github:
+  #         - "openproblems-bio/core#subdirectory=packages/python/openproblems"
 
 runners:
   # This platform allows running the component natively
@@ -98,4 +109,4 @@ runners:
   # Allows turning the component into a Nextflow module / pipeline.
   - type: nextflow
     directives:
-      label: [midtime,midmem,midcpu]
+      label: [midtime,midmem,midcpu,gpu]
diff --git a/src/metrics/bras/script.py b/src/metrics/bras/script.py
index 0ac51fb8..98229ad0 100644
--- a/src/metrics/bras/script.py
+++ b/src/metrics/bras/script.py
@@ -57,7 +57,6 @@
     labels=ct_labels_s1,
     batch=batch_labels_s1,
     metric="euclidean",
-    chunk_size=1000,
 )
 
 batch_labels_s2 = integrated_s2.obs["batch"].values

From 8e65e0327a940afea14f8d5369bf7510534d4262 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Tue, 14 Oct 2025 13:57:11 +1100
Subject: [PATCH 09/31] testing bras with cuda

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

diff --git a/src/metrics/bras/config.vsh.yaml b/src/metrics/bras/config.vsh.yaml
index 6d7434cf..84a3583f 100644
--- a/src/metrics/bras/config.vsh.yaml
+++ b/src/metrics/bras/config.vsh.yaml
@@ -76,6 +76,7 @@ engines:
           - git
       - type: python
         packages:
+          - jax[cuda_13]
           - anndata~=0.11.0
           - scanpy~=1.11.0
           - pyyaml

From 9569f6be1741b5fb454f7ad97810d879a0512eb4 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Tue, 14 Oct 2025 14:06:42 +1100
Subject: [PATCH 10/31] missing pip and sci-b metrics *facepalm

---
 src/metrics/bras/config.vsh.yaml | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/src/metrics/bras/config.vsh.yaml b/src/metrics/bras/config.vsh.yaml
index 84a3583f..dcd8303e 100644
--- a/src/metrics/bras/config.vsh.yaml
+++ b/src/metrics/bras/config.vsh.yaml
@@ -72,6 +72,7 @@ engines:
     setup:
       - type: apt
         packages:
+          - python3-pip
           - procps
           - git
       - type: python
@@ -79,6 +80,7 @@ engines:
           - jax[cuda_13]
           - anndata~=0.11.0
           - scanpy~=1.11.0
+          - scib-metrics~=0.5.7
           - pyyaml
           - requests
           - jsonschema
@@ -110,4 +112,4 @@ runners:
   # Allows turning the component into a Nextflow module / pipeline.
   - type: nextflow
     directives:
-      label: [midtime,midmem,midcpu,gpu]
+      label: [midtime,midmem,lowcpu,gpu]

From 42e8bd21b81e5cef28d1b61489ca9e9d5828afe1 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Tue, 14 Oct 2025 14:51:10 +1100
Subject: [PATCH 11/31] downgrading image

---
 scripts/run_benchmark/run_full_seqeracloud.sh | 6 +++---
 src/metrics/bras/config.vsh.yaml              | 5 ++---
 2 files changed, 5 insertions(+), 6 deletions(-)

diff --git a/scripts/run_benchmark/run_full_seqeracloud.sh b/scripts/run_benchmark/run_full_seqeracloud.sh
index 979440c9..e59eba42 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: '{"metrics_exclude": ["cms"], "methods_include": ["mnnpy", "cytovi"]}'
+settings: '{"metrics_include": ["emd", "ratio_inconsistent_peaks", "n_inconsistent_peaks"], "methods_include": ["harmonypy", "cycombine_no_controls_to_goal", "cycombine_all_controls_to_goal", "cytonorm_no_controls_to_goal", "cytonorm_all_controls_to_goal"]}'
 publish_dir: "$publish_dir"
 HERE
 
 tw launch https://github.com/openproblems-bio/task_cyto_batch_integration.git \
-  --revision build/fix_failed_stuff \
+  --revision build/update_n_inconsistent_peak \
   --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,mnnnpy
+  --labels task_cyto_batch_integration,test_subset
diff --git a/src/metrics/bras/config.vsh.yaml b/src/metrics/bras/config.vsh.yaml
index dcd8303e..cecca8e7 100644
--- a/src/metrics/bras/config.vsh.yaml
+++ b/src/metrics/bras/config.vsh.yaml
@@ -68,7 +68,7 @@ engines:
   # Specifications for the Docker image for this component.
   # testing gpu jax version
   - type: docker
-    image: nvidia/cuda:13.0.1-cudnn-devel-ubuntu24.04
+    image: nvidia/cuda:12.4-runtime-ubuntu22.04
     setup:
       - type: apt
         packages:
@@ -77,14 +77,13 @@ engines:
           - git
       - type: python
         packages:
-          - jax[cuda_13]
+          - jax[cuda_12_pip]
           - anndata~=0.11.0
           - scanpy~=1.11.0
           - scib-metrics~=0.5.7
           - pyyaml
           - requests
           - jsonschema
-          - scikit-learn
         github:
           - openproblems-bio/core#subdirectory=packages/python/openproblems
   # - type: docker

From a7cef47b1f3df4da5fde6f961c0838d2f16f2705 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Tue, 14 Oct 2025 16:06:38 +1100
Subject: [PATCH 12/31] update the image name

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

diff --git a/src/metrics/bras/config.vsh.yaml b/src/metrics/bras/config.vsh.yaml
index cecca8e7..c8d26438 100644
--- a/src/metrics/bras/config.vsh.yaml
+++ b/src/metrics/bras/config.vsh.yaml
@@ -68,7 +68,7 @@ engines:
   # Specifications for the Docker image for this component.
   # testing gpu jax version
   - type: docker
-    image: nvidia/cuda:12.4-runtime-ubuntu22.04
+    image: nvidia/cuda:12.4.0-runtime-ubuntu22.04
     setup:
       - type: apt
         packages:

From 913a70b5fcb8b7e82b37a6707da78e86a4e51227 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Tue, 14 Oct 2025 20:29:47 +1100
Subject: [PATCH 13/31] testing openproblems image

---
 src/metrics/bras/config.vsh.yaml | 14 +-------------
 1 file changed, 1 insertion(+), 13 deletions(-)

diff --git a/src/metrics/bras/config.vsh.yaml b/src/metrics/bras/config.vsh.yaml
index c8d26438..788b2fd2 100644
--- a/src/metrics/bras/config.vsh.yaml
+++ b/src/metrics/bras/config.vsh.yaml
@@ -68,24 +68,12 @@ engines:
   # Specifications for the Docker image for this component.
   # testing gpu jax version
   - type: docker
-    image: nvidia/cuda:12.4.0-runtime-ubuntu22.04
+    image: openproblems/base_pytorch_nvidia:1.1
     setup:
-      - type: apt
-        packages:
-          - python3-pip
-          - procps
-          - git
       - type: python
         packages:
           - jax[cuda_12_pip]
-          - anndata~=0.11.0
-          - scanpy~=1.11.0
           - scib-metrics~=0.5.7
-          - pyyaml
-          - requests
-          - jsonschema
-        github:
-          - openproblems-bio/core#subdirectory=packages/python/openproblems
   # - type: docker
   #   image: python:3.11
   #   setup:

From f8c9a4dbc227dff7f904a79427984fcbf6948abf Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Tue, 14 Oct 2025 20:30:07 +1100
Subject: [PATCH 14/31] undo changes to run script

---
 scripts/run_benchmark/run_full_seqeracloud.sh | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/scripts/run_benchmark/run_full_seqeracloud.sh b/scripts/run_benchmark/run_full_seqeracloud.sh
index e59eba42..baa19b65 100755
--- a/scripts/run_benchmark/run_full_seqeracloud.sh
+++ b/scripts/run_benchmark/run_full_seqeracloud.sh
@@ -17,12 +17,11 @@ 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: '{"metrics_include": ["emd", "ratio_inconsistent_peaks", "n_inconsistent_peaks"], "methods_include": ["harmonypy", "cycombine_no_controls_to_goal", "cycombine_all_controls_to_goal", "cytonorm_no_controls_to_goal", "cytonorm_all_controls_to_goal"]}'
 publish_dir: "$publish_dir"
 HERE
 
 tw launch https://github.com/openproblems-bio/task_cyto_batch_integration.git \
-  --revision build/update_n_inconsistent_peak \
+  --revision build/main \
   --pull-latest \
   --main-script target/nextflow/workflows/run_benchmark/main.nf \
   --workspace 53907369739130 \

From 320fd3ca4c1864c2ed47b79f41bcc97d364e21b2 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Tue, 14 Oct 2025 20:55:21 +1100
Subject: [PATCH 15/31] downgrade scib metrics package

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

diff --git a/src/metrics/bras/config.vsh.yaml b/src/metrics/bras/config.vsh.yaml
index 788b2fd2..a5b73c6a 100644
--- a/src/metrics/bras/config.vsh.yaml
+++ b/src/metrics/bras/config.vsh.yaml
@@ -73,7 +73,7 @@ engines:
       - type: python
         packages:
           - jax[cuda_12_pip]
-          - scib-metrics~=0.5.7
+          - scib-metrics~=0.5.6
   # - type: docker
   #   image: python:3.11
   #   setup:

From 7c98d9684a8c91eb0987b1c756ff479e26e9246d Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Wed, 15 Oct 2025 09:57:03 +1100
Subject: [PATCH 16/31] reverting as gpu doesn't work

---
 src/metrics/bras/config.vsh.yaml | 42 ++++++++++++++++----------------
 1 file changed, 21 insertions(+), 21 deletions(-)

diff --git a/src/metrics/bras/config.vsh.yaml b/src/metrics/bras/config.vsh.yaml
index a5b73c6a..b7348dc1 100644
--- a/src/metrics/bras/config.vsh.yaml
+++ b/src/metrics/bras/config.vsh.yaml
@@ -67,31 +67,31 @@ resources:
 engines:
   # Specifications for the Docker image for this component.
   # testing gpu jax version
-  - type: docker
-    image: openproblems/base_pytorch_nvidia:1.1
-    setup:
-      - type: python
-        packages:
-          - jax[cuda_12_pip]
-          - scib-metrics~=0.5.6
   # - type: docker
-  #   image: python:3.11
+  #   image: openproblems/base_pytorch_nvidia:1.1
   #   setup:
-  #     - type: apt
-  #       packages:
-  #         - procps
   #     - type: python
   #       packages:
-  #         - jax~=0.6.2
-  #         - jaxlib~=0.6.2
-  #         - anndata~=0.11.0
-  #         - scanpy~=1.11.0
+  #         - jax[cuda_12_pip]
   #         - scib-metrics~=0.5.6
-  #         - pyyaml
-  #         - requests
-  #         - jsonschema
-  #       github:
-  #         - "openproblems-bio/core#subdirectory=packages/python/openproblems"
+  - 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.6
+          - pyyaml
+          - requests
+          - jsonschema
+        github:
+          - "openproblems-bio/core#subdirectory=packages/python/openproblems"
 
 runners:
   # This platform allows running the component natively
@@ -99,4 +99,4 @@ runners:
   # Allows turning the component into a Nextflow module / pipeline.
   - type: nextflow
     directives:
-      label: [midtime,midmem,lowcpu,gpu]
+      label: [midtime,midmem,midcpu]

From d9e146f738de54f80b41483bc5781237b22b608e Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Sat, 1 Nov 2025 22:51:17 +1100
Subject: [PATCH 17/31] testing bin shifting

---
 common                                         |  2 +-
 src/metrics/ratio_inconsistent_peaks/helper.py | 17 +++++++++++++++++
 2 files changed, 18 insertions(+), 1 deletion(-)

diff --git a/common b/common
index 79b884b4..67da19a3 160000
--- a/common
+++ b/common
@@ -1 +1 @@
-Subproject commit 79b884b4c7fed300972d83a6ca025abb6116cbdc
+Subproject commit 67da19a36ae56ea068804d15ccadec88a06da920
diff --git a/src/metrics/ratio_inconsistent_peaks/helper.py b/src/metrics/ratio_inconsistent_peaks/helper.py
index 008eb211..04363f8c 100644
--- a/src/metrics/ratio_inconsistent_peaks/helper.py
+++ b/src/metrics/ratio_inconsistent_peaks/helper.py
@@ -48,6 +48,15 @@ def get_kde_density(expression_array, return_xgrid=False, plot=False):
     x_grid = np.linspace(min_val, max_val, 100)
     density = kde(x_grid)
 
+    # If the highest value is at the first bin, shift bins by one and adjust x_grid
+    if np.argmax(density) == 0 and density.size > 1:
+        # Prepend a zero so the former first bin becomes index 1
+        density = np.concatenate([[0.0], density])[: density.size]
+        # Have to use actual grid spacing to keep uniform spacing in x_grid.
+        # Can't just blindly add 1.
+        step = (max_val - min_val) / (len(x_grid) - 1) if len(x_grid) > 1 else 0.0
+        x_grid = np.concatenate([[min_val - step], x_grid])[: x_grid.size]
+
     if plot:
         fig, ax = plt.subplots()
         sns.scatterplot(x=x_grid, y=density, ax=ax)
@@ -102,6 +111,14 @@ def persistent_peak_count(ys, persistence_cutoff=0.08):
         int: number of significant peaks
     """
 
+    y = np.asarray(ys)
+    if y.size == 0:
+        return 0
+
+    # Shift if max is at the first bin
+    if y.size > 1 and np.argmax(y) == 0:
+        y = np.concatenate([[0.0], y[:-1]])
+
     # Invert to turn peaks into "holes" for 0D persistence
     Y = -ys.reshape(-1, 1)
     diagram = ripser(Y, maxdim=0)["dgms"][0]

From ea44a387ec504b97016453f8ca8ed19a0dcaab04 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Sat, 1 Nov 2025 22:58:06 +1100
Subject: [PATCH 18/31] increase chunk size for bras

---
 src/metrics/bras/script.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/metrics/bras/script.py b/src/metrics/bras/script.py
index 98229ad0..db3428a1 100644
--- a/src/metrics/bras/script.py
+++ b/src/metrics/bras/script.py
@@ -57,6 +57,7 @@
     labels=ct_labels_s1,
     batch=batch_labels_s1,
     metric="euclidean",
+    chunk_size=512,
 )
 
 batch_labels_s2 = integrated_s2.obs["batch"].values
@@ -67,6 +68,7 @@
     labels=ct_labels_s2,
     batch=batch_labels_s2,
     metric="euclidean",
+    chunk_size=512,
 )
 
 bras_score = np.mean([bras_s1, bras_s2])

From c50a0c0d52455d4e33fa01dc7ddd39ad3b2087eb Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Mon, 3 Nov 2025 23:43:04 +1100
Subject: [PATCH 19/31] small changes to the bin

---
 src/metrics/ratio_inconsistent_peaks/helper.py | 15 +++++++++++----
 src/metrics/ratio_inconsistent_peaks/script.py | 16 +++++++++++-----
 2 files changed, 22 insertions(+), 9 deletions(-)

diff --git a/src/metrics/ratio_inconsistent_peaks/helper.py b/src/metrics/ratio_inconsistent_peaks/helper.py
index 04363f8c..461f6915 100644
--- a/src/metrics/ratio_inconsistent_peaks/helper.py
+++ b/src/metrics/ratio_inconsistent_peaks/helper.py
@@ -50,12 +50,19 @@ def get_kde_density(expression_array, return_xgrid=False, plot=False):
 
     # If the highest value is at the first bin, shift bins by one and adjust x_grid
     if np.argmax(density) == 0 and density.size > 1:
-        # Prepend a zero so the former first bin becomes index 1
-        density = np.concatenate([[0.0], density])[: density.size]
+        print("Shifting KDE bins by one as the highest density is at the first bin.")
+        # orig_x_grid = x_grid.copy()
+        # recale the grid so we only have 99 bins and shift everything by one to the right..
+        x_grid = np.linspace(min_val, max_val, 99)
+        density = kde(x_grid)
+
+        # Prepend a zero so the beginning, but remove the last value to keep size consistent
+        # as otherwise we will end up with an extra bin...
+        density = np.concatenate([[0.0], density])
         # Have to use actual grid spacing to keep uniform spacing in x_grid.
         # Can't just blindly add 1.
-        step = (max_val - min_val) / (len(x_grid) - 1) if len(x_grid) > 1 else 0.0
-        x_grid = np.concatenate([[min_val - step], x_grid])[: x_grid.size]
+        step = (max_val - min_val) / (len(x_grid)) if len(x_grid) > 1 else 0.0
+        x_grid = np.concatenate([[min_val - step], x_grid])
 
     if plot:
         fig, ax = plt.subplots()
diff --git a/src/metrics/ratio_inconsistent_peaks/script.py b/src/metrics/ratio_inconsistent_peaks/script.py
index 9bc26657..2102bec4 100644
--- a/src/metrics/ratio_inconsistent_peaks/script.py
+++ b/src/metrics/ratio_inconsistent_peaks/script.py
@@ -9,6 +9,7 @@
 # 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",
+    # "input_unintegrated": "/Users/putri.g/Documents/cytobenchmark/benchmark_out_20251015/human_blood_mass_cytometry/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",
@@ -116,6 +117,10 @@
             # TODO uncomment me when done
             continue
 
+        # unintegrated for split 1
+        u_view_ct_s1 = u_view_ct[u_view_ct.obs["split"] == 1]
+        u_view_ct_s2 = u_view_ct[u_view_ct.obs["split"] == 2]
+
         for marker in s1_view_ct.var.index:
             # for testing only
             # marker = u_view_ct.var.index[0]
@@ -124,9 +129,6 @@
 
             print("--------------------------------", flush=True)
             print("Computing peaks for unintegrated", flush=True)
-            # unintegrated for split 1
-            u_view_ct_s1 = u_view_ct[u_view_ct.obs["split"] == 1]
-            u_view_ct_s2 = u_view_ct[u_view_ct.obs["split"] == 2]
 
             print("Standardising marker expression", flush=True)
             # standardise marker expression based on pooled mean and sd of
@@ -139,8 +141,12 @@
                 u_s2_unscaled,
             )
             print("Computing KDE density", flush=True)
-            density_dist_u_s1 = metric_helper.get_kde_density(u_s1_scaled)
-            density_dist_u_s2 = metric_helper.get_kde_density(u_s2_scaled)
+            density_dist_u_s1 = metric_helper.get_kde_density(
+                expression_array=u_s1_scaled
+            )
+            density_dist_u_s2 = metric_helper.get_kde_density(
+                expression_array=u_s2_scaled
+            )
 
             print("Calling peaks", flush=True)
 

From 695562a3660c8c8c7ef95e8fd6b780a2c873f559 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Wed, 5 Nov 2025 00:14:10 +1100
Subject: [PATCH 20/31] disable metrics

---
 src/metrics/bras/config.vsh.yaml                 | 9 +--------
 src/metrics/n_inconsistent_peaks/config.vsh.yaml | 2 +-
 2 files changed, 2 insertions(+), 9 deletions(-)

diff --git a/src/metrics/bras/config.vsh.yaml b/src/metrics/bras/config.vsh.yaml
index b7348dc1..a99b43f5 100644
--- a/src/metrics/bras/config.vsh.yaml
+++ b/src/metrics/bras/config.vsh.yaml
@@ -1,16 +1,9 @@
-# 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: bras
-
-
-
+status: disabled
 # Metadata for your component
 info:
   metrics:
diff --git a/src/metrics/n_inconsistent_peaks/config.vsh.yaml b/src/metrics/n_inconsistent_peaks/config.vsh.yaml
index 569be185..9c8bae0e 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:

From 48fb41c1ebcc45ac1faa641d33194bbfe59f2479 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Wed, 5 Nov 2025 10:03:28 +1100
Subject: [PATCH 21/31] disable two metrics again

---
 src/workflows/run_benchmark/config.vsh.yaml | 5 +++--
 src/workflows/run_benchmark/main.nf         | 6 +++---
 2 files changed, 6 insertions(+), 5 deletions(-)

diff --git a/src/workflows/run_benchmark/config.vsh.yaml b/src/workflows/run_benchmark/config.vsh.yaml
index c104629c..de188e22 100644
--- a/src/workflows/run_benchmark/config.vsh.yaml
+++ b/src/workflows/run_benchmark/config.vsh.yaml
@@ -109,12 +109,13 @@ dependencies:
   - name: methods/rpca_to_mid
   - name: methods/cytovi
   - name: metrics/emd
-  - name: metrics/n_inconsistent_peaks
+  # - name: metrics/bras
+  # - name: metrics/n_inconsistent_peaks
   - name: metrics/ratio_inconsistent_peaks
   - name: metrics/average_batch_r2
   - name: metrics/flowsom_mapping_similarity
   - name: metrics/lisi
-  - name: metrics/bras
+  
 
 runners:
   - type: nextflow
diff --git a/src/workflows/run_benchmark/main.nf b/src/workflows/run_benchmark/main.nf
index 0cd386eb..1bb72733 100644
--- a/src/workflows/run_benchmark/main.nf
+++ b/src/workflows/run_benchmark/main.nf
@@ -40,12 +40,12 @@ methods = [
 // construct list of metrics
 metrics = [
   emd,
-  n_inconsistent_peaks,
+  // bras,
+  // n_inconsistent_peaks,
   ratio_inconsistent_peaks,
   average_batch_r2,
   flowsom_mapping_similarity,
-  lisi,
-  bras
+  lisi
 ]
 
 workflow run_wf {

From 956325241c42f6de1c87cfcf51ad68fe6ac9b3e2 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Wed, 19 Nov 2025 00:00:42 +1100
Subject: [PATCH 22/31] update cytovi

---
 src/methods/cytovi/config.vsh.yaml |  3 +-
 src/methods/cytovi/script.py       | 82 +++++++++++++++++++-----------
 2 files changed, 54 insertions(+), 31 deletions(-)

diff --git a/src/methods/cytovi/config.vsh.yaml b/src/methods/cytovi/config.vsh.yaml
index eb6ff7d8..206df6a6 100644
--- a/src/methods/cytovi/config.vsh.yaml
+++ b/src/methods/cytovi/config.vsh.yaml
@@ -68,13 +68,14 @@ engines:
         packages:
           - anndata>=0.11.0
           - scanpy[skmisc]>=1.10
-          - scvi-tools==1.4.0
+          # - scvi-tools==1.4.0
           - pyyaml
           - requests
           - jsonschema
           - scikit-learn
         github:
           - openproblems-bio/core#subdirectory=packages/python/openproblems
+          - YosefLab/cytovi-reference-implementation
 
 runners:
   # This platform allows running the component natively
diff --git a/src/methods/cytovi/script.py b/src/methods/cytovi/script.py
index 6a9768af..9716535b 100644
--- a/src/methods/cytovi/script.py
+++ b/src/methods/cytovi/script.py
@@ -1,8 +1,11 @@
 import anndata as ad
+import cytovi
 import numpy as np
-from scvi.external import cytovi
-from sklearn.cluster import KMeans
-from threadpoolctl import threadpool_limits
+
+# from scvi.external import cytovi
+
+# from sklearn.cluster import KMeans
+# from threadpoolctl import threadpool_limits
 
 ## VIASH START
 par = {
@@ -20,6 +23,7 @@
 adata = ad.read_h5ad(par["input"])
 
 adata.obs["batch_str"] = adata.obs["batch"].astype(str)
+adata.obs["sample_key_str"] = adata.obs["sample"].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()
@@ -29,46 +33,64 @@
 print("Scaling data", flush=True)
 
 # scale data. this will add a layer "scaled" to the anndata
-cytovi.scale(
+cytovi.pp.scale(
     adata=adata_to_correct,
     transformed_layer_key="preprocessed",
     batch_key="batch_str",
     inplace=True,
 )
 
-print("Clustering using k-means with k =", par["n_clusters"], flush=True)
-# cluster data using Kmeans
-with threadpool_limits(limits=1):
-    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",
+    f"Train CytoVI on {adata_to_correct.shape[0]} cells",
     flush=True,
 )
 
-cytovi.CYTOVI.setup_anndata(adata_subsampled, layer="scaled", batch_key="batch_str")
-model = cytovi.CYTOVI(
-    adata=adata_subsampled, n_hidden=par["n_hidden"], n_layers=par["n_layers"]
+cytovi.CytoVI.setup_anndata(
+    adata_to_correct, layer="scaled", batch_key="batch_str", sample_key="sample_key_str"
+)
+
+model = cytovi.CytoVI(
+    adata_to_correct, n_hidden=par["n_hidden"], n_layers=par["n_layers"]
 )
+
+
+print("Start training CytoVI model", flush=True)
 model.train()
 
+# Todo: re-enable subsampling if needed..
+# print("Clustering using k-means with k =", par["n_clusters"], flush=True)
+# # cluster data using Kmeans
+# with threadpool_limits(limits=1):
+#     adata_to_correct.obs["clusters"] = (
+#         KMeans(n_clusters=par["n_clusters"], random_state=0)
+#         .fit_predict(adata_to_correct.layers["scaled"])
+#         .astype(str)
+#     )
+# # concatenate obs so we can use it for subsampling
+# adata_to_correct.obs["sample_cluster"] = (
+#     adata_to_correct.obs["sample"].astype(str) + "_" + adata_to_correct.obs["clusters"]
+# )
+# # subsample cells without replacement
+# print("Subsampling cells", flush=True)
+# subsampled_cells = adata_to_correct.obs.groupby("sample_cluster")[
+#     "sample_cluster"
+# ].apply(lambda x: x.sample(n=round(len(x) * par["subsample_fraction"]), replace=False))
+# # need the cell id included in the subsample
+# subsampled_cells_idx = [x[1] for x in subsampled_cells.index.to_list()]
+
+# adata_subsampled = adata_to_correct[subsampled_cells_idx, :].copy()
+
+# print(
+#     f"Train CytoVI on subsampled data containing {adata_subsampled.shape[0]} cells",
+#     flush=True,
+# )
+
+# cytovi.CYTOVI.setup_anndata(adata_subsampled, layer="scaled", batch_key="batch_str")
+# model = cytovi.CYTOVI(
+#     adata=adata_subsampled, n_hidden=par["n_hidden"], n_layers=par["n_layers"]
+# )
+# model.train()
+
 # get batch corrected data
 print("Correcting data", flush=True)
 corrected_data = model.get_normalized_expression(adata=adata_to_correct)

From d361dd0bd7ad25a7f6d9e6e1b07d740fd329f7b3 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Thu, 20 Nov 2025 10:38:53 +1100
Subject: [PATCH 23/31] switched training to TF32

---
 src/methods/cytovi/script.py | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/src/methods/cytovi/script.py b/src/methods/cytovi/script.py
index 9716535b..658f9bba 100644
--- a/src/methods/cytovi/script.py
+++ b/src/methods/cytovi/script.py
@@ -1,6 +1,7 @@
 import anndata as ad
 import cytovi
 import numpy as np
+import torch
 
 # from scvi.external import cytovi
 
@@ -19,6 +20,9 @@
 meta = {"name": "cytovi"}
 ## VIASH END
 
+# setting calculation to TF32 to speed up training
+torch.backends.cuda.matmul.allow_tf32 = True
+
 print("Reading and preparing input files", flush=True)
 adata = ad.read_h5ad(par["input"])
 

From 328929e17fad8ac0c7676d3d92c483db0219fc61 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Thu, 20 Nov 2025 22:49:48 +1100
Subject: [PATCH 24/31] remove persistent peaks

---
 .../ratio_inconsistent_peaks/config.vsh.yaml  |  8 ++-
 .../ratio_inconsistent_peaks/script.py        | 53 -------------------
 2 files changed, 3 insertions(+), 58 deletions(-)

diff --git a/src/metrics/ratio_inconsistent_peaks/config.vsh.yaml b/src/metrics/ratio_inconsistent_peaks/config.vsh.yaml
index 933c9aad..d4ca192e 100644
--- a/src/metrics/ratio_inconsistent_peaks/config.vsh.yaml
+++ b/src/metrics/ratio_inconsistent_peaks/config.vsh.yaml
@@ -16,9 +16,9 @@ info:
       # Can contain only lowercase letters or underscores.
     - name: ratio_inconsistent_peaks
       label: Ratio of inconsistent peaks
-      summary: "Ratio of the number of cell‑type marker‑expression peaks between unintegrated and batch‑normalized data."
+      summary: "Ratio of the number of cell‑type marker‑expression peaks between unintegrated and batch-integrated data."
       description: |
-        The metric compares the number of cell type specific marker expression peaks between unintegrated and batch normalized data.
+        The metric compares the number of cell type specific marker expression peaks between unintegrated and batch integrated 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 unintegrated and batch-normalized data.
         The (cell type) marker expression profiles are first smoothed using kernel density estimation (KDE) (`scipy.stats.gaussian_kde`),
@@ -27,9 +27,7 @@ info:
         Ratio of inconsistent peaks is defined as number of cases where the number of peaks differ between the two splits in the batch
         normalized data divided by the total number of cases.
         Cases where there are different number of peaks between the two splits in the unintegrated data are ignored from the denominator.
-        A lower score indicates better performance, means there are less cases with inconsistent peaks after batch correction.
-        An alternative peak counting method using persistent homology is also implemented for comparison because peak calling 
-        is sensitive to noise and parameter choices.
+        A lower score indicates better performance, means there are less cases with inconsistent peaks after batch integration.
         
       references:
         doi: 
diff --git a/src/metrics/ratio_inconsistent_peaks/script.py b/src/metrics/ratio_inconsistent_peaks/script.py
index 2102bec4..e1aeae75 100644
--- a/src/metrics/ratio_inconsistent_peaks/script.py
+++ b/src/metrics/ratio_inconsistent_peaks/script.py
@@ -74,10 +74,6 @@
 # so we can see where each cases comes from
 case_details = defaultdict(list)
 
-# for comparison only
-persistent_peaks_res = []
-
-
 for donor in donor_list:
     # for testing only
     # donor = donor_list[0]
@@ -153,14 +149,6 @@
             peaks_u_s1 = metric_helper.call_peaks(density_dist_u_s1)
             peaks_u_s2 = metric_helper.call_peaks(density_dist_u_s2)
 
-            # use persistent peak only if the peak calling method is too sensitive...
-            persistent_peak_count_u_s1 = metric_helper.persistent_peak_count(
-                density_dist_u_s1
-            )
-            persistent_peak_count_u_s2 = metric_helper.persistent_peak_count(
-                density_dist_u_s2
-            )
-
             print("--------------------------------", flush=True)
 
             print("\n", flush=True)
@@ -187,14 +175,6 @@
             peaks_s1 = metric_helper.call_peaks(density_dist_s1)
             peaks_s2 = metric_helper.call_peaks(density_dist_s2)
 
-            # use persistent peak only if the peak calling method is too sensitive...
-            persistent_peak_count_s1 = metric_helper.persistent_peak_count(
-                density_dist_s1
-            )
-            persistent_peak_count_s2 = metric_helper.persistent_peak_count(
-                density_dist_s2
-            )
-
             print("--------------------------------", flush=True)
 
             print("\n", flush=True)
@@ -223,22 +203,6 @@
                 )
                 case_details["case2or4"].append((donor, celltype, marker))
 
-            # for comparison only
-            persistent_peaks_res.append(
-                [
-                    donor,
-                    celltype,
-                    marker,
-                    peaks_u_s1,
-                    peaks_u_s2,
-                    persistent_peak_count_u_s1,
-                    persistent_peak_count_u_s2,
-                    peaks_s1,
-                    peaks_s2,
-                    persistent_peak_count_s1,
-                    persistent_peak_count_s2,
-                ]
-            )
             print("Done comparing peaks.", flush=True)
             print("\n", flush=True)
 
@@ -254,22 +218,6 @@
 else:
     metric_val = n_case3 / (n_case1 + n_case3)
 
-persistent_peaks_res = pd.DataFrame(
-    persistent_peaks_res,
-    columns=[
-        "donor",
-        "celltype",
-        "marker",
-        "peaks_u_s1",
-        "peaks_u_s2",
-        "persistent_peaks_u_s1",
-        "persistent_peaks_u_s2",
-        "peaks_s1",
-        "peaks_s2",
-        "persistent_peaks_s1",
-        "persistent_peaks_s2",
-    ],
-)
 
 print("Write output AnnData to file", flush=True)
 output = ad.AnnData(
@@ -284,7 +232,6 @@
             "case2or4": len(case_details["case2or4"]),
         },
         "case_details": dict(case_details),
-        "peak_calling_results_comparison": persistent_peaks_res,
     }
 )
 output.write_h5ad(par["output"], compression="gzip")

From 51552a4405b0b62c20226db2ec2943c0969c1e3f Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Thu, 20 Nov 2025 22:51:57 +1100
Subject: [PATCH 25/31] removed scaling from cytovi

---
 src/methods/cytovi/script.py | 16 ++++------------
 1 file changed, 4 insertions(+), 12 deletions(-)

diff --git a/src/methods/cytovi/script.py b/src/methods/cytovi/script.py
index 658f9bba..47ebb941 100644
--- a/src/methods/cytovi/script.py
+++ b/src/methods/cytovi/script.py
@@ -34,30 +34,22 @@
 
 adata_to_correct = adata[:, markers_to_correct].copy()
 
-print("Scaling data", flush=True)
-
-# scale data. this will add a layer "scaled" to the anndata
-cytovi.pp.scale(
-    adata=adata_to_correct,
-    transformed_layer_key="preprocessed",
-    batch_key="batch_str",
-    inplace=True,
-)
-
 print(
     f"Train CytoVI on {adata_to_correct.shape[0]} cells",
     flush=True,
 )
 
 cytovi.CytoVI.setup_anndata(
-    adata_to_correct, layer="scaled", batch_key="batch_str", sample_key="sample_key_str"
+    adata_to_correct,
+    layer="preprocessed",
+    batch_key="batch_str",
+    sample_key="sample_key_str",
 )
 
 model = cytovi.CytoVI(
     adata_to_correct, n_hidden=par["n_hidden"], n_layers=par["n_layers"]
 )
 
-
 print("Start training CytoVI model", flush=True)
 model.train()
 

From 165dc9907febf040e62dc5b523eb8806fc428780 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Sat, 22 Nov 2025 12:32:08 +1100
Subject: [PATCH 26/31] increase batch size

---
 src/methods/cytovi/config.vsh.yaml |  4 +--
 src/methods/cytovi/script.py       | 56 ++++++++----------------------
 2 files changed, 16 insertions(+), 44 deletions(-)

diff --git a/src/methods/cytovi/config.vsh.yaml b/src/methods/cytovi/config.vsh.yaml
index 206df6a6..8ba4f536 100644
--- a/src/methods/cytovi/config.vsh.yaml
+++ b/src/methods/cytovi/config.vsh.yaml
@@ -68,14 +68,12 @@ engines:
         packages:
           - anndata>=0.11.0
           - scanpy[skmisc]>=1.10
-          # - scvi-tools==1.4.0
+          - scvi-tools==1.4.0.post1
           - pyyaml
           - requests
           - jsonschema
-          - scikit-learn
         github:
           - openproblems-bio/core#subdirectory=packages/python/openproblems
-          - YosefLab/cytovi-reference-implementation
 
 runners:
   # This platform allows running the component natively
diff --git a/src/methods/cytovi/script.py b/src/methods/cytovi/script.py
index 47ebb941..138c8567 100644
--- a/src/methods/cytovi/script.py
+++ b/src/methods/cytovi/script.py
@@ -1,9 +1,10 @@
+import time
+
 import anndata as ad
-import cytovi
 import numpy as np
+import scvi
 import torch
-
-# from scvi.external import cytovi
+from scvi.external import cytovi
 
 # from sklearn.cluster import KMeans
 # from threadpoolctl import threadpool_limits
@@ -14,8 +15,6 @@
     "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
@@ -23,6 +22,9 @@
 # setting calculation to TF32 to speed up training
 torch.backends.cuda.matmul.allow_tf32 = True
 
+# increase num workers for data loading
+scvi.settings.num_workers = 95
+
 print("Reading and preparing input files", flush=True)
 adata = ad.read_h5ad(par["input"])
 
@@ -39,53 +41,25 @@
     flush=True,
 )
 
-cytovi.CytoVI.setup_anndata(
+cytovi.CYTOVI.setup_anndata(
     adata_to_correct,
     layer="preprocessed",
     batch_key="batch_str",
     sample_key="sample_key_str",
 )
 
-model = cytovi.CytoVI(
+model = cytovi.CYTOVI(
     adata_to_correct, n_hidden=par["n_hidden"], n_layers=par["n_layers"]
 )
 
 print("Start training CytoVI model", flush=True)
-model.train()
-
-# Todo: re-enable subsampling if needed..
-# print("Clustering using k-means with k =", par["n_clusters"], flush=True)
-# # cluster data using Kmeans
-# with threadpool_limits(limits=1):
-#     adata_to_correct.obs["clusters"] = (
-#         KMeans(n_clusters=par["n_clusters"], random_state=0)
-#         .fit_predict(adata_to_correct.layers["scaled"])
-#         .astype(str)
-#     )
-# # concatenate obs so we can use it for subsampling
-# adata_to_correct.obs["sample_cluster"] = (
-#     adata_to_correct.obs["sample"].astype(str) + "_" + adata_to_correct.obs["clusters"]
-# )
-# # subsample cells without replacement
-# print("Subsampling cells", flush=True)
-# subsampled_cells = adata_to_correct.obs.groupby("sample_cluster")[
-#     "sample_cluster"
-# ].apply(lambda x: x.sample(n=round(len(x) * par["subsample_fraction"]), replace=False))
-# # need the cell id included in the subsample
-# subsampled_cells_idx = [x[1] for x in subsampled_cells.index.to_list()]
-
-# adata_subsampled = adata_to_correct[subsampled_cells_idx, :].copy()
-
-# print(
-#     f"Train CytoVI on subsampled data containing {adata_subsampled.shape[0]} cells",
-#     flush=True,
-# )
 
-# cytovi.CYTOVI.setup_anndata(adata_subsampled, layer="scaled", batch_key="batch_str")
-# model = cytovi.CYTOVI(
-#     adata=adata_subsampled, n_hidden=par["n_hidden"], n_layers=par["n_layers"]
-# )
-# model.train()
+start = time.time()
+model.train(
+    batch_size=8192,
+)
+end = time.time()
+print(f"Training took {end - start:.2f} seconds", flush=True)
 
 # get batch corrected data
 print("Correcting data", flush=True)

From 65e48a6cea8385bf895502fc1ec6cb402559dcde Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Sat, 22 Nov 2025 12:42:01 +1100
Subject: [PATCH 27/31] reduce max epochs and train size

---
 src/methods/cytovi/config.vsh.yaml | 10 +++++-----
 src/methods/cytovi/script.py       |  4 ++++
 2 files changed, 9 insertions(+), 5 deletions(-)

diff --git a/src/methods/cytovi/config.vsh.yaml b/src/methods/cytovi/config.vsh.yaml
index 8ba4f536..ff91ddab 100644
--- a/src/methods/cytovi/config.vsh.yaml
+++ b/src/methods/cytovi/config.vsh.yaml
@@ -38,13 +38,13 @@ arguments:
     type: integer
     default: 1
     description: Number of layers.
-  - name: --n_clusters
+  - name: --max_epochs
     type: integer
-    default: 20
-    description: Number of clusters to use for subsampling.
-  - name: --subsample_fraction
+    default: 500
+    description: Number of epochs to train the model.
+  - name: --train_size
     type: double
-    default: 0.5
+    default: 0.7
     description: Fraction of cells to subsample from each cluster for training.
   
 # Resources required to run the component
diff --git a/src/methods/cytovi/script.py b/src/methods/cytovi/script.py
index 138c8567..6d8e0f28 100644
--- a/src/methods/cytovi/script.py
+++ b/src/methods/cytovi/script.py
@@ -15,6 +15,8 @@
     "output": "resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/output_cytovi_split2.h5ad",
     "n_hidden": 128,
     "n_layers": 1,
+    "max_epochs": 500,
+    "train_size": 0.7,
 }
 meta = {"name": "cytovi"}
 ## VIASH END
@@ -57,6 +59,8 @@
 start = time.time()
 model.train(
     batch_size=8192,
+    max_epochs=par["max_epochs"],
+    train_size=par["train_size"],
 )
 end = time.time()
 print(f"Training took {end - start:.2f} seconds", flush=True)

From e94a30ea4b1828fd56f55ae9f3e0599991f9b643 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Sat, 22 Nov 2025 18:51:58 +1100
Subject: [PATCH 28/31] reverting config to default values

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

diff --git a/src/methods/cytovi/config.vsh.yaml b/src/methods/cytovi/config.vsh.yaml
index ff91ddab..2d9a097f 100644
--- a/src/methods/cytovi/config.vsh.yaml
+++ b/src/methods/cytovi/config.vsh.yaml
@@ -40,11 +40,11 @@ arguments:
     description: Number of layers.
   - name: --max_epochs
     type: integer
-    default: 500
+    default: 1000
     description: Number of epochs to train the model.
   - name: --train_size
     type: double
-    default: 0.7
+    default: 0.9
     description: Fraction of cells to subsample from each cluster for training.
   
 # Resources required to run the component
diff --git a/src/methods/cytovi/script.py b/src/methods/cytovi/script.py
index 6d8e0f28..fc53cc55 100644
--- a/src/methods/cytovi/script.py
+++ b/src/methods/cytovi/script.py
@@ -15,8 +15,8 @@
     "output": "resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/output_cytovi_split2.h5ad",
     "n_hidden": 128,
     "n_layers": 1,
-    "max_epochs": 500,
-    "train_size": 0.7,
+    "max_epochs": 1000,
+    "train_size": 0.9,
 }
 meta = {"name": "cytovi"}
 ## VIASH END

From 96d688bc0a7c9d75d23941ae3f8bd8b9afd34da4 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Sun, 23 Nov 2025 21:35:00 +1100
Subject: [PATCH 29/31] update description

---
 .../shuffle_integration/config.vsh.yaml            |  7 +++----
 .../shuffle_integration_by_batch/config.vsh.yaml   |  8 +++-----
 .../config.vsh.yaml                                | 14 +++++++++-----
 3 files changed, 15 insertions(+), 14 deletions(-)

diff --git a/src/control_methods/shuffle_integration/config.vsh.yaml b/src/control_methods/shuffle_integration/config.vsh.yaml
index 23f22899..9ba34a08 100644
--- a/src/control_methods/shuffle_integration/config.vsh.yaml
+++ b/src/control_methods/shuffle_integration/config.vsh.yaml
@@ -3,10 +3,9 @@ name: shuffle_integration
 label: Shuffle Integration
 summary: Randomly shuffle cells in the whole dataset.
 description: |
-  This negative control randomly permutes cell-to-sample (hence batch)
-  assignments while keeping each cell's measured markers unchanged. 
-  This destroys any biological and batch specific structure but preserves marker expression.
-
+  This negative control randomly shuffles all cells in the input data,
+  destroying any biological structure (e.g., sample to cell mapping or batch assignments).
+  
   Purpose:
   - Provide a baseline to verify that integration methods outperform
     random assignment of cells to batches.
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 9bf6a1d1..4f8ba2ea 100644
--- a/src/control_methods/shuffle_integration_by_batch/config.vsh.yaml
+++ b/src/control_methods/shuffle_integration_by_batch/config.vsh.yaml
@@ -3,11 +3,9 @@ name: shuffle_integration_by_batch
 label: Shuffle Integration — within batches
 summary: Randomly reassign cells to any samples within the same batch.
 description: | 
-  This negative-control method randomly permutes cell-to-cell type assignments.
-  Cells remain assigned to their original batch (batch effects preserved).
-  Within each batch, cells are reassigned to random samples, destroying
-  biological/sample-specific structure (e.g., KO vs WT differences).
-
+  This negative-control method randomly shuffles cells within each batch independently, 
+  destroying cell to sample mapping while preserving batch-specific distributions. 
+  
   Purpose:
   - Evaluate whether an integration method preserves differences between samples
     and biological groups while removing batch effects.
diff --git a/src/control_methods/shuffle_integration_by_cell_type/config.vsh.yaml b/src/control_methods/shuffle_integration_by_cell_type/config.vsh.yaml
index 47bc46df..bb9862bd 100644
--- a/src/control_methods/shuffle_integration_by_cell_type/config.vsh.yaml
+++ b/src/control_methods/shuffle_integration_by_cell_type/config.vsh.yaml
@@ -3,17 +3,21 @@ name: shuffle_integration_by_cell_type
 label: Shuffle Integration — within cell type
 summary: Randomly reassign cells to any cell types
 description: | 
-  This negative-control method randomly permutes cell-to-cell type assignments.
-  Cells will be assigned to any cell types, regardless of their original cell type
-  or sample of origin or batch of origin.
+  This negative-control method randomly shuffles cells within each cell type independently, 
+  destroying batch structure while preserving cell type-specific distributions. 
+  This serves as a negative control that maintains biological groupings but 
+  eliminates batch grouping in each cell type.
 
   Purpose:
   - Evaluate whether an integration method preserves differences between cell types
   while removing batch effects.
 
   Example:
-  - A Neutrophil from a KO sample in batch 1 may be reassigned to any cell type
-    (B cell, T cell, Monocyte, etc.) from any sample in any batch.
+  - A Neutrophil in batch 1 from KO sample may be reassigned to a Neutrophil in batch 2
+  KO or WT sample or remain in a KO sample in batch 1 but assigned to different donor,
+  or moved to a WT sample in batch 1 or 2, or remain in the same sample,
+  but it will never be re-assigned to another cell type.
+
 # status: disabled
 resources:
   - type: python_script

From ca35329934029ee02b805e39ca8e5b84ca2e02d3 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Mon, 24 Nov 2025 19:08:12 +1100
Subject: [PATCH 30/31] adding scaling back into cytovi

---
 src/methods/cytovi/script.py | 13 ++++++++++++-
 1 file changed, 12 insertions(+), 1 deletion(-)

diff --git a/src/methods/cytovi/script.py b/src/methods/cytovi/script.py
index fc53cc55..dc8b6e07 100644
--- a/src/methods/cytovi/script.py
+++ b/src/methods/cytovi/script.py
@@ -38,6 +38,17 @@
 
 adata_to_correct = adata[:, markers_to_correct].copy()
 
+print("Scaling data", flush=True)
+
+# scale data. this will add a layer "scaled" to the anndata
+cytovi.scale(
+    adata=adata_to_correct,
+    transformed_layer_key="preprocessed",
+    batch_key="batch_str",
+    scaled_layer_key="scaled",
+    inplace=True,
+)
+
 print(
     f"Train CytoVI on {adata_to_correct.shape[0]} cells",
     flush=True,
@@ -45,7 +56,7 @@
 
 cytovi.CYTOVI.setup_anndata(
     adata_to_correct,
-    layer="preprocessed",
+    layer="scaled",
     batch_key="batch_str",
     sample_key="sample_key_str",
 )

From f9fe9f2fe6134c1b9d167324055b623ea540b262 Mon Sep 17 00:00:00 2001
From: Givanna Putri <givanna.h@gmail.com>
Date: Wed, 10 Dec 2025 12:14:50 +1100
Subject: [PATCH 31/31] add changelog

---
 CHANGELOG.md | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index a39cf059..0515bcd8 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -57,7 +57,6 @@
   * Added CytoNorm with aggregate of samples as controls (`methods/cytonorm_no_controls`).
   * Added parameters to tune CytoNorm.
 
-
 * Added CytoNorm correction to a goal batch (PR #92).
 * Added cyCombine correction to a reference batch (PR #90).
 * Added `metrics/bras` (PR #91).
@@ -66,6 +65,8 @@
 
 * Added processing scripts for CLL dataset (PR #106).
 
+* Added new metric `ratio_inconsistent_peaks` (PR #114).
+
 ## MAJOR CHANGES
 
 * Updated file schema (PR #18): 
@@ -100,6 +101,8 @@
 
 * Fix problems identified during a full run (PR #99).
 
+* Update CytoVI (PR #114).
+
 ## MINOR CHANGES
 
 * Enabled unit tests (PR #2).