Skills › Scientific Agent Skills › pyopenms › plot_ms_data.py
plot_ms_data.py
129 lines · skills/pyopenms/scripts/plot_ms_data.py at 5f2bc03
1#!/usr/bin/env python3
2"""
3Visualize MS Data
5Quick plots for inspecting MS data and results:
6 spectrum single MS spectrum (peak/stick plot) by index or RT
7 tic total ion chromatogram
8 featuremap 2D feature map (RT vs m/z, sized/colored by intensity)
9 map2d 2D heatmap of MS1 signal (RT vs m/z)
11Usage:
12 python plot_ms_data.py spectrum data.mzML --index 0 --out spec.png
13 python plot_ms_data.py tic data.mzML --out tic.png
14 python plot_ms_data.py featuremap features.featureXML --out fmap.png
15 python plot_ms_data.py map2d data.mzML --out map.png
16"""
18import argparse
19import os
20import sys
22try:
23 import pyopenms as ms
24 import numpy as np
25 import matplotlib
26 matplotlib.use("Agg")
27 import matplotlib.pyplot as plt
28except ImportError:
29 print("Error: pyopenms/matplotlib not installed. Install with: uv pip install pyopenms matplotlib")
30 sys.exit(1)
33def load_exp(path):
34 exp = ms.MSExperiment()
35 ms.FileHandler().loadExperiment(path, exp)
36 return exp
39def plot_spectrum(args):
40 exp = load_exp(args.input)
41 if args.rt is not None:
42 spec = min((s for s in exp), key=lambda s: abs(s.getRT() - args.rt))
43 else:
44 spec = exp.getSpectrum(args.index)
45 mz, inten = spec.get_peaks()
46 plt.figure(figsize=(10, 5))
47 plt.vlines(mz, 0, inten, linewidth=0.8)
48 plt.xlabel("m/z")
49 plt.ylabel("Intensity")
50 plt.title(f"Spectrum (MS{spec.getMSLevel()}, RT={spec.getRT():.1f}s, {len(mz)} peaks)")
51 plt.tight_layout()
54def plot_tic(args):
55 exp = load_exp(args.input)
56 rts, tic = [], []
57 for s in exp:
58 if s.getMSLevel() != 1:
59 continue
60 _, inten = s.get_peaks()
61 rts.append(s.getRT())
62 tic.append(float(inten.sum()) if len(inten) else 0.0)
63 plt.figure(figsize=(10, 5))
64 plt.plot(rts, tic, linewidth=1)
65 plt.xlabel("Retention time (s)")
66 plt.ylabel("Total ion current")
67 plt.title("Total Ion Chromatogram")
68 plt.tight_layout()
71def plot_featuremap(args):
72 fm = ms.FeatureMap()
73 ms.FeatureXMLFile().load(args.input, fm)
74 df = fm.get_df()
75 plt.figure(figsize=(10, 6))
76 sizes = 10 + 40 * (df["intensity"] / df["intensity"].max())
77 sc = plt.scatter(df["rt"], df["mz"], s=sizes, c=np.log10(df["intensity"] + 1),
78 cmap="viridis", alpha=0.7)
79 plt.colorbar(sc, label="log10(intensity)")
80 plt.xlabel("Retention time (s)")
81 plt.ylabel("m/z")
82 plt.title(f"Feature map ({fm.size()} features)")
83 plt.tight_layout()
86def plot_map2d(args):
87 exp = load_exp(args.input)
88 rts, mzs, ints = [], [], []
89 for s in exp:
90 if s.getMSLevel() != 1:
91 continue
92 mz, inten = s.get_peaks()
93 rt = s.getRT()
94 for m, it in zip(mz, inten):
95 rts.append(rt); mzs.append(m); ints.append(it)
96 plt.figure(figsize=(10, 6))
97 sc = plt.scatter(rts, mzs, c=np.log10(np.array(ints) + 1), s=2,
98 cmap="inferno", alpha=0.5)
99 plt.colorbar(sc, label="log10(intensity)")
100 plt.xlabel("Retention time (s)")
101 plt.ylabel("m/z")
102 plt.title("MS1 signal map")
103 plt.tight_layout()
106PLOTS = {"spectrum": plot_spectrum, "tic": plot_tic,
107 "featuremap": plot_featuremap, "map2d": plot_map2d}
110def main():
111 parser = argparse.ArgumentParser(description="Visualize MS data.")
112 parser.add_argument("kind", choices=list(PLOTS), help="Plot type")
113 parser.add_argument("input", help="Input file")
114 parser.add_argument("--index", type=int, default=0, help="Spectrum index (spectrum plot)")
115 parser.add_argument("--rt", type=float, help="Select spectrum nearest this RT (spectrum plot)")
116 parser.add_argument("--out", required=True, help="Output image path (PNG/PDF/SVG)")
117 args = parser.parse_args()
119 if not os.path.exists(args.input):
120 print(f"Error: file not found: {args.input}")
121 sys.exit(1)
123 PLOTS[args.kind](args)
124 plt.savefig(args.out, dpi=150)
125 print(f"Wrote {args.out}")
128if __name__ == "__main__":
129 main()