Skills › Scientific Agent Skills › pyopenms › process_spectra.py
process_spectra.py
108 lines · skills/pyopenms/scripts/process_spectra.py at 5f2bc03
1#!/usr/bin/env python3
2"""
3Signal Processing for Spectra
5Apply a configurable chain of signal-processing steps to all (or selected MS-level)
6spectra in an MS file: smoothing, centroiding (peak picking), normalization, and
7intensity/S-N thresholding. Steps run in the order listed below.
9Steps (enable with flags):
10 --smooth gauss|sgolay Smooth profile data
11 --pick Centroid profile data (PeakPickerHiRes)
12 --normalize to_one|to_TIC Normalize intensities
13 --threshold FLOAT Remove peaks below absolute intensity
14 --sn FLOAT Remove peaks below this signal-to-noise ratio
16Usage:
17 python process_spectra.py raw.mzML centroided.mzML --smooth gauss --pick
18 python process_spectra.py data.mzML out.mzML --normalize to_one --threshold 0.01
19 python process_spectra.py data.mzML out.mzML --ms-level 1 --sn 2.0
20"""
22import argparse
23import os
24import sys
26try:
27 import pyopenms as ms
28except ImportError:
29 print("Error: pyopenms not installed. Install with: uv pip install pyopenms")
30 sys.exit(1)
33def main():
34 parser = argparse.ArgumentParser(description="Apply signal processing to spectra.")
35 parser.add_argument("input", help="Input mzML/mzXML file")
36 parser.add_argument("output", help="Output mzML file")
37 parser.add_argument("--ms-level", type=int, help="Process only this MS level (others passed through)")
38 parser.add_argument("--smooth", choices=["gauss", "sgolay"], help="Smoothing filter")
39 parser.add_argument("--gaussian-width", type=float, default=0.2, help="Gaussian width in m/z (Th; default 0.2)")
40 parser.add_argument("--pick", action="store_true", help="Centroid via PeakPickerHiRes")
41 parser.add_argument("--signal-to-noise", type=float, default=0.0,
42 help="PeakPicker S/N threshold (0 = off, default)")
43 parser.add_argument("--normalize", choices=["to_one", "to_TIC"], help="Normalization method")
44 parser.add_argument("--threshold", type=float, help="Remove peaks below this absolute intensity")
45 parser.add_argument("--sn", type=float, help="Remove peaks below this S/N (SignalToNoiseEstimatorMedian)")
46 parser.add_argument("--assume-profile", action="store_true", help="Allow smoothing/picking unknown spectrum type after independent verification")
47 args = parser.parse_args()
49 if not os.path.exists(args.input):
50 print(f"Error: file not found: {args.input}")
51 sys.exit(1)
53 exp = ms.MSExperiment()
54 ms.FileHandler().loadExperiment(args.input, exp)
55 print(f"Loaded {exp.getNrSpectra()} spectra")
57 # Work on copies and explicitly replace the experiment spectra. Iterating
58 # an MSExperiment does not promise writable references to its C++ storage.
59 processed = []
60 for original in exp:
61 spec = ms.MSSpectrum(original)
62 if args.ms_level is not None and spec.getMSLevel() != args.ms_level:
63 processed.append(spec)
64 continue
65 spec.sortByPosition()
66 if args.smooth or args.pick:
67 kind = spec.getType()
68 if kind == ms.SpectrumSettings.SpectrumType.CENTROID:
69 parser.error("smoothing/peak picking requires profile spectra; select a profile MS level")
70 if kind == ms.SpectrumSettings.SpectrumType.UNKNOWN and not args.assume_profile:
71 parser.error("spectrum type unknown; verify profile acquisition before --assume-profile")
72 if args.smooth:
73 smoother = ms.GaussFilter() if args.smooth == "gauss" else ms.SavitzkyGolayFilter()
74 if args.smooth == "gauss":
75 p = smoother.getParameters()
76 p.setValue("gaussian_width", args.gaussian_width)
77 smoother.setParameters(p)
78 smoother.filter(spec)
79 if args.pick:
80 picker = ms.PeakPickerHiRes()
81 p = picker.getParameters()
82 p.setValue("signal_to_noise", args.signal_to_noise)
83 picker.setParameters(p)
84 picked = ms.MSSpectrum()
85 picker.pick(spec, picked)
86 spec = picked
87 if args.normalize:
88 normalizer = ms.Normalizer()
89 p = normalizer.getParameters()
90 p.setValue("method", args.normalize)
91 normalizer.setParameters(p)
92 normalizer.filterSpectrum(spec)
93 if args.sn is not None:
94 estimator = ms.SignalToNoiseEstimatorMedian()
95 estimator.init(spec)
96 spec.select([i for i in range(len(spec)) if estimator.getSignalToNoise(i) >= args.sn])
97 if args.threshold is not None:
98 _, intensity = spec.get_peaks()
99 spec.select((intensity >= args.threshold).nonzero()[0].tolist())
100 processed.append(spec)
101 exp.setSpectra(processed)
103 ms.MzMLFile().store(args.output, exp)
104 print(f"Wrote {args.output}")
107if __name__ == "__main__":
108 main()