Skip to main content

koth_ms/output/
report.rs

1use std::collections::BTreeMap;
2use std::path::Path;
3
4use serde::Serialize;
5
6use crate::error::KothError;
7use crate::models::{Hill, ScoredFeature};
8
9// ── helpers ──────────────────────────────────────────────────────────────────
10
11fn median(v: &mut [f64]) -> f64 {
12    if v.is_empty() {
13        return f64::NAN;
14    }
15    v.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
16    let n = v.len();
17    if n % 2 == 0 {
18        (v[n / 2 - 1] + v[n / 2]) / 2.0
19    } else {
20        v[n / 2]
21    }
22}
23
24fn stats(values: &mut [f64]) -> StatSummary {
25    StatSummary {
26        median: median(values),
27        std: crate::stats::std_dev(values),
28    }
29}
30
31// ── serialisable types ────────────────────────────────────────────────────────
32
33#[derive(Serialize)]
34pub struct StatSummary {
35    pub median: f64,
36    pub std: f64,
37}
38
39#[derive(Serialize)]
40pub struct HillsReport {
41    pub n: usize,
42    pub mz_std: StatSummary,
43    pub im_std: StatSummary,
44    pub rt_width: StatSummary,
45    pub n_scans: StatSummary,
46    pub intensity_sum: StatSummary,
47}
48
49#[derive(Serialize)]
50pub struct FeaturesReport {
51    pub n: usize,
52    /// Charge state → feature count, sorted by charge.
53    pub charge_distribution: BTreeMap<u8, usize>,
54    /// Chromatographic cosine score (mean adjacent-isotope-hill cosine).
55    pub cosine_score: StatSummary,
56    /// Bhattacharyya isotope-pattern score.
57    pub isotope_score: StatSummary,
58    /// `isotope_score × cosine_score`.
59    pub combined_score: StatSummary,
60    pub ppm_error: StatSummary,
61    pub n_isotopes: StatSummary,
62    pub rt_width: StatSummary,
63}
64
65#[derive(Serialize)]
66pub struct RunReport {
67    pub hills: HillsReport,
68    pub features: FeaturesReport,
69    /// Which 1/K0 scale the `im` columns are on. Absent from reports written
70    /// before 0.10.0, whose Bruker `im` values are on the linear scale.
71    pub mobility: MobilityReport,
72}
73
74/// Provenance of the reported ion-mobility values.
75#[derive(Serialize, Debug, Clone, PartialEq)]
76pub struct MobilityReport {
77    /// `"bruker-acquisition-calibrated-1/K0"`, `"timsrust-linear-1/K0"`, or
78    /// `"as-read-from-input"` for mzML/Thermo, where koth does not convert.
79    pub scale: String,
80    /// The `TimsCalibration` rows applied (calibrated Bruker input only).
81    #[serde(skip_serializing_if = "Vec::is_empty")]
82    pub calibration: Vec<CalibrationRow>,
83}
84
85#[derive(Serialize, Debug, Clone, PartialEq)]
86pub struct CalibrationRow {
87    pub id: i64,
88    pub model_type: i64,
89    /// `C0..C9` as stored in `analysis.tdf`.
90    pub coefficients: [f64; 10],
91}
92
93impl MobilityReport {
94    /// Non-Bruker input: ion mobility is whatever the input file carries.
95    pub fn as_read() -> Self {
96        Self {
97            scale: crate::input::SCALE_AS_READ.into(),
98            calibration: Vec::new(),
99        }
100    }
101
102    /// Bruker input read on `conv`'s scale.
103    #[cfg(feature = "tdf")]
104    pub fn bruker(conv: &crate::io::tims_calibration::ScanToMobility) -> Self {
105        Self {
106            scale: conv.scale().label().into(),
107            calibration: conv
108                .calibration_rows()
109                .iter()
110                .map(|&(id, coefficients)| CalibrationRow {
111                    id,
112                    model_type: crate::io::tims_calibration::SUPPORTED_MODEL_TYPE,
113                    coefficients,
114                })
115                .collect(),
116        }
117    }
118}
119
120// ── builders ──────────────────────────────────────────────────────────────────
121
122pub fn build_hills_report(hills: &[Hill]) -> HillsReport {
123    let mut mz_stds = hills.iter().map(|h| h.mz_std).collect::<Vec<_>>();
124    let mut im_stds = hills.iter().map(|h| h.im_std).collect::<Vec<_>>();
125    let mut rt_widths = hills.iter().map(|h| h.rt_width).collect::<Vec<_>>();
126    let mut n_scans = hills.iter().map(|h| h.n_scans as f64).collect::<Vec<_>>();
127    let mut int_sums = hills.iter().map(|h| h.intensity_sum).collect::<Vec<_>>();
128
129    HillsReport {
130        n: hills.len(),
131        mz_std: stats(&mut mz_stds),
132        im_std: stats(&mut im_stds),
133        rt_width: stats(&mut rt_widths),
134        n_scans: stats(&mut n_scans),
135        intensity_sum: stats(&mut int_sums),
136    }
137}
138
139pub fn build_features_report(features: &[ScoredFeature]) -> FeaturesReport {
140    let scored: Vec<&ScoredFeature> = features.iter().filter(|f| f.feature.charge > 0).collect();
141
142    let mut charge_dist: BTreeMap<u8, usize> = BTreeMap::new();
143    for sf in &scored {
144        *charge_dist.entry(sf.feature.charge).or_insert(0) += 1;
145    }
146
147    let mut cosine = scored.iter().map(|sf| sf.cosine_score).collect::<Vec<_>>();
148    let mut isotope = scored.iter().map(|sf| sf.isotope_score).collect::<Vec<_>>();
149    let mut combined = scored
150        .iter()
151        .map(|sf| sf.combined_score)
152        .collect::<Vec<_>>();
153    let mut ppm = scored
154        .iter()
155        .map(|sf| sf.feature.ppm_error)
156        .collect::<Vec<_>>();
157    let mut n_iso = scored
158        .iter()
159        .map(|sf| sf.feature.hills.len() as f64)
160        .collect::<Vec<_>>();
161    let mut rt_w = scored
162        .iter()
163        .map(|sf| sf.feature.rt_end() - sf.feature.rt_start())
164        .collect::<Vec<_>>();
165
166    FeaturesReport {
167        n: scored.len(),
168        charge_distribution: charge_dist,
169        cosine_score: stats(&mut cosine),
170        isotope_score: stats(&mut isotope),
171        combined_score: stats(&mut combined),
172        ppm_error: stats(&mut ppm),
173        n_isotopes: stats(&mut n_iso),
174        rt_width: stats(&mut rt_w),
175    }
176}
177
178// ── writer ────────────────────────────────────────────────────────────────────
179
180pub fn write_report(report: &RunReport, path: &Path) -> Result<(), KothError> {
181    let json = serde_json::to_string_pretty(report)?;
182    std::fs::write(path, json)?;
183    log::info!("Wrote report to {}", path.display());
184    Ok(())
185}