Too high fragment m/z reported

[MatteoLacki]

Hi,

Basically I submit data as mgf with maximal fragment m/z of 1698.1472 but in the annotated fragments I get:

         psm_id fragment_type  fragment_ordinals  fragment_charge  fragment_mz_calculated  fragment_mz_experimental  fragment_intensity
1255868  193625             y                 33                1               3694.9072                 3694.9343            202362.0

So it reported fragment_mz_experimental = 3694.9343 , whereas this is impossible.

I can make a reproducible example but I would need to share data with you, ~800MiB.
I would then put it online but obviously you might have some preferences about where that could be to avoid any problems.

I also asked copilot on your github page about the issue and it already says it seems some bug:
I have declared in config "max_fragment_charge": 1, but that this does not get respected could obviously result in the problem above.

Best wishes,

MatteoLacki

[lazear]

This is almost certainly from deisotoping. Sage treats all experimental and theoretical peaks as z=1 - if you have an observed peak that appears to be at a higher charge state (e.g. isotopes with neutron/z spacing), Sage will multiply by z and store the decharged peak internally.

max_fragment_charge is orthogonal to the deisotoping feature - during database search, it multiples the experimental peaks by increasing charge to try and find a match.

So if you re-run the search with max_fragment_charge set to 3 (or whatever) and deisotope to False, you should see the same y33 annotation at the experimentally observed charge.

The more "straightforward" approach would be storing charged theoretical fragments, but this obviously leads to an unwanted explosion of memory consumption. A potential improvement here could be to store the charge of deisotoped peaks and report that in the annotated matches instead of always reporting z=1 in your case.

[MatteoLacki]

Hi,

Claude made a reproducible minimal example for you https://github.com/MatteoLacki/sagebug220
Note that I did set max_fragment_charge to 1. The theoretical calculation for peptide NSLLFNTSLDLLLSVVGFSFAYMGKELPTNYNEAK is correct:
y33 fragment is LLFNTSLDLLLSVVGFSFAYMGKELPTNYNEAK and that should show up as [M + H]⁺ (z = 1) ≈ 3694.9076 which it does.

In reported precursors the entry is:

   psm_id                              peptide               proteins  \
0       1  NSLLFNTSLDLLLSVVGFSFAYMGKELPTNYNEAK  sp|Q8TE23|TS1R2_HUMAN   

          protein_groups  num_proteins  num_protein_groups  filename  \
0  sp|Q8TE23|TS1R2_HUMAN             1                   1  tiny.mgf   

                                              scannr  rank  label   expmass  \
0  "idx=1080501 frame=8515 scan=130 tof=268640 ii...     1      1  3895.935   

    calcmass  charge  peptide_len  missed_cleavages  semi_enzymatic  \
0  3894.9749       4           35                 1               0   

   isotope_error  precursor_ppm  fragment_ppm  hyperscore  delta_next  \
0        1.00335      -11.07714      7.379143   23.915547    1.066943   

   delta_best      rt  aligned_rt  predicted_rt  delta_rt_model  ion_mobility  \
0         0.0  8.0162    0.890689           0.0        0.890689           0.0   

   predicted_mobility  delta_mobility  matched_peaks  longest_b  longest_y  \
0                 0.0             0.0              5          1          3   

   longest_y_pct  matched_intensity_pct  scored_candidates   poisson  \
0       0.085714               8.687813                827 -2.042693   

   sage_discriminant_score  posterior_error  spectrum_q  peptide_q  protein_q  \
0                 1.141314              1.0         1.0        1.0        1.0   

   protein_group_q  ms2_intensity  
0              1.0       212286.0  

and here it makes much more sense, as result is called mass and charge is 4.

But in reported fragments I see:

   psm_id fragment_type  fragment_ordinals  fragment_charge  \
0       1             b                  4                1   
1       1             y                 33                1   
2       1             y                  7                1   
3       1             y                  6                1   
4       1             y                  5                1   

   fragment_mz_calculated  fragment_mz_experimental  fragment_intensity  
0                428.2503                  428.2510              2678.0  
1               3694.9072                 3694.9343            202362.0  
2                839.3898                  839.3860               884.0  
3                738.3421                  738.3420                 0.0  
4                624.2992                  624.2920              6362.0  

So at the very least the columns are not correctly named, but seeing also there charge 1 instead of 4 also makes little sense.
So I think currently --annotate-matches simply does not work.

Best wishes,

[MatteoLacki]

AI pointed to this https://github.com/copilot/share/08015198-0a80-84d4-a101-280a80824019

[MatteoLacki]

