Skip to content

10_features.md

Latest commit

 

History

History
829 lines (720 loc) · 49.7 KB

10_features.md

File metadata and controls

829 lines (720 loc) · 49.7 KB
Raw

features (Stage E): the feature battery + PIN

Part of the MuMDIA developer documentation (see docs/README.md).

Purpose

Stage E (mumdia features) turns one extracted PSM (apex-level identification plus its per-fragment chromatograms) into a fixed, named, versioned feature vector for the semi-supervised rescorer. It computes a config-selected feature set, a scalar prelim_score used by compete, a Percolator PIN, and a blake3 schema id that pins the ordered column list so the classifier is never trained or applied under a mismatched feature layout.

The stage is pure per-PSM: each row's features are a function of that row's own inputs only, with two cross-row exceptions computed up front (charge-state corroboration grouped by peptidoform, and a global elution half-width learned from the confident seed set). This makes the heavy per-PSM work embarrassingly parallel (rayon) while keeping byte-identical output to a serial run.

Everything is driven by FeaturesConfig (mumdia-core/src/config.rs:578). The three feature sets are Minimal (14), Rich (44), and Extended (381). The default is Minimal; the tuned --profile dia preset and the DIA-NN-library recipe use Extended.

Files

path role
rust/mumdia/crates/mumdia/src/stages/features.rs stage entry point, run, Evidence, evidence_from (and its test-only wrapper build_evidence), fragment_features, boundary detection, prelim_score, PIN, schema hash, the Minimal/Rich column lists, and the Extended family registry
rust/mumdia/crates/mumdia/src/stages/features/similarity.rs Extended family: observed-vs-library intensity agreement kernels (64 names)
rust/mumdia/crates/mumdia/src/stages/features/entropy.rs Extended family: spectral-entropy / information-divergence (18 names); also exports the gate kernel spectral_entropy_similarity_sqrt
rust/mumdia/crates/mumdia/src/stages/features/coelution.rs Extended family: fragment-vs-reference and pairwise co-elution + cross-correlation (38 names)
rust/mumdia/crates/mumdia/src/stages/features/interference.rs Extended family: co-isolation / chimera detection, interference removal, rank decomposition (26 names)
rust/mumdia/crates/mumdia/src/stages/features/chromatographic.rs Extended family: peak-shape descriptors of the reference profile and fragments (43 names)
rust/mumdia/crates/mumdia/src/stages/features/mass_accuracy.rs Extended family: fragment ppm distribution + positive mass evidence (17 names)
rust/mumdia/crates/mumdia/src/stages/features/ion_series.rs Extended family: b/y series coverage, runs, complementarity, per-series similarity (34 names)
rust/mumdia/crates/mumdia/src/stages/features/ms1.rs Extended family: precursor isotope-envelope agreement + MS1/MS2 XIC co-elution (25 names)
rust/mumdia/crates/mumdia/src/stages/features/rt.rs Extended family: RT-agreement variants (13 names)
rust/mumdia/crates/mumdia/src/stages/features/novel.rs Extended family: seed corroboration + precursor/charge metadata (12 names)
rust/mumdia/crates/mumdia/src/stages/features/nonzero.rs Extended family: zero-ignoring apex/co-elution variants (12 names)
rust/mumdia/crates/mumdia/src/stages/features/order_consistency.rs Extended family: prediction-free MS2-XIC rank stability (8 names)
rust/mumdia/crates/mumdia/src/stages/features/peak_scans.rs Extended family: peak-scan count / window-degeneracy indicator (2 names)
rust/mumdia/crates/mumdia/src/stages/features/apex_dispersion.rs Extended family: fragment apex dispersion + consensus peak shape (13 names)
rust/mumdia/crates/mumdia/src/stages/features/mass_uncertainty.rs Extended family: fragment mass-error dispersion + evidence breadth (10 names)
rust/mumdia/crates/mumdia/src/stats.rs shared pearson/cosine/spectral_angle kernels used by every family
rust/mumdia/crates/mumdia-core/src/config.rs:578 FeaturesConfig

Inputs and outputs

FeaturesParams (features.rs:537) names the paths: psms (extracted PSMs), chromatograms, optional seed (seed PSMs), out (features Parquet), out_pin (PIN), cfg, config_hash.

Consumed: psms_extracted (features.rs:550-588)

Columns read (getter -> column): candidate_id (u32), apex_rt (f64), apex_intensity (f32), n_matched_fragments (i32), n_predicted_fragments (i32, optional; defaults to 6 when absent), coelution_run (i32), rt_pred_cal (f64), charge (i32), label (str), base_peptide_id (u32), peptidoform (str), protein (str), precursor_mz (f64). Optional soft- competition columns default to 0.0 when absent: contested_frac, contested_count_frac, apportioned_frac. Optional MS1 apex isotope columns (opt_f64, default None): ms1_isom1, ms1_mono, ms1_iso1, ms1_iso2.

Consumed: chromatograms (features.rs:591-618)

candidate_id (u32), frag_name (str), frag_mz (f64), frag_obs_mz (f64, optional; falls back to frag_mz), predicted_intensity (f32), rt (list), intensity (list). Rows whose frag_name starts with ms1_ are routed to a separate ms1x map and fed to the MS1 XIC evidence; all others are the fragment chromatograms.

Either chromatogram layout is read (docs/15_data_dictionary.md, "Layout v2"). ChromStream passes every v2 row through a chromatograms::Decoder before the row is stored, so a chunk holds exactly the v1 rows and every feature is the same. The main pass starts each chunk at a candidate's first row, where the decoder needs nothing earlier. The confident-bounds pass reads sub-chunks on the absolute chunk_rows grid, which can start inside a row group and inside a candidate. There ChromStream::open_at opens the table at that row group's first row and follows the rows before the sub-chunk through the decoder without keeping them, so the samples are v1's in order, bit for bit.

Consumed: seed PSMs (optional, features.rs:620-652)

candidate_id (u32), score (f64), spectrum_q (f64), label (str). Builds the candidate_id -> seed score and -> identified flag maps, and the confident-target set (spectrum_q <= 0.01 and label == "target") used for bound_from_confident.

