Add residual design-effect correction and aggregation calibration parameters

alphaquant/cluster/cluster_utils.py

@@ -134,6 +134,22 @@ def aggregate_node_properties(node, only_use_mainclust, peptide_outlier_filterin
         effective_mode = aggregation_mode.get(node.type, DEFAULT_AGGREGATION_MODE)
     else:
         effective_mode = aggregation_mode
+
+    # median-on-collapse: pruning left one survivor where the parent had several, so the
+    # children were near-duplicates. Take the median of all eligible children's z instead of
+    # the survivor's, recovering within-group evidence without over-counting correlation.
+    if getattr(aqvariables, "MEDIAN_ON_COLLAPSE", False) and len(childs_zfiltered) == 1:
+        all_eligible = [
+            x for x in node.children
+            if x.is_included and (not only_use_mainclust or x.cluster == 0)
+        ]
+        if len(all_eligible) > 1:
+            all_z = get_feature_numpy_array_from_nodes(nodes=all_eligible, feature_name="z_val")
+            all_z = all_z[np.isfinite(all_z)]
+            if len(all_z) > 1:
+                zvals = all_z
+                effective_mode = "median_z"
+
     z_normed = combine_zvalues(zvals, rho=rho, mode=effective_mode)
 
     p_val = transform_znormed_to_pval(z_normed)
@@ -196,9 +212,11 @@ def _select_peptides_around_median_z(peptide_nodes, max_peptides=31):
 def get_selected_nodes_for_zvalcalc(childs, peptide_outlier_filtering, node):
     if peptide_outlier_filtering and node.type == "gene":
         filtered_childs = [x for x in childs if not x.is_outlier_peptide]
-        # Additional restriction: if more than 31 peptides, keep only 31 closest to median z-value
-        if len(filtered_childs) > 31:
-            filtered_childs = _select_peptides_around_median_z(filtered_childs, max_peptides=31)
+        # Additional restriction: cap the number of peptides (closest to median z-value).
+        # Cap is configurable via MAX_PEPTIDES_PER_PROTEIN (disabled by default); None disables it.
+        cap = aqvariables.MAX_PEPTIDES_PER_PROTEIN
+        if cap is not None and len(filtered_childs) > cap:
+            filtered_childs = _select_peptides_around_median_z(filtered_childs, max_peptides=cap)
         return filtered_childs
 
     if node.type == "frgion":

alphaquant/cluster/residual_decorrelation.py

@@ -96,7 +96,9 @@
     ("ms1_isotopes", "base"),
 )
 
-DEFAULT_CUTOFF_GRID = tuple(round(1.0 - 0.1 * k, 2) for k in range(10))
+# 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_MIN_KEEP = 1
 
@@ -427,6 +429,25 @@ def run_level_sweep(
     )
 
 
+def _corr_budget_cols(c1_cols, c2_cols):
+    """Columns to keep for the sibling-correlation estimate, per RESIDUAL_DECORR_CORR_MODE.
+    Returns None (use all) for mode 'off'. For mode 'cap', picks up to RESIDUAL_DECORR_CORR_CAP
+    columns from EACH condition, deterministically evenly spaced, keeping both conditions
+    represented -- so the correlation is never estimated more precisely than ~2*CAP samples."""
+    mode = aqvariables.RESIDUAL_DECORR_CORR_MODE
+    na, nb = len(c1_cols), len(c2_cols)
+    if mode == "cap":
+        cap = aqvariables.RESIDUAL_DECORR_CORR_CAP
+        ka, kb = min(na, cap), min(nb, cap)
+    else:
+        return None
+    def pick(cols, k):
+        k = max(1, min(int(k), len(cols)))
+        idx = sorted(set(np.linspace(0, len(cols) - 1, k).round().astype(int).tolist()))
+        return [cols[i] for i in idx]
+    return pick(list(c1_cols), ka) + pick(list(c2_cols), kb)
+
+
 def attach_lm_residuals(protnodes, df_c1_normed, df_c2_normed, min_n_per_cond=2):
     """Attach per-ion residuals from ``log2(intensity) ~ condition``.
 
@@ -455,6 +476,12 @@ def attach_lm_residuals(protnodes, df_c1_normed, df_c2_normed, min_n_per_cond=2)
     n2_ok = X[c2_cols].notna().sum(axis=1) >= int(min_n_per_cond)
     res.loc[~(n1_ok & n2_ok), :] = np.nan
 
+    # optional correlation-estimation budget: cap the columns used for the sibling-correlation
+    # so its precision (hence pruning aggressiveness) does not blow up at high sample count.
+    keep = _corr_budget_cols(c1_cols, c2_cols)
+    if keep is not None:
+        res = res[keep]
+
     for protnode in protnodes:
         # initialise residuals to None on every node before filling
         for node in PreOrderIter(protnode):
@@ -533,6 +560,7 @@ def apply_residual_decorrelation(
     cutoff_grid=DEFAULT_CUTOFF_GRID,
     aggregation_mode="stouffer_decorrelation",
     null_seed=42,
+    plot_dir=None,
 ):
     """Main entry point: run full residual decorrelation on a list of protein nodes.
 