It is difficult for me to judge if AI is right here or not. Depends what is kept in peak.mass at that moment of the algorithm.

[lazear]

I'm confused why you think that "--annotate-matches simply does not work.", columns are misnamed, or that AI might be correct.

[MatteoLacki]

Hi, so I was tired and mumbling might be a way to prove I am human ;)

But the core of the problem is that I did not expect (nor did @theGreatHerrLebert) for the reported fragments to contain m/z values from above the range. So today me and him spent some time and tokens on getting some fixes. We found out that the root of the problem is, as you already mentioned it, deisotope = true. Here is a summary of what was found out:

## Root Cause: Two Code Paths, One Formula

Sage pre-processes MS2 spectra in `SpectrumProcessor::process_ms2`
(`crates/sage/src/spectrum.rs`). There are two branches depending on the
`deisotope` configuration flag.

### Branch 1 — `deisotope: false` (`spectrum.rs:350`)

let mass = (mz - PROTON) * 1.0;
Peak { mass, intensity }

Every peak is stored as a **z=1 pseudo-neutral mass**: `peak.mass = observed_mz − H`.

### Branch 2 — `deisotope: true` (`spectrum.rs:336`, the active config)

let mass = (peak.mz - PROTON) * peak.charge.unwrap_or(1) as f32;
Peak { mass, intensity }

The deisotoping algorithm assigns a charge to each isotope envelope. The peak
mass is therefore the **true neutral mass**: `peak.mass = (observed_mz − H) × z`.

For a z=2 fragment observed at m/z 800:

| Branch | `peak.mass` stored |
|--------|--------------------|
| deisotope: false | `(800 − 1.007) × 1 = 798.993` |
| deisotope: true  | `(800 − 1.007) × 2 = 1597.986` |

### The reporting formula — `scoring.rs:704` (original)

let exp_mz = peak.mass + PROTON;

This formula is correct only for the `deisotope: false` branch, where
`peak.mass = observed_mz − H`, so `exp_mz = observed_mz`. When
`deisotope: true`, `peak.mass` is the full neutral mass, and the formula yields
the z=1 equivalent m/z — which can be far above the instrument's scan range:

peak.mass = 1597.986   (true neutral mass)
exp_mz    = 1597.986 + 1.007 = 1598.993   ← reported (wrong)
actual    = 800.000                         ← what's in the MGF

The same inconsistency affected `fragment_mz_calculated` (line 705) and
`fragment_charge` (line 728), which both used the **theoretical search charge**
(always 1 when `max_fragment_charge: 1`) rather than the actual observed charge
of the deisotoped peak.

### Why `max_fragment_charge: 1` did not prevent the problem

The `max_fragment_charge` setting controls the loop over **theoretical** charge
states during scoring (`1..max_fragment_charge`, line 691). At search charge 1,
the theoretical m/z equals the fragment neutral mass. With `deisotope: true`,
deisotoped peaks are also stored as neutral masses — so the z=1 theoretical
search found z=2 (and higher) deisotoped peaks perfectly, but then mis-reported
their m/z and charge.

We then wanted a fix for that. You can investigate both solutions in the fork I have here.

A quick, yet inefficient solution was to equip the Peak struct with charge and use it to correct the reported values, see commit cdabfa50e29832a7802c9f433e83221481167dec . That worked, but it would have significant impact on RAM locality in a hot-loop, so not anything we would like.

`select_most_intense_peak` is called for every fragment × every search charge ×
every candidate peptide — the innermost hot loop. A 64-byte cache line that used
to hold **8 peaks** now holds only **5** (64 ÷ 12 = 5.3). That is a **37%
reduction in peaks per cache line**, directly increasing cache-miss pressure in
the tightest loop of the engine.

### RAM impact

Peaks are stored in `Vec<Peak>` in every `ProcessedSpectrum`. With
`max_peaks: 800`:

| Spectra | Before (8 B) | After (12 B) | Extra |
|---------|-------------|--------------|-------|
| 10,000 | 61 MB | 91 MB | +30 MB |
| 100,000 | 610 MB | 915 MB | +305 MB |
| 1,000,000 | 6.1 GB | 9.2 GB | +3.1 GB |

A million-spectrum timsTOF run would consume ~3 GB of additional RAM purely from
struct padding — before accounting for any other data structures.

As alternative, we went for keeping a parallel copy of deisotope charges, commit c099b2caf00648e898176ff6735ded3beacacbdc, the last one on the fork.

It is more involved of course:

## Follow-up Fix: Parallel `Vec<u8>` for Deisotope Charges

