Feat xcorr

docs/api/features/chimeric.md

@@ -48,6 +48,7 @@ All columns from `FragmentMatchFeatures` with `chimeric_` prefix:
 | `chimeric_max_ion_gap` | Daltons (Da) | Largest gap between matched runner-up ions |
 | `chimeric_b_y_intensity_ratio` | Ratio | Ratio of b-ion to y-ion intensity for runner-up (including isotopic envelopes) |
 | `chimeric_spectral_angle` | Score (0-1) | Normalised spectral angle similarity between runner-up theoretical and observed intensities |
+| `chimeric_xcorr` | Score | SEQUEST fast cross-correlation score for the runner-up peptide. Measures overall agreement between the observed spectrum and the runner-up theoretical spectrum with local background correction. See [FragmentMatchFeatures — Cross-correlation Score](fragment_match.md#cross-correlation-score) for details on the algorithm. |
 
 ```python
 from winnow.calibration.features import ChimericFeatures

docs/api/features/fragment_match.md

@@ -57,6 +57,21 @@ For each theoretical peak, we search for the nearest observed peak using binary
 | `max_ion_gap` | Daltons (Da) | Largest m/z difference between two consecutive matched theoretical peaks when sorted by m/z. Large gaps may indicate missing fragmentation coverage. |
 | `b_y_intensity_ratio` | Ratio | Ratio of total matched b-ion intensity to total matched y-ion intensity (including isotopic envelopes). |
 | `spectral_angle` | Score (0-1) | Normalised spectral angle similarity between theoretical and matched observed intensity vectors. A value of 1 indicates perfect correlation, 0 indicates orthogonal vectors. |
+| `xcorr` | Score | SEQUEST fast cross-correlation score. Measures overall agreement between the observed and theoretical spectra with local background correction. Higher values indicate better matches. |
+
+### Cross-correlation Score
+
+The `xcorr` column implements the fast cross-correlation score function from SEQUEST (Eng et al., 2008). Unlike the other features which operate on individually matched peaks, xcorr evaluates the overall pattern of the observed spectrum against the theoretical spectrum using a background-corrected dot product.
+
+The observed spectrum is preprocessed by:
+
+1. **Binning** into near-unit-dalton bins (~1.0005 Da) with square-root intensity compression
+2. **Window normalization** — the maximum intensity within each of 10 equal m/z windows is normalized to a fixed value, making the score robust to different intensity scales
+3. **Background subtraction** — the mean intensity from ±75 neighbouring bins is subtracted from each bin, so that only signal specifically at the correct m/z positions contributes positively
+
+The score is then computed as the dot product of the preprocessed observed spectrum with a binary theoretical spectrum (1.0 at each predicted fragment ion bin position). This background correction naturally penalizes unmatched theoretical ions and discounts noise.
+
+**Reference:** Eng JK, Fischer B, Grossmann J, Maccoss MJ. A fast SEQUEST cross correlation algorithm. *J Proteome Res.* 2008 Oct;7(10):4598-602. doi: 10.1021/pr800420s.
 
 ## Usage
 

tests/calibration/features/test_chimeric.py

@@ -80,6 +80,7 @@ def test_properties(self, chimeric_features):
             "chimeric_max_ion_gap",
             "chimeric_b_y_intensity_ratio",
             "chimeric_spectral_angle",
+            "chimeric_xcorr",
             # Missing indicator (learn_from_missing=True by default)
             "is_missing_chimeric_features",
         ]
@@ -158,6 +159,7 @@ def test_learn_from_missing_false_columns_excludes_indicator(self):
             "chimeric_max_ion_gap",
             "chimeric_b_y_intensity_ratio",
             "chimeric_spectral_angle",
+            "chimeric_xcorr",
         ]
 
     # ------------------------------------------------------------------
@@ -232,6 +234,7 @@ def test_compute_preprocessing_pipeline(
         mock_max_ion_gap = 1
         mock_b_y_intensity_ratio = 0.5
         mock_spectral_angle = 0.8
+        mock_xcorr = [1.0, 2.0]
         mock_compute_ions.return_value = IonIdentificationResult(
             mock_match_rate,
             mock_match_intensity,
@@ -241,6 +244,7 @@ def test_compute_preprocessing_pipeline(
             mock_max_ion_gap,
             mock_b_y_intensity_ratio,
             mock_spectral_angle,
+            mock_xcorr,
         )
 
         # Run the compute method

tests/calibration/features/test_fragment_match.py

@@ -66,6 +66,7 @@ def test_properties(self, prosit_features):
             "max_ion_gap",
             "b_y_intensity_ratio",
             "spectral_angle",
+            "xcorr",
             # Missing indicator (learn_from_missing=True by default)
             "is_missing_fragment_match_features",
         ]
@@ -143,6 +144,7 @@ def test_learn_from_missing_false_columns_excludes_indicator(self):
             "max_ion_gap",
             "b_y_intensity_ratio",
             "spectral_angle",
+            "xcorr",
         ]
 
     def test_learn_from_missing_true_columns_includes_indicator(self):
@@ -256,6 +258,7 @@ def test_compute_with_mock(
         mock_max_ion_gap = 1
         mock_b_y_intensity_ratio = 0.5
         mock_spectral_angle = 0.8
+        mock_xcorr = [1.0, 2.0, 3.0]
         mock_compute_ions.return_value = IonIdentificationResult(
             mock_match_rate,
             mock_match_intensity,
@@ -265,6 +268,7 @@ def test_compute_with_mock(
             mock_max_ion_gap,
             mock_b_y_intensity_ratio,
             mock_spectral_angle,
+            mock_xcorr,
         )
 
         prosit_features.compute(sample_dataset_with_spectra)

tests/calibration/features/test_utils.py

@@ -17,6 +17,7 @@
     compute_max_ion_gap,
     compute_b_y_intensity_ratio,
     compute_spectral_angle,
+    compute_xcorr,
     parse_mz_tolerance_unit,
     _validate_mz_tolerance,
 )
@@ -969,3 +970,115 @@ def test_compute_spectral_angle_length_mismatch_raises(self):
         observed = [1.0, 2.0]  # Different length
         with pytest.raises(ValueError, match="must align"):
             compute_spectral_angle(theoretical, observed)
+
+
+class TestXcorr:
+    """Test the fast xcorr score function."""
+
+    def test_matching_peaks_positive_score(self):
+        """Peaks exactly at theoretical positions should produce a positive xcorr."""
+        theoretical_mz = [200.0, 400.0, 600.0]
+        observed_mz = [200.0, 400.0, 600.0]
+        observed_intensities = [1000.0, 2000.0, 1500.0]
+
+        score = compute_xcorr(observed_mz, observed_intensities, theoretical_mz)
+        assert score > 0.0
+
+    def test_no_matching_peaks_low_score(self):
+        """Peaks far from theoretical positions should yield a low or negative xcorr."""
+        theoretical_mz = [200.0, 400.0, 600.0]
+        observed_mz = [250.0, 450.0, 650.0]
+        observed_intensities = [1000.0, 2000.0, 1500.0]
+
+        matching_score = compute_xcorr(
+            [200.0, 400.0, 600.0], observed_intensities, [200.0, 400.0, 600.0]
+        )
+        non_matching_score = compute_xcorr(
+            observed_mz, observed_intensities, theoretical_mz
+        )
+        assert non_matching_score < matching_score
+
+    def test_nan_theoretical_returns_zero(self):
+        """NaN theoretical m/z (missing prediction) returns 0.0."""
+        import math
+
+        score = compute_xcorr([100.0], [1000.0], math.nan)
+        assert score == 0.0
+
+    def test_empty_theoretical_returns_zero(self):
+        """Empty theoretical m/z list returns 0.0."""
+        score = compute_xcorr([100.0], [1000.0], [])
+        assert score == 0.0
+
+    def test_empty_observed_returns_zero(self):
+        """Empty observed spectrum returns 0.0."""
+        score = compute_xcorr([], [], [100.0, 200.0])
+        assert score == 0.0
+
+    def test_uniform_spectrum_much_lower_than_matching(self):
+        """A uniform observed spectrum should score much lower than a matching one.
+
+        When every bin has the same intensity, the background subtraction
+        largely cancels the signal. Some residual remains due to edge effects
+        in the background window near spectrum boundaries, but the score
+        should be much lower than when peaks align with theoretical positions.
+        """
+        observed_mz = [float(i) for i in range(50, 500)]
+        observed_intensities = [100.0] * len(observed_mz)
+        theoretical_mz = [100.0, 200.0, 300.0]
+
+        uniform_score = compute_xcorr(observed_mz, observed_intensities, theoretical_mz)
+
+        sparse_mz = [100.0, 200.0, 300.0]
+        sparse_intensities = [1000.0, 1000.0, 1000.0]
+        matching_score = compute_xcorr(sparse_mz, sparse_intensities, theoretical_mz)
+
+        assert uniform_score < matching_score * 0.25
+
+    def test_scale_invariance(self):
+        """Xcorr should be the same regardless of observed intensity scale.
+
+        Window normalization rescales each region to a fixed value, so
+        multiplying all intensities by a constant should not change the score.
+        """
+        theoretical_mz = [200.0, 400.0, 600.0]
+        observed_mz = [200.0, 300.0, 400.0, 500.0, 600.0]
+        intensities_1x = [1000.0, 500.0, 2000.0, 300.0, 1500.0]
+        intensities_100x = [i * 100.0 for i in intensities_1x]
+
+        score_1x = compute_xcorr(observed_mz, intensities_1x, theoretical_mz)
+        score_100x = compute_xcorr(observed_mz, intensities_100x, theoretical_mz)
+        assert score_1x == pytest.approx(score_100x, rel=1e-10)
+
+    def test_more_matches_higher_score(self):
+        """More matching peaks should produce a higher xcorr than fewer matches."""
+        observed_mz = [100.0, 200.0, 300.0, 400.0, 500.0]
+        observed_intensities = [1000.0] * 5
+
+        score_2_matches = compute_xcorr(
+            observed_mz, observed_intensities, [100.0, 200.0]
+        )
+        score_4_matches = compute_xcorr(
+            observed_mz, observed_intensities, [100.0, 200.0, 300.0, 400.0]
+        )
+        assert score_4_matches > score_2_matches
+
+    def test_observed_preprocessing_independent_of_out_of_range_theoretical_ions(self):
+        """Unmatched high-m/z theoretical ions must not change observed normalization."""
+        observed_mz = [166.51, 239.26, 280.07, 391.08, 429.99]
+        observed_intensities = [416.6, 13651.9, 163.4, 806.1, 2155.4]
+        theoretical_mz = [239.26, 429.99, 391.08]
+
+        score = compute_xcorr(observed_mz, observed_intensities, theoretical_mz)
+        score_with_unmatched_ion = compute_xcorr(
+            observed_mz,
+            observed_intensities,
+            theoretical_mz + [2000.0],
+        )
+
+        assert score_with_unmatched_ion == pytest.approx(score)
+
+    def test_single_peak_positive(self):
+        """A single observed peak matching a single theoretical ion should score positive."""
+        score = compute_xcorr([300.0], [5000.0], [300.0])
+        assert score > 0.0

winnow/calibration/features/chimeric.py

@@ -131,6 +131,7 @@ def columns(self) -> List[str]:
             "chimeric_max_ion_gap",
             "chimeric_b_y_intensity_ratio",
             "chimeric_spectral_angle",
+            "chimeric_xcorr",
         ]
         if self.learn_from_missing:
             columns.append("is_missing_chimeric_features")
@@ -338,3 +339,4 @@ def compute(self, dataset: CalibrationDataset) -> None:
             ion_identifications.b_y_intensity_ratio
         )
         dataset.metadata["chimeric_spectral_angle"] = ion_identifications.spectral_angle
+        dataset.metadata["chimeric_xcorr"] = ion_identifications.xcorr

winnow/calibration/features/constants.py

@@ -7,3 +7,13 @@
 PROTON_MASS = 1.007276
 
 VALID_INTENSITY_MODEL_PROVIDERS = ["prosit", "ms2pip", "alphapeptdeep"]
+# Default bin width matching Comet's fragment_bin_tol for near-unit-dalton bins.
+XCORR_BIN_SIZE = 1.0005079
+# Default bin offset matching Comet's fragment_bin_offset.
+XCORR_BIN_OFFSET = 0.4
+# Number of equal-width m/z windows for intensity normalization (SEQUEST default).
+XCORR_NUM_WINDOWS = 10
+# Maximum m/z offset in bins for background subtraction (SEQUEST default ±75).
+XCORR_MAX_OFFSET = 75
+# Normalization target for the maximum intensity within each window.
+XCORR_WINDOW_NORM_VALUE = 50.0

winnow/calibration/features/fragment_match.py

@@ -130,6 +130,7 @@ def columns(self) -> List[str]:
             "max_ion_gap",
             "b_y_intensity_ratio",
             "spectral_angle",
+            "xcorr",
         ]
         if self.learn_from_missing:
             columns.append("is_missing_fragment_match_features")
@@ -313,3 +314,4 @@ def compute(self, dataset: CalibrationDataset) -> None:
             ion_identifications.b_y_intensity_ratio
         )
         dataset.metadata["spectral_angle"] = ion_identifications.spectral_angle
+        dataset.metadata["xcorr"] = ion_identifications.xcorr