Skip to main content

koth_ms/
input.rs

1/// Readers for hill and feature files produced by `koth_ff`.
2///
3/// Supports both TSV (default) and Parquet output formats.
4/// Auto-detection is based on the file extension.
5use std::path::Path;
6use std::sync::Arc;
7
8use serde::Deserialize;
9
10use crate::error::KothError;
11use crate::models::{Feature, Hill, IsolationWindow, ScoredFeature};
12
13const C13_NEUTRON: f64 = 1.003_354_835;
14
15// ──────────────────────────────────────────────────────────────────────────────
16// Public API
17// ──────────────────────────────────────────────────────────────────────────────
18
19/// Read hills from a file, detecting format from the extension (.tsv / .parquet).
20pub fn read_hills(path: &Path) -> Result<Vec<Hill>, KothError> {
21    match path.extension().and_then(|e| e.to_str()) {
22        Some("parquet") => read_hills_parquet(path),
23        _ => read_hills_tsv(path),
24    }
25}
26
27/// Read scored features from a file, detecting format from extension.
28pub fn read_features(path: &Path) -> Result<Vec<ScoredFeature>, KothError> {
29    match path.extension().and_then(|e| e.to_str()) {
30        Some("parquet") => read_features_parquet(path),
31        _ => read_features_tsv(path),
32    }
33}
34
35// ──────────────────────────────────────────────────────────────────────────────
36// TSV readers
37// ──────────────────────────────────────────────────────────────────────────────
38
39#[derive(Debug, Deserialize)]
40struct HillRow {
41    #[serde(default, deserialize_with = "de_opt_u64")]
42    hill_id: Option<u64>,
43    mz: f64,
44    mz_std: f64,
45    #[serde(default)]
46    mz_se: f64,
47    rt: f64,
48    rt_start: f64,
49    rt_end: f64,
50    rt_width: f64,
51    im: f64,
52    im_std: f64,
53    scan_start: usize,
54    scan_apex: usize,
55    scan_end: usize,
56    n_scans: usize,
57    skipped_scans: usize,
58    intensity_sum: f64,
59    intensity_max: f64,
60    hill_score: f64,
61    intensity_profile: String,
62    #[serde(default, deserialize_with = "de_opt_float")]
63    iso_target_mz: Option<f64>,
64    #[serde(default, deserialize_with = "de_opt_float")]
65    iso_lower_mz: Option<f64>,
66    #[serde(default, deserialize_with = "de_opt_float")]
67    iso_upper_mz: Option<f64>,
68}
69
70fn de_opt_float<'de, D: serde::Deserializer<'de>>(de: D) -> Result<Option<f64>, D::Error> {
71    let s: Option<String> = Option::deserialize(de)?;
72    match s {
73        Some(s) if !s.is_empty() => s.parse::<f64>().map(Some).map_err(serde::de::Error::custom),
74        _ => Ok(None),
75    }
76}
77
78fn de_opt_u64<'de, D: serde::Deserializer<'de>>(de: D) -> Result<Option<u64>, D::Error> {
79    let s: Option<String> = Option::deserialize(de)?;
80    match s {
81        Some(s) if !s.is_empty() => s.parse::<u64>().map(Some).map_err(serde::de::Error::custom),
82        _ => Ok(None),
83    }
84}
85
86pub fn read_hills_tsv(path: &Path) -> Result<Vec<Hill>, KothError> {
87    let mut rdr = csv::ReaderBuilder::new().delimiter(b'\t').from_path(path)?;
88    let mut hills = Vec::new();
89    for (row_idx, result) in rdr.deserialize::<HillRow>().enumerate() {
90        let row = result?;
91        let profile: Vec<f32> = serde_json::from_str(&row.intensity_profile)?;
92        let isolation_window = match (row.iso_target_mz, row.iso_lower_mz, row.iso_upper_mz) {
93            (Some(t), Some(l), Some(u)) => Some(IsolationWindow {
94                target: t,
95                lower: l,
96                upper: u,
97            }),
98            _ => None,
99        };
100        hills.push(Hill {
101            hill_id: row.hill_id.unwrap_or(row_idx as u64),
102            mz: row.mz,
103            mz_std: row.mz_std,
104            mz_se: row.mz_se,
105            rt: row.rt,
106            rt_start: row.rt_start,
107            rt_end: row.rt_end,
108            rt_width: row.rt_width,
109            im: row.im,
110            im_std: row.im_std,
111            scan_start: row.scan_start,
112            scan_apex: row.scan_apex,
113            scan_end: row.scan_end,
114            n_scans: row.n_scans,
115            skipped_scans: row.skipped_scans,
116            intensity_sum: row.intensity_sum,
117            intensity_max: row.intensity_max,
118            hill_score: row.hill_score,
119            intensity_profile: Arc::from(profile.as_slice()),
120            isolation_window,
121            faims_cv: None,
122        });
123    }
124    Ok(hills)
125}
126
127#[derive(Debug, Deserialize)]
128#[allow(dead_code)]
129struct FeatureRow {
130    #[serde(rename = "massCalib")]
131    mass_calib: String,
132    mz: f64,
133    #[serde(rename = "rtApex")]
134    rt_apex: f64,
135    #[serde(rename = "rtStart")]
136    rt_start: f64,
137    #[serde(rename = "rtEnd")]
138    rt_end: f64,
139    #[serde(rename = "FAIMS", default, deserialize_with = "de_opt_float")]
140    faims_cv: Option<f64>,
141    #[serde(rename = "intensityApex")]
142    intensity_apex: f64,
143    #[serde(rename = "intensitySum")]
144    intensity_sum: f64,
145    charge: u8,
146    #[serde(rename = "nIsotopes")]
147    n_isotopes: usize,
148    #[serde(rename = "nScans")]
149    n_scans: usize,
150    #[serde(deserialize_with = "de_optional_float")]
151    im: f64,
152    cosine_score: f64,
153    ppm_error: f64,
154    neutron_offset: i8,
155    isotope_score: f64,
156    combined_score: f64,
157    theoretical_pattern: String,
158    // remaining columns are not needed for alignment/LFQ
159    isotope_profile: String,
160    elution_profile: String,
161}
162
163pub fn read_features_tsv(path: &Path) -> Result<Vec<ScoredFeature>, KothError> {
164    let mut rdr = csv::ReaderBuilder::new().delimiter(b'\t').from_path(path)?;
165    let mut features = Vec::new();
166    for result in rdr.deserialize::<FeatureRow>() {
167        let row = result?;
168        if row.charge == 0 {
169            continue;
170        }
171        let theoretical_pattern: Vec<f64> = serde_json::from_str(&row.theoretical_pattern)?;
172        features.push(scored_feature_from_row(row, theoretical_pattern));
173    }
174    Ok(features)
175}
176
177fn scored_feature_from_row(row: FeatureRow, theoretical_pattern: Vec<f64>) -> ScoredFeature {
178    // The stored `mz` is the calibrated monoisotopic mz (neutron-offset corrected).
179    // To reconstruct the raw hill mz: hill.mz = mz + offset * C13/charge
180    let raw_mz = row.mz + row.neutron_offset as f64 * C13_NEUTRON / row.charge as f64;
181
182    let hill = Hill {
183        hill_id: 0,
184        mz: raw_mz,
185        mz_std: 0.0,
186        mz_se: 0.0,
187        rt: row.rt_apex,
188        rt_start: row.rt_start,
189        rt_end: row.rt_end,
190        rt_width: row.rt_end - row.rt_start,
191        im: row.im,
192        im_std: 0.0,
193        scan_start: 0,
194        scan_apex: 0,
195        scan_end: 0,
196        n_scans: row.n_scans,
197        skipped_scans: 0,
198        intensity_sum: row.intensity_sum,
199        intensity_max: row.intensity_apex,
200        hill_score: 1.0,
201        intensity_profile: Arc::from(&[] as &[f32]),
202        isolation_window: None,
203        faims_cv: row.faims_cv.map(|v| v as f32),
204    };
205
206    ScoredFeature {
207        feature: Feature {
208            hills: vec![hill],
209            charge: row.charge,
210            cosine_score: row.cosine_score,
211            ppm_error: row.ppm_error,
212        },
213        neutron_offset: row.neutron_offset,
214        isotope_score: row.isotope_score,
215        cosine_score: row.cosine_score,
216        combined_score: row.combined_score,
217        theoretical_pattern,
218    }
219}
220
221fn de_optional_float<'de, D: serde::Deserializer<'de>>(de: D) -> Result<f64, D::Error> {
222    let s = String::deserialize(de)?;
223    if s.is_empty() {
224        Ok(0.0)
225    } else {
226        s.parse::<f64>().map_err(serde::de::Error::custom)
227    }
228}
229
230// ──────────────────────────────────────────────────────────────────────────────
231// Parquet readers
232// ──────────────────────────────────────────────────────────────────────────────
233
234pub fn read_hills_parquet(path: &Path) -> Result<Vec<Hill>, KothError> {
235    use arrow::array::{Array, Float64Array, Int64Array, StringArray};
236    use parquet::arrow::arrow_reader::ParquetRecordBatchReaderBuilder;
237
238    let file = std::fs::File::open(path)?;
239    let reader = ParquetRecordBatchReaderBuilder::try_new(file)
240        .map_err(|e| KothError::ParquetError(e.to_string()))?
241        .build()
242        .map_err(|e| KothError::ParquetError(e.to_string()))?;
243
244    let mut hills = Vec::new();
245
246    for batch in reader {
247        let batch = batch.map_err(|e| KothError::ParquetError(e.to_string()))?;
248
249        macro_rules! f64_col {
250            ($name:expr) => {
251                batch
252                    .column_by_name($name)
253                    .ok_or_else(|| KothError::ParquetError(format!("missing column: {}", $name)))?
254                    .as_any()
255                    .downcast_ref::<Float64Array>()
256                    .ok_or_else(|| KothError::ParquetError(format!("wrong type: {}", $name)))?
257            };
258        }
259        macro_rules! i64_col {
260            ($name:expr) => {
261                batch
262                    .column_by_name($name)
263                    .ok_or_else(|| KothError::ParquetError(format!("missing column: {}", $name)))?
264                    .as_any()
265                    .downcast_ref::<Int64Array>()
266                    .ok_or_else(|| KothError::ParquetError(format!("wrong type: {}", $name)))?
267            };
268        }
269        macro_rules! str_col {
270            ($name:expr) => {
271                batch
272                    .column_by_name($name)
273                    .ok_or_else(|| KothError::ParquetError(format!("missing column: {}", $name)))?
274                    .as_any()
275                    .downcast_ref::<StringArray>()
276                    .ok_or_else(|| KothError::ParquetError(format!("wrong type: {}", $name)))?
277            };
278        }
279
280        let mz_col = f64_col!("mz");
281        let mz_std_col = f64_col!("mz_std");
282        // mz_se was added later; tolerate older parquets by defaulting to 0.
283        let mz_se_col = batch
284            .column_by_name("mz_se")
285            .and_then(|c| c.as_any().downcast_ref::<Float64Array>());
286        let rt_col = f64_col!("rt");
287        let rt_start_col = f64_col!("rt_start");
288        let rt_end_col = f64_col!("rt_end");
289        let rt_width_col = f64_col!("rt_width");
290        let im_col = f64_col!("im");
291        let im_std_col = f64_col!("im_std");
292        let scan_start_col = i64_col!("scan_start");
293        let scan_apex_col = i64_col!("scan_apex");
294        let scan_end_col = i64_col!("scan_end");
295        let n_scans_col = i64_col!("n_scans");
296        let skipped_col = i64_col!("skipped_scans");
297        let int_sum_col = f64_col!("intensity_sum");
298        let int_max_col = f64_col!("intensity_max");
299        let score_col = f64_col!("hill_score");
300        let profile_col = str_col!("intensity_profile");
301
302        let iso_target_col = batch
303            .column_by_name("iso_target_mz")
304            .and_then(|c| c.as_any().downcast_ref::<Float64Array>());
305        let iso_lower_col = batch
306            .column_by_name("iso_lower_mz")
307            .and_then(|c| c.as_any().downcast_ref::<Float64Array>());
308        let iso_upper_col = batch
309            .column_by_name("iso_upper_mz")
310            .and_then(|c| c.as_any().downcast_ref::<Float64Array>());
311        let hill_id_col = batch
312            .column_by_name("hill_id")
313            .and_then(|c| c.as_any().downcast_ref::<arrow::array::UInt64Array>());
314
315        let row_offset = hills.len() as u64;
316        for i in 0..batch.num_rows() {
317            let profile: Vec<f32> = serde_json::from_str(profile_col.value(i))?;
318            let isolation_window = match (iso_target_col, iso_lower_col, iso_upper_col) {
319                (Some(t), Some(l), Some(u)) if !t.is_null(i) && !l.is_null(i) && !u.is_null(i) => {
320                    Some(IsolationWindow {
321                        target: t.value(i),
322                        lower: l.value(i),
323                        upper: u.value(i),
324                    })
325                }
326                _ => None,
327            };
328            let hill_id = hill_id_col
329                .filter(|c| !c.is_null(i))
330                .map(|c| c.value(i))
331                .unwrap_or(row_offset + i as u64);
332            hills.push(Hill {
333                hill_id,
334                mz: mz_col.value(i),
335                mz_std: mz_std_col.value(i),
336                mz_se: mz_se_col.map(|c| c.value(i)).unwrap_or(0.0),
337                rt: rt_col.value(i),
338                rt_start: rt_start_col.value(i),
339                rt_end: rt_end_col.value(i),
340                rt_width: rt_width_col.value(i),
341                im: im_col.value(i),
342                im_std: im_std_col.value(i),
343                scan_start: scan_start_col.value(i) as usize,
344                scan_apex: scan_apex_col.value(i) as usize,
345                scan_end: scan_end_col.value(i) as usize,
346                n_scans: n_scans_col.value(i) as usize,
347                skipped_scans: skipped_col.value(i) as usize,
348                intensity_sum: int_sum_col.value(i),
349                intensity_max: int_max_col.value(i),
350                hill_score: score_col.value(i),
351                intensity_profile: Arc::from(profile.as_slice()),
352                isolation_window,
353                faims_cv: None,
354            });
355        }
356    }
357
358    Ok(hills)
359}
360
361pub fn read_features_parquet(path: &Path) -> Result<Vec<ScoredFeature>, KothError> {
362    use arrow::array::{Array, Float64Array, Int64Array, Int8Array, StringArray, UInt8Array};
363    use parquet::arrow::arrow_reader::ParquetRecordBatchReaderBuilder;
364
365    let file = std::fs::File::open(path)?;
366    let reader = ParquetRecordBatchReaderBuilder::try_new(file)
367        .map_err(|e| KothError::ParquetError(e.to_string()))?
368        .build()
369        .map_err(|e| KothError::ParquetError(e.to_string()))?;
370
371    let mut features = Vec::new();
372
373    for batch in reader {
374        let batch = batch.map_err(|e| KothError::ParquetError(e.to_string()))?;
375
376        macro_rules! f64_col {
377            ($name:expr) => {
378                batch
379                    .column_by_name($name)
380                    .ok_or_else(|| KothError::ParquetError(format!("missing column: {}", $name)))?
381                    .as_any()
382                    .downcast_ref::<Float64Array>()
383                    .ok_or_else(|| KothError::ParquetError(format!("wrong type: {}", $name)))?
384            };
385        }
386
387        let mz_col = f64_col!("mz");
388        let rt_apex_col = f64_col!("rtApex");
389        let rt_start_col = f64_col!("rtStart");
390        let rt_end_col = f64_col!("rtEnd");
391        let int_apex_col = f64_col!("intensityApex");
392        let int_sum_col = f64_col!("intensitySum");
393        let im_col = f64_col!("im");
394        let cosine_col = f64_col!("cosine_score");
395        let ppm_col = f64_col!("ppm_error");
396        let isotope_col = f64_col!("isotope_score");
397        let combined_col = f64_col!("combined_score");
398        let faims_col = batch
399            .column_by_name("FAIMS")
400            .and_then(|c| c.as_any().downcast_ref::<Float64Array>());
401
402        let charge_col = batch
403            .column_by_name("charge")
404            .ok_or_else(|| KothError::ParquetError("missing column: charge".into()))?
405            .as_any()
406            .downcast_ref::<UInt8Array>()
407            .ok_or_else(|| KothError::ParquetError("wrong type: charge".into()))?;
408
409        let n_scans_col = batch
410            .column_by_name("nScans")
411            .ok_or_else(|| KothError::ParquetError("missing column: nScans".into()))?
412            .as_any()
413            .downcast_ref::<Int64Array>()
414            .ok_or_else(|| KothError::ParquetError("wrong type: nScans".into()))?;
415
416        let offset_col = batch
417            .column_by_name("neutron_offset")
418            .ok_or_else(|| KothError::ParquetError("missing column: neutron_offset".into()))?
419            .as_any()
420            .downcast_ref::<Int8Array>()
421            .ok_or_else(|| KothError::ParquetError("wrong type: neutron_offset".into()))?;
422
423        let theo_col = batch
424            .column_by_name("theoretical_pattern")
425            .ok_or_else(|| KothError::ParquetError("missing column: theoretical_pattern".into()))?
426            .as_any()
427            .downcast_ref::<StringArray>()
428            .ok_or_else(|| KothError::ParquetError("wrong type: theoretical_pattern".into()))?;
429
430        for i in 0..batch.num_rows() {
431            let charge = charge_col.value(i);
432            if charge == 0 {
433                continue;
434            }
435            let neutron_offset = offset_col.value(i);
436            let mz = mz_col.value(i);
437            let raw_mz = mz + neutron_offset as f64 * C13_NEUTRON / charge as f64;
438            let theoretical_pattern: Vec<f64> = serde_json::from_str(theo_col.value(i))?;
439
440            let hill = Hill {
441                hill_id: 0,
442                mz: raw_mz,
443                mz_std: 0.0,
444                mz_se: 0.0,
445                rt: rt_apex_col.value(i),
446                rt_start: rt_start_col.value(i),
447                rt_end: rt_end_col.value(i),
448                rt_width: rt_end_col.value(i) - rt_start_col.value(i),
449                im: im_col.value(i),
450                im_std: 0.0,
451                scan_start: 0,
452                scan_apex: 0,
453                scan_end: 0,
454                n_scans: n_scans_col.value(i) as usize,
455                skipped_scans: 0,
456                intensity_sum: int_sum_col.value(i),
457                intensity_max: int_apex_col.value(i),
458                hill_score: 1.0,
459                intensity_profile: Arc::from(&[] as &[f32]),
460                isolation_window: None,
461                faims_cv: faims_col
462                    .filter(|c| !c.is_null(i))
463                    .map(|c| c.value(i) as f32),
464            };
465
466            features.push(ScoredFeature {
467                feature: Feature {
468                    hills: vec![hill],
469                    charge,
470                    cosine_score: cosine_col.value(i),
471                    ppm_error: ppm_col.value(i),
472                },
473                neutron_offset,
474                isotope_score: isotope_col.value(i),
475                cosine_score: cosine_col.value(i),
476                combined_score: combined_col.value(i),
477                theoretical_pattern,
478            });
479        }
480    }
481
482    Ok(features)
483}
484
485// ──────────────────────────────────────────────────────────────────────────────
486// Batch discovery
487// ──────────────────────────────────────────────────────────────────────────────
488
489/// A discovered run: the name plus paths to its hills and features files.
490pub struct RunPaths {
491    pub name: String,
492    pub hills_path: std::path::PathBuf,
493    pub features_path: std::path::PathBuf,
494}
495
496/// Scan `batch_dir` for subdirectories that contain both a hills file and a
497/// features file.  Returns one `RunPaths` per discovered run, sorted by name.
498pub fn discover_runs(batch_dir: &Path) -> Result<Vec<RunPaths>, KothError> {
499    let mut runs = Vec::new();
500
501    let entries = std::fs::read_dir(batch_dir)?;
502    for entry in entries {
503        let entry = entry?;
504        let meta = entry.metadata()?;
505        if !meta.is_dir() {
506            continue;
507        }
508
509        let run_dir = entry.path();
510        let name = run_dir
511            .file_name()
512            .and_then(|n| n.to_str())
513            .unwrap_or("unknown")
514            .to_string();
515
516        let hills_path = find_file(&run_dir, "hills")?;
517        let features_path = find_file(&run_dir, "features")?;
518
519        match (hills_path, features_path) {
520            (Some(h), Some(f)) => runs.push(RunPaths {
521                name,
522                hills_path: h,
523                features_path: f,
524            }),
525            _ => {
526                log::debug!(
527                    "Skipping '{}': no hills+features pair found",
528                    run_dir.display()
529                );
530            }
531        }
532    }
533
534    runs.sort_by(|a, b| a.name.cmp(&b.name));
535    Ok(runs)
536}
537
538/// Find `<stem>.parquet` or `<stem>.tsv` in a directory (parquet preferred).
539fn find_file(dir: &Path, stem: &str) -> Result<Option<std::path::PathBuf>, KothError> {
540    for ext in &["parquet", "tsv"] {
541        let p = dir.join(format!("{}.{}", stem, ext));
542        if p.exists() {
543            return Ok(Some(p));
544        }
545    }
546    Ok(None)
547}
548
549// ──────────────────────────────────────────────────────────────────────────────
550// Ion-mobility scale guard for cross-run consumers
551// ──────────────────────────────────────────────────────────────────────────────
552
553pub use koth_core::config::{SCALE_CALIBRATED, SCALE_LINEAR};
554/// `report.json` label for input koth did not convert (mzML, Thermo).
555pub const SCALE_AS_READ: &str = "as-read-from-input";
556
557/// What one run's outputs say about their 1/K0 scale.
558#[derive(Debug, Clone, PartialEq)]
559pub struct RunMobility {
560    pub name: String,
561    /// `mobility.scale` from the run's `report.json`; `None` when the report or
562    /// the field is absent (koth 0.9.0 and earlier, or a hand-assembled batch).
563    pub recorded: Option<String>,
564    /// Whether any feature carries a nonzero ion mobility.
565    pub has_mobility: bool,
566}
567
568/// The scale a run's `im` values are on, for the purpose of mixing runs.
569#[derive(Debug, Clone, Copy, PartialEq, Eq)]
570enum EffectiveScale {
571    /// No mobility to compare (mzML/Thermo input, or no IM values at all).
572    NotApplicable,
573    Calibrated,
574    Linear {
575        recorded: bool,
576    },
577}
578
579impl RunMobility {
580    fn effective(&self) -> Result<EffectiveScale, KothError> {
581        match self.recorded.as_deref() {
582            Some(SCALE_CALIBRATED) => Ok(EffectiveScale::Calibrated),
583            Some(SCALE_LINEAR) => Ok(EffectiveScale::Linear { recorded: true }),
584            Some(SCALE_AS_READ) => Ok(EffectiveScale::NotApplicable),
585            Some(other) => Err(KothError::ConfigError(format!(
586                "run '{}': unknown ion-mobility scale '{other}' in report.json",
587                self.name
588            ))),
589            // Written before the scale was recorded: timsTOF output from koth
590            // 0.9.0 or earlier, which only had the linear scale.
591            None if self.has_mobility => Ok(EffectiveScale::Linear { recorded: false }),
592            None => Ok(EffectiveScale::NotApplicable),
593        }
594    }
595}
596
597/// Read `mobility.scale` from the `report.json` beside a run's features file.
598/// A missing report or field is `None`; an unreadable report is an error.
599pub fn recorded_mobility_scale(features_path: &Path) -> Result<Option<String>, KothError> {
600    let Some(dir) = features_path.parent() else {
601        return Ok(None);
602    };
603    let report = dir.join("report.json");
604    if !report.exists() {
605        return Ok(None);
606    }
607    let text = std::fs::read_to_string(&report)?;
608    let value: serde_json::Value = serde_json::from_str(&text)
609        .map_err(|e| KothError::ConfigError(format!("cannot parse {}: {e}", report.display())))?;
610    Ok(value
611        .get("mobility")
612        .and_then(|m| m.get("scale"))
613        .and_then(|s| s.as_str())
614        .map(str::to_owned))
615}
616
617/// Describe one run's scale from its paths and loaded features.
618pub fn run_mobility(
619    paths: &RunPaths,
620    features: &[ScoredFeature],
621) -> Result<RunMobility, KothError> {
622    Ok(RunMobility {
623        name: paths.name.clone(),
624        recorded: recorded_mobility_scale(&paths.features_path)?,
625        has_mobility: features.iter().any(|f| f.feature.im_apex() != 0.0),
626    })
627}
628
629/// Refuse a batch whose timsTOF runs report 1/K0 on different scales.
630///
631/// Calibrated and linear 1/K0 differ by up to ~0.03 (about 2%), more than the
632/// alignment and extraction IM windows, so mixing them silently mis-matches
633/// features across runs. Runs without ion mobility are ignored. A timsTOF run
634/// with no recorded scale is treated as linear (koth ≤ 0.9.0).
635pub fn check_mobility_scales(runs: &[RunMobility]) -> Result<(), KothError> {
636    let mut calibrated = Vec::new();
637    let mut linear = Vec::new();
638    for run in runs {
639        match run.effective()? {
640            EffectiveScale::NotApplicable => {}
641            EffectiveScale::Calibrated => {
642                calibrated.push(format!("{} ({SCALE_CALIBRATED})", run.name))
643            }
644            EffectiveScale::Linear { recorded: true } => {
645                linear.push(format!("{} ({SCALE_LINEAR})", run.name))
646            }
647            EffectiveScale::Linear { recorded: false } => linear.push(format!(
648                "{} (no scale recorded: koth 0.9.0 or earlier, linear 1/K0)",
649                run.name
650            )),
651        }
652    }
653    if calibrated.is_empty() || linear.is_empty() {
654        return Ok(());
655    }
656    Err(KothError::ConfigError(format!(
657        "runs report ion mobility on different 1/K0 scales and cannot be aligned together.\n  \
658         calibrated: {}\n  linear: {}\n\
659         Re-run feature detection on the linear runs with koth 0.10.0 or later (calibrated \
660         is the default), or run every input with `[file] bruker_mobility_scale = \"linear\"`.",
661        calibrated.join(", "),
662        linear.join(", ")
663    )))
664}
665
666/// [`run_mobility`] for every run, then [`check_mobility_scales`].
667pub fn check_batch_mobility_scales<'a>(
668    runs: impl IntoIterator<Item = (&'a RunPaths, &'a [ScoredFeature])>,
669) -> Result<(), KothError> {
670    let described = runs
671        .into_iter()
672        .map(|(p, f)| run_mobility(p, f))
673        .collect::<Result<Vec<_>, _>>()?;
674    check_mobility_scales(&described)
675}
676
677#[cfg(test)]
678mod mobility_scale_tests {
679    use super::*;
680
681    fn run(name: &str, recorded: Option<&str>, has_mobility: bool) -> RunMobility {
682        RunMobility {
683            name: name.into(),
684            recorded: recorded.map(str::to_owned),
685            has_mobility,
686        }
687    }
688
689    #[test]
690    fn all_calibrated_passes() {
691        let runs = [
692            run("a", Some(SCALE_CALIBRATED), true),
693            run("b", Some(SCALE_CALIBRATED), true),
694        ];
695        assert!(check_mobility_scales(&runs).is_ok());
696    }
697
698    #[test]
699    fn all_linear_passes_recorded_or_not() {
700        let runs = [
701            run("a", Some(SCALE_LINEAR), true),
702            run("b", None, true),
703            run("c", None, true),
704        ];
705        assert!(check_mobility_scales(&runs).is_ok());
706    }
707
708    #[test]
709    fn mixed_scales_fail_and_name_the_runs() {
710        let runs = [
711            run("new_a", Some(SCALE_CALIBRATED), true),
712            run("lin_b", Some(SCALE_LINEAR), true),
713        ];
714        let msg = check_mobility_scales(&runs).unwrap_err().to_string();
715        assert!(
716            msg.contains("new_a (bruker-acquisition-calibrated-1/K0)"),
717            "{msg}"
718        );
719        assert!(msg.contains("lin_b (timsrust-linear-1/K0)"), "{msg}");
720        assert!(msg.contains("bruker_mobility_scale = \"linear\""), "{msg}");
721    }
722
723    #[test]
724    fn missing_metadata_is_linear_and_conflicts_with_calibrated() {
725        let runs = [
726            run("new_a", Some(SCALE_CALIBRATED), true),
727            run("old_b", None, true),
728        ];
729        let msg = check_mobility_scales(&runs).unwrap_err().to_string();
730        assert!(msg.contains("old_b (no scale recorded"), "{msg}");
731        assert!(msg.contains("Re-run"), "{msg}");
732    }
733
734    #[test]
735    fn runs_without_mobility_are_ignored() {
736        let runs = [
737            run("tims", Some(SCALE_CALIBRATED), true),
738            run("mzml", Some(SCALE_AS_READ), true),
739            run("orbitrap_old", None, false),
740            run("tims_lin_no_im", Some(SCALE_AS_READ), false),
741        ];
742        assert!(check_mobility_scales(&runs).is_ok());
743        let only_old = [
744            run("orbi", None, false),
745            run("tims", Some(SCALE_LINEAR), true),
746        ];
747        assert!(check_mobility_scales(&only_old).is_ok());
748    }
749
750    #[test]
751    fn unknown_label_is_an_error() {
752        let runs = [run("a", Some("furlongs"), true)];
753        assert!(check_mobility_scales(&runs).is_err());
754    }
755
756    #[test]
757    fn labels_match_the_writer() {
758        use crate::io::tims_calibration::MobilityScale;
759        assert_eq!(MobilityScale::Calibrated.label(), SCALE_CALIBRATED);
760        assert_eq!(MobilityScale::Linear.label(), SCALE_LINEAR);
761        assert_eq!(
762            crate::output::MobilityReport::as_read().scale,
763            SCALE_AS_READ
764        );
765    }
766
767    #[test]
768    fn reads_scale_from_report_json_and_tolerates_absence() {
769        let dir =
770            std::env::temp_dir().join(format!("koth_mobility_scale_test_{}", std::process::id()));
771        let _ = std::fs::remove_dir_all(&dir);
772        std::fs::create_dir_all(&dir).unwrap();
773        let features = dir.join("features.parquet");
774        assert_eq!(recorded_mobility_scale(&features).unwrap(), None);
775        std::fs::write(dir.join("report.json"), r#"{"hills":{},"features":{}}"#).unwrap();
776        assert_eq!(recorded_mobility_scale(&features).unwrap(), None);
777        std::fs::write(
778            dir.join("report.json"),
779            format!(r#"{{"mobility":{{"scale":"{SCALE_CALIBRATED}"}}}}"#),
780        )
781        .unwrap();
782        assert_eq!(
783            recorded_mobility_scale(&features).unwrap().as_deref(),
784            Some(SCALE_CALIBRATED)
785        );
786        std::fs::write(dir.join("report.json"), "not json").unwrap();
787        assert!(recorded_mobility_scale(&features).is_err());
788        std::fs::remove_dir_all(&dir).unwrap();
789    }
790}