To eliminate the `Peak` struct bloat while keeping correct output, the charge
should be carried in a parallel `Vec<u8>` on `ProcessedSpectrum` rather than
embedded in every `Peak`.

### Files to change

| File | What changes |
|------|-------------|
| `crates/sage/src/spectrum.rs` | Revert `Peak` to 8 bytes; add `peak_charges` to `ProcessedSpectrum`; change `process_ms2` return type; change `select_most_intense_peak` return type |
| `crates/sage/src/scoring.rs` | Two call sites updated (lines 613 and 695) |
| `crates/sage/src/tmt.rs` | One call site updated (line 201) |

### 1. Revert `Peak` to 8 bytes — `spectrum.rs:4-18`

Remove `charge: u8` and the manual `Default` impl; restore `#[derive(Default)]`:


#[derive(PartialEq, Copy, Clone, Default, Debug)]
pub struct Peak {
    pub intensity: f32,
    pub mass: f32,
}


Also revert the three construction sites (`spectrum.rs:337-340`, `351`, `369`)
back to their pre-fix form (drop the `charge:` field).

### 2. Add `peak_charges: Vec<u8>` to `ProcessedSpectrum` — `spectrum.rs:83`


pub struct ProcessedSpectrum<T> {
    // ... existing fields unchanged ...
    pub peaks: Vec<T>,
    /// Parallel to `peaks`. Non-empty only when deisotoping is enabled;
    /// `peak_charges[i]` is the observed charge of `peaks[i]`.
    /// Empty (all charges implicitly 1) when deisotope=false.
    pub peak_charges: Vec<u8>,
    pub total_ion_current: f32,
}


### 3. Change `process_ms2` to return `(Vec<Peak>, Vec<u8>)` — `spectrum.rs:301`

**Deisotope branch** — collect as `(Peak, u8)` pairs, sort by mass together,
then unzip:


let mut pairs: Vec<(Peak, u8)> = deisotoped
    .into_iter()
    .filter(|p| p.envelope.is_none())
    .map(|p| {
        let charge = p.charge.unwrap_or(1);
        let mass = (p.mz - PROTON) * charge as f32;
        (Peak { mass, intensity: p.intensity }, charge)
    })
    .take(self.take_top_n)
    .collect();
pairs.sort_by(|(a, _), (b, _)| a.mass.total_cmp(&b.mass));
let (peaks, charges) = pairs.into_iter().unzip();
(peaks, charges)


**Non-deisotope branch** — return empty charges vec (no change to peak logic):


(peaks, vec![])


### 4. Update `process()` to store both — `spectrum.rs:360`

Remove the standalone `peaks.sort_by(...)` (sorting now happens inside
`process_ms2`). Destructure the tuple and store both in `ProcessedSpectrum`:


let (peaks, peak_charges) = match spectrum.ms_level {
    2 => self.process_ms2(self.deisotope, &spectrum),
    _ => ( /* existing MS1 peak logic */ , vec![]),
};

ProcessedSpectrum { /* ... */ peaks, peak_charges, total_ion_current }


### 5. Change `select_most_intense_peak` to return `Option<usize>` — `spectrum.rs:163`

Return the index of the best peak instead of a reference, so callers can
independently look up the parallel charge vector:


pub fn select_most_intense_peak(
    peaks: &[Peak],
    center: f32,
    tolerance: Tolerance,
    offset: Option<f32>,
) -> Option<usize>   // ← was Option<&Peak>


Inner loop returns `Some(i + idx)` (index offset by the binary-search lower
bound `i`) instead of `Some(peak)`.

### 6. Update the three callers

