Skills › Scientific Agent Skills › pyopenms › detect_features_metabo.py
detect_features_metabo.py
113 lines · skills/pyopenms/scripts/detect_features_metabo.py at 5f2bc03
1#!/usr/bin/env python3
2"""
3Untargeted Metabolomics Feature Detection
5Run the standard OpenMS small-molecule feature-finding pipeline on centroided
6LC-MS data:
8 MassTraceDetection -> ElutionPeakDetection -> FeatureFindingMetabo
10Outputs a featureXML and (optionally) a CSV table of detected features.
11This is the recommended entry point for untargeted metabolomics preprocessing.
13Usage:
14 python detect_features_metabo.py sample.mzML
15 python detect_features_metabo.py sample.mzML --out-features feats.featureXML --out-csv feats.csv
16 python detect_features_metabo.py sample.mzML --ppm 5 --noise 5000 --charge-low 1 --charge-high 3
17"""
19import argparse
20import os
21import sys
23try:
24 import pyopenms as ms
25except ImportError:
26 print("Error: pyopenms not installed. Install with: uv pip install pyopenms")
27 sys.exit(1)
30from _common import centroided_ms1
33def detect_features(exp, ppm=10.0, noise=1000.0, charge_low=1, charge_high=3,
34 remove_single=True, iso_model="metabolites (5% RMS)", assume_centroided=False):
35 """Run MTD -> EPD -> FFM. Returns a FeatureMap."""
36 exp = centroided_ms1(exp, assume_centroided)
38 # 1. Mass trace detection
39 mtd = ms.MassTraceDetection()
40 p = mtd.getDefaults()
41 p.setValue("mass_error_ppm", float(ppm))
42 p.setValue("noise_threshold_int", float(noise))
43 mtd.setParameters(p)
44 mass_traces = mtd.run(exp, 0)
46 # 2. Elution peak detection
47 epd = ms.ElutionPeakDetection()
48 p = epd.getDefaults()
49 p.setValue("width_filtering", "fixed")
50 epd.setParameters(p)
51 mt_split = epd.detectPeaks(mass_traces)
53 # 3. Feature assembly (isotope/charge resolution)
54 ffm = ms.FeatureFindingMetabo()
55 p = ffm.getDefaults()
56 p.setValue("isotope_filtering_model", iso_model)
57 p.setValue("remove_single_traces", "true" if remove_single else "false")
58 p.setValue("charge_lower_bound", int(charge_low))
59 p.setValue("charge_upper_bound", int(charge_high))
60 p.setValue("report_convex_hulls", "true")
61 ffm.setParameters(p)
62 fm = ms.FeatureMap()
63 ffm.run(mt_split, fm)
64 fm.setUniqueIds()
65 return fm, len(mass_traces)
68def main():
69 parser = argparse.ArgumentParser(description="Untargeted metabolomics feature detection.")
70 parser.add_argument("input", help="Centroided mzML file")
71 parser.add_argument("--out-features", help="Output featureXML (default: <input>.featureXML)")
72 parser.add_argument("--out-csv", help="Optional CSV table of features")
73 parser.add_argument("--ppm", type=float, default=10.0, help="Mass trace m/z tolerance in ppm (default 10)")
74 parser.add_argument("--noise", type=float, default=1000.0, help="Noise intensity threshold (default 1000)")
75 parser.add_argument("--charge-low", type=int, default=1, help="Lower charge bound (default 1)")
76 parser.add_argument("--charge-high", type=int, default=3, help="Upper charge bound (default 3)")
77 parser.add_argument("--keep-singletons", action="store_true",
78 help="Keep single-trace features (default: remove)")
79 parser.add_argument("--iso-model", default="metabolites (5% RMS)",
80 help="Isotope filtering model (e.g. 'none', 'metabolites (5%% RMS)')")
81 parser.add_argument("--assume-centroided", action="store_true", help="Accept unknown MS1 spectrum type after independent verification")
82 args = parser.parse_args()
84 if not os.path.exists(args.input):
85 print(f"Error: file not found: {args.input}")
86 sys.exit(1)
88 exp = ms.MSExperiment()
89 ms.FileHandler().loadExperiment(args.input, exp)
90 print(f"Loaded {exp.getNrSpectra()} spectra from {args.input}")
92 fm, n_traces = detect_features(
93 exp, ppm=args.ppm, noise=args.noise,
94 charge_low=args.charge_low, charge_high=args.charge_high,
95 remove_single=not args.keep_singletons, iso_model=args.iso_model,
96 assume_centroided=args.assume_centroided,
97 )
98 print(f"Mass traces: {n_traces}")
99 fm.setPrimaryMSRunPath([os.path.abspath(args.input)])
100 print(f"Features detected: {fm.size()}")
102 out_features = args.out_features or os.path.splitext(args.input)[0] + ".featureXML"
103 ms.FeatureXMLFile().store(out_features, fm)
104 print(f"Wrote {out_features}")
106 if args.out_csv:
107 df = fm.get_df()
108 df.to_csv(args.out_csv, index=False)
109 print(f"Wrote {args.out_csv} ({len(df)} rows)")
112if __name__ == "__main__":
113 main()