rustdf/sim/
dda.rs

1use mscore::algorithm::isotope::{
2    calculate_transmission_dependent_fragment_ion_isotope_distribution,
3    calculate_precursor_transmission_factor,
4};
5use crate::sim::containers::IsotopeTransmissionMode;
6use mscore::data::peptide::{PeptideIon, PeptideProductIonSeriesCollection};
7use mscore::data::spectrum::{IndexedMzSpectrum, MsType, MzSpectrum};
8use mscore::simulation::annotation::{
9    MzSpectrumAnnotated, PeakAnnotation, TimsFrameAnnotated, TimsSpectrumAnnotated,
10};
11use mscore::timstof::frame::TimsFrame;
12use mscore::timstof::quadrupole::{IonTransmission, TimsTransmissionDDA};
13use mscore::timstof::spectrum::TimsSpectrum;
14use std::collections::{BTreeMap, HashSet};
15use std::path::Path;
16use std::sync::Arc;
17
18use rand::Rng;
19use rayon::prelude::*;
20use crate::sim::containers::IsotopeTransmissionConfig;
21use crate::sim::handle::{FragmentIonsWithComplementary, TimsTofSyntheticsDataHandle};
22use crate::sim::precursor::TimsTofSyntheticsPrecursorFrameBuilder;
23
24pub struct TimsTofSyntheticsFrameBuilderDDA {
25    pub path: String,
26    pub precursor_frame_builder: TimsTofSyntheticsPrecursorFrameBuilder,
27    pub transmission_settings: TimsTransmissionDDA,
28    pub fragment_ions:
29        Option<BTreeMap<(u32, i8, i32), (PeptideProductIonSeriesCollection, Vec<MzSpectrum>)>>,
30    pub fragment_ions_annotated: Option<
31        BTreeMap<(u32, i8, i32), (PeptideProductIonSeriesCollection, Vec<MzSpectrumAnnotated>)>,
32    >,
33    /// Configuration for quad-selection dependent isotope transmission
34    pub isotope_transmission_config: IsotopeTransmissionConfig,
35    /// Fragment ions with complementary distribution data (only populated when isotope_transmission_config.enabled)
36    pub fragment_ions_with_complementary: Option<BTreeMap<(u32, i8, i32), FragmentIonsWithComplementary>>,
37}
38
39impl TimsTofSyntheticsFrameBuilderDDA {
40    /// Create a new DDA frame builder.
41    ///
42    /// # Arguments
43    ///
44    /// * `path` - Path to the simulation database
45    /// * `with_annotations` - Whether to include annotations
46    /// * `num_threads` - Number of threads for parallel processing
47    /// * `isotope_config` - Optional configuration for quad-dependent isotope transmission
48    pub fn new(
49        path: &Path,
50        with_annotations: bool,
51        num_threads: usize,
52        isotope_config: Option<IsotopeTransmissionConfig>,
53    ) -> Self {
54        let handle = TimsTofSyntheticsDataHandle::new(path).unwrap();
55        let fragment_ions_raw = handle.read_fragment_ions().unwrap();
56        let transmission_settings = handle.get_transmission_dda();
57
58        let synthetics = TimsTofSyntheticsPrecursorFrameBuilder::new(path).unwrap();
59        let config = isotope_config.unwrap_or_default();
60
61        // Build fragment ions with complementary data if any transmission mode is enabled
62        let fragment_ions_with_complementary = if config.is_enabled() {
63            Some(TimsTofSyntheticsDataHandle::build_fragment_ions_with_transmission_data(
64                &synthetics.peptides,
65                &fragment_ions_raw,
66                num_threads,
67            ))
68        } else {
69            None
70        };
71
72        match with_annotations {
73            true => {
74                let fragment_ions =
75                    Some(TimsTofSyntheticsDataHandle::build_fragment_ions_annotated(
76                        &synthetics.peptides,
77                        &fragment_ions_raw,
78                        num_threads,
79                    ));
80                Self {
81                    path: path.to_str().unwrap().to_string(),
82                    precursor_frame_builder: synthetics,
83                    transmission_settings,
84                    fragment_ions: None,
85                    fragment_ions_annotated: fragment_ions,
86                    isotope_transmission_config: config,
87                    fragment_ions_with_complementary,
88                }
89            }
90            false => {
91                let fragment_ions = Some(TimsTofSyntheticsDataHandle::build_fragment_ions(
92                    &synthetics.peptides,
93                    &fragment_ions_raw,
94                    num_threads,
95                ));
96                Self {
97                    path: path.to_str().unwrap().to_string(),
98                    precursor_frame_builder: synthetics,
99                    transmission_settings,
100                    fragment_ions,
101                    fragment_ions_annotated: None,
102                    isotope_transmission_config: config,
103                    fragment_ions_with_complementary,
104                }
105            }
106        }
107    }
108    /// Build a frame for DDA synthetic experiment
109    ///
110    /// # Arguments
111    ///
112    /// * `frame_id` - The frame id
113    /// * `fragmentation` - A boolean indicating if fragmentation is enabled, if false, the frame has same mz distribution as the precursor frame but will be quadrupole filtered
114    ///
115    /// # Returns
116    ///
117    /// A TimsFrame
118    ///
119    pub fn build_frame(
120        &self,
121        frame_id: u32,
122        fragmentation: bool,
123        mz_noise_precursor: bool,
124        uniform: bool,
125        precursor_noise_ppm: f64,
126        mz_noise_fragment: bool,
127        fragment_noise_ppm: f64,
128        right_drag: bool,
129    ) -> TimsFrame {
130        // determine if the frame is a precursor frame
131        match self
132            .precursor_frame_builder
133            .precursor_frame_id_set
134            .contains(&frame_id)
135        {
136            true => self.build_ms1_frame(
137                frame_id,
138                mz_noise_precursor,
139                uniform,
140                precursor_noise_ppm,
141                right_drag,
142            ),
143            false => self.build_ms2_frame(
144                frame_id,
145                fragmentation,
146                mz_noise_fragment,
147                uniform,
148                fragment_noise_ppm,
149                right_drag,
150            ),
151        }
152    }
153
154    pub fn build_frame_annotated(
155        &self,
156        frame_id: u32,
157        fragmentation: bool,
158        mz_noise_precursor: bool,
159        uniform: bool,
160        precursor_noise_ppm: f64,
161        mz_noise_fragment: bool,
162        fragment_noise_ppm: f64,
163        right_drag: bool,
164    ) -> TimsFrameAnnotated {
165        match self
166            .precursor_frame_builder
167            .precursor_frame_id_set
168            .contains(&frame_id)
169        {
170            true => self.build_ms1_frame_annotated(
171                frame_id,
172                mz_noise_precursor,
173                uniform,
174                precursor_noise_ppm,
175                right_drag,
176            ),
177            false => self.build_ms2_frame_annotated(
178                frame_id,
179                fragmentation,
180                mz_noise_fragment,
181                uniform,
182                fragment_noise_ppm,
183                right_drag,
184            ),
185        }
186    }
187
188    pub fn get_fragment_ion_ids(&self, precursor_frame_ids: Vec<u32>) -> Vec<u32> {
189        let mut peptide_ids: HashSet<u32> = HashSet::new();
190        // get all peptide ids for the precursor frame ids
191        for frame_id in precursor_frame_ids {
192            for (peptide_id, peptide) in self.precursor_frame_builder.peptides.iter() {
193                if peptide.frame_start <= frame_id && peptide.frame_end >= frame_id {
194                    peptide_ids.insert(*peptide_id);
195                }
196            }
197        }
198        // get all ion ids for the peptide ids
199        let mut result: Vec<u32> = Vec::new();
200        for item in peptide_ids {
201            let ions = self.precursor_frame_builder.ions.get(&item).unwrap();
202            for ion in ions.iter() {
203                result.push(ion.ion_id);
204            }
205        }
206        result
207    }
208
209    pub fn build_frames(
210        &self,
211        frame_ids: Vec<u32>,
212        fragmentation: bool,
213        mz_noise_precursor: bool,
214        uniform: bool,
215        precursor_noise_ppm: f64,
216        mz_noise_fragment: bool,
217        fragment_noise_ppm: f64,
218        right_drag: bool,
219        num_threads: usize,
220    ) -> Vec<TimsFrame> {
221        // Use thread pool with custom parallelism
222        let pool = rayon::ThreadPoolBuilder::new()
223            .num_threads(num_threads)
224            .build()
225            .unwrap();
226
227        pool.install(|| {
228            // Use indexed parallel iteration to maintain order, avoiding post-sort
229            let mut tims_frames: Vec<TimsFrame> = Vec::with_capacity(frame_ids.len());
230            unsafe { tims_frames.set_len(frame_ids.len()); }
231
232            frame_ids.par_iter().enumerate().for_each(|(idx, frame_id)| {
233                let frame = self.build_frame(
234                    *frame_id,
235                    fragmentation,
236                    mz_noise_precursor,
237                    uniform,
238                    precursor_noise_ppm,
239                    mz_noise_fragment,
240                    fragment_noise_ppm,
241                    right_drag,
242                );
243                unsafe {
244                    let ptr = tims_frames.as_ptr() as *mut TimsFrame;
245                    std::ptr::write(ptr.add(idx), frame);
246                }
247            });
248
249            tims_frames
250        })
251    }
252    pub fn build_frames_annotated(
253        &self,
254        frame_ids: Vec<u32>,
255        fragmentation: bool,
256        mz_noise_precursor: bool,
257        uniform: bool,
258        precursor_noise_ppm: f64,
259        mz_noise_fragment: bool,
260        fragment_noise_ppm: f64,
261        right_drag: bool,
262        num_threads: usize,
263    ) -> Vec<TimsFrameAnnotated> {
264        // Use thread pool with custom parallelism
265        let pool = rayon::ThreadPoolBuilder::new()
266            .num_threads(num_threads)
267            .build()
268            .unwrap();
269
270        pool.install(|| {
271            // Use indexed parallel iteration to maintain order, avoiding post-sort
272            let mut tims_frames: Vec<TimsFrameAnnotated> = Vec::with_capacity(frame_ids.len());
273            unsafe { tims_frames.set_len(frame_ids.len()); }
274
275            frame_ids.par_iter().enumerate().for_each(|(idx, frame_id)| {
276                let frame = self.build_frame_annotated(
277                    *frame_id,
278                    fragmentation,
279                    mz_noise_precursor,
280                    uniform,
281                    precursor_noise_ppm,
282                    mz_noise_fragment,
283                    fragment_noise_ppm,
284                    right_drag,
285                );
286                unsafe {
287                    let ptr = tims_frames.as_ptr() as *mut TimsFrameAnnotated;
288                    std::ptr::write(ptr.add(idx), frame);
289                }
290            });
291
292            tims_frames
293        })
294    }
295
296    fn build_ms1_frame(
297        &self,
298        frame_id: u32,
299        mz_noise_precursor: bool,
300        uniform: bool,
301        precursor_ppm: f64,
302        right_drag: bool,
303    ) -> TimsFrame {
304        let mut tims_frame = self.precursor_frame_builder.build_precursor_frame(
305            frame_id,
306            mz_noise_precursor,
307            uniform,
308            precursor_ppm,
309            right_drag,
310        );
311        let intensities_rounded = tims_frame
312            .ims_frame
313            .intensity
314            .iter()
315            .map(|x| x.round())
316            .collect::<Vec<_>>();
317        tims_frame.ims_frame.intensity = Arc::new(intensities_rounded);
318        tims_frame
319    }
320
321    fn build_ms1_frame_annotated(
322        &self,
323        frame_id: u32,
324        mz_noise_precursor: bool,
325        uniform: bool,
326        precursor_ppm: f64,
327        right_drag: bool,
328    ) -> TimsFrameAnnotated {
329        let mut tims_frame = self
330            .precursor_frame_builder
331            .build_precursor_frame_annotated(
332                frame_id,
333                mz_noise_precursor,
334                uniform,
335                precursor_ppm,
336                right_drag,
337            );
338        let intensities_rounded = tims_frame
339            .intensity
340            .iter()
341            .map(|x| x.round())
342            .collect::<Vec<_>>();
343        tims_frame.intensity = intensities_rounded;
344        tims_frame
345    }
346
347    fn build_ms2_frame(
348        &self,
349        frame_id: u32,
350        fragmentation: bool,
351        mz_noise_fragment: bool,
352        uniform: bool,
353        fragment_ppm: f64,
354        right_drag: bool,
355    ) -> TimsFrame {
356        match fragmentation {
357            false => {
358                let mut frame = self.transmission_settings.transmit_tims_frame(
359                    &self.build_ms1_frame(
360                        frame_id,
361                        mz_noise_fragment,
362                        uniform,
363                        fragment_ppm,
364                        right_drag,
365                    ),
366                    None,
367                );
368                let intensities_rounded = frame
369                    .ims_frame
370                    .intensity
371                    .iter()
372                    .map(|x| x.round())
373                    .collect::<Vec<_>>();
374                frame.ims_frame.intensity = Arc::new(intensities_rounded);
375                frame.ms_type = MsType::FragmentDia;
376                frame
377            }
378            true => {
379                let mut frame = self.build_fragment_frame(
380                    frame_id,
381                    &self.fragment_ions.as_ref().unwrap(),
382                    mz_noise_fragment,
383                    uniform,
384                    fragment_ppm,
385                    None,
386                    None,
387                    None,
388                    Some(right_drag),
389                );
390                let intensities_rounded = frame
391                    .ims_frame
392                    .intensity
393                    .iter()
394                    .map(|x| x.round())
395                    .collect::<Vec<_>>();
396                frame.ims_frame.intensity = Arc::new(intensities_rounded);
397                frame
398            }
399        }
400    }
401
402    fn build_ms2_frame_annotated(
403        &self,
404        frame_id: u32,
405        fragmentation: bool,
406        mz_noise_fragment: bool,
407        uniform: bool,
408        fragment_ppm: f64,
409        right_drag: bool,
410    ) -> TimsFrameAnnotated {
411        match fragmentation {
412            false => {
413                let mut frame = self.transmission_settings.transmit_tims_frame_annotated(
414                    &self.build_ms1_frame_annotated(
415                        frame_id,
416                        mz_noise_fragment,
417                        uniform,
418                        fragment_ppm,
419                        right_drag,
420                    ),
421                    None,
422                );
423                let intensities_rounded = frame
424                    .intensity
425                    .iter()
426                    .map(|x| x.round())
427                    .collect::<Vec<_>>();
428                frame.intensity = intensities_rounded;
429                frame.ms_type = MsType::FragmentDia;
430                frame
431            }
432            true => {
433                let mut frame = self.build_fragment_frame_annotated(
434                    frame_id,
435                    &self.fragment_ions_annotated.as_ref().unwrap(),
436                    mz_noise_fragment,
437                    uniform,
438                    fragment_ppm,
439                    None,
440                    None,
441                    None,
442                    Some(right_drag),
443                );
444                let intensities_rounded = frame
445                    .intensity
446                    .iter()
447                    .map(|x| x.round())
448                    .collect::<Vec<_>>();
449                frame.intensity = intensities_rounded;
450                frame
451            }
452        }
453    }
454
455    /// Build a fragment frame
456    ///
457    /// # Arguments
458    ///
459    /// * `frame_id` - The frame id
460    /// * `mz_min` - The minimum m/z value in fragment spectrum
461    /// * `mz_max` - The maximum m/z value in fragment spectrum
462    /// * `intensity_min` - The minimum intensity value in fragment spectrum
463    ///
464    /// # Returns
465    ///
466    /// A TimsFrame
467    ///
468    fn build_fragment_frame(
469        &self,
470        frame_id: u32,
471        fragment_ions: &BTreeMap<
472            (u32, i8, i32),
473            (PeptideProductIonSeriesCollection, Vec<MzSpectrum>),
474        >,
475        mz_noise_fragment: bool,
476        uniform: bool,
477        fragment_ppm: f64,
478        mz_min: Option<f64>,
479        mz_max: Option<f64>,
480        intensity_min: Option<f64>,
481        right_drag: Option<bool>,
482    ) -> TimsFrame {
483        // Cache frame-level lookups
484        let ms_type = if self.precursor_frame_builder.precursor_frame_id_set.contains(&frame_id) {
485            MsType::Unknown
486        } else {
487            MsType::FragmentDda
488        };
489
490        let rt = *self.precursor_frame_builder.frame_to_rt.get(&frame_id)
491            .expect("frame_to_rt should always have this frame") as f64;
492        let right_drag_val = right_drag.unwrap_or(false);
493        let mz_min_val = mz_min.unwrap_or(100.0);
494        let mz_max_val = mz_max.unwrap_or(1700.0);
495        let intensity_min_val = intensity_min.unwrap_or(1.0);
496
497        // Get PASEF meta for this frame - these define which precursors were SELECTED
498        let Some(pasef_meta) = self.transmission_settings.pasef_meta.get(&(frame_id as i32)) else {
499            return TimsFrame::new(frame_id as i32, ms_type, rt, vec![], vec![], vec![], vec![], vec![]);
500        };
501
502        // Preallocate with estimated capacity
503        let estimated_capacity = pasef_meta.len() * 4;
504        let mut tims_spectra: Vec<TimsSpectrum> = Vec::with_capacity(estimated_capacity);
505
506        // Iterate over PASEF meta entries - each represents a SELECTED precursor
507        for meta in pasef_meta.iter() {
508            // Get the selected ion_id from the PASEF precursor field
509            let ion_id = meta.precursor as u32;
510
511            // Look up peptide_id and charge from the ion_id
512            let Some(&(peptide_id, charge_state)) = self.precursor_frame_builder.ion_id_to_peptide_charge.get(&ion_id) else {
513                continue;
514            };
515
516            // Get peptide info
517            let Some(peptide) = self.precursor_frame_builder.peptides.get(&peptide_id) else {
518                continue;
519            };
520
521            // Check if this peptide is active in this frame
522            if frame_id < peptide.frame_start || frame_id > peptide.frame_end {
523                continue;
524            }
525
526            // Get ion info from peptide_to_ions
527            let Some((ion_abundances, scan_occurrences, scan_abundances, charges, spectra)) = self
528                .precursor_frame_builder
529                .peptide_to_ions
530                .get(&peptide_id)
531            else {
532                continue;
533            };
534
535            // Find the index of this specific charge state
536            let Some(ion_index) = charges.iter().position(|&c| c == charge_state) else {
537                continue;
538            };
539
540            let ion_abundance = ion_abundances[ion_index];
541            let all_scan_occurrence = &scan_occurrences[ion_index];
542            let all_scan_abundance = &scan_abundances[ion_index];
543            let spectrum = &spectra[ion_index];
544            let total_events = *self.precursor_frame_builder.peptide_to_events.get(&peptide_id).unwrap();
545
546            // Get frame abundance for this peptide
547            let frame_abundance = self.precursor_frame_builder.frame_to_abundances
548                .get(&frame_id)
549                .and_then(|(pep_ids, abundances)| {
550                    pep_ids.iter().position(|&p| p == peptide_id)
551                        .map(|idx| abundances[idx])
552                })
553                .unwrap_or(0.0);
554
555            if frame_abundance == 0.0 {
556                continue;
557            }
558
559            // Collision energy from the PASEF meta
560            let collision_energy = meta.collision_energy;
561            let collision_energy_quantized = (collision_energy * 1e1).round() as i32;
562
563            // Look up fragment ions for this (peptide_id, charge, CE)
564            let Some((_, fragment_series_vec)) = fragment_ions.get(&(peptide_id, charge_state, collision_energy_quantized)) else {
565                continue;
566            };
567
568            // Process scans within the PASEF selection window
569            for (scan, scan_abundance) in all_scan_occurrence.iter().zip(all_scan_abundance.iter()) {
570                // Check if scan is within the PASEF selection window
571                let scan_i32 = *scan as i32;
572                if scan_i32 < meta.scan_start || scan_i32 > meta.scan_end {
573                    continue;
574                }
575
576                // Get transmitted isotope indices based on config mode
577                let transmitted_indices = match self.isotope_transmission_config.mode {
578                    IsotopeTransmissionMode::None => HashSet::new(),
579                    IsotopeTransmissionMode::PrecursorScaling | IsotopeTransmissionMode::PerFragment => {
580                        self.transmission_settings.get_transmission_set(
581                            frame_id as i32,
582                            scan_i32,
583                            &spectrum.mz,
584                            Some(self.isotope_transmission_config.min_probability),
585                        )
586                    },
587                };
588
589                // Calculate abundance factor
590                let fraction_events = frame_abundance * scan_abundance * ion_abundance * total_events;
591
592                // Cache scan mobility
593                let scan_mobility = *self.precursor_frame_builder.scan_to_mobility.get(scan).unwrap() as f64;
594
595                // Calculate transmission factor for PrecursorScaling mode
596                let transmission_factor = if self.isotope_transmission_config.mode == IsotopeTransmissionMode::PrecursorScaling {
597                    if let Some(comp_data) = self.fragment_ions_with_complementary.as_ref() {
598                        if let Some(frag_data) = comp_data.get(&(peptide_id, charge_state, collision_energy_quantized)) {
599                            calculate_precursor_transmission_factor(
600                                &frag_data.precursor_isotope_distribution,
601                                &transmitted_indices,
602                            )
603                        } else {
604                            1.0
605                        }
606                    } else {
607                        1.0
608                    }
609                } else {
610                    1.0
611                };
612
613                for (series_idx, fragment_ion_series) in fragment_series_vec.iter().enumerate() {
614                    let final_spectrum = match self.isotope_transmission_config.mode {
615                        IsotopeTransmissionMode::None => {
616                            // Standard spectrum scaling
617                            fragment_ion_series.clone() * fraction_events as f64
618                        },
619                        IsotopeTransmissionMode::PrecursorScaling => {
620                            // Apply precursor-based transmission factor
621                            fragment_ion_series.clone() * (fraction_events as f64 * transmission_factor)
622                        },
623                        IsotopeTransmissionMode::PerFragment => {
624                            // Per-fragment transmission-dependent calculation
625                            if let Some(comp_data) = self.fragment_ions_with_complementary.as_ref() {
626                                if let Some(frag_data) = comp_data.get(&(peptide_id, charge_state, collision_energy_quantized)) {
627                                    if series_idx < frag_data.per_fragment_data.len() {
628                                        // Aggregate adjusted spectra from all fragment ions in this series
629                                        let series_data = &frag_data.per_fragment_data[series_idx];
630                                        let mut aggregated_mz: Vec<f64> = Vec::new();
631                                        let mut aggregated_intensity: Vec<f64> = Vec::new();
632
633                                        for frag_ion_data in series_data {
634                                            // Apply transmission-dependent calculation for this fragment
635                                            let adjusted_dist = calculate_transmission_dependent_fragment_ion_isotope_distribution(
636                                                &frag_ion_data.fragment_distribution,
637                                                &frag_ion_data.complementary_distribution,
638                                                &transmitted_indices,
639                                                self.isotope_transmission_config.max_isotopes,
640                                            );
641
642                                            // Scale by predicted intensity and fraction_events
643                                            for (mz, abundance) in adjusted_dist {
644                                                aggregated_mz.push(mz);
645                                                aggregated_intensity.push(abundance * frag_ion_data.predicted_intensity * fraction_events as f64);
646                                            }
647                                        }
648
649                                        if !aggregated_mz.is_empty() {
650                                            MzSpectrum::new(aggregated_mz, aggregated_intensity)
651                                        } else {
652                                            fragment_ion_series.clone() * fraction_events as f64
653                                        }
654                                    } else {
655                                        fragment_ion_series.clone() * fraction_events as f64
656                                    }
657                                } else {
658                                    fragment_ion_series.clone() * fraction_events as f64
659                                }
660                            } else {
661                                fragment_ion_series.clone() * fraction_events as f64
662                            }
663                        },
664                    };
665
666                    let mz_spectrum = if mz_noise_fragment {
667                        if uniform {
668                            final_spectrum.add_mz_noise_uniform(fragment_ppm, right_drag_val)
669                        } else {
670                            final_spectrum.add_mz_noise_normal(fragment_ppm)
671                        }
672                    } else {
673                        final_spectrum
674                    };
675
676                    let spectrum_len = mz_spectrum.mz.len();
677                    tims_spectra.push(TimsSpectrum::new(
678                        frame_id as i32,
679                        *scan as i32,
680                        rt,
681                        scan_mobility,
682                        ms_type.clone(),
683                        IndexedMzSpectrum::from_mz_spectrum(
684                            vec![0; spectrum_len],
685                            mz_spectrum,
686                        ).filter_ranged(100.0, 1700.0, 1.0, 1e9),
687                    ));
688                }
689
690                // Add unfragmented precursor ions (survival) if configured
691                if self.isotope_transmission_config.has_precursor_survival() {
692                    let mut rng = rand::thread_rng();
693                    let survival_fraction = rng.gen_range(
694                        self.isotope_transmission_config.precursor_survival_min
695                        ..=self.isotope_transmission_config.precursor_survival_max
696                    );
697
698                    if survival_fraction > 0.0 {
699                        // Transmit the precursor spectrum through the quadrupole
700                        let precursor_transmitted = self.transmission_settings.transmit_spectrum(
701                            frame_id as i32,
702                            *scan as i32,
703                            spectrum.clone(),
704                            Some(self.isotope_transmission_config.min_probability),
705                        );
706
707                        if !precursor_transmitted.mz.is_empty() {
708                            // Scale by survival fraction and event count
709                            let precursor_scaled = precursor_transmitted * (fraction_events as f64 * survival_fraction);
710
711                            let precursor_mz_spectrum = if mz_noise_fragment {
712                                if uniform {
713                                    precursor_scaled.add_mz_noise_uniform(fragment_ppm, right_drag_val)
714                                } else {
715                                    precursor_scaled.add_mz_noise_normal(fragment_ppm)
716                                }
717                            } else {
718                                precursor_scaled
719                            };
720
721                            let precursor_len = precursor_mz_spectrum.mz.len();
722                            tims_spectra.push(TimsSpectrum::new(
723                                frame_id as i32,
724                                *scan as i32,
725                                rt,
726                                scan_mobility,
727                                ms_type.clone(),
728                                IndexedMzSpectrum::from_mz_spectrum(
729                                    vec![0; precursor_len],
730                                    precursor_mz_spectrum,
731                                ).filter_ranged(100.0, 1700.0, 1.0, 1e9),
732                            ));
733                        }
734                    }
735                }
736            }
737        }
738
739        if tims_spectra.is_empty() {
740            return TimsFrame::new(frame_id as i32, ms_type, rt, vec![], vec![], vec![], vec![], vec![]);
741        }
742
743        let tims_frame = TimsFrame::from_tims_spectra(tims_spectra);
744        tims_frame.filter_ranged(
745            mz_min_val,
746            mz_max_val,
747            0,
748            1000,
749            0.0,
750            10.0,
751            intensity_min_val,
752            1e9,
753            0,
754            i32::MAX,
755        )
756    }
757
758    pub fn build_fragment_frame_annotated(
759        &self,
760        frame_id: u32,
761        fragment_ions: &BTreeMap<
762            (u32, i8, i32),
763            (PeptideProductIonSeriesCollection, Vec<MzSpectrumAnnotated>),
764        >,
765        mz_noise_fragment: bool,
766        uniform: bool,
767        fragment_ppm: f64,
768        mz_min: Option<f64>,
769        mz_max: Option<f64>,
770        intensity_min: Option<f64>,
771        right_drag: Option<bool>,
772    ) -> TimsFrameAnnotated {
773        // Cache frame-level lookups
774        let ms_type = if self.precursor_frame_builder.precursor_frame_id_set.contains(&frame_id) {
775            MsType::Unknown
776        } else {
777            MsType::FragmentDda
778        };
779
780        let rt = *self.precursor_frame_builder.frame_to_rt.get(&frame_id).unwrap() as f64;
781        let right_drag_val = right_drag.unwrap_or(false);
782        let mz_min_val = mz_min.unwrap_or(100.0);
783        let mz_max_val = mz_max.unwrap_or(1700.0);
784        let intensity_min_val = intensity_min.unwrap_or(1.0);
785
786        // Get PASEF meta for this frame - these define which precursors were SELECTED
787        let Some(pasef_meta) = self.transmission_settings.pasef_meta.get(&(frame_id as i32)) else {
788            return TimsFrameAnnotated::new(frame_id as i32, rt, ms_type, vec![], vec![], vec![], vec![], vec![], vec![]);
789        };
790
791        // Preallocate with estimated capacity
792        let estimated_capacity = pasef_meta.len() * 4;
793        let mut tims_spectra: Vec<TimsSpectrumAnnotated> = Vec::with_capacity(estimated_capacity);
794
795        // Iterate over PASEF meta entries - each represents a SELECTED precursor
796        for meta in pasef_meta.iter() {
797            // Get the selected ion_id from the PASEF precursor field
798            let ion_id = meta.precursor as u32;
799
800            // Look up peptide_id and charge from the ion_id
801            let Some(&(peptide_id, charge_state)) = self.precursor_frame_builder.ion_id_to_peptide_charge.get(&ion_id) else {
802                continue;
803            };
804
805            // Get peptide info
806            let Some(peptide) = self.precursor_frame_builder.peptides.get(&peptide_id) else {
807                continue;
808            };
809
810            // Check if this peptide is active in this frame
811            if frame_id < peptide.frame_start || frame_id > peptide.frame_end {
812                continue;
813            }
814
815            // Get ion info from peptide_to_ions
816            let Some((ion_abundances, scan_occurrences, scan_abundances, charges, _)) = self
817                .precursor_frame_builder
818                .peptide_to_ions
819                .get(&peptide_id)
820            else {
821                continue;
822            };
823
824            // Find the index of this specific charge state
825            let Some(ion_index) = charges.iter().position(|&c| c == charge_state) else {
826                continue;
827            };
828
829            let ion_abundance = ion_abundances[ion_index];
830            let all_scan_occurrence = &scan_occurrences[ion_index];
831            let all_scan_abundance = &scan_abundances[ion_index];
832            let total_events = *self.precursor_frame_builder.peptide_to_events.get(&peptide_id).unwrap();
833
834            // Get frame abundance for this peptide
835            let frame_abundance = self.precursor_frame_builder.frame_to_abundances
836                .get(&frame_id)
837                .and_then(|(pep_ids, abundances)| {
838                    pep_ids.iter().position(|&p| p == peptide_id)
839                        .map(|idx| abundances[idx])
840                })
841                .unwrap_or(0.0);
842
843            if frame_abundance == 0.0 {
844                continue;
845            }
846
847            // Collision energy from the PASEF meta
848            let collision_energy = meta.collision_energy;
849            let collision_energy_quantized = (collision_energy * 1e1).round() as i32;
850
851            // Look up fragment ions for this (peptide_id, charge, CE)
852            let Some((_, fragment_series_vec)) = fragment_ions.get(&(peptide_id, charge_state, collision_energy_quantized)) else {
853                continue;
854            };
855
856            // Create ion for annotation and precursor spectrum calculation
857            let ion = PeptideIon::new(
858                peptide.sequence.sequence.clone(),
859                charge_state as i32,
860                ion_abundance as f64,
861                Some(peptide_id as i32),
862            );
863            // Calculate isotopic spectrum for precursor survival
864            let precursor_spectrum_annotated = ion.calculate_isotopic_spectrum_annotated(1e-3, 1e-8, 200, 1e-4);
865
866            // Process scans within the PASEF selection window
867            for (scan, scan_abundance) in all_scan_occurrence.iter().zip(all_scan_abundance.iter()) {
868                // Check if scan is within the PASEF selection window
869                let scan_i32 = *scan as i32;
870                if scan_i32 < meta.scan_start || scan_i32 > meta.scan_end {
871                    continue;
872                }
873
874                // Calculate abundance factor
875                let fraction_events = frame_abundance * scan_abundance * ion_abundance * total_events;
876
877                // Cache scan mobility
878                let scan_mobility = *self.precursor_frame_builder.scan_to_mobility.get(scan).unwrap() as f64;
879
880                for fragment_ion_series in fragment_series_vec.iter() {
881                    let scaled_spec = fragment_ion_series.clone() * fraction_events as f64;
882
883                    let mz_spectrum = if mz_noise_fragment {
884                        if uniform {
885                            scaled_spec.add_mz_noise_uniform(fragment_ppm, right_drag_val)
886                        } else {
887                            scaled_spec.add_mz_noise_normal(fragment_ppm)
888                        }
889                    } else {
890                        scaled_spec
891                    };
892
893                    let spectrum_len = mz_spectrum.mz.len();
894                    tims_spectra.push(TimsSpectrumAnnotated::new(
895                        frame_id as i32,
896                        *scan,
897                        rt,
898                        scan_mobility,
899                        ms_type.clone(),
900                        vec![0; spectrum_len],
901                        mz_spectrum,
902                    ));
903                }
904
905                // Add unfragmented precursor ions (survival) if configured
906                if self.isotope_transmission_config.has_precursor_survival() {
907                    let mut rng = rand::thread_rng();
908                    let survival_fraction = rng.gen_range(
909                        self.isotope_transmission_config.precursor_survival_min
910                        ..=self.isotope_transmission_config.precursor_survival_max
911                    );
912
913                    if survival_fraction > 0.0 {
914                        // Create a non-annotated spectrum for transmission
915                        let precursor_mz_spectrum = MzSpectrum::new(
916                            precursor_spectrum_annotated.mz.clone(),
917                            precursor_spectrum_annotated.intensity.clone(),
918                        );
919
920                        // Transmit through the quadrupole
921                        let precursor_transmitted = self.transmission_settings.transmit_spectrum(
922                            frame_id as i32,
923                            scan_i32,
924                            precursor_mz_spectrum,
925                            Some(self.isotope_transmission_config.min_probability),
926                        );
927
928                        if !precursor_transmitted.mz.is_empty() {
929                            // Scale by survival fraction and event count
930                            let precursor_scaled = precursor_transmitted * (fraction_events as f64 * survival_fraction);
931
932                            let precursor_final = if mz_noise_fragment {
933                                if uniform {
934                                    precursor_scaled.add_mz_noise_uniform(fragment_ppm, right_drag_val)
935                                } else {
936                                    precursor_scaled.add_mz_noise_normal(fragment_ppm)
937                                }
938                            } else {
939                                precursor_scaled
940                            };
941
942                            // Convert to annotated spectrum (with precursor annotations)
943                            let annotations: Vec<PeakAnnotation> = precursor_final.mz.iter()
944                                .map(|_| PeakAnnotation { contributions: vec![] })
945                                .collect();
946                            let precursor_annotated = MzSpectrumAnnotated::new(
947                                precursor_final.mz.to_vec(),
948                                precursor_final.intensity.to_vec(),
949                                annotations,
950                            );
951
952                            let precursor_len = precursor_annotated.mz.len();
953                            tims_spectra.push(TimsSpectrumAnnotated::new(
954                                frame_id as i32,
955                                *scan,
956                                rt,
957                                scan_mobility,
958                                ms_type.clone(),
959                                vec![0; precursor_len],
960                                precursor_annotated,
961                            ));
962                        }
963                    }
964                }
965            }
966        }
967
968        if tims_spectra.is_empty() {
969            return TimsFrameAnnotated::new(frame_id as i32, rt, ms_type, vec![], vec![], vec![], vec![], vec![], vec![]);
970        }
971
972        TimsFrameAnnotated::from_tims_spectra_annotated(tims_spectra).filter_ranged(
973            mz_min_val, mz_max_val, 0.0, 10.0, 0, 1000, intensity_min_val, 1e9,
974        )
975    }
976
977    pub fn get_collision_energy(&self, frame_id: i32, scan_id: i32) -> f64 {
978        self.transmission_settings.get_collision_energy(frame_id, scan_id).unwrap_or(0.0)
979    }
980
981    pub fn get_collision_energies(&self, frame_ids: Vec<i32>, scan_ids: Vec<i32>) -> Vec<f64> {
982        let mut collision_energies: Vec<f64> = Vec::new();
983        for frame_id in frame_ids {
984            for scan_id in &scan_ids {
985                collision_energies.push(self.get_collision_energy(frame_id, *scan_id));
986            }
987        }
988        collision_energies
989    }
990}
rustdf/sim/dda.rs

rustdf/sim/
dda.rs
