diff --git a/CHANGELOG.md b/CHANGELOG.md index 0bf2525..363bcd6 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,24 @@ This project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.htm ## [Unreleased] +### Added + +- mzML export now emits a `` via + `SpectrumSource::iter_chromatograms` on `TdfSource` / `OwnedTdfSource`, + which previously used the trait's default empty implementation. Two + whole-run traces are produced, one point per frame in acquisition order, + built directly from columns `Reader::frames` already SELECTs (no extra + SQL, no peak decode): a TIC ("total ion current chromatogram", + `MS:1000235`) from `Frames.SummedIntensities`, and a basepeak + chromatogram ("basepeak chromatogram", `MS:1000628`) from + `Frames.MaxIntensity`. `MaxIntensity` is a new additive column on the + `Frames` SELECT and a new `Frame::max_intensity` field (also exposed on + the Python `Frame`). A trace whose source column is absent on every frame + is omitted rather than emitted empty. SRM/PRM transition chromatograms + are intentionally left out: they need genuine per-target cross-frame + aggregation of the decoded peak stream, not just column wiring. Closes + #25. Contributed by @Nabejo. + ### Fixed - `openmassspec-core`'s declared minimum version in `Cargo.toml` was still diff --git a/crates/opentimstdf-py/src/lib.rs b/crates/opentimstdf-py/src/lib.rs index f3ff10c..3810dc0 100644 --- a/crates/opentimstdf-py/src/lib.rs +++ b/crates/opentimstdf-py/src/lib.rs @@ -80,6 +80,8 @@ struct Frame { accumulation_time: Option, #[pyo3(get)] summed_intensities: Option, + #[pyo3(get)] + max_intensity: Option, } #[pymethods] @@ -117,6 +119,7 @@ impl From for Frame { mz_calibration_id: f.mz_calibration_id, accumulation_time: f.accumulation_time, summed_intensities: f.summed_intensities, + max_intensity: f.max_intensity, } } } @@ -535,6 +538,7 @@ impl Reader { mz_calibration_id: frame.mz_calibration_id, accumulation_time: frame.accumulation_time, summed_intensities: frame.summed_intensities, + max_intensity: frame.max_intensity, }; Ok(self .locked_inner()? @@ -565,6 +569,7 @@ impl Reader { mz_calibration_id: frame.mz_calibration_id, accumulation_time: frame.accumulation_time, summed_intensities: frame.summed_intensities, + max_intensity: frame.max_intensity, }; let peaks = guard.decode_peaks(&rs_frame).map_err(to_py_err)?; let n = peaks.len(); diff --git a/crates/opentimstdf/src/codec.rs b/crates/opentimstdf/src/codec.rs index 792376c..c94c9ec 100644 --- a/crates/opentimstdf/src/codec.rs +++ b/crates/opentimstdf/src/codec.rs @@ -4,7 +4,7 @@ use crate::types::{Frame, Peak}; /// /// Column order must match the SELECT used in `Reader::frame()` and `Reader::frames()`: /// Id, Time, NumScans, NumPeaks, TimsId, ScanMode, MsMsType, -/// MzCalibration, AccumulationTime, SummedIntensities +/// MzCalibration, AccumulationTime, SummedIntensities, MaxIntensity pub(crate) fn frame_from_row(row: &rusqlite::Row<'_>) -> rusqlite::Result { Ok(Frame { id: row.get(0)?, @@ -17,6 +17,7 @@ pub(crate) fn frame_from_row(row: &rusqlite::Row<'_>) -> rusqlite::Result mz_calibration_id: row.get::<_, Option>(7)?.unwrap_or(1) as u32, accumulation_time: row.get(8)?, summed_intensities: row.get::<_, Option>(9)?.map(|v| v as u64), + max_intensity: row.get::<_, Option>(10)?.map(|v| v as u64), }) } diff --git a/crates/opentimstdf/src/mzml.rs b/crates/opentimstdf/src/mzml.rs index a7cc5c8..ff3cff8 100644 --- a/crates/opentimstdf/src/mzml.rs +++ b/crates/opentimstdf/src/mzml.rs @@ -9,6 +9,11 @@ //! * [`write_mzml`] / [`write_indexed_mzml`] - convenience entry points //! that wrap the canonical writer in `openmassspec-core`. //! +//! Chromatograms: [`iter_chromatograms`](msc::SpectrumSource::iter_chromatograms) +//! emits whole-run TIC and basepeak (BPC) traces built directly from the +//! per-frame `Frames.SummedIntensities` / `Frames.MaxIntensity` columns (one +//! point per frame). See [`chromatogram_records_for`]. +//! //! Frame -> spectrum projection: //! //! * **MS1 frames** (`msms_type == 0`): all mobility scans pooled into one @@ -569,6 +574,80 @@ fn run_metadata_for(meta: &Metadata, bundle_name: &str) -> msc::RunMetadata { } } +/// Build the whole-run chromatogram traces from the already-decoded per-frame +/// `Frames` columns carried on `frames`. +/// +/// Two traces are emitted, each a time series with one point per frame (in +/// acquisition order) built purely from columns `Reader::frames` already +/// SELECTs - no extra SQL, no peak decode: +/// +/// * **TIC** - "total ion current chromatogram" (`MS:1000235`) from +/// [`Frame::summed_intensities`] (`Frames.SummedIntensities`, the +/// accumulation-normalized per-frame intensity sum documented in +/// `docs/docs/format/01-tdf-sqlite-schema.md` and cross-checked against the +/// decoded peak sum in `tests/roundtrip.rs`). +/// * **BPC** - "basepeak chromatogram" (`MS:1000628`) from +/// [`Frame::max_intensity`] (`Frames.MaxIntensity`, the most intense peak in +/// the frame). +/// +/// `Frame::time` is retention time in seconds, which is what +/// [`msc::ChromatogramRecord::time_sec`] wants (the canonical writer converts +/// to minutes for the mzML `time array` itself). A frame whose source column +/// is absent (`None`) contributes no point; a trace with no points at all is +/// omitted rather than emitted empty, so a bundle missing a column yields no +/// misleading zero-length chromatogram - matching the writer's contract of +/// only emitting `` when the source yields something. +/// +/// SRM/PRM transition chromatograms are intentionally not produced here: that +/// needs genuine per-target cross-frame aggregation of the decoded peak stream +/// (see Sigilweaver/OpenTimsTDF#25), not just wiring of existing columns. +fn chromatogram_records_for(frames: &[Frame]) -> Vec { + fn trace(frames: &[Frame], value_of: impl Fn(&Frame) -> Option) -> (Vec, Vec) { + let mut time_sec = Vec::new(); + let mut intensity = Vec::new(); + for f in frames { + if let Some(v) = value_of(f) { + time_sec.push(f.time as f32); + intensity.push(v as f32); + } + } + (time_sec, intensity) + } + + let mut out: Vec = Vec::new(); + + let (tic_time, tic_int) = trace(frames, |f| f.summed_intensities); + if !tic_time.is_empty() { + out.push(msc::ChromatogramRecord { + index: out.len(), + id: "TIC".to_string(), + chromatogram_type: Some(msc::CvTerm::new( + "MS:1000235", + "total ion current chromatogram", + )), + precursor_mz: None, + product_mz: None, + time_sec: tic_time, + intensity: tic_int, + }); + } + + let (bpc_time, bpc_int) = trace(frames, |f| f.max_intensity); + if !bpc_time.is_empty() { + out.push(msc::ChromatogramRecord { + index: out.len(), + id: "BPC".to_string(), + chromatogram_type: Some(msc::CvTerm::new("MS:1000628", "basepeak chromatogram")), + precursor_mz: None, + product_mz: None, + time_sec: bpc_time, + intensity: bpc_int, + }); + } + + out +} + impl<'a> msc::SpectrumSource for TdfSource<'a> { fn run_metadata(&self) -> msc::RunMetadata { run_metadata_for(&self.metadata, &self.bundle_name) @@ -576,6 +655,11 @@ impl<'a> msc::SpectrumSource for TdfSource<'a> { fn iter_spectra<'s>(&'s mut self) -> Box + 's> { Box::new(frame_iter(self.reader, &self.frames, &self.calibration)) } + fn iter_chromatograms<'s>( + &'s mut self, + ) -> Box + 's> { + Box::new(chromatogram_records_for(&self.frames).into_iter()) + } } impl msc::SpectrumSource for OwnedTdfSource { @@ -585,6 +669,11 @@ impl msc::SpectrumSource for OwnedTdfSource { fn iter_spectra<'s>(&'s mut self) -> Box + 's> { Box::new(frame_iter(&self.reader, &self.frames, &self.calibration)) } + fn iter_chromatograms<'s>( + &'s mut self, + ) -> Box + 's> { + Box::new(chromatogram_records_for(&self.frames).into_iter()) + } } /// Write a `.d/` bundle to mzML in one call. @@ -630,6 +719,7 @@ mod tests { mz_calibration_id: 1, accumulation_time: None, summed_intensities: None, + max_intensity: None, } } @@ -878,4 +968,98 @@ mod tests { assert_eq!(cv.name, term_name, "wrong CV name for {name:?}"); } } + + fn frame_with_intensities( + id: u32, + time: f64, + summed_intensities: Option, + max_intensity: Option, + ) -> Frame { + Frame { + id, + time, + num_scans: 10, + num_peaks: 3, + tims_id: 0, + scan_mode: 0, + msms_type: 0, + mz_calibration_id: 1, + accumulation_time: None, + summed_intensities, + max_intensity, + } + } + + #[test] + fn chromatogram_records_for_builds_tic_and_bpc_from_frame_columns() { + let frames = [ + frame_with_intensities(1, 1.5, Some(100), Some(40)), + frame_with_intensities(2, 3.0, Some(250), Some(90)), + ]; + let recs = chromatogram_records_for(&frames); + assert_eq!(recs.len(), 2, "expected a TIC and a BPC record"); + + let tic = &recs[0]; + assert_eq!(tic.index, 0); + assert_eq!(tic.id, "TIC"); + let tic_cv = tic.chromatogram_type.as_ref().expect("TIC CV term"); + assert_eq!(tic_cv.accession, "MS:1000235"); + assert_eq!(tic_cv.name, "total ion current chromatogram"); + assert!(tic.precursor_mz.is_none()); + assert!(tic.product_mz.is_none()); + // Retention time carried through in seconds (not minutes). + assert_eq!(tic.time_sec, vec![1.5, 3.0]); + assert_eq!(tic.intensity, vec![100.0, 250.0]); + + let bpc = &recs[1]; + assert_eq!(bpc.index, 1); + assert_eq!(bpc.id, "BPC"); + let bpc_cv = bpc.chromatogram_type.as_ref().expect("BPC CV term"); + assert_eq!(bpc_cv.accession, "MS:1000628"); + assert_eq!(bpc_cv.name, "basepeak chromatogram"); + assert_eq!(bpc.time_sec, vec![1.5, 3.0]); + assert_eq!(bpc.intensity, vec![40.0, 90.0]); + } + + #[test] + fn chromatogram_records_for_skips_frames_missing_a_column_keeping_arrays_parallel() { + // Middle frame has no SummedIntensities and no MaxIntensity: it must + // contribute no point to either trace, and the two arrays stay + // parallel (one time per intensity). + let frames = [ + frame_with_intensities(1, 1.0, Some(10), Some(5)), + frame_with_intensities(2, 2.0, None, None), + frame_with_intensities(3, 3.0, Some(30), Some(15)), + ]; + let recs = chromatogram_records_for(&frames); + assert_eq!(recs.len(), 2); + assert_eq!(recs[0].time_sec, vec![1.0, 3.0]); + assert_eq!(recs[0].intensity, vec![10.0, 30.0]); + assert_eq!(recs[1].time_sec, vec![1.0, 3.0]); + assert_eq!(recs[1].intensity, vec![5.0, 15.0]); + } + + #[test] + fn chromatogram_records_for_omits_a_trace_with_no_points() { + // SummedIntensities present on every frame but MaxIntensity absent + // (e.g. an older schema without the column): only the TIC is emitted, + // and it keeps index 0 - no empty BPC record is produced. + let frames = [ + frame_with_intensities(1, 1.0, Some(10), None), + frame_with_intensities(2, 2.0, Some(20), None), + ]; + let recs = chromatogram_records_for(&frames); + assert_eq!(recs.len(), 1); + assert_eq!(recs[0].id, "TIC"); + assert_eq!(recs[0].index, 0); + } + + #[test] + fn chromatogram_records_for_empty_when_no_columns_populated() { + let frames = [ + frame_with_intensities(1, 1.0, None, None), + frame_with_intensities(2, 2.0, None, None), + ]; + assert!(chromatogram_records_for(&frames).is_empty()); + } } diff --git a/crates/opentimstdf/src/reader.rs b/crates/opentimstdf/src/reader.rs index 2fa7634..1df4753 100644 --- a/crates/opentimstdf/src/reader.rs +++ b/crates/opentimstdf/src/reader.rs @@ -236,7 +236,7 @@ impl Reader { let conn = self.conn.lock().map_err(|_| Error::LockPoisoned)?; let frame = conn.query_row( "SELECT Id, Time, NumScans, NumPeaks, TimsId, ScanMode, MsMsType, - MzCalibration, AccumulationTime, SummedIntensities + MzCalibration, AccumulationTime, SummedIntensities, MaxIntensity FROM Frames WHERE Id = ?1", [frame_id], frame_from_row, @@ -249,7 +249,7 @@ impl Reader { let conn = self.conn.lock().map_err(|_| Error::LockPoisoned)?; let mut stmt = conn.prepare( "SELECT Id, Time, NumScans, NumPeaks, TimsId, ScanMode, MsMsType, - MzCalibration, AccumulationTime, SummedIntensities + MzCalibration, AccumulationTime, SummedIntensities, MaxIntensity FROM Frames ORDER BY Id ASC", )?; let rows = stmt.query_map([], frame_from_row)?; diff --git a/crates/opentimstdf/src/types.rs b/crates/opentimstdf/src/types.rs index 8b0928a..ea6d479 100644 --- a/crates/opentimstdf/src/types.rs +++ b/crates/opentimstdf/src/types.rs @@ -18,6 +18,9 @@ pub struct Frame { /// Raw sum of decoded intensities divided by (AccumulationTime/100). /// See SPEC ยง2.2 for the normalization formula. pub summed_intensities: Option, + /// Most intense peak in the frame (`Frames.MaxIntensity`). `None` when the + /// column is absent. Used to build a basepeak chromatogram. + pub max_intensity: Option, } /// One decoded peak. diff --git a/docs/docs/format/01-tdf-sqlite-schema.md b/docs/docs/format/01-tdf-sqlite-schema.md index f6ccaa4..e3aa40d 100644 --- a/docs/docs/format/01-tdf-sqlite-schema.md +++ b/docs/docs/format/01-tdf-sqlite-schema.md @@ -50,7 +50,8 @@ Columns that matter for the binary stream: | `Time` | REAL | Retention time (seconds from start of acquisition). | | `AccumulationTime` | REAL | Ion accumulation time in milliseconds (sv >= 3.5; absent in older schemas). | | `RampTime` | REAL | TIMS ramp (elution) time in milliseconds. Equal to `AccumulationTime` in all observed bundles. | -| `SummedIntensities` | INT | Sum of all peak intensities in the frame, **normalised to 100 ms accumulation**: `SummedIntensities = sum(raw_intensity) * 100.0 / AccumulationTime_ms`. For bundles where `AccumulationTime ~ 100` the factor is ~ 1.0. Verified on PXD022216 (AT=108.46 ms): ratio decoded_sum / SummedIntensities = 1.08460 +/- 0.00001 across all 57,886 frames. | +| `SummedIntensities` | INT | Sum of all peak intensities in the frame, **normalised to 100 ms accumulation**: `SummedIntensities = sum(raw_intensity) * 100.0 / AccumulationTime_ms`. For bundles where `AccumulationTime ~ 100` the factor is ~ 1.0. Verified on PXD022216 (AT=108.46 ms): ratio decoded_sum / SummedIntensities = 1.08460 +/- 0.00001 across all 57,886 frames. Consumed as the per-frame point of the whole-run TIC chromatogram in mzML export. | +| `MaxIntensity` | INT | Most intense single peak in the frame. Consumed as the per-frame point of the whole-run basepeak chromatogram (BPC) in mzML export. | | `MzCalibration` | INT | Foreign key into `MzCalibration` table. | | `TimsCalibration` | INT | Foreign key into `TimsCalibration` table. | diff --git a/docs/docs/reference/python-api.md b/docs/docs/reference/python-api.md index 13e47c7..e8e4c55 100644 --- a/docs/docs/reference/python-api.md +++ b/docs/docs/reference/python-api.md @@ -76,6 +76,7 @@ stream. | `mz_calibration_id` | `int` | Calibration id used for this frame | | `accumulation_time` | `float \| None` | TIMS accumulation time (ms) | | `summed_intensities` | `int \| None` | Sum of all peak intensities in the frame | +| `max_intensity` | `int \| None` | Most intense peak in the frame (`Frames.MaxIntensity`) | | `polarity` | `str` | `"positive"` or `"negative"`, derived from `mz_calibration_id` | ### Metadata