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

Author: LuLeom

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

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: