1use mscore::simulation::noise_rng::noise_rng;
2use mscore::algorithm::isotope::{
3 calculate_precursor_transmission_factor,
4 calculate_transmission_dependent_fragment_ion_isotope_distribution,
5};
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::collision::{TimsTofCollisionEnergy, TimsTofCollisionEnergyDIA};
12use mscore::timstof::frame::TimsFrame;
13use mscore::timstof::quadrupole::{IonTransmission, TimsTransmissionDIA, WindowTransmission};
14use mscore::timstof::spectrum::TimsSpectrum;
15use std::collections::{BTreeMap, HashSet};
16use std::path::Path;
17use std::sync::Arc;
18
19use rand::Rng;
20use rayon::prelude::*;
21use rayon::ThreadPoolBuilder;
22
23use crate::sim::containers::{IsotopeTransmissionConfig, IsotopeTransmissionMode};
24use crate::sim::projector::{IntensityStage, MzCoordSpace, RenderedEvent, RenderedSpectrum};
25use crate::sim::scheme::{DataMode, InstrumentCapabilities, IsolationWindow};
26use crate::sim::handle::{FragmentIonsWithComplementary, TimsTofSyntheticsDataHandle};
27use crate::sim::precursor::TimsTofSyntheticsPrecursorFrameBuilder;
28
29#[allow(clippy::too_many_arguments)]
39pub(crate) fn fragment_series_spectrum(
40 mode: IsotopeTransmissionMode,
41 fragment_ion_series: &MzSpectrum,
42 series_idx: usize,
43 fraction_events: f32,
44 transmission_factor: f64,
45 frag_data: Option<&FragmentIonsWithComplementary>,
46 transmitted_indices: &HashSet<usize>,
47 max_isotopes: usize,
48) -> MzSpectrum {
49 match mode {
50 IsotopeTransmissionMode::None => fragment_ion_series.clone() * fraction_events as f64,
51 IsotopeTransmissionMode::PrecursorScaling => {
52 fragment_ion_series.clone() * (fraction_events as f64 * transmission_factor)
53 }
54 IsotopeTransmissionMode::PerFragment => {
55 if let Some(frag_data) = frag_data {
56 if series_idx < frag_data.per_fragment_data.len() {
57 let series_data = &frag_data.per_fragment_data[series_idx];
58 let mut aggregated_mz: Vec<f64> = Vec::new();
59 let mut aggregated_intensity: Vec<f64> = Vec::new();
60 for frag_ion_data in series_data {
61 let adjusted_dist =
62 calculate_transmission_dependent_fragment_ion_isotope_distribution(
63 &frag_ion_data.fragment_distribution,
64 &frag_ion_data.complementary_distribution,
65 transmitted_indices,
66 max_isotopes,
67 );
68 for (mz, abundance) in adjusted_dist {
69 aggregated_mz.push(mz);
70 aggregated_intensity
71 .push(abundance * frag_ion_data.predicted_intensity * fraction_events as f64);
72 }
73 }
74 if !aggregated_mz.is_empty() {
75 MzSpectrum::new(aggregated_mz, aggregated_intensity)
76 } else {
77 fragment_ion_series.clone() * fraction_events as f64
78 }
79 } else {
80 fragment_ion_series.clone() * fraction_events as f64
81 }
82 } else {
83 fragment_ion_series.clone() * fraction_events as f64
84 }
85 }
86 }
87}
88
89pub struct TimsTofSyntheticsFrameBuilderDIA {
90 pub path: String,
91 pub precursor_frame_builder: TimsTofSyntheticsPrecursorFrameBuilder,
92 pub transmission_settings: TimsTransmissionDIA,
93 pub fragmentation_settings: TimsTofCollisionEnergyDIA,
94 pub fragment_ions:
95 Option<BTreeMap<(u32, i8, i32), (PeptideProductIonSeriesCollection, Vec<MzSpectrum>)>>,
96 pub fragment_ions_annotated: Option<
97 BTreeMap<(u32, i8, i32), (PeptideProductIonSeriesCollection, Vec<MzSpectrumAnnotated>)>,
98 >,
99 pub isotope_config: IsotopeTransmissionConfig,
102 pub fragment_ions_with_transmission:
104 Option<BTreeMap<(u32, i8, i32), FragmentIonsWithComplementary>>,
105 pub capabilities: InstrumentCapabilities,
107 pub noise_seed: u64,
110}
111
112impl TimsTofSyntheticsFrameBuilderDIA {
113
114 pub fn set_noise_seed(&mut self, seed: u64) {
116 self.noise_seed = seed;
117 }
118 pub fn new(path: &Path, with_annotations: bool, num_threads: usize) -> rusqlite::Result<Self> {
119 Self::new_with_config(path, with_annotations, num_threads, IsotopeTransmissionConfig::default())
120 }
121
122 pub fn new_with_config(
123 path: &Path,
124 with_annotations: bool,
125 num_threads: usize,
126 isotope_config: IsotopeTransmissionConfig,
127 ) -> rusqlite::Result<Self> {
128 Self::new_with_config_and_source(
129 path,
130 with_annotations,
131 num_threads,
132 isotope_config,
133 &crate::sim::projector::DistributionSource::Columns,
134 )
135 }
136
137 pub fn new_with_config_and_source(
140 path: &Path,
141 with_annotations: bool,
142 num_threads: usize,
143 isotope_config: IsotopeTransmissionConfig,
144 source: &crate::sim::projector::DistributionSource,
145 ) -> rusqlite::Result<Self> {
146 let synthetics = TimsTofSyntheticsPrecursorFrameBuilder::from_source(path, source)?;
147 let handle = TimsTofSyntheticsDataHandle::new(path)?;
148
149 handle
152 .read_prediction_set()?
153 .assert_render_compatible()
154 .map_err(|_| rusqlite::Error::InvalidQuery)?;
155
156 let fragment_ions = handle.read_fragment_ions()?;
157
158 let capabilities = InstrumentCapabilities::default();
161 let isotope_config = isotope_config.gated_by(capabilities);
162
163 let fragmentation_settings = handle.get_collision_energy_dia();
165 let transmission_settings = handle.get_transmission_dia();
167
168 let fragment_ions_with_transmission = if isotope_config.is_enabled() {
170 Some(TimsTofSyntheticsDataHandle::build_fragment_ions_with_transmission_data(
171 &synthetics.peptides,
172 &fragment_ions,
173 num_threads,
174 ))
175 } else {
176 None
177 };
178
179 match with_annotations {
180 true => {
181 let fragment_ions_annotated =
182 Some(TimsTofSyntheticsDataHandle::build_fragment_ions_annotated(
183 &synthetics.peptides,
184 &fragment_ions,
185 num_threads,
186 ));
187 Ok(Self {
188 noise_seed: 0,
189 path: path.to_str().unwrap().to_string(),
190 precursor_frame_builder: synthetics,
191 transmission_settings,
192 fragmentation_settings,
193 fragment_ions: None,
194 fragment_ions_annotated,
195 isotope_config,
196 fragment_ions_with_transmission,
197 capabilities,
198 })
199 }
200
201 false => {
202 let fragment_ions = Some(TimsTofSyntheticsDataHandle::build_fragment_ions(
203 &synthetics.peptides,
204 &fragment_ions,
205 num_threads,
206 ));
207 Ok(Self {
208 noise_seed: 0,
209 path: path.to_str().unwrap().to_string(),
210 precursor_frame_builder: synthetics,
211 transmission_settings,
212 fragmentation_settings,
213 fragment_ions,
214 fragment_ions_annotated: None,
215 isotope_config,
216 fragment_ions_with_transmission,
217 capabilities,
218 })
219 }
220 }
221 }
222
223 pub fn build_frame(
235 &self,
236 frame_id: u32,
237 fragmentation: bool,
238 mz_noise_precursor: bool,
239 uniform: bool,
240 precursor_noise_ppm: f64,
241 mz_noise_fragment: bool,
242 fragment_noise_ppm: f64,
243 right_drag: bool,
244 ) -> TimsFrame {
245 match self
247 .precursor_frame_builder
248 .precursor_frame_id_set
249 .contains(&frame_id)
250 {
251 true => self.build_ms1_frame(
252 frame_id,
253 mz_noise_precursor,
254 uniform,
255 precursor_noise_ppm,
256 right_drag,
257 ),
258 false => self.build_ms2_frame(
259 frame_id,
260 fragmentation,
261 mz_noise_fragment,
262 uniform,
263 fragment_noise_ppm,
264 right_drag,
265 ),
266 }
267 }
268
269 pub fn build_frame_annotated(
270 &self,
271 frame_id: u32,
272 fragmentation: bool,
273 mz_noise_precursor: bool,
274 uniform: bool,
275 precursor_noise_ppm: f64,
276 mz_noise_fragment: bool,
277 fragment_noise_ppm: f64,
278 right_drag: bool,
279 ) -> TimsFrameAnnotated {
280 match self
281 .precursor_frame_builder
282 .precursor_frame_id_set
283 .contains(&frame_id)
284 {
285 true => self.build_ms1_frame_annotated(
286 frame_id,
287 mz_noise_precursor,
288 uniform,
289 precursor_noise_ppm,
290 right_drag,
291 ),
292 false => self.build_ms2_frame_annotated(
293 frame_id,
294 fragmentation,
295 mz_noise_fragment,
296 uniform,
297 fragment_noise_ppm,
298 right_drag,
299 ),
300 }
301 }
302
303 pub fn get_fragment_ion_ids(&self, precursor_frame_ids: Vec<u32>) -> Vec<u32> {
304 let mut peptide_ids: HashSet<u32> = HashSet::new();
305 for frame_id in precursor_frame_ids {
307 for (peptide_id, peptide) in self.precursor_frame_builder.peptides.iter() {
308 if peptide.frame_start <= frame_id && peptide.frame_end >= frame_id {
309 peptide_ids.insert(*peptide_id);
310 }
311 }
312 }
313 let mut result: Vec<u32> = Vec::new();
315 for item in peptide_ids {
316 let ions = self.precursor_frame_builder.ions.get(&item).unwrap();
317 for ion in ions.iter() {
318 result.push(ion.ion_id);
319 }
320 }
321 result
322 }
323
324 pub fn build_frames(
325 &self,
326 frame_ids: Vec<u32>,
327 fragmentation: bool,
328 mz_noise_precursor: bool,
329 uniform: bool,
330 precursor_noise_ppm: f64,
331 mz_noise_fragment: bool,
332 fragment_noise_ppm: f64,
333 right_drag: bool,
334 num_threads: usize,
335 ) -> Vec<TimsFrame> {
336 let pool = rayon::ThreadPoolBuilder::new()
338 .num_threads(num_threads)
339 .build()
340 .unwrap();
341
342 pool.install(|| {
343 let mut tims_frames: Vec<TimsFrame> = Vec::with_capacity(frame_ids.len());
345 unsafe { tims_frames.set_len(frame_ids.len()); }
347
348 frame_ids.par_iter().enumerate().for_each(|(idx, frame_id)| {
349 let frame = self.build_frame(
350 *frame_id,
351 fragmentation,
352 mz_noise_precursor,
353 uniform,
354 precursor_noise_ppm,
355 mz_noise_fragment,
356 fragment_noise_ppm,
357 right_drag,
358 );
359 unsafe {
361 let ptr = tims_frames.as_ptr() as *mut TimsFrame;
362 std::ptr::write(ptr.add(idx), frame);
363 }
364 });
365
366 tims_frames
367 })
368 }
369
370 pub fn build_frames_annotated(
371 &self,
372 frame_ids: Vec<u32>,
373 fragmentation: bool,
374 mz_noise_precursor: bool,
375 uniform: bool,
376 precursor_noise_ppm: f64,
377 mz_noise_fragment: bool,
378 fragment_noise_ppm: f64,
379 right_drag: bool,
380 num_threads: usize,
381 ) -> Vec<TimsFrameAnnotated> {
382 let pool = rayon::ThreadPoolBuilder::new()
384 .num_threads(num_threads)
385 .build()
386 .unwrap();
387
388 pool.install(|| {
389 let mut tims_frames: Vec<TimsFrameAnnotated> = Vec::with_capacity(frame_ids.len());
391 unsafe { tims_frames.set_len(frame_ids.len()); }
392
393 frame_ids.par_iter().enumerate().for_each(|(idx, frame_id)| {
394 let frame = self.build_frame_annotated(
395 *frame_id,
396 fragmentation,
397 mz_noise_precursor,
398 uniform,
399 precursor_noise_ppm,
400 mz_noise_fragment,
401 fragment_noise_ppm,
402 right_drag,
403 );
404 unsafe {
405 let ptr = tims_frames.as_ptr() as *mut TimsFrameAnnotated;
406 std::ptr::write(ptr.add(idx), frame);
407 }
408 });
409
410 tims_frames
411 })
412 }
413
414 fn build_ms1_frame(
415 &self,
416 frame_id: u32,
417 mz_noise_precursor: bool,
418 uniform: bool,
419 precursor_ppm: f64,
420 right_drag: bool,
421 ) -> TimsFrame {
422 let mut tims_frame = self.precursor_frame_builder.build_precursor_frame(
423 frame_id,
424 mz_noise_precursor,
425 uniform,
426 precursor_ppm,
427 right_drag,
428 );
429 let intensities_rounded = tims_frame
430 .ims_frame
431 .intensity
432 .iter()
433 .map(|x| x.round())
434 .collect::<Vec<_>>();
435 tims_frame.ims_frame.intensity = Arc::new(intensities_rounded);
436 tims_frame
437 }
438
439 fn build_ms1_frame_annotated(
440 &self,
441 frame_id: u32,
442 mz_noise_precursor: bool,
443 uniform: bool,
444 precursor_ppm: f64,
445 right_drag: bool,
446 ) -> TimsFrameAnnotated {
447 let mut tims_frame = self
448 .precursor_frame_builder
449 .build_precursor_frame_annotated(
450 frame_id,
451 mz_noise_precursor,
452 uniform,
453 precursor_ppm,
454 right_drag,
455 );
456 let intensities_rounded = tims_frame
457 .intensity
458 .iter()
459 .map(|x| x.round())
460 .collect::<Vec<_>>();
461 tims_frame.intensity = intensities_rounded;
462 tims_frame
463 }
464
465 fn build_ms2_frame(
466 &self,
467 frame_id: u32,
468 fragmentation: bool,
469 mz_noise_fragment: bool,
470 uniform: bool,
471 fragment_ppm: f64,
472 right_drag: bool,
473 ) -> TimsFrame {
474 match fragmentation {
475 false => {
476 let mut frame = self.transmission_settings.transmit_tims_frame(
477 &self.build_ms1_frame(
478 frame_id,
479 mz_noise_fragment,
480 uniform,
481 fragment_ppm,
482 right_drag,
483 ),
484 None,
485 );
486 let intensities_rounded = frame
487 .ims_frame
488 .intensity
489 .iter()
490 .map(|x| x.round())
491 .collect::<Vec<_>>();
492 frame.ims_frame.intensity = Arc::new(intensities_rounded);
493 frame.ms_type = MsType::FragmentDia;
494 frame
495 }
496 true => {
497 let mut frame = self.build_fragment_frame(
498 frame_id,
499 &self.fragment_ions.as_ref().unwrap(),
500 mz_noise_fragment,
501 uniform,
502 fragment_ppm,
503 None,
504 None,
505 None,
506 Some(right_drag),
507 );
508 let intensities_rounded = frame
509 .ims_frame
510 .intensity
511 .iter()
512 .map(|x| x.round())
513 .collect::<Vec<_>>();
514 frame.ims_frame.intensity = Arc::new(intensities_rounded);
515 frame
516 }
517 }
518 }
519
520 fn build_ms2_frame_annotated(
521 &self,
522 frame_id: u32,
523 fragmentation: bool,
524 mz_noise_fragment: bool,
525 uniform: bool,
526 fragment_ppm: f64,
527 right_drag: bool,
528 ) -> TimsFrameAnnotated {
529 match fragmentation {
530 false => {
531 let mut frame = self.transmission_settings.transmit_tims_frame_annotated(
532 &self.build_ms1_frame_annotated(
533 frame_id,
534 mz_noise_fragment,
535 uniform,
536 fragment_ppm,
537 right_drag,
538 ),
539 None,
540 );
541 let intensities_rounded = frame
542 .intensity
543 .iter()
544 .map(|x| x.round())
545 .collect::<Vec<_>>();
546 frame.intensity = intensities_rounded;
547 frame.ms_type = MsType::FragmentDia;
548 frame
549 }
550 true => {
551 let mut frame = self.build_fragment_frame_annotated(
552 frame_id,
553 &self.fragment_ions_annotated.as_ref().unwrap(),
554 mz_noise_fragment,
555 uniform,
556 fragment_ppm,
557 None,
558 None,
559 None,
560 Some(right_drag),
561 );
562 let intensities_rounded = frame
563 .intensity
564 .iter()
565 .map(|x| x.round())
566 .collect::<Vec<_>>();
567 frame.intensity = intensities_rounded;
568 frame
569 }
570 }
571 }
572
573 fn build_fragment_frame(
587 &self,
588 frame_id: u32,
589 fragment_ions: &BTreeMap<
590 (u32, i8, i32),
591 (PeptideProductIonSeriesCollection, Vec<MzSpectrum>),
592 >,
593 mz_noise_fragment: bool,
594 uniform: bool,
595 fragment_ppm: f64,
596 mz_min: Option<f64>,
597 mz_max: Option<f64>,
598 intensity_min: Option<f64>,
599 right_drag: Option<bool>,
600 ) -> TimsFrame {
601 let mut rng = noise_rng(self.noise_seed, &[frame_id as u64]);
604 let ms_type = if self.precursor_frame_builder.precursor_frame_id_set.contains(&frame_id) {
606 MsType::Unknown
607 } else {
608 MsType::FragmentDia
609 };
610
611 let rt = *self.precursor_frame_builder.frame_to_rt.get(&frame_id).unwrap() as f64;
612 let right_drag_val = right_drag.unwrap_or(false);
613 let mz_min_val = mz_min.unwrap_or(100.0);
614 let mz_max_val = mz_max.unwrap_or(1700.0);
615 let intensity_min_val = intensity_min.unwrap_or(1.0);
616
617 let Some((peptide_ids, frame_abundances)) = self
619 .precursor_frame_builder
620 .frame_to_abundances
621 .get(&frame_id)
622 else {
623 return TimsFrame::new(frame_id as i32, ms_type, rt, vec![], vec![], vec![], vec![], vec![]);
624 };
625
626 let estimated_capacity = peptide_ids.len() * 4;
628 let mut tims_spectra: Vec<TimsSpectrum> = Vec::with_capacity(estimated_capacity);
629
630 for (peptide_id, frame_abundance) in peptide_ids.iter().zip(frame_abundances.iter()) {
632 let Some((ion_abundances, scan_occurrences, scan_abundances, charges, spectra)) = self
634 .precursor_frame_builder
635 .peptide_to_ions
636 .get(peptide_id)
637 else {
638 continue;
639 };
640
641 let total_events = *self.precursor_frame_builder.peptide_to_events.get(peptide_id).unwrap();
643
644 for (index, ion_abundance) in ion_abundances.iter().enumerate() {
645 let all_scan_occurrence = &scan_occurrences[index];
646 let all_scan_abundance = &scan_abundances[index];
647 let spectrum = &spectra[index];
648 let charge_state = charges[index];
649
650 for (scan, scan_abundance) in all_scan_occurrence.iter().zip(all_scan_abundance.iter()) {
651 let (is_transmitted, transmitted_indices) = match self.isotope_config.mode {
653 IsotopeTransmissionMode::None => {
654 let any = self.transmission_settings.any_transmitted(
656 frame_id as i32,
657 *scan as i32,
658 &spectrum.mz,
659 None,
660 );
661 (any, HashSet::new())
662 },
663 IsotopeTransmissionMode::PrecursorScaling | IsotopeTransmissionMode::PerFragment => {
664 let indices = self.transmission_settings.get_transmission_set(
665 frame_id as i32,
666 *scan as i32,
667 &spectrum.mz,
668 Some(self.isotope_config.min_probability),
669 );
670 (!indices.is_empty(), indices)
671 },
672 };
673
674 if !is_transmitted {
675 continue;
676 }
677
678 let fraction_events = frame_abundance * scan_abundance * ion_abundance * total_events;
680
681 let collision_energy = self.fragmentation_settings.get_collision_energy(frame_id as i32, *scan as i32);
683 let Some(collision_energy_quantized) = crate::sim::handle::resolve_fragment_ce_key(
686 fragment_ions, *peptide_id, charge_state, collision_energy,
687 ) else {
688 if crate::sim::handle::fragment_prefix_exists(fragment_ions, *peptide_id, charge_state) {
692 panic!(
693 "DIA fragment lookup miss: peptide {} charge {} applied CE {:.4} eV \
694 has predicted fragments, but none within 0.1 eV — the prediction \
695 set does not cover this instrument's collision energy",
696 *peptide_id, charge_state, collision_energy,
697 );
698 }
699 continue;
700 };
701 let (_, fragment_series_vec) = fragment_ions
702 .get(&(*peptide_id, charge_state, collision_energy_quantized))
703 .expect("resolve_fragment_ce_key returned a present key");
704
705 let scan_mobility = *self.precursor_frame_builder.scan_to_mobility.get(scan).unwrap() as f64;
707
708 let transmission_factor = if self.isotope_config.mode == IsotopeTransmissionMode::PrecursorScaling {
710 if let Some(comp_data) = self.fragment_ions_with_transmission.as_ref() {
711 if let Some(frag_data) = comp_data.get(&(*peptide_id, charge_state, collision_energy_quantized)) {
712 calculate_precursor_transmission_factor(
713 &frag_data.precursor_isotope_distribution,
714 &transmitted_indices,
715 )
716 } else {
717 1.0
718 }
719 } else {
720 1.0
721 }
722 } else {
723 1.0
724 };
725
726 let frag_data = self
729 .fragment_ions_with_transmission
730 .as_ref()
731 .and_then(|c| c.get(&(*peptide_id, charge_state, collision_energy_quantized)));
732
733 for (series_idx, fragment_ion_series) in fragment_series_vec.iter().enumerate() {
734 let final_spectrum = fragment_series_spectrum(
735 self.isotope_config.mode,
736 fragment_ion_series,
737 series_idx,
738 fraction_events,
739 transmission_factor,
740 frag_data,
741 &transmitted_indices,
742 self.isotope_config.max_isotopes,
743 );
744
745 let mz_spectrum = if mz_noise_fragment {
746 if uniform {
747 final_spectrum.add_mz_noise_uniform_with_rng(fragment_ppm, right_drag_val, &mut rng)
748 } else {
749 final_spectrum.add_mz_noise_normal_with_rng(fragment_ppm, &mut rng)
750 }
751 } else {
752 final_spectrum
753 };
754
755 let spectrum_len = mz_spectrum.mz.len();
756 tims_spectra.push(TimsSpectrum::new(
757 frame_id as i32,
758 *scan as i32,
759 rt,
760 scan_mobility,
761 ms_type.clone(),
762 IndexedMzSpectrum::from_mz_spectrum(
763 vec![0; spectrum_len],
764 mz_spectrum,
765 ).filter_ranged(100.0, 1700.0, 1.0, 1e9),
766 ));
767 }
768
769 if self.isotope_config.has_precursor_survival() {
771 let mut rng = rand::thread_rng();
772 let survival_fraction = rng.gen_range(
773 self.isotope_config.precursor_survival_min
774 ..=self.isotope_config.precursor_survival_max
775 );
776
777 if survival_fraction > 0.0 {
778 let precursor_transmitted = self.transmission_settings.transmit_spectrum(
780 frame_id as i32,
781 *scan as i32,
782 spectrum.clone(),
783 Some(self.isotope_config.min_probability),
784 );
785
786 if !precursor_transmitted.mz.is_empty() {
787 let precursor_scaled = precursor_transmitted * (fraction_events as f64 * survival_fraction);
789
790 let precursor_mz_spectrum = if mz_noise_fragment {
791 if uniform {
792 precursor_scaled.add_mz_noise_uniform_with_rng(fragment_ppm, right_drag_val, &mut rng)
793 } else {
794 precursor_scaled.add_mz_noise_normal_with_rng(fragment_ppm, &mut rng)
795 }
796 } else {
797 precursor_scaled
798 };
799
800 let precursor_len = precursor_mz_spectrum.mz.len();
801 tims_spectra.push(TimsSpectrum::new(
802 frame_id as i32,
803 *scan as i32,
804 rt,
805 scan_mobility,
806 ms_type.clone(),
807 IndexedMzSpectrum::from_mz_spectrum(
808 vec![0; precursor_len],
809 precursor_mz_spectrum,
810 ).filter_ranged(100.0, 1700.0, 1.0, 1e9),
811 ));
812 }
813 }
814 }
815 }
816 }
817 }
818
819 if tims_spectra.is_empty() {
820 return TimsFrame::new(frame_id as i32, ms_type, rt, vec![], vec![], vec![], vec![], vec![]);
821 }
822
823 let tims_frame = TimsFrame::from_tims_spectra(tims_spectra);
824 tims_frame.filter_ranged(
825 mz_min_val,
826 mz_max_val,
827 0,
828 1000,
829 0.0,
830 10.0,
831 intensity_min_val,
832 1e9,
833 0,
834 i32::MAX,
835 )
836 }
837
838 pub fn render_fragment_scan(
865 &self,
866 frame_id: u32,
867 window: &WindowTransmission,
868 nce: f64,
869 data_mode: DataMode,
870 ) -> RenderedEvent {
871 let rt = *self
872 .precursor_frame_builder
873 .frame_to_rt
874 .get(&frame_id)
875 .unwrap_or(&0.0) as f64;
876 let isolation = Some(IsolationWindow {
877 center_mz: window.center_mz,
878 width_mz: window.width_mz,
879 });
880 let make_scan = |peaks: MzSpectrum| RenderedEvent::Scan {
881 ms_level: 2,
882 retention_time_s: rt,
883 isolation,
884 spectrum: RenderedSpectrum {
885 mz: (*peaks.mz).clone(),
886 intensity: (*peaks.intensity).clone(),
887 coords: MzCoordSpace::Physical,
888 mode: data_mode,
889 detector_applied: false,
890 stage: IntensityStage::Transmitted,
891 },
892 };
893
894 let (Some(fragment_ions), Some((peptide_ids, frame_abundances))) = (
895 self.fragment_ions.as_ref(),
896 self.precursor_frame_builder.frame_to_abundances.get(&frame_id),
897 ) else {
898 return make_scan(MzSpectrum::from_collection(vec![]));
899 };
900
901 let mut specs: Vec<MzSpectrum> = Vec::new();
902 for (peptide_id, frame_abundance) in peptide_ids.iter().zip(frame_abundances.iter()) {
903 let Some((ion_abundances, _scan_occ, _scan_abu, charges, spectra)) =
904 self.precursor_frame_builder.peptide_to_ions.get(peptide_id)
905 else {
906 continue;
907 };
908 let total_events = *self
909 .precursor_frame_builder
910 .peptide_to_events
911 .get(peptide_id)
912 .unwrap();
913 for (index, ion_abundance) in ion_abundances.iter().enumerate() {
914 let spectrum = &spectra[index];
915 let charge_state = charges[index];
916 if !window.any_transmitted(&spectrum.mz, None) {
918 continue;
919 }
920 let Some(ce_key) = crate::sim::handle::resolve_fragment_ce_key(
928 fragment_ions,
929 *peptide_id,
930 charge_state,
931 nce,
932 ) else {
933 continue;
934 };
935 let (_, fragment_series_vec) = fragment_ions
936 .get(&(*peptide_id, charge_state, ce_key))
937 .expect("resolve_fragment_ce_key returned a present key");
938
939 let fraction_events = frame_abundance * ion_abundance * total_events;
941 for (series_idx, fragment_ion_series) in fragment_series_vec.iter().enumerate() {
942 specs.push(fragment_series_spectrum(
943 IsotopeTransmissionMode::None,
944 fragment_ion_series,
945 series_idx,
946 fraction_events,
947 1.0,
948 None,
949 &HashSet::new(),
950 self.isotope_config.max_isotopes,
951 ));
952 }
953 }
954 }
955 make_scan(MzSpectrum::from_collection(specs))
956 }
957
958 pub fn build_fragment_frame_annotated(
959 &self,
960 frame_id: u32,
961 fragment_ions: &BTreeMap<
962 (u32, i8, i32),
963 (PeptideProductIonSeriesCollection, Vec<MzSpectrumAnnotated>),
964 >,
965 mz_noise_fragment: bool,
966 uniform: bool,
967 fragment_ppm: f64,
968 mz_min: Option<f64>,
969 mz_max: Option<f64>,
970 intensity_min: Option<f64>,
971 right_drag: Option<bool>,
972 ) -> TimsFrameAnnotated {
973 let mut rng = noise_rng(self.noise_seed, &[frame_id as u64]);
976 let ms_type = if self.precursor_frame_builder.precursor_frame_id_set.contains(&frame_id) {
978 MsType::Unknown
979 } else {
980 MsType::FragmentDia
981 };
982
983 let rt = *self.precursor_frame_builder.frame_to_rt.get(&frame_id).unwrap() as f64;
984 let right_drag_val = right_drag.unwrap_or(false);
985 let mz_min_val = mz_min.unwrap_or(100.0);
986 let mz_max_val = mz_max.unwrap_or(1700.0);
987 let intensity_min_val = intensity_min.unwrap_or(1.0);
988
989 let Some((peptide_ids, frame_abundances)) = self
991 .precursor_frame_builder
992 .frame_to_abundances
993 .get(&frame_id)
994 else {
995 return TimsFrameAnnotated::new(frame_id as i32, rt, ms_type, vec![], vec![], vec![], vec![], vec![], vec![]);
996 };
997
998 let estimated_capacity = peptide_ids.len() * 4;
1000 let mut tims_spectra: Vec<TimsSpectrumAnnotated> = Vec::with_capacity(estimated_capacity);
1001
1002 for (peptide_id, frame_abundance) in peptide_ids.iter().zip(frame_abundances.iter()) {
1003 let Some((ion_abundances, scan_occurrences, scan_abundances, charges, _)) = self
1005 .precursor_frame_builder
1006 .peptide_to_ions
1007 .get(peptide_id)
1008 else {
1009 continue;
1010 };
1011
1012 let total_events = *self.precursor_frame_builder.peptide_to_events.get(peptide_id).unwrap();
1014 let peptide = self.precursor_frame_builder.peptides.get(peptide_id).unwrap();
1015
1016 for (index, ion_abundance) in ion_abundances.iter().enumerate() {
1017 let all_scan_occurrence = &scan_occurrences[index];
1018 let all_scan_abundance = &scan_abundances[index];
1019 let charge_state = charges[index];
1020
1021 let ion = PeptideIon::new(
1022 peptide.sequence.sequence.clone(),
1023 charge_state as i32,
1024 *ion_abundance as f64,
1025 Some(*peptide_id as i32),
1026 );
1027 let spectrum = ion.calculate_isotopic_spectrum_annotated(1e-3, 1e-8, 200, 1e-4);
1029
1030 for (scan, scan_abundance) in all_scan_occurrence.iter().zip(all_scan_abundance.iter()) {
1031 if !self.transmission_settings.any_transmitted(
1032 frame_id as i32,
1033 *scan as i32,
1034 &spectrum.mz,
1035 None,
1036 ) {
1037 continue;
1038 }
1039
1040 let fraction_events = frame_abundance * scan_abundance * ion_abundance * total_events;
1041
1042 let collision_energy = self.fragmentation_settings.get_collision_energy(frame_id as i32, *scan as i32);
1043 let Some(collision_energy_quantized) = crate::sim::handle::resolve_fragment_ce_key(
1045 fragment_ions, *peptide_id, charge_state, collision_energy,
1046 ) else {
1047 if crate::sim::handle::fragment_prefix_exists(fragment_ions, *peptide_id, charge_state) {
1048 panic!(
1049 "DIA (annotated) fragment lookup miss: peptide {} charge {} applied CE \
1050 {:.4} eV has predicted fragments, but none within 0.1 eV — the \
1051 prediction set does not cover this instrument's collision energy",
1052 *peptide_id, charge_state, collision_energy,
1053 );
1054 }
1055 continue;
1056 };
1057 let (_, fragment_series_vec) = fragment_ions
1058 .get(&(*peptide_id, charge_state, collision_energy_quantized))
1059 .expect("resolve_fragment_ce_key returned a present key");
1060
1061 let scan_mobility = *self.precursor_frame_builder.scan_to_mobility.get(scan).unwrap() as f64;
1063
1064 for fragment_ion_series in fragment_series_vec.iter() {
1065 let scaled_spec = fragment_ion_series.clone() * fraction_events as f64;
1066
1067 let mz_spectrum = if mz_noise_fragment {
1068 if uniform {
1069 scaled_spec.add_mz_noise_uniform_with_rng(fragment_ppm, right_drag_val, &mut rng)
1070 } else {
1071 scaled_spec.add_mz_noise_normal_with_rng(fragment_ppm, &mut rng)
1072 }
1073 } else {
1074 scaled_spec
1075 };
1076
1077 let spectrum_len = mz_spectrum.mz.len();
1078 tims_spectra.push(TimsSpectrumAnnotated::new(
1079 frame_id as i32,
1080 *scan,
1081 rt,
1082 scan_mobility,
1083 ms_type.clone(),
1084 vec![0; spectrum_len],
1085 mz_spectrum,
1086 ));
1087 }
1088
1089 if self.isotope_config.has_precursor_survival() {
1091 let mut rng = rand::thread_rng();
1092 let survival_fraction = rng.gen_range(
1093 self.isotope_config.precursor_survival_min
1094 ..=self.isotope_config.precursor_survival_max
1095 );
1096
1097 if survival_fraction > 0.0 {
1098 let precursor_mz_spectrum = MzSpectrum::new(spectrum.mz.clone(), spectrum.intensity.clone());
1100
1101 let precursor_transmitted = self.transmission_settings.transmit_spectrum(
1103 frame_id as i32,
1104 *scan as i32,
1105 precursor_mz_spectrum,
1106 Some(self.isotope_config.min_probability),
1107 );
1108
1109 if !precursor_transmitted.mz.is_empty() {
1110 let precursor_scaled = precursor_transmitted * (fraction_events as f64 * survival_fraction);
1112
1113 let precursor_final = if mz_noise_fragment {
1114 if uniform {
1115 precursor_scaled.add_mz_noise_uniform_with_rng(fragment_ppm, right_drag_val, &mut rng)
1116 } else {
1117 precursor_scaled.add_mz_noise_normal_with_rng(fragment_ppm, &mut rng)
1118 }
1119 } else {
1120 precursor_scaled
1121 };
1122
1123 let annotations: Vec<PeakAnnotation> = precursor_final.mz.iter()
1125 .map(|_| PeakAnnotation { contributions: vec![] })
1126 .collect();
1127 let precursor_annotated = MzSpectrumAnnotated::new(
1128 precursor_final.mz.to_vec(),
1129 precursor_final.intensity.to_vec(),
1130 annotations,
1131 );
1132
1133 let precursor_len = precursor_annotated.mz.len();
1134 tims_spectra.push(TimsSpectrumAnnotated::new(
1135 frame_id as i32,
1136 *scan,
1137 rt,
1138 scan_mobility,
1139 ms_type.clone(),
1140 vec![0; precursor_len],
1141 precursor_annotated,
1142 ));
1143 }
1144 }
1145 }
1146 }
1147 }
1148 }
1149
1150 if tims_spectra.is_empty() {
1151 return TimsFrameAnnotated::new(frame_id as i32, rt, ms_type, vec![], vec![], vec![], vec![], vec![], vec![]);
1152 }
1153
1154 TimsFrameAnnotated::from_tims_spectra_annotated(tims_spectra).filter_ranged(
1155 mz_min_val, mz_max_val, 0.0, 10.0, 0, 1000, intensity_min_val, 1e9,
1156 )
1157 }
1158
1159 pub fn get_ion_transmission_matrix(
1160 &self,
1161 peptide_id: u32,
1162 charge: i8,
1163 include_precursor_frames: bool,
1164 ) -> Vec<Vec<f32>> {
1165 let maybe_peptide_sim = self.precursor_frame_builder.peptides.get(&peptide_id);
1166
1167 let mut frame_ids = match maybe_peptide_sim {
1168 Some(maybe_peptide_sim) => maybe_peptide_sim.frame_distribution.occurrence.clone(),
1169 _ => vec![],
1170 };
1171
1172 if !include_precursor_frames {
1173 frame_ids = frame_ids
1174 .iter()
1175 .filter(|frame_id| {
1176 !self
1177 .precursor_frame_builder
1178 .precursor_frame_id_set
1179 .contains(frame_id)
1180 })
1181 .cloned()
1182 .collect();
1183 }
1184
1185 let ion = self
1186 .precursor_frame_builder
1187 .ions
1188 .get(&peptide_id)
1189 .unwrap()
1190 .iter()
1191 .find(|ion| ion.charge == charge)
1192 .unwrap();
1193 let spectrum = ion.simulated_spectrum.clone();
1194 let scan_distribution = &ion.scan_distribution;
1195
1196 let mut transmission_matrix =
1197 vec![vec![0.0; frame_ids.len()]; scan_distribution.occurrence.len()];
1198
1199 for (frame_index, frame) in frame_ids.iter().enumerate() {
1200 for (scan_index, scan) in scan_distribution.occurrence.iter().enumerate() {
1201 if self.transmission_settings.all_transmitted(
1202 *frame as i32,
1203 *scan as i32,
1204 &spectrum.mz,
1205 None,
1206 ) {
1207 transmission_matrix[scan_index][frame_index] = 1.0;
1208 } else if self.transmission_settings.any_transmitted(
1209 *frame as i32,
1210 *scan as i32,
1211 &spectrum.mz,
1212 None,
1213 ) {
1214 let transmitted_spectrum = self.transmission_settings.transmit_spectrum(
1215 *frame as i32,
1216 *scan as i32,
1217 spectrum.clone(),
1218 None,
1219 );
1220 let percentage_transmitted = transmitted_spectrum.intensity.iter().sum::<f64>()
1221 / spectrum.intensity.iter().sum::<f64>();
1222 transmission_matrix[scan_index][frame_index] = percentage_transmitted as f32;
1223 }
1224 }
1225 }
1226
1227 transmission_matrix
1228 }
1229
1230 pub fn count_number_transmissions(&self, peptide_id: u32, charge: i8) -> (usize, usize) {
1231 let frame_ids: Vec<_> = self
1232 .precursor_frame_builder
1233 .peptides
1234 .get(&peptide_id)
1235 .unwrap()
1236 .frame_distribution
1237 .occurrence
1238 .clone()
1239 .iter()
1240 .filter(|frame_id| {
1241 !self
1242 .precursor_frame_builder
1243 .precursor_frame_id_set
1244 .contains(frame_id)
1245 })
1246 .cloned()
1247 .collect();
1248 let ion = self
1249 .precursor_frame_builder
1250 .ions
1251 .get(&peptide_id)
1252 .unwrap()
1253 .iter()
1254 .find(|ion| ion.charge == charge)
1255 .unwrap();
1256 let spectrum = ion.simulated_spectrum.clone();
1257 let scan_distribution = &ion.scan_distribution;
1258
1259 let mut frame_count = 0;
1260 let mut scan_count = 0;
1261
1262 for frame in frame_ids.iter() {
1263 let mut frame_transmitted = false;
1264 for scan in scan_distribution.occurrence.iter() {
1265 if self.transmission_settings.any_transmitted(
1266 *frame as i32,
1267 *scan as i32,
1268 &spectrum.mz,
1269 None,
1270 ) {
1271 frame_transmitted = true;
1272 scan_count += 1;
1273 }
1274 }
1275 if frame_transmitted {
1276 frame_count += 1;
1277 }
1278 }
1279
1280 (frame_count, scan_count)
1281 }
1282
1283 pub fn count_number_transmissions_parallel(
1284 &self,
1285 peptide_ids: Vec<u32>,
1286 charge: Vec<i8>,
1287 num_threads: usize,
1288 ) -> Vec<(usize, usize)> {
1289 let thread_pool = ThreadPoolBuilder::new()
1290 .num_threads(num_threads)
1291 .build()
1292 .unwrap();
1293 let result: Vec<(usize, usize)> = thread_pool.install(|| {
1294 peptide_ids
1295 .par_iter()
1296 .zip(charge.par_iter())
1297 .map(|(peptide_id, charge)| self.count_number_transmissions(*peptide_id, *charge))
1298 .collect()
1299 });
1300
1301 result
1302 }
1303}
1304
1305impl TimsTofCollisionEnergy for TimsTofSyntheticsFrameBuilderDIA {
1306 fn get_collision_energy(&self, frame_id: i32, scan_id: i32) -> f64 {
1307 self.fragmentation_settings
1308 .get_collision_energy(frame_id, scan_id)
1309 }
1310}
1311
1312#[cfg(test)]
1313mod p6a_fragment_kernel_tests {
1314 use super::*;
1315 use std::collections::HashSet;
1316
1317 #[test]
1322 fn fragment_series_spectrum_scaling_branches() {
1323 let series = MzSpectrum::new(vec![200.0, 300.0], vec![10.0, 20.0]);
1324 let empty: HashSet<usize> = HashSet::new();
1325
1326 let none = fragment_series_spectrum(
1328 IsotopeTransmissionMode::None, &series, 0, 2.0, 99.0 , None, &empty, 10,
1329 );
1330 assert_eq!(*none.mz, vec![200.0, 300.0]);
1331 assert_eq!(*none.intensity, vec![20.0, 40.0]);
1332
1333 let ps = fragment_series_spectrum(
1335 IsotopeTransmissionMode::PrecursorScaling, &series, 0, 2.0, 0.5, None, &empty, 10,
1336 );
1337 assert_eq!(*ps.intensity, vec![10.0, 20.0]); let pf = fragment_series_spectrum(
1341 IsotopeTransmissionMode::PerFragment, &series, 0, 2.0, 0.5, None, &empty, 10,
1342 );
1343 assert_eq!(*pf.intensity, vec![20.0, 40.0]);
1344 }
1345
1346 use crate::sim::precursor::TimsTofSyntheticsPrecursorFrameBuilder;
1347
1348 fn one_fragment_frame_dia() -> TimsTofSyntheticsFrameBuilderDIA {
1353 let mut frame_to_abundances = BTreeMap::new();
1354 frame_to_abundances.insert(2u32, (vec![10u32], vec![0.8f32]));
1355 let mut peptide_to_ions = BTreeMap::new();
1356 peptide_to_ions.insert(
1357 10u32,
1358 (
1359 vec![1.0f32],
1360 vec![vec![100u32, 101u32]],
1361 vec![vec![0.3f32, 0.2f32]], vec![2i8],
1363 vec![MzSpectrum::new(vec![600.0], vec![1.0])], ),
1365 );
1366 let mut peptide_to_events = BTreeMap::new();
1367 peptide_to_events.insert(10u32, 1000.0f32);
1368 let mut frame_to_rt = BTreeMap::new();
1369 frame_to_rt.insert(2u32, 30.0f32);
1370
1371 let precursor = TimsTofSyntheticsPrecursorFrameBuilder {
1372 noise_seed: 0,
1373 ions: BTreeMap::new(),
1374 peptides: BTreeMap::new(),
1375 scans: Vec::new(),
1376 frames: Vec::new(),
1377 precursor_frame_id_set: HashSet::new(),
1378 frame_to_abundances,
1379 peptide_to_ions,
1380 frame_to_rt,
1381 scan_to_mobility: BTreeMap::new(),
1382 peptide_to_events,
1383 ion_id_to_peptide_charge: BTreeMap::new(),
1384 };
1385
1386 let mut fragment_ions = BTreeMap::new();
1387 fragment_ions.insert(
1388 (10u32, 2i8, 270i32), (
1390 PeptideProductIonSeriesCollection::new(vec![]),
1391 vec![MzSpectrum::new(vec![200.0], vec![1.0])],
1392 ),
1393 );
1394
1395 TimsTofSyntheticsFrameBuilderDIA {
1396 noise_seed: 0,
1397 path: String::new(),
1398 precursor_frame_builder: precursor,
1399 transmission_settings: TimsTransmissionDIA::new(
1400 vec![], vec![], vec![], vec![], vec![], vec![], vec![], None,
1401 ),
1402 fragmentation_settings: TimsTofCollisionEnergyDIA::new(
1403 vec![], vec![], vec![], vec![], vec![], vec![],
1404 ),
1405 fragment_ions: Some(fragment_ions),
1406 fragment_ions_annotated: None,
1407 isotope_config: IsotopeTransmissionConfig::default(),
1408 fragment_ions_with_transmission: None,
1409 capabilities: InstrumentCapabilities::astral(),
1410 }
1411 }
1412
1413 #[test]
1414 fn astral_ms2_render_full_marginal_and_window_gating() {
1415 let b = one_fragment_frame_dia();
1416
1417 let win = WindowTransmission::new(600.0, 50.0, 15.0);
1421 let RenderedEvent::Scan { ms_level, isolation, spectrum, .. } =
1422 b.render_fragment_scan(2, &win, 27.0, DataMode::Centroid)
1423 else {
1424 panic!("MS2 render must be a Scan");
1425 };
1426 assert_eq!(ms_level, 2);
1427 let iso = isolation.expect("MS2 scan carries an isolation window");
1428 assert!((iso.center_mz - 600.0).abs() < 1e-9 && (iso.width_mz - 50.0).abs() < 1e-9);
1429 assert_eq!(spectrum.mz, vec![200.0]);
1430 assert!((spectrum.intensity[0] - 800.0).abs() < 1e-2, "full marginal, got {}", spectrum.intensity[0]);
1431
1432 let off = WindowTransmission::new(900.0, 50.0, 15.0);
1435 let RenderedEvent::Scan { spectrum: empty, .. } =
1436 b.render_fragment_scan(2, &off, 27.0, DataMode::Centroid)
1437 else {
1438 panic!("must be a Scan");
1439 };
1440 assert!(empty.mz.is_empty(), "no precursor transmitted -> empty MS2 scan");
1441
1442 let RenderedEvent::Scan { spectrum: s2, .. } =
1444 b.render_fragment_scan(2, &win, 27.0, DataMode::Centroid)
1445 else { panic!() };
1446 assert_eq!(s2.mz, spectrum.mz);
1447 assert_eq!(s2.intensity, spectrum.intensity);
1448 }
1449}