Describe the peak calling in ratio_consistent_peaks as implemented (#132)

* describe the peak calling in ratio_consistent_peaks as implemented

* Document the standardisation step and the actual height/prominence thresholds
* Drop persistent_peak_count() and the scikit-tda dependency, it was never called
* Only import matplotlib and seaborn when actually plotting
* Say case 2 where the code means case 2

* update changelog

---------

Co-authored-by: Givanna Putri <givanna.h@gmail.com>

Author: rcannood

CHANGELOG.md

@@ -241,6 +241,9 @@
 * Update flowsom mapping similarity so we subset to just markers to correct, and lisi to remove control samples
   and unlabelled cells (PR #119).
 
+* Update the `ratio_consistent_peaks` description to match the implemented peak calling
+  (standardisation step, `height` 0.1 and `prominence` 0.01), and remove the unused
+  `persistent_peak_count()` helper together with the `scikit-tda` dependency (PR #132).
 * Point `scripts/run_benchmark/wehi_hpc/run_full_hpc.sh` at `build/main` instead of
   `build/update_ilisi`, and label the seqera full run as `full` instead of `test_subset` (PR #129).
 
@@ -249,7 +252,6 @@
 
 * Clean up stale mock parameters and dead code (PR #133).
 
-
 * Fix bug in `average_batch_r2` where the R2 was computed on all cell types of a donor at once
   instead of on each cell type separately (PR #127).
 

src/metrics/ratio_consistent_peaks/config.vsh.yaml

@@ -20,9 +20,13 @@ info:
       description: |
         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 (cell type) marker expression profiles are first smoothed using kernel density estimation (KDE) (`scipy.stats.gaussian_kde`),
-        and then peaks are 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 (cell type) marker expression values of the two splits are first standardised together,
+        using the pooled mean and standard deviation of the two splits, so that the peak calling
+        thresholds mean the same thing for every marker.
+        They are then smoothed using kernel density estimation (KDE) (`scipy.stats.gaussian_kde`,
+        evaluated on a grid of 100 points), and peaks are identified using the
+        `scipy.signal.find_peaks` function.
+        For peak calling, the `height` parameter is set to 0.1 and the `prominence` parameter is set to 0.01.
         
         Case Definitions:
         - Case 1: Consistent peaks in both unintegrated AND integrated data (ideal outcome).
@@ -75,10 +79,6 @@ engines:
   # Specifications for the Docker image for this component.
   - type: docker
     image: openproblems/base_python:1
-    setup:
-      - type: python
-        packages:
-          - scikit-tda
 
 runners:
   # This platform allows running the component natively

src/metrics/ratio_consistent_peaks/helper.py

@@ -1,7 +1,4 @@
-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
 
@@ -65,6 +62,10 @@ def get_kde_density(expression_array, return_xgrid=False, plot=False):
         x_grid = np.concatenate([[min_val - step], x_grid])
 
     if plot:
+        # only needed when debugging locally, so don't import at module level
+        import matplotlib.pyplot as plt
+        import seaborn as sns
+
         fig, ax = plt.subplots()
         sns.scatterplot(x=x_grid, y=density, ax=ax)
         ax.set_title("KDE Density Estimation")
@@ -98,40 +99,3 @@ def call_peaks(density):
     num_peaks = len(peaks)
 
     return num_peaks
-
-
-def persistent_peak_count(ys, persistence_cutoff=0.08):
-    """
-    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
-            Default to 0.08 to biased towards strong peaks but not overly.
-
-    Returns:
-        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]
-    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

src/metrics/ratio_consistent_peaks/script.py

@@ -64,8 +64,8 @@
 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
+# case 2 = consistent peaks in unintegrated but inconsistent in integrated
+# case NGT = unintegrated is inconsistent, recorded but excluded from the ratio
 n_case1 = 0
 n_case2 = 0
 
@@ -180,7 +180,7 @@
                 np.std(s2_unscaled) == 0 and np.std(u_s2_unscaled) != 0
             ):
                 print(
-                    f"WARNING: Marker {marker}, donor {donor}, cell type {celltype}: has no variance either before or after integration. Automatic to case 3.",
+                    f"WARNING: Marker {marker}, donor {donor}, cell type {celltype}: has no variance after integration, but does before. Automatic to case 2.",
                     flush=True,
                 )
                 n_case2 += 1