Peptide outlier guard and decoupled deff gate

alphaquant/cluster/outlier_filtering.py

@@ -6,8 +6,34 @@
 def apply_peptide_outlier_filtering(protnodes: list[anytree.Node], aggregation_mode="stouffer_decorrelation"):
     regulation_score = calculate_regulation_score(protnodes)
     for protnode in protnodes:
-        _determine_and_annotate_outlier_status_of_peptides(protnode, regulation_score)
+        _filter_and_aggregate_protein(protnode, regulation_score, aggregation_mode)
 
+
+def _filter_and_aggregate_protein(protnode, regulation_score, aggregation_mode):
+    """Flags the outlier peptides of one protein and re-aggregates it.
+
+    Peptides are flagged on their two-sided p-value, which carries no direction, while the
+    protein-level aggregation (``cluster_utils.sum_and_re_scale_zvalues``) sums *signed*
+    z-values. Discarding a peptide that opposes the protein's direction therefore raises
+    ``|sum(z)|`` and lowers ``n``, which would make the protein *more* significant - the
+    opposite of what this filter is for. Whenever that happens the flags are dropped and
+    the unfiltered aggregate is restored, so removing peptides can only ever lower a
+    protein's significance.
+
+    Args:
+        protnode: Protein node with peptide children. On entry its ``z_val``/``p_val``
+            hold the aggregate over all peptides, i.e. the unfiltered state.
+        regulation_score: Float between 0 and 1 representing overall regulation context
+        aggregation_mode: Strategy for combining child z-values, see
+            ``cluster_utils.aggregate_node_properties``
+    """
+    z_val_unfiltered = protnode.z_val
+
+    _determine_and_annotate_outlier_status_of_peptides(protnode, regulation_score)
+    aqcluster_utils.aggregate_node_properties(protnode, only_use_mainclust=True, peptide_outlier_filtering=True, aggregation_mode=aggregation_mode)
+
+    if abs(protnode.z_val) > abs(z_val_unfiltered):
+        _annotate_peptides(protnode.children, is_outlier=False)
         aqcluster_utils.aggregate_node_properties(protnode, only_use_mainclust=True, peptide_outlier_filtering=True, aggregation_mode=aggregation_mode)
 
 
