Skip to main content

mscore/algorithm/
peptide.rs

1use crate::chemistry::amino_acid::{amino_acid_composition, amino_acid_masses};
2use crate::chemistry::constants::{MASS_CO, MASS_NH3, MASS_PROTON, MASS_WATER};
3use crate::chemistry::formulas::calculate_mz;
4use crate::chemistry::unimod::{
5    modification_atomic_composition, unimod_modifications_mass_numerical,
6};
7use crate::chemistry::utility::{find_unimod_patterns, unimod_sequence_to_tokens};
8use crate::data::peptide::{FragmentType, PeptideProductIon, PeptideSequence};
9use rayon::prelude::*;
10use rayon::ThreadPoolBuilder;
11use regex::Regex;
12use statrs::distribution::{Binomial, Discrete};
13use std::collections::HashMap;
14
15/// calculate the monoisotopic mass of a peptide sequence
16///
17/// Arguments:
18///
19/// * `sequence` - peptide sequence
20///
21/// Returns:
22///
23/// * `mass` - monoisotopic mass of the peptide
24///
25/// # Examples
26///
27/// ```
28/// use mscore::algorithm::peptide::calculate_peptide_mono_isotopic_mass;
29/// use mscore::data::peptide::PeptideSequence;
30///
31/// let peptide_sequence = PeptideSequence::new("PEPTIDEH".to_string(), Some(1));
32/// let mass = calculate_peptide_mono_isotopic_mass(&peptide_sequence);
33/// let mass_quantized = (mass * 1e6).round() as i32;
34/// assert_eq!(mass_quantized, 936418877);
35/// ```
36pub fn calculate_peptide_mono_isotopic_mass(peptide_sequence: &PeptideSequence) -> f64 {
37    let amino_acid_masses = amino_acid_masses();
38    let modifications_mz_numerical = unimod_modifications_mass_numerical();
39    let pattern = Regex::new(r"\[UNIMOD:(\d+)]").unwrap();
40
41    let sequence = peptide_sequence.sequence.as_str();
42
43    // Find all occurrences of the pattern
44    let modifications: Vec<u32> = pattern
45        .find_iter(sequence)
46        .filter_map(|mat| mat.as_str()[8..mat.as_str().len() - 1].parse().ok())
47        .collect();
48
49    // Remove the modifications from the sequence
50    let sequence = pattern.replace_all(sequence, "");
51
52    // Count occurrences of each amino acid
53    let mut aa_counts = HashMap::new();
54    for char in sequence.chars() {
55        *aa_counts.entry(char).or_insert(0) += 1;
56    }
57
58    // Mass of amino acids and modifications. HashMap iteration order is
59    // randomized, while floating-point addition is not associative. Sort the
60    // residue counts to make monoisotopic mass calculation bitwise reproducible.
61    let mut aa_counts_sorted: Vec<_> = aa_counts.into_iter().collect();
62    aa_counts_sorted.sort_unstable_by_key(|(aa, _)| *aa);
63
64    let mass_sequence: f64 = aa_counts_sorted
65        .into_iter()
66        .map(|(aa, count)| {
67            amino_acid_masses.get(&aa.to_string()[..]).unwrap_or(&0.0) * count as f64
68        })
69        .sum();
70    let mass_modifications: f64 = modifications
71        .iter()
72        .map(|&mod_id| modifications_mz_numerical.get(&mod_id).unwrap_or(&0.0))
73        .sum();
74
75    mass_sequence + mass_modifications + MASS_WATER
76}
77
78/// calculate the monoisotopic mass of a peptide product ion for a given fragment type
79///
80/// Arguments:
81///
82/// * `sequence` - peptide sequence
83/// * `kind` - fragment type
84///
85/// Returns:
86///
87/// * `mass` - monoisotopic mass of the peptide
88///
89/// # Examples
90/// ```
91/// use mscore::algorithm::peptide::calculate_peptide_product_ion_mono_isotopic_mass;
92/// use mscore::data::peptide::FragmentType;
93/// let sequence = "PEPTIDEH";
94/// let mass = calculate_peptide_product_ion_mono_isotopic_mass(sequence, FragmentType::Y);
95/// assert_eq!(mass, 936.4188766862999);
96/// ```
97pub fn calculate_peptide_product_ion_mono_isotopic_mass(sequence: &str, kind: FragmentType) -> f64 {
98    let (sequence, modifications) = find_unimod_patterns(sequence);
99
100    // Return mz of empty sequence
101    if sequence.is_empty() {
102        return 0.0;
103    }
104
105    let amino_acid_masses = amino_acid_masses();
106
107    // Add up raw amino acid masses and potential modifications
108    let mass_sequence: f64 = sequence
109        .chars()
110        .map(|aa| amino_acid_masses.get(&aa.to_string()[..]).unwrap_or(&0.0))
111        .sum();
112
113    let mass_modifications: f64 = modifications.iter().sum();
114
115    // Calculate total mass
116    let mass = mass_sequence + mass_modifications + MASS_WATER;
117
118    let mass = match kind {
119        FragmentType::A => mass - MASS_CO - MASS_WATER,
120        FragmentType::B => mass - MASS_WATER,
121        FragmentType::C => mass + MASS_NH3 - MASS_WATER,
122        FragmentType::X => mass + MASS_CO - 2.0 * MASS_PROTON,
123        FragmentType::Y => mass,
124        FragmentType::Z => mass - MASS_NH3,
125    };
126
127    mass
128}
129
130/// calculate the monoisotopic m/z of a peptide product ion for a given fragment type and charge
131///
132/// Arguments:
133///
134/// * `sequence` - peptide sequence
135/// * `kind` - fragment type
136/// * `charge` - charge
137///
138/// Returns:
139///
140/// * `mz` - monoisotopic mass of the peptide
141///
142/// # Examples
143/// ```
144/// use mscore::algorithm::peptide::calculate_product_ion_mz;
145/// use mscore::chemistry::constants::MASS_PROTON;
146/// use mscore::data::peptide::FragmentType;
147/// let sequence = "PEPTIDEH";
148/// let mz = calculate_product_ion_mz(sequence, FragmentType::Y, Some(1));
149/// assert_eq!(mz, 936.4188766862999 + MASS_PROTON);
150/// ```
151pub fn calculate_product_ion_mz(sequence: &str, kind: FragmentType, charge: Option<i32>) -> f64 {
152    let mass = calculate_peptide_product_ion_mono_isotopic_mass(sequence, kind);
153    calculate_mz(mass, charge.unwrap_or(1))
154}
155
156/// get a count dictionary of the amino acid composition of a peptide sequence
157///
158/// Arguments:
159///
160/// * `sequence` - peptide sequence
161///
162/// Returns:
163///
164/// * `composition` - a dictionary of amino acid composition
165///
166/// # Examples
167///
168/// ```
169/// use mscore::algorithm::peptide::calculate_amino_acid_composition;
170///
171/// let sequence = "PEPTIDEH";
172/// let composition = calculate_amino_acid_composition(sequence);
173/// assert_eq!(composition.get("P"), Some(&2));
174/// assert_eq!(composition.get("E"), Some(&2));
175/// assert_eq!(composition.get("T"), Some(&1));
176/// assert_eq!(composition.get("I"), Some(&1));
177/// assert_eq!(composition.get("D"), Some(&1));
178/// assert_eq!(composition.get("H"), Some(&1));
179/// ```
180pub fn calculate_amino_acid_composition(sequence: &str) -> HashMap<String, i32> {
181    let mut composition = HashMap::new();
182    for char in sequence.chars() {
183        *composition.entry(char.to_string()).or_insert(0) += 1;
184    }
185    composition
186}
187
188/// calculate the atomic composition of a peptide sequence
189pub fn peptide_sequence_to_atomic_composition(
190    peptide_sequence: &PeptideSequence,
191) -> HashMap<&'static str, i32> {
192    let token_sequence = unimod_sequence_to_tokens(peptide_sequence.sequence.as_str(), false);
193    let mut collection: HashMap<&'static str, i32> = HashMap::new();
194
195    // Assuming amino_acid_composition and modification_composition return appropriate mappings...
196    let aa_compositions = amino_acid_composition();
197    let mod_compositions = modification_atomic_composition();
198
199    // No need for conversion to HashMap<String, ...> as long as you're directly accessing
200    // the HashMap provided by modification_composition() if it uses String keys.
201    for token in token_sequence {
202        if token.len() == 1 {
203            let char = token.chars().next().unwrap();
204            if let Some(composition) = aa_compositions.get(&char) {
205                for (key, value) in composition.iter() {
206                    *collection.entry(key).or_insert(0) += *value;
207                }
208            }
209        } else {
210            // Directly use &token without .as_str() conversion
211            if let Some(composition) = mod_compositions.get(&token) {
212                for (key, value) in composition.iter() {
213                    *collection.entry(key).or_insert(0) += *value;
214                }
215            }
216        }
217    }
218
219    // Add water
220    *collection.entry("H").or_insert(0) += 2; //
221    *collection.entry("O").or_insert(0) += 1; //
222
223    collection
224}
225
226/// calculate the atomic composition of a product ion
227///
228/// Arguments:
229///
230/// * `product_ion` - a PeptideProductIon instance
231///
232/// Returns:
233///
234/// * `Vec<(&str, i32)>` - a vector of tuples representing the atomic composition of the product ion
235pub fn atomic_product_ion_composition(product_ion: &PeptideProductIon) -> Vec<(&str, i32)> {
236    let mut composition = peptide_sequence_to_atomic_composition(&product_ion.ion.sequence);
237
238    match product_ion.kind {
239        FragmentType::A => {
240            // A: peptide_mass - CO - Water
241            *composition.entry("H").or_insert(0) -= 2;
242            *composition.entry("O").or_insert(0) -= 2;
243            *composition.entry("C").or_insert(0) -= 1;
244        }
245        FragmentType::B => {
246            // B: peptide_mass - Water
247            *composition.entry("H").or_insert(0) -= 2;
248            *composition.entry("O").or_insert(0) -= 1;
249        }
250        FragmentType::C => {
251            // C: peptide_mass + NH3 - Water
252            *composition.entry("H").or_insert(0) += 1;
253            *composition.entry("N").or_insert(0) += 1;
254            *composition.entry("O").or_insert(0) -= 1;
255        }
256        FragmentType::X => {
257            // X: peptide_mass + CO
258            *composition.entry("C").or_insert(0) += 1; // Add 1 for CO
259            *composition.entry("O").or_insert(0) += 1; // Add 1 for CO
260            *composition.entry("H").or_insert(0) -= 2; // Subtract 2 for 2 protons
261        }
262        FragmentType::Y => (),
263        FragmentType::Z => {
264            *composition.entry("H").or_insert(0) -= 3;
265            *composition.entry("N").or_insert(0) -= 1;
266        }
267    }
268
269    composition.iter().map(|(k, v)| (*k, *v)).collect()
270}
271
272/// Calculate the atomic composition of the complementary fragment.
273///
274/// When a peptide is fragmented, the resulting fragment ion has a complementary
275/// portion (the rest of the precursor). This function calculates the atomic
276/// composition of that complementary fragment, which is needed for quad-selection
277/// dependent isotope distribution calculations.
278///
279/// # Arguments
280///
281/// * `precursor_composition` - atomic composition of the full precursor
282/// * `fragment_composition` - atomic composition of the fragment ion (from atomic_product_ion_composition)
283///
284/// # Returns
285///
286/// * `HashMap<String, i32>` - atomic composition of the complementary fragment
287///
288/// # Examples
289///
290/// ```
291/// use mscore::algorithm::peptide::{peptide_sequence_to_atomic_composition, atomic_product_ion_composition, calculate_complementary_fragment_composition};
292/// use mscore::data::peptide::{PeptideSequence, PeptideProductIon, PeptideIon, FragmentType};
293///
294/// let precursor = PeptideSequence::new("PEPTIDE".to_string(), Some(1));
295/// let precursor_comp = peptide_sequence_to_atomic_composition(&precursor);
296///
297/// // Create a b3 ion (PEP)
298/// let b3_ion = PeptideProductIon {
299///     kind: FragmentType::B,
300///     ion: PeptideIon {
301///         sequence: PeptideSequence::new("PEP".to_string(), Some(1)),
302///         charge: 1,
303///         intensity: 1.0,
304///     },
305/// };
306/// let fragment_comp: Vec<(&str, i32)> = atomic_product_ion_composition(&b3_ion);
307/// let fragment_map: std::collections::HashMap<&str, i32> = fragment_comp.into_iter().collect();
308///
309/// let complementary = calculate_complementary_fragment_composition(&precursor_comp, &fragment_map);
310/// // Complementary should be TIDE as y-ion (precursor - b-ion)
311/// ```
312pub fn calculate_complementary_fragment_composition(
313    precursor_composition: &HashMap<&str, i32>,
314    fragment_composition: &HashMap<&str, i32>,
315) -> HashMap<String, i32> {
316    let mut complementary: HashMap<String, i32> = HashMap::new();
317
318    // Start with precursor - fragment
319    for (element, &prec_count) in precursor_composition.iter() {
320        let frag_count = fragment_composition.get(element).copied().unwrap_or(0);
321        let diff = prec_count - frag_count;
322        if diff != 0 {
323            complementary.insert(element.to_string(), diff);
324        }
325    }
326
327    // Check for any elements in fragment that might not be in precursor (edge case)
328    for (element, &frag_count) in fragment_composition.iter() {
329        if !precursor_composition.contains_key(element) {
330            complementary.insert(element.to_string(), -frag_count);
331        }
332    }
333
334    complementary
335}
336
337/// calculate the atomic composition of a peptide product ion series
338/// Arguments:
339///
340/// * `product_ions` - a vector of PeptideProductIon instances
341/// * `num_threads` - an usize representing the number of threads to use
342/// Returns:
343///
344/// * `Vec<Vec<(String, i32)>>` - a vector of vectors of tuples representing the atomic composition of each product ion
345///
346pub fn fragments_to_composition(
347    product_ions: Vec<PeptideProductIon>,
348    num_threads: usize,
349) -> Vec<Vec<(String, i32)>> {
350    let thread_pool = ThreadPoolBuilder::new()
351        .num_threads(num_threads)
352        .build()
353        .unwrap();
354    let result = thread_pool.install(|| {
355        product_ions
356            .par_iter()
357            .map(|ion| atomic_product_ion_composition(ion))
358            .map(|composition| {
359                composition
360                    .iter()
361                    .map(|(k, v)| (k.to_string(), *v))
362                    .collect()
363            })
364            .collect()
365    });
366    result
367}
368
369/// count the number of protonizable sites in a peptide sequence
370///
371/// # Arguments
372///
373/// * `sequence` - a string representing the peptide sequence
374///
375/// # Returns
376///
377/// * `usize` - the number of protonizable sites in the peptide sequence
378///
379/// # Example
380///
381/// ```
382/// use mscore::algorithm::peptide::get_num_protonizable_sites;
383///
384/// let sequence = "PEPTIDEH";
385/// let num_protonizable_sites = get_num_protonizable_sites(sequence);
386/// assert_eq!(num_protonizable_sites, 2);
387/// ```
388pub fn get_num_protonizable_sites(sequence: &str) -> usize {
389    let mut sites = 1; // n-terminus
390    for s in sequence.chars() {
391        match s {
392            'H' | 'R' | 'K' => sites += 1,
393            _ => {}
394        }
395    }
396    sites
397}
398
399/// simulate the charge state distribution for a peptide sequence
400///
401/// # Arguments
402///
403/// * `sequence` - a string representing the peptide sequence
404/// * `max_charge` - an optional usize representing the maximum charge state to simulate
405/// * `charged_probability` - an optional f64 representing the probability of a site being charged
406///
407/// # Returns
408///
409/// * `Vec<f64>` - a vector of f64 representing the probability of each charge state
410///
411/// # Example
412///
413/// ```
414/// use mscore::algorithm::peptide::simulate_charge_state_for_sequence;
415///
416/// let sequence = "PEPTIDEH";
417/// let charge_state_probs = simulate_charge_state_for_sequence(sequence, None, None);
418/// assert_eq!(charge_state_probs, vec![0.03999999999999999, 0.32, 0.64, 0.0, 0.0]);
419pub fn simulate_charge_state_for_sequence(
420    sequence: &str,
421    max_charge: Option<usize>,
422    charged_probability: Option<f64>,
423) -> Vec<f64> {
424    let charged_prob = charged_probability.unwrap_or(0.8);
425    let max_charge = max_charge.unwrap_or(4)+1;
426    let num_protonizable_sites = get_num_protonizable_sites(sequence);
427    let mut charge_state_probs = vec![0.0; max_charge];
428    let binom = Binomial::new(charged_prob, num_protonizable_sites as u64).unwrap();
429    
430    for charge in 0..max_charge {
431        charge_state_probs[charge] = binom.pmf(charge as u64);
432    }
433    charge_state_probs
434}
435
436/// simulate the charge state distribution for a list of peptide sequences
437///
438/// # Arguments
439///
440/// * `sequences` - a vector of strings representing the peptide sequences
441/// * `num_threads` - an usize representing the number of threads to use
442/// * `max_charge` - an optional usize representing the maximum charge state to simulate
443/// * `charged_probability` - an optional f64 representing the probability of a site being charged
444///
445/// # Returns
446///
447/// * `Vec<Vec<f64>>` - a vector of vectors of f64 representing the probability of each charge state for each sequence
448///
449/// # Example
450///
451/// ```
452/// use mscore::algorithm::peptide::simulate_charge_states_for_sequences;
453///
454/// let sequences = vec!["PEPTIDEH", "PEPTIDEH", "PEPTIDEH"];
455/// let charge_state_probs = simulate_charge_states_for_sequences(sequences, 4, None, None);
456/// assert_eq!(charge_state_probs, vec![vec![0.03999999999999999, 0.32, 0.64, 0.0, 0.0], vec![0.03999999999999999, 0.32, 0.64, 0.0, 0.0], vec![0.03999999999999999, 0.32, 0.64, 0.0, 0.0]]);
457/// ```
458pub fn simulate_charge_states_for_sequences(
459    sequences: Vec<&str>,
460    num_threads: usize,
461    max_charge: Option<usize>,
462    charged_probability: Option<f64>,
463) -> Vec<Vec<f64>> {
464    let pool = ThreadPoolBuilder::new()
465        .num_threads(num_threads)
466        .build()
467        .unwrap();
468    pool.install(|| {
469        sequences
470            .par_iter()
471            .map(|sequence| {
472                simulate_charge_state_for_sequence(sequence, max_charge, charged_probability)
473            })
474            .collect()
475    })
476}
477
478#[cfg(test)]
479mod tests {
480    use super::*;
481    use std::collections::BTreeSet;
482
483    #[test]
484    fn peptide_monoisotopic_mass_is_bitwise_stable() {
485        let peptide = PeptideSequence::new("CNHHDGPSHADGK".to_string(), Some(1));
486
487        let mut observed_bits = BTreeSet::new();
488        for _ in 0..1024 {
489            observed_bits.insert(calculate_peptide_mono_isotopic_mass(&peptide).to_bits());
490        }
491
492        assert_eq!(
493            observed_bits.len(),
494            1,
495            "identical peptide sequence produced multiple exact monoisotopic masses: {observed_bits:?}"
496        );
497    }
498}