Skills › Scientific Agent Skills › pyopenms › accurate_mass_search.py
accurate_mass_search.py
103 lines · skills/pyopenms/scripts/accurate_mass_search.py at 5f2bc03
1#!/usr/bin/env python3
2"""
3Accurate Mass Search (Metabolite Annotation)
5Annotate detected features with putative metabolite identities by matching
6accurate masses against the bundled HMDB databases using AccurateMassSearchEngine.
7Outputs an annotated mzTab and a flat CSV of hits. By default uses the HMDB
8mapping/structure files shipped with pyOpenMS.
10Usage:
11 python accurate_mass_search.py features.featureXML --out-mztab hits.mzTab
12 python accurate_mass_search.py features.featureXML --negative --ppm 5 --csv hits.csv
13"""
15import argparse
16import os
17import sys
19try:
20 import pyopenms as ms
21except ImportError:
22 print("Error: pyopenms not installed. Install with: uv pip install pyopenms")
23 sys.exit(1)
26def main():
27 parser = argparse.ArgumentParser(description="Accurate-mass metabolite annotation.")
28 parser.add_argument("input", help="Input featureXML (or consensusXML)")
29 parser.add_argument("--out-mztab", help="Output mzTab file")
30 parser.add_argument("--csv", help="Output CSV of annotation hits")
31 parser.add_argument("--negative", action="store_true", help="Negative ionization mode")
32 parser.add_argument("--ppm", type=float, default=5.0, help="Mass error tolerance in ppm (default 5)")
33 parser.add_argument("--db-mapping", help="Custom HMDB mapping TSV (default: bundled)")
34 parser.add_argument("--db-struct", help="Custom HMDB structure TSV (default: bundled)")
35 args = parser.parse_args()
37 if not os.path.exists(args.input):
38 print(f"Error: file not found: {args.input}")
39 sys.exit(1)
41 is_consensus = args.input.lower().endswith(".consensusxml")
42 if is_consensus:
43 fmap = ms.ConsensusMap()
44 ms.ConsensusXMLFile().load(args.input, fmap)
45 else:
46 fmap = ms.FeatureMap()
47 ms.FeatureXMLFile().load(args.input, fmap)
48 print(f"Loaded {fmap.size()} features")
50 if args.ppm <= 0:
51 parser.error("--ppm must be positive")
52 for supplied in (args.db_mapping, args.db_struct):
53 if supplied and not os.path.isfile(supplied):
54 parser.error(f"database file not found: {supplied}")
56 engine = ms.AccurateMassSearchEngine()
57 p = engine.getDefaults()
58 p.setValue("mass_error_value", args.ppm)
59 p.setValue("mass_error_unit", "ppm")
60 p.setValue("ionization_mode", "negative" if args.negative else "positive")
61 if args.db_mapping:
62 p.setValue("db:mapping", [os.path.abspath(args.db_mapping)])
63 if args.db_struct:
64 p.setValue("db:struct", [os.path.abspath(args.db_struct)])
65 engine.setParameters(p)
66 engine.init()
68 mztab = ms.MzTab()
69 if is_consensus:
70 engine.run(fmap, mztab)
71 else:
72 engine.run(fmap, mztab)
74 out_mztab = args.out_mztab or os.path.splitext(args.input)[0] + ".mzTab"
75 ms.MzTabFile().store(out_mztab, mztab)
76 print(f"Wrote {out_mztab}")
78 if args.csv:
79 # mzTab small-molecule section -> CSV via pandas
80 import pandas as pd
81 try:
82 sm = mztab.getSmallMoleculeSectionRows()
83 print(f"Small-molecule annotation rows: {len(sm)}")
84 except Exception:
85 pass
86 # Re-read the mzTab text for a robust flat dump of the SML table
87 rows = []
88 with open(out_mztab) as fh:
89 header = None
90 for line in fh:
91 if line.startswith("SMH"):
92 header = line.rstrip("\n").split("\t")
93 elif line.startswith("SML") and header:
94 rows.append(dict(zip(header, line.rstrip("\n").split("\t"))))
95 if rows:
96 pd.DataFrame(rows).to_csv(args.csv, index=False)
97 print(f"Wrote {args.csv} ({len(rows)} annotation rows)")
98 else:
99 print("No small-molecule annotations to write.")
102if __name__ == "__main__":
103 main()