Skills › Scientific Agent Skills › pyopenms › theoretical_spectrum.py
theoretical_spectrum.py
82 lines · skills/pyopenms/scripts/theoretical_spectrum.py at 5f2bc03
1#!/usr/bin/env python3
2"""
3Theoretical Fragment Spectrum Generator
5Generate a theoretical fragment-ion spectrum (b/y, optionally a/c/x/z and
6losses) for a peptide using TheoreticalSpectrumGenerator. Prints annotated
7fragment ions and optionally writes the spectrum to mzML and/or a CSV peak list.
9Usage:
10 python theoretical_spectrum.py DFPIANGER
11 python theoretical_spectrum.py PEPTIDEK --charge 2 --ions b y a --losses
12 python theoretical_spectrum.py DFPIANGER --out-mzml theo.mzML --out-csv peaks.csv
13"""
15import argparse
16import sys
18try:
19 import pyopenms as ms
20except ImportError:
21 print("Error: pyopenms not installed. Install with: uv pip install pyopenms")
22 sys.exit(1)
25def main():
26 parser = argparse.ArgumentParser(description="Generate a theoretical peptide spectrum.")
27 parser.add_argument("peptide", help="Amino-acid sequence (OpenMS mod syntax allowed)")
28 parser.add_argument("--charge", type=int, default=1, help="Max fragment charge (default 1)")
29 parser.add_argument("--ions", nargs="+", default=["b", "y"],
30 choices=["a", "b", "c", "x", "y", "z"], help="Ion series (default b y)")
31 parser.add_argument("--losses", action="store_true", help="Include neutral losses")
32 parser.add_argument("--precursor", action="store_true", help="Include precursor peaks")
33 parser.add_argument("--out-mzml", help="Write spectrum to mzML")
34 parser.add_argument("--out-csv", help="Write annotated peak list to CSV")
35 args = parser.parse_args()
36 if args.charge < 1:
37 parser.error("--charge must be positive")
39 seq = ms.AASequence.fromString(args.peptide)
40 tsg = ms.TheoreticalSpectrumGenerator()
41 p = tsg.getParameters()
42 for ion in ["a", "b", "c", "x", "y", "z"]:
43 p.setValue(f"add_{ion}_ions", "true" if ion in args.ions else "false")
44 p.setValue("add_losses", "true" if args.losses else "false")
45 p.setValue("add_precursor_peaks", "true" if args.precursor else "false")
46 p.setValue("add_metainfo", "true")
47 tsg.setParameters(p)
49 spec = ms.MSSpectrum()
50 spec.setMSLevel(2)
51 spec.setType(ms.SpectrumSettings.SpectrumType.CENTROID)
52 tsg.getSpectrum(spec, seq, 1, args.charge)
53 print(f"Peptide: {seq.toString()} ({len(spec)} fragment peaks)")
55 mz, inten = spec.get_peaks()
56 annotations = spec.getStringDataArrays()
57 names = []
58 for i in range(len(mz)):
59 label = annotations[0][i] if annotations else ""
60 names.append(label.decode("utf-8") if isinstance(label, bytes) else str(label))
61 rows = sorted(zip(mz, inten, names), key=lambda r: r[0])
62 for m, _, name in rows:
63 print(f" {name:12s} m/z {m:.4f}")
65 if args.out_mzml:
66 exp = ms.MSExperiment()
67 exp.addSpectrum(spec)
68 ms.MzMLFile().store(args.out_mzml, exp)
69 print(f"Wrote {args.out_mzml}")
71 if args.out_csv:
72 import csv
73 with open(args.out_csv, "w", newline="") as fh:
74 w = csv.writer(fh)
75 w.writerow(["ion", "mz", "intensity"])
76 for m, it, name in rows:
77 w.writerow([name, f"{m:.5f}", f"{it:.3f}"])
78 print(f"Wrote {args.out_csv}")
81if __name__ == "__main__":
82 main()