Speed up open/chimeric searching 🚀

- Rework isotope errors, precursor filtering interally

- Use u8 for internal residue representation, instead of char

Author: lazear

.github/workflows/build.yml

@@ -79,7 +79,7 @@ jobs:
           sage_version=${{ github.ref_name }} 
           staging="sage-$sage_version-${{ matrix.target }}"
           mkdir -p "$staging"
-          cp {README.md,LICENSE} "$staging/"
+          cp {README.md,DOCS.md,LICENSE} "$staging/"
           cp "target/${{ matrix.target }}/release/sage.exe" "$staging/"
           7z a "$staging.zip" "$staging"
           echo "ASSET=$staging.zip" >> $GITHUB_ENV

CHANGELOG.md

@@ -4,6 +4,12 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [0.12.0]
+### Added
+- Add `wide_window` option to configuration file. This option turns off `precursor_tol`, instead using the isolation window written in the mzML file.
+### Changed
+- Changed internal calculation of precursor tolerances when searching with `isotope_errors`. The new version should be more accurate. This change also enables a significant boost to search speed for open searches.
+
 ## [0.11.2]
 ### Added
 - Add rank & charge features to LDA

Cargo.lock

@@ -6,6 +6,6 @@ members = [
 ]
 
 [profile.release]
-#lto = "fat"
-#codegen-units = 1
+lto = "fat"
+codegen-units = 1
 panic = "abort"

Cargo.toml

