From a686d4cc999cd81c6d0f0327b93a77bcbf2d58ba Mon Sep 17 00:00:00 2001 From: ammarcsj <70114795+ammarcsj@users.noreply.github.com> Date: Fri, 21 Aug 2026 14:44:11 +0200 Subject: [PATCH 1/2] Never let peptide outlier filtering increase protein significance Peptide outlier candidates are flagged on their two-sided p-value, which carries no direction, while the protein-level aggregation in sum_and_re_scale_zvalues sums *signed* z-values. Discarding a peptide whose z-value opposes the protein's direction therefore raises |sum(z)| and lowers n, making the protein more significant - the opposite of what the filter is for, and contrary to what its docstring claimed. apply_peptide_outlier_filtering now delegates each protein to _filter_and_aggregate_protein, which re-aggregates and, if |z| grew, un-flags every peptide and restores the unfiltered aggregate. Removing peptides can now only ever lower a protein's significance. Verified by replaying both Spectronaut species-mixture benchmarks from their exported iontrees (the replay reproduces the pipeline's stored per-peptide is_outlier_peptide flags and gene p-values exactly): - LargeFC: 1005 of 8595 proteins had their significance inflated, worst case 46x (p 0.945 -> 0.021). All were unregulated background. The sets of significant proteins are identical before and after this change at 5%, 1% and 0.1% FDR, and the empirical FDR is unchanged at 3.30%. - SmallFC: 189 of 10634 affected, worst case 37222x. 7 of 1337 significant proteins are lost (0.5% at 5% FDR), all spiked-in true positives that reached significance through the removal of a single peptide. Empirical FDR unchanged at 1.50%. - No protein gains significance at any threshold in either dataset, and the mechanism never rendered an unregulated protein significant. A sign-aware variant (drop only candidates concordant with the protein) was tested and rejected: it is worse than no fix at all, because the discordant candidate that is then retained keeps inflating the aggregate once its concordant partner has been removed. Co-Authored-By: Claude Opus 5 (1M context) --- alphaquant/cluster/outlier_filtering.py | 29 ++++++++++++++++++++++++- 1 file changed, 28 insertions(+), 1 deletion(-) diff --git a/alphaquant/cluster/outlier_filtering.py b/alphaquant/cluster/outlier_filtering.py index d6e8f1d8..12cc6fd1 100644 --- a/alphaquant/cluster/outlier_filtering.py +++ b/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 From aa50ad92e5f7c4b1b3afd5a2a7f6f9c6ee3b66e2 Mon Sep 17 00:00:00 2001 From: ammarcsj <70114795+ammarcsj@users.noreply.github.com> Date: Thu, 1 Oct 2026 07:57:26 +0200 Subject: [PATCH 2/2] Decouple deff gate from pruning tolerance; default tolerance 0.08 residual_decorrelation_tolerance default 0.10 -> 0.08. New residual_deff_gate_tolerance (default 0.05): the deff gate no longer uses the pruning tolerance. None restores the old coupled behaviour. Co-Authored-By: Claude Opus 5.5 (1M context) --- alphaquant/cluster/residual_decorrelation.py | 13 +++++++++---- alphaquant/config/variables.py | 9 +++++++++ alphaquant/run_pipeline.py | 11 ++++++++--- 3 files changed, 26 insertions(+), 7 deletions(-) diff --git a/alphaquant/cluster/residual_decorrelation.py b/alphaquant/cluster/residual_decorrelation.py index 0b2dee6d..aed27968 100644 --- a/alphaquant/cluster/residual_decorrelation.py +++ b/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: diff --git a/alphaquant/config/variables.py b/alphaquant/config/variables.py index d0249dfd..52e0c67c 100644 --- a/alphaquant/config/variables.py +++ b/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 diff --git a/alphaquant/run_pipeline.py b/alphaquant/run_pipeline.py index ff99242d..71f5c59f 100644 --- a/alphaquant/run_pipeline.py +++ b/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)