Update mzTab protein parsing

Author: douweschulte

.github/workflows/rust.yml

@@ -33,30 +33,30 @@ jobs:
         cargo build -p rustyms --no-default-features --features imgt,rayon
         cargo build -p rustyms --no-default-features --features imgt,rayon,align
         cargo build -p rustyms --no-default-features --features imgt,align
-  
+
   fmt:
     runs-on: ubuntu-latest
-    
+
     needs: build
     steps:
     - uses: actions/checkout@v4
     - name: Run fmt
       run: cargo fmt --check
-      
+
   clippy:
     runs-on: ubuntu-latest
-    
+
     needs: build
     steps:
     - uses: actions/checkout@v4
     - name: Run clippy
       run: cargo clippy
-      
+
   test:
     runs-on: ubuntu-latest
-    env: 
+    env:
       CARGO_ENCODED_RUSTFLAGS: --cfg=github_action
-    
+
     needs: build
     steps:
     - uses: actions/checkout@v4

CITATION.cff

@@ -17,7 +17,7 @@ cff-version: 1.2.0
 identifiers:
   - description: Main paper (preprint)
     type: doi
-    value: 10.1101/2025.01.18.633732 
+    value: 10.1101/2025.01.18.633732
   - description: Mass alignment algorithm paper
     type: doi
     value: 10.1021/acs.jproteome.4c00188

examples/ion-explorer/figure.r

@@ -9,7 +9,7 @@ library(dplyr)
 #         geom_point() +
 #         scale_size_continuous(range = c(0.25, 3)) +
 #         xlab("Difference to theoretical ion (Da)") +
-#         theme_bw() 
+#         theme_bw()
 
 #     if (file == "fragment_y") {
 #         plot = plot + geom_vline(xintercept=0, linetype="dashed") + geom_vline(xintercept=+18.011-0.984, linetype="dashed") + geom_vline(xintercept=-25.979+18.011, linetype="dashed")
@@ -48,7 +48,7 @@ make_plot = function(fragment, element=NA, mode=NA) {
         geom_point() +
         scale_size_continuous(range = c(0.25, 3)) +
         xlab("Difference to theoretical ion (Da)") +
-        theme_bw() 
+        theme_bw()
 
     if (fragment == "y") {
         plot = plot + geom_vline(xintercept=0, linetype="dashed", alpha=0.25) + geom_vline(xintercept=0.984-18.011, linetype="dashed", alpha=0.25) + geom_vline(xintercept=25.979-18.011, linetype="dashed", alpha=0.25) + xlim(-50, 50)
@@ -74,16 +74,16 @@ compare_plot = function(fragment1, fragment2, element1, element2, mode1, mode2)
         geom_point(data=data2, aes(x=mass, y=-count / max2, size=total_intensity / count), colour="red") +
         scale_size_continuous(range = c(0.125, 3)) +
         xlab("Difference to theoretical ion (Da)") +
-        theme_bw() 
+        theme_bw()
 
     if (fragment1 == "y" && fragment2 == "y") {
-        plot = plot + 
-            geom_vline(xintercept=0, linetype="dashed", alpha=0.25) + 
-            geom_vline(xintercept=0.984-18.011, linetype="dashed", alpha=0.25) + 
-            geom_vline(xintercept=0.984-18.011-15.023, linetype="dotted", alpha=0.25) + 
-            geom_vline(xintercept=0.984-18.011-29.039, linetype="dotted", alpha=0.25) + 
-            geom_vline(xintercept=0.984-18.011-43.088, linetype="dotted", alpha=0.25) + 
-            geom_vline(xintercept=25.979-18.011, linetype="dashed", alpha=0.25) + 
+        plot = plot +
+            geom_vline(xintercept=0, linetype="dashed", alpha=0.25) +
+            geom_vline(xintercept=0.984-18.011, linetype="dashed", alpha=0.25) +
+            geom_vline(xintercept=0.984-18.011-15.023, linetype="dotted", alpha=0.25) +
+            geom_vline(xintercept=0.984-18.011-29.039, linetype="dotted", alpha=0.25) +
+            geom_vline(xintercept=0.984-18.011-43.088, linetype="dotted", alpha=0.25) +
+            geom_vline(xintercept=25.979-18.011, linetype="dashed", alpha=0.25) +
             xlim(-75, 50)
     }
     if (fragment1 == "b" && fragment2 == "b") {

examples/multi-annotator/src/main.rs

@@ -21,17 +21,28 @@ use fragment::FragmentType;
 use itertools::Itertools;
 use mzdata::{
     io::{MZFileReader, SpectrumSource},
-    mzpeaks::{peak_set::PeakSetVec, PeakCollection},
+    mzpeaks::{PeakCollection, peak_set::PeakSetVec},
     mzsignal::PeakPicker,
     spectrum::{SignalContinuity, SpectrumLike},
 };
 use rayon::prelude::*;
 use rustyms::{
     annotation::{
-        model::{parse_custom_models, FragmentationModel, MatchingParameters}, AnnotatableSpectrum, AnnotatedPeak, Score, Scores
-    }, chemistry::MassMode, fragment::{DiagnosticPosition, Fragment}, glycan::MonoSaccharide, identification::{csv::write_csv, BasicCSVData, IdentifiedPeptidoformSource}, quantities::Tolerance, sequence::{
-        parse_custom_modifications, AminoAcid, GnoComposition, SequencePosition, SimpleModificationInner
-    }, spectrum::PeakSpectrum, system::MassOverCharge, *
+        AnnotatableSpectrum, AnnotatedPeak, Score, Scores,
+        model::{FragmentationModel, MatchingParameters, parse_custom_models},
+    },
+    chemistry::MassMode,
+    fragment::{DiagnosticPosition, Fragment},
+    glycan::MonoSaccharide,
+    identification::{BasicCSVData, IdentifiedPeptidoformSource, csv::write_csv},
+    quantities::Tolerance,
+    sequence::{
+        AminoAcid, GnoComposition, SequencePosition, SimpleModificationInner,
+        parse_custom_modifications,
+    },
+    spectrum::PeakSpectrum,
+    system::MassOverCharge,
+    *,
 };
 
 /// The command line interface arguments
@@ -224,16 +235,16 @@ fn main() {
                         // USI spectra are mostly loaded as the binary array maps instead of peaks regardless of the signal continuity level
                         spectrum.peaks = spectrum.arrays.as_ref().map(Into::into);
                     }
-                    if let Some(threshold) = args.tic_noise_threshold && let Some(peaks) = spectrum.peaks.as_mut() { 
+                    if let Some(threshold) = args.tic_noise_threshold && let Some(peaks) = spectrum.peaks.as_mut() {
                         let threshold = peaks.total_ion_current() * threshold;
-                        peaks.peaks.retain(|p: &mzdata::mzpeaks::CentroidPeak| p.intensity > threshold); 
+                        peaks.peaks.retain(|p: &mzdata::mzpeaks::CentroidPeak| p.intensity > threshold);
                     }
-                    if let Some(threshold) = args.basepeak_noise_threshold && let Some(peaks) = spectrum.peaks.as_mut() { 
+                    if let Some(threshold) = args.basepeak_noise_threshold && let Some(peaks) = spectrum.peaks.as_mut() {
                         let threshold = peaks.base_peak().map_or(0.0, |v| v.intensity  * threshold);
-                        peaks.peaks.retain(|p: &mzdata::mzpeaks::CentroidPeak| p.intensity > threshold); 
+                        peaks.peaks.retain(|p: &mzdata::mzpeaks::CentroidPeak| p.intensity > threshold);
                     }
-                    if let Some(threshold) = args.absolute_noise_threshold && let Some(peaks) = spectrum.peaks.as_mut() { 
-                            peaks.peaks.retain(|p: &mzdata::mzpeaks::CentroidPeak| p.intensity > threshold); 
+                    if let Some(threshold) = args.absolute_noise_threshold && let Some(peaks) = spectrum.peaks.as_mut() {
+                            peaks.peaks.retain(|p: &mzdata::mzpeaks::CentroidPeak| p.intensity > threshold);
                     }
                     let fragments = line.sequence.generate_theoretical_fragments(line.z, selected_model);
                     let annotated = spectrum.annotate(

rustyms/README.md

@@ -1,8 +1,8 @@
 # Match those fragments!
 
 Handle mass spectrometry data in Rust. This crate is set up to handle very complex peptides with
-loads of ambiguity and complexity. It pivots around the [`CompoundPeptidoformIon`](crate::sequence::CompoundPeptidoformIon), 
-[`PeptidoformIon`](crate::sequence::PeptidoformIon) and [`Peptidoform`](crate::sequence::Peptidoform) which encode the 
+loads of ambiguity and complexity. It pivots around the [`CompoundPeptidoformIon`](crate::sequence::CompoundPeptidoformIon),
+[`PeptidoformIon`](crate::sequence::PeptidoformIon) and [`Peptidoform`](crate::sequence::Peptidoform) which encode the
 [ProForma](https://github.com/HUPO-PSI/ProForma) specification. Additionally, this crate enables the
 reading of [mgf](spectrum::mgf), doing [spectrum annotation](crate::annotation::AnnotatableSpectrum::annotate)
 (BU/MD/TD), finding [isobaric sequences](crate::prelude::find_isobaric_sets), doing [alignments of peptides](crate::align::align)
@@ -41,7 +41,7 @@ let fragments = peptide.generate_theoretical_fragments(Charge::new::<e>(2), mode
 let parameters = MatchingParameters::default();
 // Annotate the raw data with the theoretical fragments
 let annotated = spectrum[0].annotate(peptide, &fragments, &parameters, MassMode::Monoisotopic);
-// Calculate a peak false discovery rate for this annotation 
+// Calculate a peak false discovery rate for this annotation
 let (fdr, _) = annotated.fdr(&fragments, &parameters, MassMode::Monoisotopic);
 // This is the incorrect sequence for this spectrum so the peak FDR will indicate this
 # dbg!(&fdr, fdr.peaks_sigma(), fdr.peaks_fdr(), fdr.peaks_score());
@@ -60,9 +60,9 @@ let second_peptide = Peptidoform::pro_forma("LVQVET", None)?.into_simple_linear(
 // LVQVET B
 // ─  ╶╴
 let alignment = align::<4, &Peptidoform<SimpleLinear>, &Peptidoform<SimpleLinear>>(
-                  &first_peptide, 
+                  &first_peptide,
                   &second_peptide,
-                  AlignScoring::default(), 
+                  AlignScoring::default(),
                   AlignType::GLOBAL);
 # dbg!(&alignment);
 // Calculate some more statistics on this alignment

rustyms/data/real_examples.mgf

@@ -3353,7 +3353,7 @@ HIGHEST_OBSERVED_M/Z=25743.64907663
 BASE_PEAK_INTENSITY=62787.5703125
 [THERMO_TRAILER_EXTRA]MICRO_SCAN_COUNT=5
 [THERMO_TRAILER_EXTRA]MONOISOTOPIC_M/Z=-1
-[THERMO_TRAILER_EXTRA]HCD_ENERGY=                         
+[THERMO_TRAILER_EXTRA]HCD_ENERGY=
 SEQUENCE=EVQLVESGGGLVQPGGSLRLSC[C:Cystine#XL1]AASGFNIKDTYIHWVRQAPGKGLEWVARIYPTNGYTRYADSVKGRFTISADTSKNTAYLQMNSLRAEDTAVYYC[#XL1]SRWGGDGFYAMDYWGQGTLVTVSSASTKGPSVFPLAPSSKSTSGGTAALGC[C:Cystine#XL2]LVKDYFPEPVTVSWNSGALTSGVHTFPAVLQSSGLYSLSSVVTVPSSSLGTQTYIC[#XL2]NVNHKPSNTKVDKKVEPKSC[C:Cystine#XL3]DKT//DIQMTQSPSSLSASVGDRVTITC[C:Cystine#XL4]RASQDVNTAVAWYQQKPGKAPKLLIYSASFLYSGVPSRFSGSRSGTDFTLTISSLQPEDFATYYC[#XL4]QQHYTTPPTFGQGTKVEIKRTVAAPSVFIFPPSDEQLKSGTASVVC[C:Cystine#XL5]LLNNFYPREAKVQWKVDNALQSGNSQESVTEQDSKDSTYSLSSTLTLSKADYEKHKVYAC[#XL5]EVTHQGLSSPVTKSFNRGEC[#XL3]
 2740.35347663 1054.28
 2745.36057663 897.24

rustyms/src/align/mass_alignment.rs

@@ -214,21 +214,22 @@ pub(super) fn determine_final_score<A, B>(
     path: &[Piece],
     scoring: AlignScoring<'_>,
 ) -> Score {
-    let maximal_score = isize::midpoint(seq_a.sequence()
-        [start_a..start_a + path.iter().map(|p| p.step_a as usize).sum::<usize>()]
-        .iter()
-        .map(|a| {
-            scoring.matrix[a.aminoacid.aminoacid() as usize][a.aminoacid.aminoacid() as usize]
-                as isize
-        })
-        .sum::<isize>(), seq_b.sequence()
-            [start_b..start_b + path.iter().map(|p| p.step_b as usize).sum::<usize>()]
+    let maximal_score = isize::midpoint(
+        seq_a.sequence()[start_a..start_a + path.iter().map(|p| p.step_a as usize).sum::<usize>()]
             .iter()
             .map(|a| {
                 scoring.matrix[a.aminoacid.aminoacid() as usize][a.aminoacid.aminoacid() as usize]
                     as isize
             })
-            .sum::<isize>());
+            .sum::<isize>(),
+        seq_b.sequence()[start_b..start_b + path.iter().map(|p| p.step_b as usize).sum::<usize>()]
+            .iter()
+            .map(|a| {
+                scoring.matrix[a.aminoacid.aminoacid() as usize][a.aminoacid.aminoacid() as usize]
+                    as isize
+            })
+            .sum::<isize>(),
+    );
     let absolute_score = path.last().map(|p| p.score).unwrap_or_default();
     Score {
         absolute: absolute_score,

rustyms/src/chemistry/formula/unimod.rs

@@ -38,7 +38,7 @@ fn parse_unimod_composition_brick(
                 Ok(Brick::Formula(ms.formula_inner(SequencePosition::default(),0)))
             } else {
                 Err(BoxedError::new(BasicKind::Error,
-                    "Invalid Unimod chemical formula", 
+                    "Invalid Unimod chemical formula",
                     "Unknown Unimod composition brick, use an element or one of the unimod shorthands. Eg: 'H(13) C(12) N O(3)'.",
                      Context::line_range(None, text, range)))
             }, |el| Ok(Brick::Element(el)))

rustyms/src/fragment/mzpaf/read.rs

@@ -477,8 +477,8 @@ fn parse_ion<'a>(
                         interpretation.and_then(|i| {
                             i.into_semi_ambiguous()
                                 .ok_or_else(|| BoxedError::new(BasicKind::Error,
-                                    "Invalid mzPAF interpretation", 
-                                    "An mzPAF interpretation should be limited to `base-ProForma compliant` without any labile modifications", 
+                                    "Invalid mzPAF interpretation",
+                                    "An mzPAF interpretation should be limited to `base-ProForma compliant` without any labile modifications",
                                     Context::line_range(None, line, range.start_index()..location)))
                                 .map(|i| (location + 1, Some(i)))
                         })?

rustyms/src/fragment/neutral_loss.rs

@@ -37,8 +37,8 @@ impl ParseJson for NeutralLoss {
                     if let Some(n) = arr[0].as_u64() {
                         Ok((
                             u16::try_from(n).map_err(|_| BoxedError::new(BasicKind::Error,
-                                "Invalid NeutralLoss", 
-                                format!("The {context} amount is too big, the number has to be below {}", u16::MAX), 
+                                "Invalid NeutralLoss",
+                                format!("The {context} amount is too big, the number has to be below {}", u16::MAX),
                                 Context::show(n.to_string())))?,
                             MolecularFormula::from_json_value(arr.pop().unwrap())?,
                         ))