Skills › Scientific Agent Skills › pyopenms › extract_chromatograms.py
extract_chromatograms.py
109 lines · skills/pyopenms/scripts/extract_chromatograms.py at 5f2bc03
1#!/usr/bin/env python3
2"""
3Extract Ion Chromatograms (XIC/EIC) and TIC/BPC
5Build chromatograms from MS1 data: total ion chromatogram (TIC), base peak
6chromatogram (BPC), and extracted ion chromatograms (XIC) for target m/z values
7within a ppm tolerance. Writes a tidy CSV (rt, trace, intensity) and optionally
8a PNG plot.
10Usage:
11 python extract_chromatograms.py data.mzML --tic --bpc --out chrom.csv
12 python extract_chromatograms.py data.mzML --mz 300.15 450.22 --ppm 10 --out xic.csv --plot xic.png
13"""
15import argparse
16import os
17import sys
19try:
20 import pyopenms as ms
21 import numpy as np
22except ImportError:
23 print("Error: pyopenms/numpy not installed. Install with: uv pip install pyopenms")
24 sys.exit(1)
27def collect_ms1(exp):
28 rts, mz_arrays, int_arrays = [], [], []
29 for spec in exp:
30 if spec.getMSLevel() != 1:
31 continue
32 mz, inten = spec.get_peaks()
33 rts.append(spec.getRT())
34 mz_arrays.append(mz)
35 int_arrays.append(inten)
36 return rts, mz_arrays, int_arrays
39def main():
40 parser = argparse.ArgumentParser(description="Extract TIC/BPC/XIC chromatograms.")
41 parser.add_argument("input", help="Input mzML/mzXML")
42 parser.add_argument("--tic", action="store_true", help="Compute total ion chromatogram")
43 parser.add_argument("--bpc", action="store_true", help="Compute base peak chromatogram")
44 parser.add_argument("--mz", type=float, nargs="+", help="Target m/z values for XIC")
45 parser.add_argument("--ppm", type=float, default=10.0, help="XIC m/z tolerance in ppm (default 10)")
46 parser.add_argument("--out", help="Output CSV (tidy: rt, trace, intensity)")
47 parser.add_argument("--plot", help="Output PNG plot")
48 args = parser.parse_args()
50 if not (args.tic or args.bpc or args.mz):
51 parser.error("request at least one of --tic, --bpc, or --mz")
52 if not os.path.exists(args.input):
53 print(f"Error: file not found: {args.input}")
54 sys.exit(1)
56 exp = ms.MSExperiment()
57 ms.FileHandler().loadExperiment(args.input, exp)
58 rts, mz_arrays, int_arrays = collect_ms1(exp)
59 if not rts:
60 parser.error("no MS1 scans available for TIC/BPC/XIC extraction")
61 if args.ppm <= 0 or (args.mz and any(m <= 0 for m in args.mz)):
62 parser.error("ppm and target m/z values must be positive")
63 rts = np.array(rts)
64 print(f"Collected {len(rts)} MS1 scans")
66 traces = {}
67 if args.tic:
68 traces["TIC"] = np.array([a.sum() if len(a) else 0.0 for a in int_arrays])
69 if args.bpc:
70 traces["BPC"] = np.array([a.max() if len(a) else 0.0 for a in int_arrays])
71 if args.mz:
72 for target in args.mz:
73 tol = target * args.ppm / 1e6
74 trace = np.empty(len(rts))
75 for i, (mz, inten) in enumerate(zip(mz_arrays, int_arrays)):
76 mask = np.abs(mz - target) <= tol
77 trace[i] = inten[mask].sum() if mask.any() else 0.0
78 traces[f"XIC_{target:.4f}"] = trace
79 apex = trace.argmax() if len(trace) else 0
80 print(f" XIC m/z {target:.4f}: apex RT {rts[apex]:.1f}s, max {trace[apex]:.3e}")
82 if args.out:
83 import csv
84 with open(args.out, "w", newline="") as fh:
85 w = csv.writer(fh)
86 w.writerow(["rt", "trace", "intensity"])
87 for name, trace in traces.items():
88 for rt, val in zip(rts, trace):
89 w.writerow([f"{rt:.3f}", name, f"{val:.4f}"])
90 print(f"Wrote {args.out}")
92 if args.plot:
93 import matplotlib
94 matplotlib.use("Agg")
95 import matplotlib.pyplot as plt
96 plt.figure(figsize=(10, 5))
97 for name, trace in traces.items():
98 plt.plot(rts, trace, label=name, linewidth=1)
99 plt.xlabel("Retention time (s)")
100 plt.ylabel("Intensity")
101 plt.legend()
102 plt.title("Chromatograms")
103 plt.tight_layout()
104 plt.savefig(args.plot, dpi=150)
105 print(f"Wrote {args.plot}")
108if __name__ == "__main__":
109 main()