Skills › Scientific Agent Skills › pyopenms › detect_adducts.py
detect_adducts.py
92 lines · skills/pyopenms/scripts/detect_adducts.py at 5f2bc03
1#!/usr/bin/env python3
2"""
3Adduct Detection / Feature Decharging
5Group features that are different ionization forms (adducts/charge variants) of
6the same neutral compound using MetaboliteFeatureDeconvolution. Annotates each
7feature with its inferred adduct and neutral mass, and writes the decharged
8feature map plus a consensus map of grouped adduct families.
10Adducts use OpenMS deconvolution syntax 'Elements:Charge:Probability', where
11the charge is given as '+'/'-' signs (e.g. 'Ca:++:0.5' is +2) and losses are
12written like 'H-2O-1'. Example: 'H:+:0.4'. This differs from the bracket
13notation (e.g. [M+H]+) used elsewhere.
15Usage:
16 python detect_adducts.py features.featureXML --out-features decharged.featureXML
17 python detect_adducts.py features.featureXML --negative
18 python detect_adducts.py features.featureXML --adducts "H:+:0.6,Na:+:0.2,K:+:0.2"
19"""
21import argparse
22import os
23import sys
25try:
26 import pyopenms as ms
27except ImportError:
28 print("Error: pyopenms not installed. Install with: uv pip install pyopenms")
29 sys.exit(1)
31# 'Elements:Charge:Probability' (OpenMS deconvolution syntax)
32DEFAULT_POS = "H:+:0.4,Na:+:0.25,NH4:+:0.25,K:+:0.1,H-2O-1:0:0.05"
33DEFAULT_NEG = "H-1:-:0.6,Cl:-:0.2,H-3O-1:-:0.1,CH2O2H-1:-:0.1"
36def main():
37 parser = argparse.ArgumentParser(description="Detect adducts / decharge a feature map.")
38 parser.add_argument("input", help="Input featureXML")
39 parser.add_argument("--out-features", help="Output decharged featureXML (default: <input>_decharged.featureXML)")
40 parser.add_argument("--out-consensus", help="Optional consensusXML of adduct groups")
41 parser.add_argument("--negative", action="store_true", help="Negative ionization mode")
42 parser.add_argument("--adducts", help="Comma-separated potential adducts (overrides defaults)")
43 parser.add_argument("--charge-min", type=int, default=1, help="Minimum charge magnitude (positive integer)")
44 parser.add_argument("--charge-max", type=int, default=1, help="Maximum charge magnitude (positive integer)")
45 parser.add_argument("--mass-max-diff", type=float, default=0.05, help="Max mass difference (Da)")
46 parser.add_argument("--rt-max-diff", type=float, default=10.0, help="Max RT difference (s)")
47 args = parser.parse_args()
48 if not 1 <= args.charge_min <= args.charge_max:
49 parser.error("charge magnitudes must satisfy 1 <= min <= max")
51 if not os.path.exists(args.input):
52 print(f"Error: file not found: {args.input}")
53 sys.exit(1)
55 fm = ms.FeatureMap()
56 ms.FeatureXMLFile().load(args.input, fm)
57 print(f"Loaded {fm.size()} features")
59 adducts = args.adducts
60 if adducts is None:
61 adducts = DEFAULT_NEG if args.negative else DEFAULT_POS
63 mfd = ms.MetaboliteFeatureDeconvolution()
64 p = mfd.getDefaults()
65 p.setValue("potential_adducts", [a.strip().encode() for a in adducts.split(",")])
66 p.setValue("charge_min", -args.charge_max if args.negative else args.charge_min)
67 p.setValue("charge_max", -args.charge_min if args.negative else args.charge_max)
68 p.setValue("mass_max_diff", args.mass_max_diff)
69 p.setValue("retention_max_diff", args.rt_max_diff)
70 p.setValue("negative_mode", "true" if args.negative else "false")
71 mfd.setParameters(p)
73 fm_out = ms.FeatureMap()
74 groups = ms.ConsensusMap()
75 edges = ms.ConsensusMap()
76 mfd.compute(fm, fm_out, groups, edges)
78 annotated = sum(1 for f in fm_out if f.metaValueExists("dc_charge_adducts"))
79 print(f"Decharged feature map: {fm_out.size()} features ({annotated} adduct-annotated)")
80 print(f"Adduct groups: {groups.size()}")
82 out_features = args.out_features or os.path.splitext(args.input)[0] + "_decharged.featureXML"
83 ms.FeatureXMLFile().store(out_features, fm_out)
84 print(f"Wrote {out_features}")
86 if args.out_consensus:
87 ms.ConsensusXMLFile().store(args.out_consensus, groups)
88 print(f"Wrote {args.out_consensus}")
91if __name__ == "__main__":
92 main()