**`scoring.rs:613`** (chimera peak removal — doesn't need charge):


if let Some(i) = crate::spectrum::select_most_intense_peak(...) {
    to_remove.push(query.peaks[i]);
}


**`scoring.rs:695`** (`score_candidate` — needs charge):


if let Some(i) = crate::spectrum::select_most_intense_peak(
    &query.peaks, mz, self.fragment_tol, None,
) {
    let peak = &query.peaks[i];
    let peak_charge = query.peak_charges.get(i).copied().unwrap_or(1);

    let exp_mz  = peak.mass / peak_charge as f32 + PROTON;
    let calc_mz = frag.monoisotopic_mass / peak_charge as f32 + PROTON;
    // ...
    fragments_details.charges.push(peak_charge as i32);
}

**`tmt.rs:201`** (reporter ion matching — doesn't need charge):


.map(|&label| {
    spectrum::select_most_intense_peak(peaks, label, label_tolerance, Some(-PROTON))
        .map(|i| &peaks[i])
})


### Result and memory/runtime implications

`Peak` is restored to **8 bytes** (two `f32`, no padding). Cache-line density
in `select_most_intense_peak` — the innermost hot loop — is back to 8 peaks per
64-byte line, identical to the original unfixed code.

The `peak_charges` field adds a fixed **24-byte** `Vec` header to every
`ProcessedSpectrum<T>` regardless of path. When non-empty (deisotope path only)
it additionally allocates `n_peaks` bytes on the heap:

| Path | Peaks | `peak_charges` | Total per spectrum |
|------|-------|----------------|--------------------|
| deisotope: false | 800 × 8 B = 6.4 KB | 24 B, no heap alloc | **6.4 KB** |
| deisotope: true | 800 × 8 B = 6.4 KB | 24 B + 800 B | **7.2 KB** |
| previous fix (charge in Peak) | 800 × 12 B = 9.6 KB | — | **9.6 KB** |

The parallel-vec design is **25% leaner** than embedding the charge in `Peak`
for the deisotope case, and has zero extra memory cost for the non-deisotope case.

The only extra work per matched fragment in `score_candidate` is one array index
(`query.peaks[i]`) and one byte read (`query.peak_charges.get(i)`), both on
already-hot cache lines. These execute only on matched peaks, which are rare
relative to the total `(fragment, charge)` pairs tried, so their cost is
negligible in practice.

One minor structural note: `peak_charges` is a field on `ProcessedSpectrum<T>`
for all `T`, including `IMPeak` (MS1 mobility spectra). Those always carry an
empty `Vec` (24 bytes in the struct, no heap allocation). This is harmless; if
it ever becomes a concern, `peak_charges` could be moved to a
`ProcessedSpectrum<Peak>`-specific wrapper.

---

## Files Changed

| File | Change |
|------|--------|
| `crates/sage/src/spectrum.rs` | `Peak` stays at 8 bytes; `peak_charges: Vec<u8>` added to `ProcessedSpectrum`; `process_ms2` returns `(Vec<Peak>, Vec<u8>)`; deisotope branch collects pairs and unzips after sorting; `select_most_intense_peak` returns `Option<usize>` instead of `Option<&Peak>` |
| `crates/sage/src/scoring.rs` | `score_candidate`: index lookup replaces reference; `peak_charge` from `query.peak_charges`; `exp_mz`, `calc_mz`, and reported charge all use `peak_charge`. Chimera scorer updated to index `query.peaks[i]` |
| `crates/sage/src/tmt.rs` | `find_reporter_ions`: maps index back to `&peaks[i]` |

Both solutions have resolved the problem to me:

main matteo@pingu$ csvlook out/matched_fragments.sage.tsv 
| psm_id | fragment_type | fragment_ordinals | fragment_charge | fragment_mz_calculated | fragment_mz_experimental | fragment_intensity |
| ------ | ------------- | ----------------- | --------------- | ---------------------- | ------------------------ | ------------------ |
|   True | b             |                 4 |               1 |               428,250… |                  428,251 |              2 678 |
|   True | y             |                33 |               4 |               924,482… |                  924,489 |            202 362 |
|   True | y             |                 7 |               1 |               839,390… |                  839,386 |                884 |
|   True | y             |                 6 |               1 |               738,342… |                  738,342 |                  0 |
|   True | y             |                 5 |               1 |               624,299… |                  624,292 |              6 362 |

If you would be interested in any of that solutions, just drop me a message and I will make a PR. If you have something better than this on top of your head - well, the better of course. Here I am but a simple consumer of AI and of course AI giveth and AI taketh away. In the very least it might be some inspiration.

Very best wishes and thanks for all the help and a truly nice open-source software to work with!

MatteoLacki

[lazear]

Second option is definitely preferable - see #178 - would be great to move all of the Peak representations to struct-of-array instead, since it's more performant and saves memory for cases like this.

So we would want a ProcessedSpectrum to ultimately be something like:

struct ProcessedSpectrum {
   mass: Vec<f32>,
   intensity: Vec<f32>,
   charge: Vec<u8>,
   im: Vec<f32>, // or maybe Option<Vec<f32>>
...
}

Preserves cache locality better, and opens the door for stuff like in-memory compression

[lazear]

@MatteoLacki wanna check out this PR? I believe this should address the concerns. Outstanding question is how to handle cases where users turn on deisotope and max_fragment_charge > 1. Should we report a set of columns like assumed_charge (deisotoped) and searched_charge (= assumed_charge unless max_fragment_charge > 1 and we get a match at z>1)