@@ -39,6 +65,7 @@ def _determine_and_annotate_outlier_status_of_peptides(protnode, regulation_scor
     """
     We look at the distribution of p-values of the peptides of a protein and focus on a particular class of proteins that are potentially dominated by outliers. For these protein 1) the majority of peptides is not significant 2) there are a few peptides are more strongly significant than the majority.
     Depending on the the overall context of the experiment which is quantified by the regulation score (low regulation score means few weakly regulated proteins, high regulation score means many strongly regulated proteins), we are more or less tolerant to the outliers.
+    Flagging is based on the two-sided p-value only and is thus blind to the direction of the effect; _filter_and_aggregate_protein discards the result for proteins where that would have increased significance.
 
     Args:
         protnode: Protein node with peptide children

alphaquant/cluster/residual_decorrelation.py

@@ -99,7 +99,7 @@
 # 1.0 down to -1.0 in steps of 0.1. The negative part is only reached when no cutoff
 # meets the tolerance, in which case the tightest one prunes down to min_keep.
 DEFAULT_CUTOFF_GRID = tuple(round(1.0 - 0.1 * k, 2) for k in range(21))
-DEFAULT_TOLERANCE = 0.10
+DEFAULT_TOLERANCE = 0.08
 DEFAULT_MIN_KEEP = 1
 
 
@@ -654,14 +654,19 @@ def apply_residual_decorrelation(
         # gene->seq only: that level carries the protein-level random effect shared across a
         # protein's peptides, which the ion-variance model does not capture. Gated on d_before
         # so that peptides no more correlated than the null stay a no-op.
-        gate_open = sweep.d_before > tolerance
+        # independent of the pruning tolerance: sharing one value meant a loose tolerance
+        # switched the correction off where it was still needed. None = old coupled behaviour.
+        gate_tol = aqvariables.RESIDUAL_DEFF_GATE_TOLERANCE
+        if gate_tol is None:
+            gate_tol = tolerance
+        gate_open = sweep.d_before > gate_tol
         # below this sample count the per-dataset correlation is unmeasurable, so survivor rho
         # reads ~0 while the correlation still leaks into the Stouffer sum: use raw rho instead.
         smalln = aqvariables.RESIDUAL_DEFF_SMALLN_TOTAL
         use_raw = bool(smalln) and n_total <= smalln
         if aqvariables.RESIDUAL_DEFF_CORRECTION and parent_level == "gene":
-            LOGGER.info("deff gate gene->seq: d_before=%.4f tolerance=%.4f n_total=%d -> %s%s",
-                        sweep.d_before, tolerance, n_total,
+            LOGGER.info("deff gate gene->seq: d_before=%.4f gate_tol=%.4f n_total=%d -> %s%s",
+                        sweep.d_before, gate_tol, n_total,
                         "OPEN" if gate_open else "CLOSED (deff off)",
                         " [raw small-n ICC]" if use_raw else "")
         if aqvariables.RESIDUAL_DEFF_CORRECTION and parent_level == "gene" and gate_open:

alphaquant/config/variables.py

@@ -20,6 +20,8 @@
 RESIDUAL_DEFF_CORRECTION = True
 
 RESIDUAL_DEFF_SMALLN_TOTAL = 7
+# deff gate threshold, deliberately independent of the pruning tolerance
+RESIDUAL_DEFF_GATE_TOLERANCE = 0.05
 
 RESIDUAL_DECORR_CORR_MODE = "cap"
 RESIDUAL_DECORR_CORR_CAP = 10
@@ -72,6 +74,13 @@ def set_residual_deff_correction(residual_deff_correction):
     RESIDUAL_DEFF_CORRECTION = bool(residual_deff_correction)
 
 
+def set_residual_deff_gate_tolerance(residual_deff_gate_tolerance):
+    """None means: fall back to the pruning tolerance (the pre-decoupling behaviour)."""
+    global RESIDUAL_DEFF_GATE_TOLERANCE
+    RESIDUAL_DEFF_GATE_TOLERANCE = (None if residual_deff_gate_tolerance is None
+                                    else float(residual_deff_gate_tolerance))
+
+
 def set_residual_deff_smalln_total(residual_deff_smalln_total):
     global RESIDUAL_DEFF_SMALLN_TOTAL
     RESIDUAL_DEFF_SMALLN_TOTAL = int(residual_deff_smalln_total) if residual_deff_smalln_total else 0

alphaquant/run_pipeline.py

@@ -57,11 +57,12 @@ def run_pipeline(input_file: str,
                 cluster_threshold_fcfc: float = 0,
                 fcdiff_cutoff_clustermerge = 0.5,
                 use_ml: bool = True,
-                residual_decorrelation_tolerance: float = 0.10,
+                residual_decorrelation_tolerance: float = 0.08,
                 residual_decorrelation_min_keep: int = 1,
                 residual_decorrelation_cutoff_grid: Optional[List[float]] = None,
                 median_on_collapse: bool = True,
                 residual_deff_correction: bool = True,
+                residual_deff_gate_tolerance: Optional[float] = 0.05,
                 residual_deff_smalln_total: int = 7,
                 residual_decorr_corr_mode: str = "cap",
                 residual_decorr_corr_cap: int = 10,
@@ -131,7 +132,7 @@ def run_pipeline(input_file: str,
     cluster_threshold_fcfc (float): Fold change threshold for clustering. Defaults to 0.
     fcdiff_cutoff_clustermerge (float): Fold change difference cutoff for merging peptide clusters. Defaults to 0.5.
     use_ml (bool): Enable machine learning analysis. Defaults to True.
-    residual_decorrelation_tolerance (float): Maximum allowed one-sided excess-CDF distance between corrected and null sibling-correlation distributions. Defaults to 0.10.
+    residual_decorrelation_tolerance (float): Maximum allowed one-sided excess-CDF distance between corrected and null sibling-correlation distributions. Defaults to 0.08.
     residual_decorrelation_min_keep (int): Minimum number of children to retain per parent during residual decorrelation pruning. Defaults to 1.
     residual_decorrelation_cutoff_grid (list[float] | None): Correlation cutoffs
         scanned (loose->tight) during residual-decorrelation pruning. None (default)
@@ -157,13 +158,16 @@ def run_pipeline(input_file: str,
         (no droppable subset). Two restrictions make it a no-op on well-calibrated
         data: it is applied ONLY at the between-peptide (gene->seq) level, which is
         where the shared protein-level random effect lives, and only when that
-        level's pre-pruning distance exceeded residual_decorrelation_tolerance. If
+        level's pre-pruning distance exceeded residual_deff_gate_tolerance. If
         peptides are no more correlated than the shuffle null to begin with, the
         gate stays closed, rho remains 0.0 and aggregation is unchanged. When the
         gate opens, a single level-wide excess rho (mean survivor correlation minus
         the shuffle-null mean, clipped at zero once on the level mean rather than
         per parent, to avoid rectifying per-parent noise into a positive bias) is
         applied to every parent at that level. Defaults to True.
+    residual_deff_gate_tolerance (float | None): Excess-CDF distance above which the deff
+        correction switches on, independent of residual_decorrelation_tolerance. None ties
+        it back to that tolerance (the old coupled behaviour). Defaults to 0.05.
     residual_deff_smalln_total (int): Total-sample threshold below which
         residual_deff_correction sources its rho from the RAW, pre-pruning peptide
         correlations (cutoff 1.0) pooled across the dataset instead of from the
@@ -339,6 +343,7 @@ def run_pipeline(input_file: str,
     aqvariables.set_peptide_outlier_filtering(peptide_outlier_filtering)
     aqvariables.set_median_on_collapse(median_on_collapse)
     aqvariables.set_residual_deff_correction(residual_deff_correction)
+    aqvariables.set_residual_deff_gate_tolerance(residual_deff_gate_tolerance)
     aqvariables.set_residual_deff_smalln_total(residual_deff_smalln_total)
     aqvariables.set_residual_decorr_corr_mode(residual_decorr_corr_mode)
     aqvariables.set_residual_decorr_corr_cap(residual_decorr_corr_cap)