Produced: features Parquet (features.rs:912-931)

Bookkeeping columns, in order: candidate_id (u32), label (str), base_peptide_id (u32), peptidoform (str), protein (str), apex_rt (f64), elution_lo (f64), elution_hi (f64), precursor_mz (f64), prelim_score (f64). elution_lo/elution_hi are the RT bounds the stage actually used (emitted so downstream and plotting read them rather than re-derive). Then one column per name in active_features(cfg.set), in order.

Each feature column is Float32, holding the f64 value the stage computed narrowed by v as f32, except the five F64_FEATURE_COLUMNS (charge, n_matched_fragments, unique_fragment_count, peak_contested_frac, contested_frac), which stay Float64 because compete or rescore reads them as f64 before narrowing. Every classifier narrows every feature to f32 the same way, so the scored outputs do not change. docs/15_data_dictionary.md ("features.parquet") has the storage contract and the version-1 layout.

Produced: companion + PIN + report

  • <out>.schema.json (features.rs:936): FeatureSchema { feature_columns, schema_id }, read back by FeatureSchema::read (features.rs:245).
  • <out_pin>: Percolator PIN, header SpecId Label ScanNr ExpMass CalcMass <features...> Peptide Proteins (write_pin, features.rs:1537).
  • <out>.report.json (ArtifactReport, features.rs:969): logical name features, schema version 2 (version 1 stored every feature column as Float64; compete and rescore still read it), stats carrying feature_schema_id, n_features, and set; params records set and coelution_corr_threshold (features.rs:976); content_hash is the blake3 of the features Parquet; model_identity is None.

How it works

Control flow of run (features.rs:548):

  1. Read the PSM table and pull the scalar columns (features.rs:550-588).
  2. Read chromatograms and group ChromRows by candidate_id into chrom, splitting off ms1_* rows into ms1x (features.rs:591-618). A ChromRow (features.rs:279) holds frag_name, frag_mz, frag_obs_mz, pred_int, and the rt/inten vectors.
  3. Build the seed maps and confident-target set from the optional seed table (features.rs:620-652).
  4. If bound_from_confident is set, learn a global elution half-width (see "bound_from_confident" below), producing Option<(L, R)> in seconds (features.rs:658-702).
  5. Compute the gradient as max(apex_rt).max(1.0) (features.rs:704); this is the run length used to normalize RT errors.
  6. Compute the three cross-charge corroboration columns by grouping rows by peptidoform (features.rs:714-734).
  7. In parallel over rows (features.rs:757-808), compute the two expensive per-PSM pieces: fragment_features (the Minimal/Rich fragment battery) and, when Extended, the Evidence (peak_window, PeakTraces, evidence_from) + extended_values. Results go into a per-row FragFeatures vector and one flat rows x n_ext value buffer, both indexed by row, preserving order.
  8. Serially assemble fmap (name -> per-row value vector), prelim, and the elution bounds (features.rs:810-905). The serial loop reads per[i] and pushes each named value with the push closure (features.rs:738).
  9. Insert the three cross-charge columns (features.rs:908-910).
  10. Build the output columns (bookkeeping + active_features), write Parquet, write the schema JSON, build the feature matrix, write the PIN, and emit the report (features.rs:912-983).

Fragment features (fragment_features, features.rs:1193)

For each fragment it takes the observed intensity at the scan nearest the apex RT (obs), the predicted intensity (pred), and |ppm|. It then computes: frag_corr=pearson(obs,pred), frag_cosine=cosine(obs,pred), spectral_angle; sum-normalized L1 (norm_manhattan) and rmsd; intensity- weighted and unweighted mean |ppm|; b/y intensity sums and counts.

It then aligns all fragment traces on the union RT axis (axis_full), restricts to the elution peak (lo_i..=hi_i) so co-elution/profile features are not diluted over the whole +/- w_rt window, and computes pairwise Pearson statistics (coelution_mean, coelution_best, n_coelution_above), lag- optimized cross-correlation (xcorr_coelution mean absolute lag, xcorr_shape), the DIA-NN-style profile block (profile_cos = elution^2- weighted spectral cosine, ref_corr = mean fragment-vs-reference Pearson, best_ref_corr, low_frag_coel), and the interference-correction block (evidence = summed fragment-vs-reference correlations, contrast_min, resid_corr, coel_clean, shadow_frac from the 1.5*r*ref cap). Finally log_sn from apex vs median trace point. elution_lo/elution_hi/ base_width_rt/n_observations come from the peak-bounded axis.

Elution-boundary detection

peak_bounds (features.rs:1049) descends from the apex-nearest scan while the smoothed profile stays >= frac * apex_height, bridging up to grace consecutive sub-threshold scans. If the supplied apex sits at zero height it is relocated to the global maximum first (features.rs:1057), so a zero-height apex does not collapse the window. The reference profile for boundary finding is the smooth3 (features.rs:1028) of the summed top-3-predicted-intensity fragment XICs.

prelim_score (features.rs:901)

prelim = n_matched * (0.5 + max(0, frag_corr))
       + max(0, coelution_mean)
       + 0.1 * ln(1 + apex_intensity)
       - rt_err / gradient

A cheap heuristic (not a trained score) that rewards matched-fragment count scaled by spectral correlation, co-elution, and log intensity, penalized by gradient-normalized RT error. compete uses it as the within-group ranking key.

Directly-computed Minimal/Rich scalars (serial loop, features.rs:810-905)

Several Minimal/Rich columns are assembled inline in the serial loop rather than inside fragment_features:

  • rt_error_abs = calibrated_rt_error(apex_rt, rt_pred_cal) (features.rs:271), i.e. |apex_rt - rt_pred_cal| when both are finite and 0.0 when either is non-finite; rt_error_rel = that over gradient (features.rs:820-821). Stage B marks an unavailable RT calibration (fewer than two anchors) as NaN, so this guard keeps the sentinel from leaking NaN into the feature matrix or prelim (tested at features.rs:1608).
  • log_apex_intensity = ln(1 + apex_intensity) (features.rs:824).
  • n_matched_fragments, coelution_run, charge pass through from the PSM columns; peptide_length calls peptide_length(peptidoform).
  • n_proteins = protein.matches(';').count() + 1 (features.rs:842): the protein string is semicolon-delimited, so this counts group membership.
  • diff_by_intensity = sum_b_intensity - sum_y_intensity (features.rs:852).
  • matched_fraction = n_matched / max(1, n_predicted) (features.rs:874- 878).

