Skip to content

Commit 2387fe4

Browse files
authored
feat: implement SRM/PRM chromatogram traces for mzML export (#30)
Wire the remaining piece of iter_chromatograms scoped out of #26: one SRM-shaped ChromatogramRecord ("selected reaction monitoring chromatogram", MS:1001473) per scheduled PrmTarget, built on top of the TIC/BPC traces #26 already wired up. For each PrmTarget, walk every msms_type == 10 (prm-PASEF) frame, decode its peak payload via the same per-frame decode path iter_spectra already uses (Reader::decode_peaks), sum the decoded intensity restricted to that target's PrmFrameMsMsInfo scan-number range, and append one point per frame. A target with a genuine zero-intensity frame keeps that point (real data, unlike a frame missing a Frames column); a target absent from every PrmFrameMsMsInfo row is omitted rather than emitted with an empty trace, matching the TIC/BPC no-empty-chromatogram contract. precursor_mz is the target's scheduled PrmTargets.MonoisotopicMz. product_mz is left unset: prm-PASEF records a full fragment-ion scan per precursor rather than a discrete triple-quadrupole product-ion transition, and the TDF schema has no product-ion m/z field to report - populating one would mean inventing a value with no schema backing, which the clean-room policy in CONTRIBUTING.md rules out. New pure helpers (sum_intensity_in_range, build_srm_record) are unit tested directly; the frame-walking orchestration (srm_chromatograms_for) follows the same shape as the existing spectra_for_frame and isn't unit tested in isolation for the same reason that one isn't - it's thin glue over already-tested pieces and a live Reader/bundle. The existing cache-gated prm_pasef_pxd028279 test in tests/roundtrip.rs already exercises the metadata this relies on (PrmFrameMsMsInfo/PrmTargets shape) and continues to pass. Distinct from #13 (prm-PASEF frames skipped in the mzML *spectrum* projection): that issue is about spectra, this is about the chromatogram trace. Addresses the remaining SRM/PRM scope of #25 (TIC/BPC portion done in #26). Claude-Session: https://claude.ai/code/session_018ijUinjq3SKNZXFNWqTTmV
1 parent 47ca1bb commit 2387fe4

2 files changed

Lines changed: 246 additions & 11 deletions

File tree

CHANGELOG.md

Lines changed: 18 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -19,10 +19,24 @@ This project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.htm
1919
`Frames.MaxIntensity`. `MaxIntensity` is a new additive column on the
2020
`Frames` SELECT and a new `Frame::max_intensity` field (also exposed on
2121
the Python `Frame`). A trace whose source column is absent on every frame
22-
is omitted rather than emitted empty. SRM/PRM transition chromatograms
23-
are intentionally left out: they need genuine per-target cross-frame
24-
aggregation of the decoded peak stream, not just column wiring; tracked
25-
as a rescoped #25. Contributed by @Nabejo.
22+
is omitted rather than emitted empty. Contributed by @Nabejo.
23+
- `iter_chromatograms` now also emits one SRM/PRM chromatogram trace
24+
("selected reaction monitoring chromatogram", `MS:1001473`) per scheduled
25+
`PrmTarget`, completing the chromatogram work scoped out of the entry
26+
above (rescoped #25; distinct from #13, which is about the prm-PASEF
27+
spectrum projection, not the chromatogram trace). For each target, every
28+
`msms_type == 10` (prm-PASEF) frame is decoded via the same per-frame
29+
peak decode `iter_spectra` already uses, the decoded intensity is summed
30+
restricted to that target's `PrmFrameMsMsInfo` scan-number range, and one
31+
point is appended per frame - a target with a genuine zero-intensity
32+
frame keeps that point (it is real data, unlike a frame missing a
33+
`Frames` column), while a target absent from every `PrmFrameMsMsInfo` row
34+
is omitted rather than emitted with an empty trace. `precursor_mz` is the
35+
target's scheduled `PrmTargets.MonoisotopicMz`; `product_mz` is left
36+
unset because prm-PASEF records a full fragment-ion scan per precursor
37+
rather than a discrete triple-quadrupole product-ion transition, and the
38+
TDF schema has no product-ion m/z field to report. Contributed by
39+
@Nabejo.
2640

2741
## [1.3.1] - 2026-07-25
2842

crates/opentimstdf/src/mzml.rs

Lines changed: 228 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -12,7 +12,10 @@
1212
//! Chromatograms: [`iter_chromatograms`](msc::SpectrumSource::iter_chromatograms)
1313
//! emits whole-run TIC and basepeak (BPC) traces built directly from the
1414
//! per-frame `Frames.SummedIntensities` / `Frames.MaxIntensity` columns (one
15-
//! point per frame). See [`chromatogram_records_for`].
15+
//! point per frame), plus one SRM/PRM trace per scheduled `PrmTarget` built by
16+
//! decoding and summing prm-PASEF (`msms_type == 10`) frames restricted to
17+
//! that target's scan-number range. See [`chromatogram_records_for`] and
18+
//! [`srm_chromatograms_for`].
1619
//!
1720
//! Frame -> spectrum projection:
1821
//!
@@ -44,7 +47,8 @@ use openmassspec_core as msc;
4447

4548
use crate::error::Result;
4649
use crate::{
47-
Calibration, DiaWindow, Frame, Metadata, PasefMsMsInfo, Peak, Precursor as TdfPrecursor, Reader,
50+
Calibration, DiaWindow, Frame, Metadata, PasefMsMsInfo, Peak, Precursor as TdfPrecursor,
51+
PrmTarget, Reader,
4852
};
4953

5054
const SOFTWARE_NAME: &str = "opentimstdf";
@@ -604,9 +608,9 @@ fn run_metadata_for(meta: &Metadata, bundle_name: &str) -> msc::RunMetadata {
604608
/// misleading zero-length chromatogram - matching the writer's contract of
605609
/// only emitting `<chromatogramList>` when the source yields something.
606610
///
607-
/// SRM/PRM transition chromatograms are intentionally not produced here: that
608-
/// needs genuine per-target cross-frame aggregation of the decoded peak stream
609-
/// (see Sigilweaver/OpenTimsTDF#25), not just wiring of existing columns.
611+
/// SRM/PRM transition chromatograms need genuine per-target cross-frame
612+
/// aggregation of the decoded peak stream rather than wiring of existing
613+
/// columns, and are built separately by [`srm_chromatograms_for`].
610614
fn chromatogram_records_for(frames: &[Frame]) -> Vec<msc::ChromatogramRecord> {
611615
fn trace(frames: &[Frame], value_of: impl Fn(&Frame) -> Option<u64>) -> (Vec<f32>, Vec<f32>) {
612616
let mut time_sec = Vec::new();
@@ -654,6 +658,135 @@ fn chromatogram_records_for(frames: &[Frame]) -> Vec<msc::ChromatogramRecord> {
654658
out
655659
}
656660

661+
/// Sum the intensity of `peaks` whose mobility scan number falls in
662+
/// `[scan_lo, scan_hi)`. Returns `0.0` (not `None`) when nothing is in
663+
/// range - unlike a missing `Frames` column, a target genuinely having no
664+
/// signal on a given prm-PASEF frame is real data, not absence of it, so
665+
/// [`srm_chromatograms_for`] keeps the point rather than dropping it.
666+
fn sum_intensity_in_range(peaks: &[Peak], scan_lo: u32, scan_hi: u32) -> f64 {
667+
peaks
668+
.iter()
669+
.filter(|p| p.scan >= scan_lo && p.scan < scan_hi)
670+
.map(|p| p.intensity as f64)
671+
.sum()
672+
}
673+
674+
/// Build one SRM-shaped [`msc::ChromatogramRecord`] for `target` from an
675+
/// already-accumulated `(time_sec, intensity)` series. `index` is left at
676+
/// `0`; callers renumber after concatenating with the TIC/BPC traces so the
677+
/// whole `<chromatogramList>` gets contiguous indices (see
678+
/// [`chromatogram_records_for_source`]).
679+
///
680+
/// `precursor_mz` is [`PrmTarget::monoisotopic_mz`] - the single scheduled
681+
/// per-target value from `PrmTargets`, stable across every frame the target
682+
/// appears in (as opposed to `PrmFrameMsMsInfo::isolation_mz`, which is
683+
/// recorded per frame and is redundant with it in the observed corpus).
684+
///
685+
/// `product_mz` is left `None`. prm-PASEF isolates one precursor per target
686+
/// and records a full fragment-ion MS2 scan over the downstream m/z range
687+
/// (`docs/docs/format/05-instrument-tables.md`'s `PrmFrameMsMsInfo` section
688+
/// documents `IsolationMz`/`IsolationWidth` as describing the *precursor*
689+
/// isolation window only) - unlike a triple-quadrupole SRM/MRM transition,
690+
/// the TDF schema records no discrete product-ion m/z to report here, and
691+
/// clean-room correctness means not inventing one.
692+
fn build_srm_record(
693+
target: &PrmTarget,
694+
time_sec: Vec<f32>,
695+
intensity: Vec<f32>,
696+
) -> msc::ChromatogramRecord {
697+
msc::ChromatogramRecord {
698+
index: 0,
699+
id: format!("SRM_{}", target.id),
700+
chromatogram_type: Some(msc::CvTerm::new(
701+
"MS:1001473",
702+
"selected reaction monitoring chromatogram",
703+
)),
704+
precursor_mz: Some(target.monoisotopic_mz),
705+
product_mz: None,
706+
time_sec,
707+
intensity,
708+
}
709+
}
710+
711+
/// Build one SRM/PRM chromatogram trace per scheduled `PrmTarget`.
712+
///
713+
/// Walks every `msms_type == 10` (prm-PASEF) frame in `frames`, in order,
714+
/// and for each one:
715+
///
716+
/// 1. Looks up that frame's `PrmFrameMsMsInfo` rows (one per active target on
717+
/// that frame) via [`Reader::prm_msms_info_for_frame`].
718+
/// 2. Decodes the frame's peaks via [`Reader::decode_peaks`] - the same
719+
/// per-frame decode path `iter_spectra` uses for every other frame kind.
720+
/// 3. For each row, sums the decoded intensity restricted to that row's
721+
/// `[scan_num_begin, scan_num_end)` mobility range
722+
/// ([`sum_intensity_in_range`]) and appends `(frame.time, sum)` to that
723+
/// target's running series.
724+
///
725+
/// A frame that fails to decode, or whose `PrmFrameMsMsInfo` lookup fails or
726+
/// comes back empty, is skipped - matching [`spectra_for_frame`]'s
727+
/// silent-skip contract for a single bad frame. A target that never appears
728+
/// in any `PrmFrameMsMsInfo` row (e.g. an unused entry in a shared target
729+
/// list) is omitted from the output entirely, rather than emitted with an
730+
/// empty trace - the same no-empty-chromatogram contract
731+
/// [`chromatogram_records_for`] applies to TIC/BPC.
732+
///
733+
/// Target iteration order is by `PrmTargets.Id` ascending (a `BTreeMap` key),
734+
/// which is deterministic but otherwise arbitrary; per-target point order is
735+
/// frame-acquisition order, since `frames` itself is walked in that order.
736+
fn srm_chromatograms_for(reader: &Reader, frames: &[Frame]) -> Vec<msc::ChromatogramRecord> {
737+
let mut per_target: std::collections::BTreeMap<u32, (Vec<f32>, Vec<f32>)> =
738+
std::collections::BTreeMap::new();
739+
740+
for frame in frames {
741+
if frame.msms_type != 10 {
742+
continue;
743+
}
744+
let Ok(infos) = reader.prm_msms_info_for_frame(frame.id) else {
745+
continue;
746+
};
747+
if infos.is_empty() {
748+
continue;
749+
}
750+
let Ok(peaks) = reader.decode_peaks(frame) else {
751+
continue;
752+
};
753+
for info in &infos {
754+
let sum = sum_intensity_in_range(&peaks, info.scan_num_begin, info.scan_num_end);
755+
let series = per_target.entry(info.target_id).or_default();
756+
series.0.push(frame.time as f32);
757+
series.1.push(sum as f32);
758+
}
759+
}
760+
761+
let mut out = Vec::with_capacity(per_target.len());
762+
for (target_id, (time_sec, intensity)) in per_target {
763+
let Ok(Some(target)) = reader.prm_target(target_id) else {
764+
continue;
765+
};
766+
out.push(build_srm_record(&target, time_sec, intensity));
767+
}
768+
out
769+
}
770+
771+
/// Full chromatogram list for [`TdfSource`]/[`OwnedTdfSource`]: TIC and BPC
772+
/// (from [`chromatogram_records_for`]) followed by one SRM/PRM trace per
773+
/// scheduled target (from [`srm_chromatograms_for`]), reindexed into one
774+
/// contiguous 0-based sequence across the combined list - mzML's
775+
/// `<chromatogramList>` expects each `<chromatogram>`'s `index` attribute to
776+
/// be unique and contiguous across the whole document, not just within one
777+
/// trace family.
778+
fn chromatogram_records_for_source(
779+
reader: &Reader,
780+
frames: &[Frame],
781+
) -> Vec<msc::ChromatogramRecord> {
782+
let mut out = chromatogram_records_for(frames);
783+
out.extend(srm_chromatograms_for(reader, frames));
784+
for (i, rec) in out.iter_mut().enumerate() {
785+
rec.index = i;
786+
}
787+
out
788+
}
789+
657790
impl<'a> msc::SpectrumSource for TdfSource<'a> {
658791
fn run_metadata(&self) -> msc::RunMetadata {
659792
run_metadata_for(&self.metadata, &self.bundle_name)
@@ -664,7 +797,7 @@ impl<'a> msc::SpectrumSource for TdfSource<'a> {
664797
fn iter_chromatograms<'s>(
665798
&'s mut self,
666799
) -> Box<dyn Iterator<Item = msc::ChromatogramRecord> + 's> {
667-
Box::new(chromatogram_records_for(&self.frames).into_iter())
800+
Box::new(chromatogram_records_for_source(self.reader, &self.frames).into_iter())
668801
}
669802
}
670803

@@ -678,7 +811,7 @@ impl msc::SpectrumSource for OwnedTdfSource {
678811
fn iter_chromatograms<'s>(
679812
&'s mut self,
680813
) -> Box<dyn Iterator<Item = msc::ChromatogramRecord> + 's> {
681-
Box::new(chromatogram_records_for(&self.frames).into_iter())
814+
Box::new(chromatogram_records_for_source(&self.reader, &self.frames).into_iter())
682815
}
683816
}
684817

@@ -1068,4 +1201,92 @@ mod tests {
10681201
];
10691202
assert!(chromatogram_records_for(&frames).is_empty());
10701203
}
1204+
1205+
#[test]
1206+
fn sum_intensity_in_range_includes_low_bound_excludes_high_bound() {
1207+
let peaks = [
1208+
Peak {
1209+
scan: 4,
1210+
tof: 100,
1211+
intensity: 10,
1212+
},
1213+
Peak {
1214+
scan: 5,
1215+
tof: 100,
1216+
intensity: 20,
1217+
},
1218+
Peak {
1219+
scan: 9,
1220+
tof: 100,
1221+
intensity: 40,
1222+
},
1223+
Peak {
1224+
scan: 10,
1225+
tof: 100,
1226+
intensity: 80,
1227+
},
1228+
];
1229+
// [5, 10): includes scan 5 and 9, excludes 4 and 10.
1230+
assert_eq!(sum_intensity_in_range(&peaks, 5, 10), 60.0);
1231+
}
1232+
1233+
#[test]
1234+
fn sum_intensity_in_range_zero_when_nothing_in_range() {
1235+
let peaks = [Peak {
1236+
scan: 1,
1237+
tof: 100,
1238+
intensity: 10,
1239+
}];
1240+
// A target with genuinely no signal on a frame gets a real 0.0, not
1241+
// a dropped point - distinct from chromatogram_records_for's
1242+
// missing-column contract.
1243+
assert_eq!(sum_intensity_in_range(&peaks, 5, 10), 0.0);
1244+
}
1245+
1246+
#[test]
1247+
fn sum_intensity_in_range_empty_peaks_is_zero() {
1248+
assert_eq!(sum_intensity_in_range(&[], 0, 100), 0.0);
1249+
}
1250+
1251+
fn sample_prm_target() -> PrmTarget {
1252+
PrmTarget {
1253+
id: 42,
1254+
external_id: "PEPTIDE_A".into(),
1255+
time: 12.3,
1256+
one_over_k0: 0.85,
1257+
monoisotopic_mz: 524.3,
1258+
charge: 2,
1259+
description: String::new(),
1260+
}
1261+
}
1262+
1263+
#[test]
1264+
fn build_srm_record_wires_up_cv_term_and_precursor_mz_no_product_mz() {
1265+
let target = sample_prm_target();
1266+
let rec = build_srm_record(&target, vec![1.0, 2.0], vec![100.0, 200.0]);
1267+
1268+
assert_eq!(rec.id, "SRM_42");
1269+
let cv = rec.chromatogram_type.as_ref().expect("SRM CV term present");
1270+
assert_eq!(cv.accession, "MS:1001473");
1271+
assert_eq!(cv.name, "selected reaction monitoring chromatogram");
1272+
assert_eq!(rec.precursor_mz, Some(524.3));
1273+
// No discrete product ion m/z exists in the TDF schema for
1274+
// prm-PASEF (full fragment-ion scan, not a triple-quad transition) -
1275+
// see build_srm_record's doc comment.
1276+
assert!(rec.product_mz.is_none());
1277+
assert_eq!(rec.time_sec, vec![1.0, 2.0]);
1278+
assert_eq!(rec.intensity, vec![100.0, 200.0]);
1279+
}
1280+
1281+
#[test]
1282+
fn build_srm_record_id_is_unique_per_target_id() {
1283+
let mut a = sample_prm_target();
1284+
a.id = 1;
1285+
let mut b = sample_prm_target();
1286+
b.id = 2;
1287+
assert_ne!(
1288+
build_srm_record(&a, vec![], vec![]).id,
1289+
build_srm_record(&b, vec![], vec![]).id
1290+
);
1291+
}
10711292
}

0 commit comments

Comments
 (0)