@@ -156,6 +156,8 @@ Example:
     "isotope_errors": [-1, 3]
     ```
 
+**NOTE**: Searching with isotope errors is slower than searching with a wider precursor tolerance that encompasses the isotope errors, e.g. `"da": [-3.5, 1.25]`. Using the wider precursor tolerance will generally increase the number of confidently identified PSMs as well.
+
 ## Other Settings
 
 Note on the settings below:

DOCS.md

@@ -245,7 +245,7 @@ For additional information about configuration options and output file formats,
   ],
   "deisotope": false,       // Optional[bool] {default=false}: perform deisotoping and charge state deconvolution
   "chimera": false,         // Optional[bool] {default=false}: search for chimeric/co-fragmenting PSMS
-  "wide_window": false,     // Optional[bool] {default=false}: ignore `precursor_tol` and search in wide-window/DIA mode
+  "wide_window": false,     // Optional[bool] {default=false}: _ignore_ `precursor_tol` and search in wide-window/DIA mode
   "predict_rt": false,    // Optional[bool] {default=true}: use retention time prediction model as an feature for LDA
   "min_peaks": 15,          // Optional[int] {default=15}: only process MS2 spectra with at least N peaks
   "max_peaks": 150,         // Optional[int] {default=150}: take the top N most intense MS2 peaks to search,

README.md

@@ -1,6 +1,6 @@
 [package]
 name = "sage-cli"
-version = "0.11.2"
+version = "0.12.0"
 authors = ["Michael Lazear <[email protected]"]
 edition = "2021"
 rust-version = "1.62"

crates/sage-cli/Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "sage-core"
-version = "0.11.2"
+version = "0.12.0"
 authors = ["Michael Lazear <[email protected]"]
 edition = "2021"
 rust-version = "1.62"

crates/sage/Cargo.toml

@@ -1,7 +1,7 @@
 use crate::enzyme::{Enzyme, EnzymeParameters};
 use crate::fasta::Fasta;
 use crate::ion_series::{IonSeries, Kind};
-use crate::mass::{Tolerance, NEUTRON};
+use crate::mass::Tolerance;
 use crate::peptide::Peptide;
 use dashmap::DashMap;
 use fnv::FnvBuildHasher;
@@ -92,7 +92,7 @@ impl Builder {
         let mut output = HashMap::new();
         if let Some(input) = input {
             for (ch, mass) in input {
-                if crate::mass::VALID_AA.contains(&ch) || "^$[]".contains(ch) {
+                if crate::mass::VALID_AA.contains(&(ch as u8)) || "^$[]".contains(ch) {
                     output.insert(ch, mass);
                 } else {
                     error!(
@@ -359,25 +359,21 @@ impl IndexedDatabase {
         precursor_mass: f32,
         precursor_tol: Tolerance,
         fragment_tol: Tolerance,
-        min_isotope_err: i8,
-        max_isotope_err: i8,
     ) -> IndexedQuery<'_> {
         let (precursor_lo, precursor_hi) = precursor_tol.bounds(precursor_mass);
 
         let (pre_idx_lo, pre_idx_hi) = binary_search_slice(
             &self.peptides,
             |p, bounds| p.monoisotopic.total_cmp(bounds),
-            precursor_lo - (NEUTRON * max_isotope_err as f32).abs(),
-            precursor_hi + (NEUTRON * min_isotope_err as f32).abs(),
+            precursor_lo,
+            precursor_hi,
         );
 
         IndexedQuery {
             db: self,
             precursor_mass,
             precursor_tol,
             fragment_tol,
-            min_isotope_err,
-            max_isotope_err,
             pre_idx_lo,
             pre_idx_hi,
         }
@@ -421,8 +417,6 @@ pub struct IndexedQuery<'d> {
     precursor_mass: f32,
     precursor_tol: Tolerance,
     fragment_tol: Tolerance,
-    min_isotope_err: i8,
-    max_isotope_err: i8,
     pub pre_idx_lo: usize,
     pub pre_idx_hi: usize,
 }
@@ -464,11 +458,24 @@ impl<'d> IndexedQuery<'d> {
 
             // Finally, filter down our slice into exact matches only
             slice[inner_left..inner_right].iter().filter(move |frag| {
-                let neutral = self.db[frag.peptide_index].monoisotopic;
-                (self.min_isotope_err..=self.max_isotope_err).any(|isotope_err| {
-                    let delta = isotope_err as f32 * NEUTRON;
-                    (neutral >= precursor_lo - delta) && (neutral <= precursor_hi - delta)
-                }) && frag.fragment_mz >= fragment_lo
+                // This looks somewhat complicated, but it's a consequence of
+                // how the `binary_search_slice` function works - it will return
+                // the set of indices that maximally cover the desired range - the exact
+                // `left` and `right` indices may be valid, or just outside of the range.
+                // Anything interior of `left` and `right` is guaranteed to be within the
+                // precursor tolerance, so we just need to check the edge cases
+                //
+                // Previously, a direct lookup to check the mass of the current fragment was
+                // performed, but the pointer indirection + float comparison can slow down
+                // open searches by as much as 2x!!
+                // e.g. used to be `self.db[frag.peptide_index].monoisotopic >= precursor_lo`
+                (frag.peptide_index.0 > self.pre_idx_lo as u32
+                    || (frag.peptide_index.0 == self.pre_idx_lo as u32
+                        && self.db[frag.peptide_index].monoisotopic >= precursor_lo))
+                    && (frag.peptide_index.0 < self.pre_idx_hi as u32
+                        || (frag.peptide_index.0 == self.pre_idx_hi as u32
+                            && self.db[frag.peptide_index].monoisotopic <= precursor_hi))
+                    && frag.fragment_mz >= fragment_lo
                     && frag.fragment_mz <= fragment_hi
             })
         })

crates/sage/src/database.rs

@@ -94,13 +94,15 @@ pub struct Enzyme {
 impl Enzyme {
     pub fn new(cleave: &str, skip_suffix: Option<char>, c_terminal: bool) -> Option<Self> {
         assert!(
-            cleave.chars().all(|x| VALID_AA.contains(&x)) || cleave == "$",
+            cleave.chars().all(|x| VALID_AA.contains(&(x as u8))) || cleave == "$",
             "Enzyme cleavage sequence contains non-amino acid characters: {}",
             cleave
         );
 
         assert!(
-            skip_suffix.map(|x| VALID_AA.contains(&x)).unwrap_or(true),
+            skip_suffix
+                .map(|x| VALID_AA.contains(&(x as u8)))
+                .unwrap_or(true),
             "Enzyme cleavage restriction is non-amino acid character: {}",
             skip_suffix.unwrap(),
         );

crates/sage/src/enzyme.rs

@@ -1,5 +1,5 @@
 use crate::database::{binary_search_slice, IndexedDatabase, PeptideIx};
-use crate::mass::{Tolerance, NEUTRON};
+use crate::mass::{Composition, Mass, Tolerance, NEUTRON};
 use crate::ml::{matrix::Matrix, retention_alignment::Alignment};
 use crate::scoring::Feature;
 use crate::spectrum::ProcessedSpectrum;
@@ -218,8 +218,15 @@ impl FeatureMap {
                                 .entry((entry.peptide, entry.decoy))
                                 .or_insert_with(|| {
                                     let p = &db[entry.peptide];
-                                    let dist =
-                                        crate::isotopes::peptide_isotopes(p.carbons, p.sulfurs);
+                                    let composition = p
+                                        .sequence
+                                        .iter()
+                                        .map(|r| r.composition())
+                                        .sum::<Composition>();
+                                    let dist = crate::isotopes::peptide_isotopes(
+                                        composition.carbon,
+                                        composition.sulfur,
+                                    );
                                     Grid::new(entry, RT_TOL, dist, alignments.len(), GRID_SIZE)
                                 });
 

crates/sage/src/lfq.rs

@@ -57,9 +57,9 @@ pub trait Mass {
 #[derive(Clone, Debug, PartialEq, PartialOrd, Serialize)]
 pub enum Residue {
     // Standard amino acid residue
-    Just(char),
+    Just(u8),
     // Amino acid residue with a mass modification
-    Mod(char, f32),
+    Mod(u8, f32),
 }
 
 impl Mass for Residue {
@@ -78,65 +78,65 @@ impl Mass for Residue {
     }
 }
 
-pub const VALID_AA: [char; 22] = [
-    'A', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'K', 'L', 'M', 'N', 'P', 'Q', 'R', 'S', 'T', 'V', 'W',
-    'Y', 'U', 'O',
+pub const VALID_AA: [u8; 22] = [
+    b'A', b'C', b'D', b'E', b'F', b'G', b'H', b'I', b'K', b'L', b'M', b'N', b'P', b'Q', b'R', b'S',
+    b'T', b'V', b'W', b'Y', b'U', b'O',
 ];
 
-impl Mass for char {
+impl Mass for u8 {
     fn monoisotopic(&self) -> f32 {
         match self {
-            'A' => 71.03711,
-            'R' => 156.1011,
-            'N' => 114.04293,
-            'D' => 115.02694,
-            'C' => 103.00919,
-            'E' => 129.04259,
-            'Q' => 128.05858,
-            'G' => 57.02146,
-            'H' => 137.05891,
-            'I' => 113.08406,
-            'L' => 113.08406,
-            'K' => 128.09496,
-            'M' => 131.0405,
-            'F' => 147.0684,
-            'P' => 97.05276,
-            'S' => 87.03203,
-            'T' => 101.04768,
-            'W' => 186.07931,
-            'Y' => 163.06333,
-            'V' => 99.06841,
-            'U' => 150.95363,
-            'O' => 237.14773,
-            _ => unreachable!("BUG: invalid amino acid {}", self),
+            b'A' => 71.03711,
+            b'R' => 156.1011,
+            b'N' => 114.04293,
+            b'D' => 115.02694,
+            b'C' => 103.00919,
+            b'E' => 129.04259,
+            b'Q' => 128.05858,
+            b'G' => 57.02146,
+            b'H' => 137.05891,
+            b'I' => 113.08406,
+            b'L' => 113.08406,
+            b'K' => 128.09496,
+            b'M' => 131.0405,
+            b'F' => 147.0684,
+            b'P' => 97.05276,
+            b'S' => 87.03203,
+            b'T' => 101.04768,
+            b'W' => 186.07931,
+            b'Y' => 163.06333,
+            b'V' => 99.06841,
+            b'U' => 150.95363,
+            b'O' => 237.14773,
+            _ => unreachable!("BUG: invalid amino acid {}", *self as char),
         }
     }
 
     fn composition(&self) -> Composition {
         match self {
-            'A' => Composition::new(3, 2, 0),
-            'R' => Composition::new(6, 2, 0),
-            'N' => Composition::new(4, 3, 0),
-            'D' => Composition::new(4, 4, 0),
-            'C' => Composition::new(3, 2, 1),
-            'E' => Composition::new(5, 4, 0),
-            'Q' => Composition::new(5, 3, 0),
-            'G' => Composition::new(2, 2, 0),
-            'H' => Composition::new(6, 2, 0),
-            'I' => Composition::new(6, 2, 0),
-            'L' => Composition::new(6, 2, 0),
-            'K' => Composition::new(6, 2, 0),
-            'M' => Composition::new(5, 2, 1),
-            'F' => Composition::new(9, 2, 0),
-            'P' => Composition::new(5, 2, 0),
-            'S' => Composition::new(3, 3, 0),
-            'T' => Composition::new(4, 3, 0),
-            'W' => Composition::new(11, 2, 0),
-            'Y' => Composition::new(9, 3, 0),
-            'V' => Composition::new(5, 2, 0),
-            'U' => Composition::new(3, 2, 0),
-            'O' => Composition::new(12, 3, 0),
-            _ => unreachable!("BUG: invalid amino acid {}", self),
+            b'A' => Composition::new(3, 2, 0),
+            b'R' => Composition::new(6, 2, 0),
+            b'N' => Composition::new(4, 3, 0),
+            b'D' => Composition::new(4, 4, 0),
+            b'C' => Composition::new(3, 2, 1),
+            b'E' => Composition::new(5, 4, 0),
+            b'Q' => Composition::new(5, 3, 0),
+            b'G' => Composition::new(2, 2, 0),
+            b'H' => Composition::new(6, 2, 0),
+            b'I' => Composition::new(6, 2, 0),
+            b'L' => Composition::new(6, 2, 0),
+            b'K' => Composition::new(6, 2, 0),
+            b'M' => Composition::new(5, 2, 1),
+            b'F' => Composition::new(9, 2, 0),
+            b'P' => Composition::new(5, 2, 0),
+            b'S' => Composition::new(3, 3, 0),
+            b'T' => Composition::new(4, 3, 0),
+            b'W' => Composition::new(11, 2, 0),
+            b'Y' => Composition::new(9, 3, 0),
+            b'V' => Composition::new(5, 2, 0),
+            b'U' => Composition::new(3, 2, 0),
+            b'O' => Composition::new(12, 3, 0),
+            _ => unreachable!("BUG: invalid amino acid {}", *self as char),
         }
     }
 }
@@ -167,12 +167,12 @@ impl Sum for Composition {
 impl std::fmt::Display for Residue {
     fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
         match self {
-            Residue::Just(c) => f.write_char(*c),
+            Residue::Just(c) => f.write_char(*c as char),
             Residue::Mod(c, m) => {
                 if m.is_sign_positive() {
-                    write!(f, "{}[+{}]", c, m)
+                    write!(f, "{}[+{}]", *c as char, m)
                 } else {
-                    write!(f, "{}[{}]", c, m)
+                    write!(f, "{}[{}]", *c as char, m)
                 }
             }
         }

crates/sage/src/mass.rs

@@ -140,9 +140,9 @@ pub fn score_psms(scores: &mut [Feature]) -> Option<()> {
                 (perc.delta_next).ln_1p(),
                 (perc.delta_best).ln_1p(),
                 (perc.delta_mass as f64).ln_1p(),
-                perc.isotope_error as f64,
-                perc.average_ppm as f64,
-                poisson,
+                (perc.isotope_error as f64),
+                (perc.average_ppm as f64),
+                (poisson),
                 (perc.matched_intensity_pct as f64).ln_1p(),
                 (perc.matched_peaks as f64),
                 (perc.longest_b as f64).ln_1p(),
@@ -162,7 +162,6 @@ pub fn score_psms(scores: &mut [Feature]) -> Option<()> {
         .map(|sc| sc.label == -1)
         .collect::<Vec<_>>();
     let features = Matrix::new(features, scores.len(), FEATURES);
-
     let lda = LinearDiscriminantAnalysis::train(&features, &decoys)?;
     if !lda.eigenvector.iter().all(|f| f.is_finite()) {
         log::error!(

crates/sage/src/ml/linear_discriminant.rs

@@ -47,7 +47,7 @@ impl RetentionModel {
                 crate::mass::Residue::Just(c) => c,
                 crate::mass::Residue::Mod(c, _) => c,
             };
-            let idx = map[((*c) as u8 - b'A') as usize];
+            let idx = map[(c - b'A') as usize];
             embedding[idx] += 1.0;
             // Embed N- and C-terminal AA's (2 on each end, excluding K/R)
             match aa_idx {
@@ -67,7 +67,7 @@ impl RetentionModel {
         // Create a mapping from amino acid character to vector embedding
         let mut map = [0; 26];
         for (idx, aa) in VALID_AA.iter().enumerate() {
-            map[((*aa as u8) - b'A') as usize] = idx;
+            map[(aa - b'A') as usize] = idx;
         }
 
         let rt = training_set