MS1 isotope features (Rich set, isotope_features, features.rs:1515)

The four Rich MS1 columns (isotope_corr, ms1_isom1_ratio, log_mono_ms1, has_ms1) come from isotope_features (features.rs:1515) over the apex isotope intensities carried on the PSM rows. It fits a Poisson-averagine envelope [1, lambda, lambda^2/2] with lambda = 0.00052 * neutral_mass (features.rs:1525), the neutral mass being precursor_mz * charge - charge * PROTON (features.rs:816). This 0.00052 averagine differs from the 0.000594 used by the Extended ms1 family (ms1.rs:85); the two MS1 code paths are independent and both are kept. isotope_corr = pearson([mono, +1, +2], theo), ms1_isom1_ratio = isom1 / (mono + 1), log_mono_ms1 = ln(1 + mono), has_ms1 = 1.0 when mono/+1/+2 are all present, else all four return 0.0.

The Evidence struct (features.rs:293)

evidence_from constructs the per-PSM Evidence handed to every Extended family, from parts the chunked pass builds once per PSM (see the end of this section). It mirrors the alignment and peak-bounding of fragment_features so families see the same elution peak. Fields:

  • Time series: axis (RT seconds over the detected elution peak), traces (per-fragment intensity over axis, zero-filled, fragment order), axis_full / traces_full (whole extracted window), apex_idx (index of apex in axis), ref_profile (predicted-intensity-weighted sum of the peak traces).
  • Fragment-indexed arrays (shared order): pred (library intensity), obs_apex (intensity at the apex scan; > 0 defines "matched"), is_b, ordinal, frag_charge, frag_mz (theoretical), frag_obs_mz (intensity- weighted observed), mass_err_ppm (signed ppm).
  • apex_rt is set inside evidence_from itself from its apex_rt argument, not by the caller.
  • Scalars filled by the caller after build (features.rs:784-798): rt_pred_cal, rt_err (via calibrated_rt_error), gradient, precursor_mz, charge, seq_len, n_matched, n_predicted, seed_score, seed_identified, apex_intensity (plus the MS1 apex isotopes below, features.rs:795-798). All start at a zero/None default set by evidence_from.
  • MS1: ms1_mono/ms1_iso1/ms1_iso2/ms1_isom1 (apex isotope intensities, None when no MS1) and ms1_xic (the [mono,+1,+2] XICs resampled onto axis; empty unless the extract stage persisted ms1_* chromatogram rows). A row that samples axis_full itself (what extract writes: the fragments' window grid) is sliced directly; any other grid goes through the RT-keyed map.
  • ref_profile_full (the raw-weighted reference over the whole window) and pair_stats (see "Per-PSM shared statistics" below).

In the chunked pass the parts are built once per PSM and shared with fragment_features: apex_intensities (the nearest-apex intensity per row), peak_window (the elution-peak window; global half-widths or the walk down the smoothed top-3 profile), and PeakTraces (the window, the sliced traces, both reference profiles and a PairStats). evidence_from then assembles the Evidence from them. build_evidence is a test-only (#[cfg(test)]) wrapper that builds those parts and calls evidence_from in one step; a release build has no build_evidence.

Per-PSM shared statistics (PairStats)

Several readers computed the same correlations from the same peak traces: fragment_features and coelution both built the pair Pearson matrix and every pair's lag-optimised cross-correlation, ion_series and nonzero read subsets of that matrix, four readers correlated every fragment with ref_profile, and two correlated every full-window trace with ref_profile_full. PairStats holds each of these once: the k x k Pearson matrix (1.0 on the diagonal), the cross-correlation for a < b, the per-trace norms, and the reference and full-window reference correlations. Every entry comes from the kernel call the reader made, on the same arguments in the same order, so a read equals the computation it replaces bit for bit. A reader uses the cache only when PairStats::fits matches its evidence (and, for the full-window values, only when it reads ref_profile_full itself rather than a rebuild); otherwise it computes the value, which is what the test fixtures (pair_stats: None) exercise. fragment_features borrows the shared window only when bound_features is set, because without it it scores the whole window.

parse_ion (features.rs:346) parses b3, y7, b3^2 into (is_b, ordinal, charge).

The Extended battery

FAMILIES (features.rs:52) is an ordered array of (NAMES, values) pairs; the order is part of the frozen schema and is append-only. Each family exposes NAMES: &[&str] and values(&Evidence) -> Vec<f64> of identical length and matching order. extended_values (features.rs:138) calls each family and applies the precomputed dedup plan (extended_value_plan, features.rs:112), keeping names and values in lockstep, and coerces any non-finite value to 0.0.

Deduplication (extended_name_refs, features.rs:84): a name that already appears in MINIMAL_FEATURES or RICH_EXTRA (reserved_names, features.rs:72), or that repeats across families, is kept only on first appearance. In the current tree four Extended names are dropped as reserved collisions: spectral_angle (similarity vs Minimal), rt_error_abs (rt vs Minimal), peptide_length and seed_identified (novel vs Minimal/Rich). So the 335 raw family names reduce to 331 unique Extended names.

The three feature sets

active_features(set) (features.rs:210) returns the ordered active column list:

  • Minimal (14, MINIMAL_FEATURES features.rs:158): rt_error_abs, rt_error_rel, n_matched_fragments, coelution_run, log_apex_intensity, frag_corr, frag_cosine, spectral_angle, coelution_mean, coelution_best, n_coelution_above, charge, peptide_length, n_proteins.
  • Rich (44 = Minimal + 30, RICH_EXTRA features.rs:176): adds library_norm_manhattan, library_rmsd, xcorr_coelution, xcorr_shape, sum_b_intensity, sum_y_intensity, diff_by_intensity, n_b_ions, n_y_ions, weighted_mass_error, mean_mass_error, isotope_corr, ms1_isom1_ratio, log_mono_ms1, has_ms1, log_sn, n_observations, base_width_rt, seed_score, seed_identified, matched_fraction, profile_cos, ref_corr, best_ref_corr, low_frag_coel, evidence, contrast_min, resid_corr, coel_clean, shadow_frac.
  • Extended (381 = Minimal + Rich + 331 family names + 6 psms-derived): appends extended_names() then six columns computed outside the family registry. Three co-elution peak-contest metrics (peak_contested_frac, peak_contested_count_frac, peak_apportioned_frac) come from the PSM columns (default 0.0 when absent) and separate peak-borrowing decoys from genuine IDs. Three charge-corroboration columns (n_charge_states, charge_multi_flag, cross_charge_intensity_log) aggregate across the charge states of one peptidoform, an axis invisible to the per-PSM Evidence families.

The size invariant is asserted in feature_sets_sized (features.rs:1590): Extended.len() == 14 + 30 + extended_names().len() + 6, and all Extended names are unique. FeatureSet (config.rs:65) has exactly Minimal, Rich, Extended; the earlier Custom variant is gone (do not document it).

Extended family reference

Each family below lists its name count and what it measures. All values are finite (NaN/Inf coerced to 0.0), and every family returns a stable-length vector even on degenerate evidence.

similarity (64, similarity.rs)

Observed-vs-library fragment-intensity agreement under many kernels. o is the per-fragment apex intensity, l the predicted intensity; _matched restricts to o_i > 0, _area replaces o with the per-fragment peak-XIC trapezoid area, on/ln are sum-normalized. Names: spectrum_cosine_matched, spectrum_cosine_sqrt, spectrum_cosine_log, spectral_angle (dropped as reserved), spectral_angle_sqrt, spectral_angle_matched, pearson_intensity_matched, pearson_intensity_log, spearman_intensity, spearman_intensity_matched, kendall_tau_intensity, dot_product_raw, dot_product_norm, library_recall_intensity, manhattan_sim, manhattan_sqrt, rmsd_norm, mae_norm, mse_log, mae_weighted_pred, abs_diff_q3, max_positive_residual, chebyshev_dist, minkowski_p3, bray_curtis, bray_curtis_sqrt, canberra, canberra_matched, wave_hedges, chi_square_pearson, chi_square_symmetric, divergence_distance, bhattacharyya_coef, hellinger, squared_chord, harmonic_mean_sim, jaccard_presence, dice_presence, intensity_weighted_pearson, regression_slope, gini_diff, wasserstein_mz, footrule_norm, rank_overlap_top3, top1_frag_match, top1_predicted_observed, frac_top3_predicted_observed, count_strong_predicted_absent, frac_predicted_absent, cosine_area, pearson_area, spectral_angle_area, cosine_fullwindow, stein_scott_weighted_dot, and the unbounded/granular scores that address cosine saturation (log_dot_product, spectral_log_evidence, scribe_score, plus their _area twins), cosine_high_ordinal (drop b1/b2/y1/y2), and the robust trimmed-cosine trajectory (cosine_robust_trim1/2/3). Local helpers include tie-corrected ranks, kendall_tau_b, interpolated quantile, gini, weighted_pearson, and trapz.

entropy (18, entropy.rs)

Li spectral-entropy similarity and information divergences between sum- normalized o and l. entropy_sim (entropy.rs:107) is 1 - (2 H(m) - H(o) - H(l)) / ln 4 with m = (o+l)/2, clamped to [0,1]. Names: spectral_entropy_similarity, weighted_spectral_entropy_similarity (Li per-spectrum weighting), spectral_entropy_similarity_sqrt, spectral_entropy_similarity_topk (top-6 by predicted intensity), spectral_entropy_similarity_area, jensen_shannon_divergence, jeffreys_divergence, kl_obs_pred, kl_pred_obs, cross_entropy_obs_pred, obs_spectrum_entropy, pred_spectrum_entropy, entropy_diff, entropy_ratio, obs_normalized_entropy (Pielou evenness), normalized_entropy_diff, residual_spectrum_entropy, entropy_weight_obs. The public spectral_entropy_similarity_sqrt(obs, pred) (entropy.rs:92) is reused by the extraction gate GateMode::SpectralEntropy so the kernel is not duplicated.

coelution (38, coelution.rs)

Fragment co-elution against the predicted-weighted reference profile R and pairwise between fragments. Three groups: per-fragment vs R (peak/full windows, leave-one-out), pairwise Pearson statistics and lag cross-correlation (best_xcorr, MAXLAG=5), and structured cross-correlations (b vs y ions, charge-1 vs multiply-charged). Names include frag_ref_corr_mean, frag_ref_corr_obsweighted, frag_ref_corr_min, frag_ref_corr_std, frag_ref_corr_sq_mean, frag_ref_corr_topk_weighted, n_frag_ref_corr_above_0_9, frac_frag_ref_corr_above_0_8, frag_ref_corr_mean_full, full_vs_peak_corr_gain, pairwise_coelution_weighted, pairwise_coelution_min/median/std/frac_negative, pairwise_coelution_hi/lo, coelution_hi_lo_contrast, coelution_corr_entropy (Shannon entropy of a 10-bin NBINS_CORR pairwise-Pearson histogram over [-1, 1]), the xcorr_shape_* and xcorr_lag_* statistics (_mean/_min/_std/_mean_abs/_iqr/_frac_zero/_max_abs/_entropy), ref_xcorr_lag_mean, ref_xcorr_shape_mean, observed_sum_vs_template_corr, frag_loo_ref_corr_mean/_min, frac_frags_apex_aligned, top3_frag_ref_corr, by_cross_coelution, by_cross_lag_mean, charge_cross_coelution. Note the JSON coelution_weighted_mean is an exact alias of pairwise_coelution_weighted and is emitted once under the latter.

interference (26, interference.rs)

Co-isolation/chimera detection via the least-squares scale r_f = <x_f, R> / <R, R> projecting each fragment onto R. Two gating constants: IFS_MIN_CORR = 0.6 (interference.rs:48) prunes the least-coherent matched fragment in the iterative remove_ifs loop until its leave-one-out ref-corr stays above it or only 3 fragments remain, and COHERENT_THR = 0.7 (interference.rs:47) gates explained_apex_intensity_frac and apex_purity (only fragments with ref-corr at or above it count as explained); n_interfered_fragments flags a fragment whose apex intensity exceeds 2 * r_f * R_apex. Names: explained_variance_ref, profile_residual_fraction, n_interfered_fragments, corrected_vs_raw_cos, corrected_vs_raw_ratio, the iterative interference- removal block (ifs_removed_count, ifs_removed_intensity_frac, ifs_corr_gain, ifs_retained_frac, matched_frac_after_ifs), the peak-vs- full area ratios (peak_to_full_area_ratio_profile/_frag_mean/_weighted, out_of_peak_intensity_frac), profile_corr_full_vs_peak_delta, frac_frag_ref_corr_below_0_5, explained_apex_intensity_frac, apex_purity, interference_apex_residual_fraction, dominant_frag_ref_corr, the rank decomposition of the matched fragment x time Gram matrix by power iteration (explained_variance_ratio, second_component_fraction), competing-peak descriptors on the full-window profile (profile_second_peak_ratio, n_competing_peaks_in_window), matched_pred_intensity_fraction, and top_pred_frag_matched.

chromatographic (43, chromatographic.rs)

Peak-shape quality of R (Gaussian moment-match fit, EMG grid-fit comparison via erfc) and of individual fragments. Names: gaussian_fit_r2, gaussian_cosine, emg_fit_improvement, apex_prominence, profile_peak_snr (MAD-based over out-of-peak scans), width descriptors (fwhm_seconds, fwhm_to_window_ratio, width_at_10pct, width_ratio_10_50), asymmetry/tailing (hwhm_asymmetry, tailing_factor_usp, asymmetry_factor_10pct), apex shape (apex_sharpness, apex_curvature, apex_to_boundary_ratio, apex_dominance), roughness (zigzag_index, jaggedness, roughness_2nd_deriv), multimodality (n_local_maxima, modality), RT moments (rt_skewness, rt_excess_kurtosis, rt_std_seconds, mean_mode_offset), area descriptors (fraction_area_within_fwhm, triangle_area_similarity), baseline_fraction, peak_completeness, apex_centering_offset, intensity_score, total_xic_log, per-fragment descriptors (frag_fwhm_cv/_mean, frag_apex_rt_dispersion/_weighted, frag_apex_offset_from_profile_mean, frag_gaussianity_mean/_weighted, frag_zigzag_mean), sumtrace_unweighted_gaussian_r2, and reference_profile_rt_entropy_peak/_ratio. Has unit tests (chromatographic.rs:979).

mass_accuracy (17, mass_accuracy.rs)

Fragment ppm-error distribution over matched fragments (uses a fixed FRAG_TOL_PPM = 20.0; the config tolerance is not carried in Evidence). Names: median_abs_frag_ppm, signed_mean_frag_ppm, ppm_std, ppm_iqr, ppm_range, max_abs_frag_ppm, intensity_weighted_abs_ppm, intensity_weighted_signed_ppm, intensity_weighted_ppm_std, lib_weighted_abs_ppm, frac_frag_within_half_tol, high_ppm_intensity_frac, ppm_intensity_anticorr, mass_error_mz_trend, mean_abs_mz_error_da, and the positive DIA-NN-style evidence mass_evidence_gauss (predicted-weighted Gaussian concentration, sigma 10 ppm) and mass_log_evidence. precursor_mass_error_ppm is deliberately skipped (no theoretical precursor m/z in Evidence).

ion_series (34, ion_series.rs)

b/y series coverage, ladder contiguity, complementarity, per-series similarity. Names: n_matched_b/_y, frac_matched_b/_y, by_count_balance, by_intensity_ratio, by_ratio_agreement (tanh-squashed log-odds), by_ratio_consistency, longest_b_run/_y_run, longest_run_max, longest_run_frac_length, series_coverage_b/_y, sequence_coverage (cleavage-site union), series_gap_fraction, by_complement_count, by_complement_mz_consistency (b+y obs m/z vs M + 2 proton), by_complement_coelution, ordinal_intensity_concordance_y/_b, series_coelution_y/_b, spectral_angle_b/_y, pearson_b/_y, cosine_charge1/_charge2, charge_corr_balance, mean_matched_ordinal_norm, by_ion_contiguous_intensity, by_ion_contiguous_lib_frac, both_series_present. Uses PROTON from mumdia-core::constants.

ms1 (25, ms1.rs)

Precursor isotope-envelope agreement against a Poisson-averagine model (lambda = 0.000594 * M) plus MS1/MS2 XIC co-elution. Apex names: ms1_isotope_cosine_apex, ms1_isotope_spectral_angle_apex, ms1_isotope_chi2_apex, ms1_isotope_manhattan_apex, iso_ratio_1_0, iso_ratio_2_0, iso_plus1_ratio_dev, iso_plus2_ratio_dev, iso_minus_one_fraction, iso_overlap_flag, log_ms1_mono, ms1_total_isotope_log, has_ms1_signal, ms1_isotope_apex_entropy_3, ms1_m1_entropy_contribution. The XIC block (ms1_ms2_time_corr, ms1_ms2_envelope_time_corr, ms1_iso_coelution, ms1_ms2_apex_rt_delta, ms1_iso_ratio_stability, ms1_mono_gaussianity, ms1_ms2_fwhm_ratio, ms1_isotope_corr_xic, ms1_envelope_over_time_corr, ms1_isotope_xic_shape_consistency) reads Evidence.ms1_xic. Extract now persists ms1_mono, ms1_iso1, and ms1_iso2 chromatogram rows when MS1 input and a scan grid are available, so these features are populated in normal orchestrated runs. They remain 0.0 for older artifacts, standalone extraction without --ms1, or candidates without usable MS1/grid evidence.

rt (13, rt.rs)

RT agreement between the observed apex and the calibrated predicted RT. Names: rt_error_signed, rt_error_abs (dropped as reserved), rt_error_squared, rt_error_signed_norm_gradient, rt_error_abs_norm_gradient, observed_rt_raw, predicted_rt_raw, observed_rt_fraction, predicted_rt_fraction, rt_error_over_peak_width (base width at 10%), rt_error_over_fwhm, rt_diff_profile_apex (vs full-window profile argmax RT), predicted_rt_in_gradient.

novel (12, novel.rs)

Seed corroboration and precursor/charge metadata. Names: log_seed_hyperscore, seed_hyperscore_per_matched, seed_identified (dropped as reserved), peptide_length (dropped as reserved), precursor_charge, charge_is_2/_is_3/_is_4plus, precursor_mass, log_total_matched_intensity, n_matched_frags, n_predicted_frags. Sequence-dependent features (missed cleavages, C-terminal residue, modification count) are skipped because Evidence carries only seq_len.

nonzero (12, nonzero.rs)

Zero-ignoring variants of the apex spectral and co-elution features. The apex scan samples one grid point, so a fragment peaking one scan off reads 0 there (~8% of fragments); these recompute over the per-fragment peak-max and over present-only scans. Names: frag_corr_peakmax, frag_cosine_peakmax, spectral_angle_peakmax, frag_corr_matched_nz, frag_cosine_matched_nz, peakmax_apex_gain, n_frag_present_inpeak, frac_frag_present_inpeak, coelution_mean_bothpos, coelution_mean_summpos, ref_corr_nz, profile_cos_nz.

order_consistency (8, order_consistency.rs)

Prediction-free MS2-XIC rank stability: at each scan the fragments are ranked by observed intensity, and the features measure whether that ranking persists across the peak (orthogonal to library-agreement families). Names: rank_corr_vs_apex_mean, rank_corr_vs_apex_std, rank_corr_adjacent_mean, kendall_vs_apex_mean, top1_frag_persistence, top2_order_persistence, argmax_frag_entropy, self_cosine_vs_apex_mean. Degenerate below 3 fragments or 3 non-empty scans. Has unit tests (order_consistency.rs:300).

peak_scans (2, peak_scans.rs)

Label-blind window-degeneracy indicator, emitted for every PSM so the rescorer can tell an undefined zero from a measured zero when the window-based families collapse. Names: n_peak_scans, peak_window_degenerate (1 when fewer than 3 non-empty scans, mirroring order_consistency::MIN_SCANS).

apex_dispersion (13, apex_dispersion.rs)

Fragment apex dispersion and consensus peak shape, intensity-independent (breadth of co-elution rather than height). Names: frag_apex_rt_std, frag_apex_rt_mad, frag_apex_max_dev, frag_apex_mean_dev, frag_apex_agree_frac, precursor_frag_apex_delta, peak_symmetry, peak_tailing, peak_n_local_maxima, peak_shoulder_score, peak_fwhm_scans, peak_truncation, apex_frac_of_window. precursor_frag_apex_delta reads the mono MS1 XIC (ms1_xic[0]) and is populated under the same MS1/grid conditions as the ms1 XIC block; otherwise it is 0.0. Has unit tests (apex_dispersion.rs:240).

mass_uncertainty (10, mass_uncertainty.rs)

Fragment mass-error distribution over matched fragments plus evidence-breadth. Names: frag_mass_err_median, frag_mass_err_abs_median, frag_mass_err_std, frag_mass_err_iqr, frag_mass_err_max_abs, frag_mass_err_range, effective_frag_count (inverse participation ratio), evidence_concentration (fraction in the strongest fragment), frac_top3_pred_observed, frac_top5_pred_observed. Has unit tests (mass_uncertainty.rs:142).

The Percolator PIN (write_pin, features.rs:1537)

Streamed row-by-row through a BufWriter (not materialized as one String). Header: SpecId\tLabel\tScanNr\tExpMass\tCalcMass\t<features joined by tab>\t Peptide\tProteins. Per row: SpecId = cand_<candidate_id>, Label = -1 when label == "decoy" else 1, ScanNr = candidate_id, ExpMass = CalcMass = precursor_mz ({:.5}), each feature at {:.6}, then Peptide = -.<peptidoform>.- and Proteins = <protein>. The feature matrix passed to the PIN is built column-parallel from fmap in active_features order (features.rs:945), so the PIN and the Parquet share the same ordered feature list.

The feature-schema hash (feature_schema_id, features.rs:233)

blake3_str(cols.join(",")) of the ordered active column list. Written to <out>.schema.json and recorded in the report stats. It is a content hash of the exact ordered names, so any addition, removal, or reordering changes the id and a classifier trained under one schema is never silently applied under another. Because the family registry order is frozen and append-only, appending a new family or feature at the end changes the id predictably while leaving all prior positions stable.

bound_from_confident (elution-boundary calibration, features.rs:658)

When bound_from_confident is true (the default), the stage learns one pair of elution half-widths (L, R) in seconds from the confident-target seed set (spectrum_q <= 0.01, label == "target"; the same anchor set used for RT calibration and DeepLC fine-tune). For each confident candidate it detects the per-candidate peak with elution_peak_rt_bounds (features.rs:1141, which returns None for a candidate with fewer than 3 distinct scans, so it does not contribute an anchor) and records apex - lo and hi - apex. If at least 20 anchors resolve, it takes the bound_confident_pct percentile of the left and right half-widths (median by default) and returns Some((L, R)); every candidate is then bounded on [apex - L, apex + R] via global_bound_indices (features.rs:1114), which falls back to the single apex-nearest scan (ai, ai) when the mapped window collapses between grid points (sparse grid or a half-width below one cycle). This removes per-candidate boundary manipulation, so a chimeric decoy is scored over a real-peptide-width window centred on its apex rather than one it can widen or narrow. With fewer than 20 anchors the stage logs a warning and falls back to per-candidate boundary detection for that run. When the flag is false, every candidate detects its own peak boundary from its top-3-predicted-fragment profile (the legacy path). The same global_bounds argument threads into peak_window, which sets the window of both fragment_features (under bound_features) and the Extended Evidence, so Minimal/Rich and Extended see the identical window.

The shared stats kernel (stats.rs)

One implementation of pearson, cosine, spectral_angle, used by fragment_features and every family (do not reimplement). Two faster forms of pearson return its value bit for bit: pearson_vs(a, r, &Centered::new(r)) correlates many vectors with one reference centred once, and pearson_pairs(rows, out) produces every pair a < b from rows centred once, four covariances at a time. Both reproduce pearson's operations exactly (the mean by the same iter().sum() over the same count, deviations and squares accumulated from 0.0 in ascending order) and hand a length mismatch to pearson itself. pearson (stats.rs:6) is population Pearson with a zero-variance guard returning 0.0 for n < 2 or zero variance. cosine (stats.rs:29) returns 0.0 if either vector is all-zero. spectral_angle (stats.rs:44) is 1 - 2 * acos(clamp(cosine, -1, 1)) / pi, in [0,1] with 1 = identical. Families add local specializations that do not belong in the shared kernel (weighted Pearson, Spearman via pearson on average ranks, Kendall tau, windowed cross- correlation via super::best_xcorr, features.rs:1484, which returns the best normalized correlation and its integer lag over [-maxlag, maxlag] as (lag_of_max, value.max(0.0))), but the base Pearson/cosine call the shared functions. best_xcorr_normed takes the two norms (xcorr_norm) so a pair matrix computes each once; at the lag window every caller uses (5) it accumulates the 11 lags side by side (lag_dots), each lag's terms still in ascending order, so the result is the per-lag loop's bit for bit.

The per-PSM order statistics (peak_snr's median and MAD, the log_sn noise median, coelution's median and IQR) are taken by selection (order_stat, median_select, quantile_select), not by sorting a copy; under f64::total_cmp every order statistic is unique bit for bit, so the values are those of the sorted copy. The interference power iteration stops as soon as its 100th iterate is known exactly (a bitwise fixed point, or a bitwise period-2 cycle whose parity gives it), and order_consistency ranks each scan column once.

Measured single-threaded on one pinned P-core over a 60,000-PSM HYE-shaped fixture (write_kernel_bench_fixture, minimum of 5-7 interleaved runs, a desktop shared with other builds): the per-PSM compute went from 6.76 s to 5.24 s over the seven kernel changes (MS1 slice -2.3%, lag-parallel cross-correlation -6.2%, centred Pearson -3 to -8%, selection -4.4%, power iteration -1.9%, rank cache within noise, shared statistics -6.1%). Every change leaves the features table byte-identical; the golden digests extended_features_match_the_pre_permutation_build and production_shaped_features_match_the_pre_kernel_build pin the values.

Configuration

FeaturesConfig (config.rs:578) is #[serde(default, deny_unknown_fields)], so every field defaults independently and an unknown config key is a hard load error. The set field's default is t() (config.rs:614, helper at config.rs:163), a generic Default-forwarding helper, so it resolves to FeatureSet::default() = Minimal. The config was pruned of dead fields (the FeatureSet::Custom variant no longer exists).

field default effect
set Minimal which set active_features returns (Minimal 14 / Rich 44 / Extended 381)
coelution_corr_threshold 0.9 threshold for n_coelution_above (count of pairwise fragment correlations at or above it)
prec_tol_ppm 20.0 precursor tolerance carried for feature bookkeeping
bound_features true restrict trace-based features to the elution peak instead of the whole extracted window; gates only the Minimal/Rich fragment_features path (the Extended Evidence always peak-bounds, see gotchas)
bound_peak_fraction 1/3 peak-boundary threshold as a fraction of apex height (DIA-NN-style; matched DIA-NN RT bounds best)
bound_peak_grace 0 consecutive sub-threshold scans to bridge before stopping (0 = stop at first miss; 1 bridges a single-scan dip)
bound_from_confident true learn one global left/right half-width from the confident seed set and apply it to every candidate; false = per-candidate detection
bound_confident_pct 50.0 percentile of the confident-set half-widths taken as the global half-width (50 = median)
chrom_loaders 3 chromatogram decode threads in the main pass, an upper bound that --threads and the chunk count also cap (see "The chunked pass" below); changes time and memory only, never a value or a byte of the features table

Note that the fragment tolerance used inside mass_accuracy is a hardcoded FRAG_TOL_PPM = 20.0 (mass_accuracy.rs:44), not prec_tol_ppm; Evidence does not carry the configured fragment tolerance.

The chunked pass: loaders, computation, writer

run_chunked processes the run in chunks planned up front from the candidate_id columns (plan_chunks; a chunk closes at CHUNK_CHROM_ROWS chromatogram rows or CHUNK_PSM_ROWS PSM rows and never cuts a candidate). Each chunk passes three stages on their own threads:

  1. Loaders. features.chrom_loaders threads (default 3) each claim the next unclaimed chunk, open that chunk's row span of the chromatogram table (TableFile::span) and decode it into a ChromChunk with a fragment-name table of its own (load_chunk). A loader may claim chunk j only while j < taken + loaders, where taken counts the chunks the computation has taken, so at most chrom_loaders + 1 chunks are resident; one loader is the previous single-loader bound of two. The setting is an upper bound: a pass runs no more loaders than it has chunks or than the engine's thread pool has threads (main_loaders_wanted), so --threads 1 decodes on one loader as it did before the loaders were parallel. Loaders beyond the first come from a process-wide pool of four (MAIN_LOADER_EXTRAS), so concurrent bands or runs share it rather than multiply it.
  2. Computation. The calling thread takes the chunks strictly in table order (ChunkLoader::take), runs the per-PSM kernels in parallel (rayon), then assembles the value matrix and the output columns serially.
  3. Writer. One thread encodes the previous chunk's columns (TableWriter::write_cols) and publishes the table only after the commit marker.

Because a chunk's rows are fixed by the plan and the computation takes chunks in order, every write_cols call receives the same columns at every loader count: the features table is byte-identical. the_loader_count_moves_no_byte_of_the_features_table checks 1, 2, 3 and 7 loaders exactly (bypassing the shared pool, so no arm can quietly run one) and the production leased path, on a chromatogram table with 64-row groups so that chunk spans start inside a row group, cover several and share a boundary group with the next chunk, and with one PSM per chunk for the zero-row chunks. A failed chunk is handed over in order like any other, so the error surfaces at the same chunk as with one loader; a loader panic stops the pass, and the computation reports that no loader is left rather than waiting (the thread scope then re-raises the panic, as it did for the single loader). Both are tested with an injected failure at chunk 3 while chunk 2 is still decoding on another loader: chunks 0 to 2 arrive intact and in order, then the failure, and no table is published (a_loader_failure_reaches_the_computation_in_chunk_order_and_never_hangs, a_loader_failure_mid_pass_publishes_nothing).

Every pass logs one features: pass timers line: loader busy and blocked time (summed over loaders), the computation's wait for the loaders, per-PSM compute, serial assembly, the wait for the writer, writer busy and idle time, and a binding hint naming the busiest stage per thread. A large wait_for_loader_ms with a small loader_blocked_ms means the pass is decode-bound; a large wait_for_writer_ms means the encoder binds.

Measured on the local desktop (32 threads, shared with other builds, so single runs carry noise):

input loaders wall binding notes
Astral 15-min run, 522,237 PSMs, ~8 points per trace, 8 chunks 1 7.6 s writer loader busy 5.3 s, compute 2.0 s, assembly 3.1 s, writer 6.8 s
same 3 8.2-8.5 s writer no gain: the encoder binds
HYE-shaped fixture (write_kernel_bench_fixture, 400,000 PSMs, 150-240-point traces, 7 chunks) 1 25.0-26.2 s loader loader busy 24 s
same 3 13.3-13.6 s compute / writer peak working set 1.83 -> 3.23 GiB

So the loaders pay where traces are long (a 2 h gradient, the immunopeptidomics windows). Where the encoder binds they buy nothing and cost memory plus some contention with the writer (the Astral run was 0.6-0.9 s slower at 3 loaders); set chrom_loaders: 1 there if that matters.

Invariants, determinism, gotchas

  • Every family must return exactly NAMES.len() values in NAMES order; a debug_assert_eq! in extended_values (features.rs:143) and in most families catches a mismatch in debug builds. Non-finite values are coerced to 0.0 at family boundaries and again in extended_values.
  • The FAMILIES registry order and each family's NAMES order are the frozen schema; they are append-only. Reordering or renaming changes schema_id and invalidates any trained classifier.
  • Deduplication is stable and precomputed once (extended_value_plan, features.rs:112), reproducing the same survivors and order as extended_name_refs, so names and values stay in lockstep across runs.
  • Determinism: the parallel per-PSM pass collects into a Vec indexed by row, so the serial assembly is byte-identical to a serial run regardless of thread count. RT-axis alignment maps intensities keyed by f32::to_bits (features.rs:407), which is exact-equality safe because the same rt values are reused, not recomputed.
  • "Matched" throughout the families means obs_apex[i] > 0.0 (observed at the apex scan), which differs subtly from "present in the peak" used by the nonzero family (per-fragment peak-max > 0).
  • The apex scan samples a single grid point; a fragment peaking one scan off reads 0.0 at the apex. The nonzero family exists specifically to give the classifier zero-tolerant variants alongside the originals.
  • peptide_length (features.rs:250) strips a leading DECOY_ prefix before counting residues, and ignores bracketed modifications, so the decoy marker is not a length-based target/decoy label leak (tested at features.rs:1581).
  • MS1 XIC features require ms1_* chromatogram rows. Normal run supplies converted MS1 spectra to extract, which writes those rows when a usable grid is present; older artifacts or standalone extraction without --ms1 legitimately leave the features at 0.0.
  • bound_features gates only the Minimal/Rich fragment_features path (features.rs:1286): when false, that path scores over the whole extracted window. The Extended Evidence build (peak_window, PeakTraces, evidence_from) takes no such flag and always peak-bounds axis/traces while still retaining axis_full/traces_full, so Extended families read whichever window they name regardless of bound_features, and global_bounds from bound_from_confident always applies to them.
  • The peptidoform grouping for cross-charge features uses the ProForma string, which is charge-independent and keeps DECOY_ peptidoforms grouped among themselves, so it is not a target/decoy label leak (features.rs:714).

How to extend / modify

  • To add a new Extended family: create a module under stages/features/ exposing pub const NAMES: &[&str] and pub fn values(&Evidence) -> Vec<f64> of matching length and order, mod-declare it (features.rs:32), and append (name, values) to FAMILIES (features.rs:52). Appending at the end keeps all prior schema positions stable. Reuse crate::stats and the parent helpers (mean, normalize_sum, best_xcorr, smooth3, peak_bounds) rather than reimplementing kernels. A family that needs a pair correlation, a pair cross-correlation or a fragment-vs-reference correlation should read it from e.pair_stats when PairStats::fits holds, and compute it otherwise. Add an arity unit test (values(&e).len() == NAMES.len()) and a degenerate-evidence finiteness test.
  • To add a Minimal/Rich feature: append the name to MINIMAL_FEATURES or RICH_EXTRA and push its value in the serial loop (features.rs:820); make sure no Extended family already uses the name, or it will be dropped as a reserved collision.
  • New names must be globally unique across Minimal, Rich, and every family; a collision is silently dropped by the dedup filter, so run feature_sets_sized (features.rs:1590) after any change to confirm the size and uniqueness invariants.
  • Do not reach into vendor formats or duplicate the mass model / stats kernel; Evidence is the sole per-PSM interface for families, and any new scalar a family needs must be added to Evidence (features.rs:293) and filled by the caller (features.rs:784).
  • Prefer adding a config field (backed by a default) over hardcoding a threshold, consistent with the project convention; the current hardcoded fragment tolerance in mass_accuracy is a documented exception awaiting a tolerance carried on Evidence.