feat: ion mobility model

added mobility to outputs and model details

simplified model

removed redundant logging

changed feature engineering

removed old comments

removed resolved comment

Author: jspaezp

crates/sage-cli/src/main.rs

@@ -275,6 +275,8 @@ impl Runner {
                 self.parameters.mzml_paths.len(),
             );
             let _ = sage_core::ml::retention_model::predict(&self.database, &mut outputs.features);
+            let _ =
+                sage_core::ml::mobility_model::predict(&self.database, &mut outputs.features);
             Some(alignments)
         } else {
             None

crates/sage-cli/src/output.rs

@@ -58,6 +58,13 @@ impl Runner {
         record.push_field(ryu::Buffer::new().format(feature.aligned_rt).as_bytes());
         record.push_field(ryu::Buffer::new().format(feature.predicted_rt).as_bytes());
         record.push_field(ryu::Buffer::new().format(feature.delta_rt_model).as_bytes());
+        record.push_field(ryu::Buffer::new().format(feature.ims).as_bytes());
+        record.push_field(ryu::Buffer::new().format(feature.predicted_ims).as_bytes());
+        record.push_field(
+            ryu::Buffer::new()
+                .format(feature.delta_ims_model)
+                .as_bytes(),
+        );
         record.push_field(itoa::Buffer::new().format(feature.matched_peaks).as_bytes());
         record.push_field(itoa::Buffer::new().format(feature.longest_b).as_bytes());
         record.push_field(itoa::Buffer::new().format(feature.longest_y).as_bytes());
@@ -174,6 +181,9 @@ impl Runner {
             "aligned_rt",
             "predicted_rt",
             "delta_rt_model",
+            "ion_mobility",
+            "predicted_mobility",
+            "delta_mobility",
             "matched_peaks",
             "longest_b",
             "longest_y",
@@ -261,6 +271,7 @@ impl Runner {
         record.push_field(ryu::Buffer::new().format(feature.calcmass).as_bytes());
         record.push_field(filenames[feature.file_id].as_bytes());
         record.push_field(ryu::Buffer::new().format(feature.rt).as_bytes());
+        record.push_field(ryu::Buffer::new().format(feature.ims).as_bytes());
         record.push_field(itoa::Buffer::new().format(feature.rank).as_bytes());
         record.push_field(
             itoa::Buffer::new()
@@ -336,6 +347,12 @@ impl Runner {
                 .format(feature.delta_rt_model.clamp(0.001, 1.0).sqrt())
                 .as_bytes(),
         );
+        record.push_field(ryu::Buffer::new().format(feature.predicted_ims).as_bytes());
+        record.push_field(
+            ryu::Buffer::new()
+                .format(feature.delta_ims_model)
+                .as_bytes(),
+        );
         record.push_field(itoa::Buffer::new().format(feature.matched_peaks).as_bytes());
         record.push_field(itoa::Buffer::new().format(feature.longest_b).as_bytes());
         record.push_field(itoa::Buffer::new().format(feature.longest_y).as_bytes());
@@ -384,6 +401,7 @@ impl Runner {
             "CalcMass",
             "FileName",
             "retentiontime",
+            "ion_mobility",
             "rank",
             "z=2",
             "z=3",
@@ -403,6 +421,8 @@ impl Runner {
             "aligned_rt",
             "predicted_rt",
             "sqrt(delta_rt_model)",
+            "predicted_mobility",
+            "sqrt(delta_mobility)",
             "matched_peaks",
             "longest_b",
             "longest_y",

crates/sage-cloudpath/src/mzml.rs

@@ -59,6 +59,8 @@ const SELECTED_ION_CHARGE: &[u8] = b"MS:1000041";
 const ISO_WINDOW_LOWER: &[u8] = b"MS:1000828";
 const ISO_WINDOW_UPPER: &[u8] = b"MS:1000829";
 
+const INVERSE_ION_MOBILITY: &[u8] = b"MS:1002815";
+
 pub struct MzMLReader {
     ms_level: Option<u8>,
     // If set to Some(level) and noise intensities are present in the MzML file,
@@ -236,6 +238,9 @@ impl MzMLReader {
                             SELECTED_ION_INT => {
                                 precursor.intensity = Some(extract_value!(ev));
                             }
+                            INVERSE_ION_MOBILITY => {
+                                precursor.inverse_ion_mobility = Some(extract_value!(ev));
+                            }
                             _ => {}
                         }
                     }
@@ -448,6 +453,7 @@ mod test {
                         <selectedIon>
                             <cvParam cvRef="MS" accession="MS:1000744" name="selected ion m/z" value="457.723968505859" unitAccession="MS:1000040" unitName="m/z" unitCvRef="MS" />
                             <cvParam cvRef="MS" accession="MS:1000041" name="charge state" value="2" />
+                            <cvParam cvRef="MS" accession="MS:1002815" name="inverse reduced ion mobility" value="1.078628" unitAccession="MS:1002814" unitName="volt-second per square centimeter"/>
                         </selectedIon>
                     </selectedIonList>
                     <activation>
@@ -483,6 +489,10 @@ mod test {
         assert_eq!(s.precursors.len(), 1);
         assert_eq!(s.precursors[0].charge, Some(2));
         assert!((s.precursors[0].mz - 457.723968) < 0.0001);
+        assert!(match s.precursors[0].inverse_ion_mobility {
+            Some(x) => (x - 1.0786) < 0.0001,
+            None => false,
+        });
         assert_eq!(
             s.precursors[0].isolation_window,
             Some(Tolerance::Da(-1.5, 0.75))

crates/sage-cloudpath/src/tdf.rs

@@ -23,6 +23,7 @@ impl TdfReader {
                 // precursor.ion_mobility = Option::from(dda_precursor.im as f32);
                 precursor.intensity = Option::from(dda_precursor.intensity as f32);
                 precursor.spectrum_ref = Option::from(dda_precursor.frame_index.to_string());
+                precursor.inverse_ion_mobility = Option::from(dda_precursor.im as f32);
                 let spectrum: RawSpectrum = RawSpectrum {
                     file_id,
                     precursors: vec![precursor],

crates/sage/src/ml/linear_discriminant.rs

@@ -16,7 +16,7 @@ use crate::mass::Tolerance;
 use crate::scoring::Feature;
 
 // Declare, so that we have compile time checking of matrix dimensions
-const FEATURES: usize = 18;
+const FEATURES: usize = 20;
 const FEATURE_NAMES: [&str; FEATURES] = [
     "rank",
     "charge",
@@ -35,7 +35,9 @@ const FEATURE_NAMES: [&str; FEATURES] = [
     "ln1p(peptide_len)",
     "missed_cleavages",
     "rt",
+    "ims",
     "sqrt(delta_rt_model)",
+    "sqrt(delta_ims_model)",
 ];
 
 struct Features<'a>(&'a [f64]);
@@ -174,7 +176,9 @@ pub fn score_psms(scores: &mut [Feature], precursor_tol: Tolerance) -> Option<()
                 (perc.peptide_len as f64).ln_1p(),
                 (perc.missed_cleavages as f64),
                 (perc.aligned_rt as f64),
+                (perc.ims as f64),
                 (perc.delta_rt_model as f64).clamp(0.001, 0.999).sqrt(),
+                (perc.delta_ims_model as f64).clamp(0.001, 0.999).sqrt(),
             ];
             x
         })