Skills › Scientific Agent Skills › pyopenms › inspect_ms_data.py
inspect_ms_data.py
167 lines · skills/pyopenms/scripts/inspect_ms_data.py at 5f2bc03
1#!/usr/bin/env python3
2"""
3Inspect Mass Spectrometry Data
5Load any supported MS file (mzML, mzXML, featureXML, consensusXML, idXML) and
6print a structured summary: spectrum counts by MS level, RT/m/z ranges, TIC,
7chromatograms, precursor info, and instrument metadata. Optionally dump a
8per-spectrum table to CSV.
10Usage:
11 python inspect_ms_data.py data.mzML
12 python inspect_ms_data.py data.mzML --spectra-csv spectra.csv
13 python inspect_ms_data.py features.featureXML
14 python inspect_ms_data.py ids.idXML
15"""
17import argparse
18import os
19import sys
21try:
22 import pyopenms as ms
23except ImportError:
24 print("Error: pyopenms not installed. Install with: uv pip install pyopenms")
25 sys.exit(1)
28def summarize_experiment(path):
29 exp = ms.MSExperiment()
30 ms.FileHandler().loadExperiment(path, exp)
32 n_spec = exp.getNrSpectra()
33 n_chrom = exp.getNrChromatograms()
34 print(f"File: {path}")
35 print(f"Spectra: {n_spec}")
36 print(f"Chromatograms: {n_chrom}")
38 if n_spec:
39 levels = {}
40 rts, total_tic, total_peaks = [], 0.0, 0
41 mz_min, mz_max = float("inf"), float("-inf")
42 n_precursors = 0
43 for spec in exp:
44 lvl = spec.getMSLevel()
45 levels[lvl] = levels.get(lvl, 0) + 1
46 rts.append(spec.getRT())
47 mz, inten = spec.get_peaks()
48 total_peaks += len(mz)
49 if len(mz):
50 total_tic += float(inten.sum())
51 mz_min = min(mz_min, float(mz.min()))
52 mz_max = max(mz_max, float(mz.max()))
53 if spec.getPrecursors():
54 n_precursors += 1
55 print("\nSpectra by MS level:")
56 for lvl in sorted(levels):
57 print(f" MS{lvl}: {levels[lvl]}")
58 if rts:
59 print(f"\nRT range: {min(rts):.1f} - {max(rts):.1f} s "
60 f"({min(rts)/60:.2f} - {max(rts)/60:.2f} min)")
61 if mz_min != float("inf"):
62 print(f"m/z range: {mz_min:.4f} - {mz_max:.4f}")
63 print(f"Total peaks: {total_peaks}")
64 print(f"Total ion current (sum): {total_tic:.3e}")
65 print(f"Spectra with precursors (MSn): {n_precursors}")
67 # Instrument / metadata
68 instr = exp.getInstrument()
69 if instr.getName():
70 print(f"\nInstrument: {instr.getName()}")
71 if instr.getVendor():
72 print(f"Vendor: {instr.getVendor()}")
73 return exp
76def dump_spectra_csv(exp, out_csv):
77 import csv
78 with open(out_csv, "w", newline="") as fh:
79 w = csv.writer(fh)
80 w.writerow(["index", "ms_level", "rt", "n_peaks", "tic",
81 "base_peak_mz", "base_peak_int", "precursor_mz", "precursor_charge"])
82 for i, spec in enumerate(exp):
83 mz, inten = spec.get_peaks()
84 tic = float(inten.sum()) if len(inten) else 0.0
85 bp_mz, bp_int = ("", "")
86 if len(inten):
87 j = int(inten.argmax())
88 bp_mz, bp_int = f"{mz[j]:.4f}", f"{inten[j]:.1f}"
89 prec_mz, prec_z = ("", "")
90 if spec.getPrecursors():
91 p = spec.getPrecursors()[0]
92 prec_mz, prec_z = f"{p.getMZ():.4f}", p.getCharge()
93 w.writerow([i, spec.getMSLevel(), f"{spec.getRT():.3f}", len(mz),
94 f"{tic:.1f}", bp_mz, bp_int, prec_mz, prec_z])
95 print(f"\nWrote per-spectrum table: {out_csv}")
98def summarize_feature_map(path):
99 fm = ms.FeatureMap()
100 ms.FeatureXMLFile().load(path, fm)
101 print(f"File: {path}")
102 print(f"Features: {fm.size()}")
103 if fm.size():
104 df = fm.get_df()
105 print(f"\nColumns: {list(df.columns)}")
106 print(f"RT range: {df['rt'].min():.1f} - {df['rt'].max():.1f} s")
107 print(f"m/z range: {df['mz'].min():.4f} - {df['mz'].max():.4f}")
108 print(f"Intensity: median={df['intensity'].median():.3e} max={df['intensity'].max():.3e}")
109 if "charge" in df:
110 print(f"Charge states: {sorted(df['charge'].unique().tolist())}")
113def summarize_consensus_map(path):
114 cm = ms.ConsensusMap()
115 ms.ConsensusXMLFile().load(path, cm)
116 print(f"File: {path}")
117 print(f"Consensus features: {cm.size()}")
118 headers = cm.getColumnHeaders()
119 print(f"Samples/maps: {len(headers)}")
120 for idx, h in headers.items():
121 print(f" map {idx}: {h.filename} (size={h.size}, label={h.label})")
124def summarize_idxml(path):
125 prot_ids = []
126 pep_ids = ms.PeptideIdentificationList() # required type in pyOpenMS 3.5+
127 ms.IdXMLFile().load(path, prot_ids, pep_ids)
128 print(f"File: {path}")
129 print(f"Protein ID runs: {len(prot_ids)}")
130 total_prot = sum(len(p.getHits()) for p in prot_ids)
131 print(f"Protein hits: {total_prot}")
132 print(f"Peptide identifications (spectra): {len(pep_ids)}")
133 total_pep = sum(len(p.getHits()) for p in pep_ids)
134 print(f"Peptide hits: {total_pep}")
135 if pep_ids:
136 scored = [p for p in pep_ids if p.getHits()]
137 if scored:
138 st = scored[0].getScoreType()
139 print(f"Score type: {st}")
140 print(f"Higher score better: {scored[0].isHigherScoreBetter()}")
143def main():
144 parser = argparse.ArgumentParser(description="Inspect a mass spectrometry data file.")
145 parser.add_argument("input", help="Input file (mzML, mzXML, featureXML, consensusXML, idXML)")
146 parser.add_argument("--spectra-csv", help="Write a per-spectrum summary table to this CSV path")
147 args = parser.parse_args()
149 if not os.path.exists(args.input):
150 print(f"Error: file not found: {args.input}")
151 sys.exit(1)
153 ext = args.input.lower()
154 if ext.endswith(".featurexml"):
155 summarize_feature_map(args.input)
156 elif ext.endswith(".consensusxml"):
157 summarize_consensus_map(args.input)
158 elif ext.endswith(".idxml"):
159 summarize_idxml(args.input)
160 else:
161 exp = summarize_experiment(args.input)
162 if args.spectra_csv:
163 dump_spectra_csv(exp, args.spectra_csv)
166if __name__ == "__main__":
167 main()