@@ -566,6 +594,8 @@ def apply_residual_decorrelation(
         for node in PreOrderIter(protnode):
             node.exclude_residual_decorrelation = False
             node.exclude_ptm_fragment_selection = False
+            if aqvariables.RESIDUAL_DEFF_CORRECTION:
+                node.icc_correction = 0.0
 
     # step 1: compute within-condition residuals and attach them to every node
     attach_lm_residuals(protnodes, df_c1_normed, df_c2_normed)
@@ -582,6 +612,7 @@ def apply_residual_decorrelation(
             for pp in parents
         ]
         null_sorted = np.sort(_cross_parent_shuffle_null(mats, rng))
+        n_total = next((m.shape[1] for m in mats if getattr(m, "size", 0)), 0)
         sweep = run_level_sweep(
             parents,
             null_sorted,
@@ -601,12 +632,75 @@ def apply_residual_decorrelation(
         LOGGER.info(msg)
         print(msg, flush=True)
 
+        # design effect on the residual correlation, recorded on the parent so that
+        # aggregation applies deff=1+(n-1)*rho. Excess over the null mean rather than raw
+        # rho: the null absorbs the finite-sample floor of short residual vectors.
+        null_mean = float(np.mean(null_sorted)) if null_sorted.size else 0.0
+        survivor_means = []
+        parents_with_pairs = []
         # mark children that did not survive the chosen cutoff
         for pp in parents:
             survivors = pp.survivors_at(sweep.cutoff, min_keep)
             for keep, child in zip(survivors, pp.child_nodes):
                 if not keep:
                     child.exclude_residual_decorrelation = True
+            if aqvariables.RESIDUAL_DEFF_CORRECTION:
+                resid_rs = _pair_rs_from_C(pp.C, survivors)
+                if resid_rs.size:
+                    survivor_means.append(float(np.mean(resid_rs)))
+                    parents_with_pairs.append(pp)
+                else:
+                    pp.parent_node.icc_correction = 0.0
+        # 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
+        # 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,
+                        "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:
+            if use_raw:
+                # raw pooled ICC: mean pairwise correlation among ALL children (no pruning)
+                source_means = []
+                for pp in parents:
+                    raw_rs = _pair_rs_from_C(pp.C, pp.survivors_at(1.0, min_keep))
+                    if raw_rs.size:
+                        source_means.append(float(np.mean(raw_rs)))
+            else:
+                source_means = survivor_means
+            # clip once on the level mean, not per parent: per-parent means are noisy and
+            # clipping each at zero would rectify that noise into a positive bias.
+            level_excess = (max(0.0, float(np.mean(source_means)) - null_mean)
+                            if source_means else 0.0)
+            for pp in parents_with_pairs:
+                pp.parent_node.icc_correction = level_excess
+            LOGGER.info(
+                "deff %s->%s: mean %s rho=%.4f  null mean=%.4f  "
+                "LEVEL excess rho=%.4f (parents=%d)",
+                parent_level, child_level, "RAW" if use_raw else "survivor",
+                float(np.mean(source_means)) if source_means else 0.0,
+                null_mean, level_excess, len(source_means),
+            )
+
+    # optional: save the per-level distribution diagnostics (before/after/null CDFs
+    # + cutoff sweep trace) using AlphaQuant's own plotting.
+    if plot_dir is not None:
+        import os
+        os.makedirs(plot_dir, exist_ok=True)
+        for sweep in level_results:
+            try:
+                fig = plot_level_sweep_diagnostics(sweep)
+                fig.savefig(os.path.join(plot_dir, f"decorr_{sweep.level[0]}__{sweep.level[1]}.png"),
+                            dpi=120)
+                plt.close(fig)
+            except Exception as exc:
+                LOGGER.warning("could not save decorrelation plot for %s: %s", sweep.level, exc)
 
     # step 3 (optional): apply PTM fragment selection on top of decorrelation exclusions
     if aqvariables.PTM_FRAGMENT_SELECTION:

alphaquant/config/variables.py

@@ -13,6 +13,18 @@
 CLASSIC_FRAGMENT_OUTLIER_FILTERING = False
 ICC_NULL_PVAL_THRESHOLD = 0.1
 NUM_BG_CONTEXTS = 10
+
+MEDIAN_ON_COLLAPSE = True
+
+
+RESIDUAL_DEFF_CORRECTION = True
+
+RESIDUAL_DEFF_SMALLN_TOTAL = 7
+
+RESIDUAL_DECORR_CORR_MODE = "cap"
+RESIDUAL_DECORR_CORR_CAP = 10
+
+MAX_PEPTIDES_PER_PROTEIN = None
 CONDITION_PAIR_SEPARATOR = "_VS_"
 
 #prefixes for the different ion types
@@ -51,6 +63,34 @@ def set_peptide_outlier_filtering(peptide_outlier_filtering):
     global PEPTIDE_OUTLIER_FILTERING
     PEPTIDE_OUTLIER_FILTERING = peptide_outlier_filtering
 
+def set_median_on_collapse(median_on_collapse):
+    global MEDIAN_ON_COLLAPSE
+    MEDIAN_ON_COLLAPSE = bool(median_on_collapse)
+
+def set_residual_deff_correction(residual_deff_correction):
+    global RESIDUAL_DEFF_CORRECTION
+    RESIDUAL_DEFF_CORRECTION = bool(residual_deff_correction)
+
+
+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
+
+
+def set_residual_decorr_corr_mode(residual_decorr_corr_mode):
+    global RESIDUAL_DECORR_CORR_MODE
+    RESIDUAL_DECORR_CORR_MODE = str(residual_decorr_corr_mode) if residual_decorr_corr_mode else "cap"
+
+
+def set_residual_decorr_corr_cap(residual_decorr_corr_cap):
+    global RESIDUAL_DECORR_CORR_CAP
+    RESIDUAL_DECORR_CORR_CAP = int(residual_decorr_corr_cap) if residual_decorr_corr_cap else 10
+
+def set_max_peptides_per_protein(max_peptides_per_protein):
+    global MAX_PEPTIDES_PER_PROTEIN
+    MAX_PEPTIDES_PER_PROTEIN = (int(max_peptides_per_protein)
+                                if max_peptides_per_protein is not None else None)
+
 def set_outlier_correction_factor(outlier_correction_factor):
     global OUTLIER_CORRECTION_FACTOR
     OUTLIER_CORRECTION_FACTOR = float(outlier_correction_factor)

alphaquant/diffquant/condpair_analysis.py

@@ -67,7 +67,8 @@ def analyze_condpair(*,runconfig, condpair):
         return
 
     df_c1_normed, df_c2_normed = aqnorm.normalize_if_specified(df_c1 = df_c1, df_c2 = df_c2, c1_samples = c1_samples, c2_samples = c2_samples, normalize_within_conds = runconfig.normalize, normalize_between_conds = runconfig.normalize,
-    runtime_plots = runconfig.runtime_plots, protein_subset_for_normalization_file=runconfig.protein_subset_for_normalization_file, pep2prot = pep2prot)#, "./test_data/normed_intensities.tsv")
+    runtime_plots = runconfig.runtime_plots, protein_subset_for_normalization_file=runconfig.protein_subset_for_normalization_file, pep2prot = pep2prot,
+    median_normalization = getattr(runconfig, 'median_normalization', False))#, "./test_data/normed_intensities.tsv")
 
     summarization_nodes = getattr(runconfig, 'summarization_nodes', [])
     if summarization_nodes:
@@ -158,7 +159,12 @@ def analyze_condpair(*,runconfig, condpair):
             df_c2_normed,
             tolerance=getattr(runconfig, "residual_decorrelation_tolerance", 0.10),
             min_keep=getattr(runconfig, "residual_decorrelation_min_keep", 1),
+            cutoff_grid=(getattr(runconfig, "residual_decorrelation_cutoff_grid", None)
+                         or aq_clust_resid.DEFAULT_CUTOFF_GRID),
             aggregation_mode=runconfig.aggregation_mode,
+            plot_dir=(os.path.join(runconfig.results_dir,
+                                   f"{aqutils.get_condpairname(condpair)}_residual_decorrelation_plots")
+                      if getattr(runconfig, "runtime_plots", False) else None),
         )
     if len(prot2missingval_diffions.keys())>0:
         LOGGER.info(f"start analysis of proteins w. completely missing values")

alphaquant/norm/normalization.py

@@ -257,7 +257,7 @@ def mode_normalization(x):
 import numpy as np
 from scipy import stats
 
-def get_betweencond_shift(df_c1_normed, df_c2_normed, enfore_median = False):
+def get_betweencond_shift(df_c1_normed, df_c2_normed, median_normalization = False):
 
     both_idx = df_c1_normed.index.intersection(df_c2_normed.index)
     df1 = df_c1_normed.loc[both_idx]
@@ -269,7 +269,7 @@ def get_betweencond_shift(df_c1_normed, df_c2_normed, enfore_median = False):
 
     diff_fcs = df1[col1].to_numpy() - df2[col2].to_numpy()
     median = np.nanmedian(diff_fcs)
-    if enfore_median:
+    if median_normalization:
         return -median
 
     if len(diff_fcs)<100:
@@ -288,7 +288,7 @@ def get_betweencond_shift(df_c1_normed, df_c2_normed, enfore_median = False):
 # Cell
 import pandas as pd
 
-def normalize_if_specified(df_c1, df_c2, c1_samples, c2_samples, normalize_within_conds = True, normalize_between_conds = True, runtime_plots = True, protein_subset_for_normalization_file = None, pep2prot =None):
+def normalize_if_specified(df_c1, df_c2, c1_samples, c2_samples, normalize_within_conds = True, normalize_between_conds = True, runtime_plots = True, protein_subset_for_normalization_file = None, pep2prot =None, median_normalization = False):
 
     if normalize_within_conds:
         df_c1 = normalize_within_cond(df_c=df_c1, samples_c= c1_samples)
@@ -299,15 +299,15 @@ def normalize_if_specified(df_c1, df_c2, c1_samples, c2_samples, normalize_withi
         aq_plot_pairwise.plot_withincond_normalization(df_c1, df_c2)
 
     if normalize_between_conds:
-        df_c1, df_c2 = get_normalized_dfs_between_conditions(df_c1, df_c2, protein_subset_for_normalization_file, pep2prot,runtime_plots = runtime_plots)
+        df_c1, df_c2 = get_normalized_dfs_between_conditions(df_c1, df_c2, protein_subset_for_normalization_file, pep2prot,runtime_plots = runtime_plots, median_normalization = median_normalization)
         LOGGER.info("normalized between conditions")
 
     return df_c1, df_c2
 
 
 
-def get_normalized_dfs_between_conditions(df_c1, df_c2, protein_subset_for_normalization_file, pep2prot,runtime_plots):
-    shift_between_cond = prepare_tables_and_get_betweencond_shift(df_c1, df_c2, protein_subset_for_normalization_file, pep2prot)
+def get_normalized_dfs_between_conditions(df_c1, df_c2, protein_subset_for_normalization_file, pep2prot,runtime_plots, median_normalization = False):
+    shift_between_cond = prepare_tables_and_get_betweencond_shift(df_c1, df_c2, protein_subset_for_normalization_file, pep2prot, median_normalization = median_normalization)
 
     LOGGER.info(f"shift comparison by {shift_between_cond}")
     df_c2 = df_c2-shift_between_cond
@@ -322,12 +322,12 @@ def normalize_within_cond(df_c, samples_c):
     df_c_normed = pd.DataFrame(apply_sampleshifts(df_c.to_numpy().T, sample2shift).T, index = df_c.index, columns = samples_c)
     return df_c_normed
 
-def prepare_tables_and_get_betweencond_shift(df_c1, df_c2, protein_subset_for_normalization_file, pep2prot):
+def prepare_tables_and_get_betweencond_shift(df_c1, df_c2, protein_subset_for_normalization_file, pep2prot, median_normalization = False):
     specified_protein_subset = read_specified_protein_subset_if_given(protein_subset_for_normalization_file)
     prepared1 = prepare_table_for_betweencond_shift(df_c1, specified_protein_subset, pep2prot)
     prepared2 = prepare_table_for_betweencond_shift(df_c2, specified_protein_subset, pep2prot)
-    enforce_median = protein_subset_for_normalization_file is not None
-    return get_betweencond_shift(prepared1, prepared2, enforce_median)
+    median_normalization = median_normalization or (protein_subset_for_normalization_file is not None)
+    return get_betweencond_shift(prepared1, prepared2, median_normalization)
 
 def read_specified_protein_subset_if_given(specified_protein_subset_file):
     if specified_protein_subset_file is not None: