Skills › Scientific Agent Skills › pyopenms › align_link_quantify.py
align_link_quantify.py
147 lines · skills/pyopenms/scripts/align_link_quantify.py at 5f2bc03
1#!/usr/bin/env python3
2"""
3Multi-Sample Align, Link, and Quantify
5End-to-end untargeted quantification across multiple LC-MS samples:
7 1. Detect features in each mzML (FeatureFindingMetabo) [or load .featureXML]
8 2. Align retention times (MapAlignmentAlgorithmPoseClustering)
9 3. Link features into a consensus map (FeatureGroupingAlgorithmQT)
10 4. Export consensusXML + a wide quantification matrix (CSV)
12Accepts either mzML files (features detected here) or pre-computed .featureXML files.
14Usage:
15 python align_link_quantify.py s1.mzML s2.mzML s3.mzML --out-prefix study
16 python align_link_quantify.py *.featureXML --out-prefix study --rt-tol 20 --mz-tol 10
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)
29# Reuse the metabolomics detector if available alongside this script
30sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
31try:
32 from detect_features_metabo import detect_features
33except ImportError:
34 detect_features = None
37def load_or_detect(path, ppm, noise, assume_centroided=False):
38 if path.lower().endswith(".featurexml"):
39 fm = ms.FeatureMap()
40 ms.FeatureXMLFile().load(path, fm)
41 return fm
42 exp = ms.MSExperiment()
43 ms.FileHandler().loadExperiment(path, exp)
44 if detect_features is None:
45 raise RuntimeError("detect_features_metabo.py not importable; pass .featureXML inputs instead.")
46 fm, _ = detect_features(exp, ppm=ppm, noise=noise, assume_centroided=assume_centroided)
47 fm.setPrimaryMSRunPath([os.path.abspath(path)])
48 return fm
51def align(feature_maps, allow_unaligned=False):
52 """Align all maps to the one with the most features (in place)."""
53 ref_idx = max(range(len(feature_maps)), key=lambda i: feature_maps[i].size())
54 aligner = ms.MapAlignmentAlgorithmPoseClustering()
55 aligner.setReference(feature_maps[ref_idx])
56 transformer = ms.MapAlignmentTransformer()
57 for i, fm in enumerate(feature_maps):
58 if i == ref_idx:
59 continue
60 trafo = ms.TransformationDescription()
61 try:
62 aligner.align(fm, trafo)
63 transformer.transformRetentionTimes(fm, trafo, True)
64 except Exception as e:
65 if not allow_unaligned:
66 raise RuntimeError(f"Alignment failed for map {i}; inspect anchors/tolerances before linking") from e
67 print(f" warning: map {i} remains unaligned: {e}")
68 return ref_idx
71def link(feature_maps, filenames, rt_tol, mz_tol, mz_unit):
72 if len(filenames) != len(set(filenames)):
73 raise ValueError("Sample filenames must be unique for consensus export")
74 grouper = ms.FeatureGroupingAlgorithmQT()
75 p = grouper.getParameters()
76 p.setValue("distance_RT:max_difference", float(rt_tol))
77 p.setValue("distance_MZ:max_difference", float(mz_tol))
78 p.setValue("distance_MZ:unit", mz_unit)
79 grouper.setParameters(p)
81 consensus = ms.ConsensusMap()
82 headers = consensus.getColumnHeaders()
83 for i, fm in enumerate(feature_maps):
84 h = ms.ColumnHeader()
85 h.filename = filenames[i]
86 h.size = fm.size()
87 h.unique_id = fm.getUniqueId()
88 headers[i] = h
89 consensus.setColumnHeaders(headers)
90 grouper.group(feature_maps, consensus)
91 consensus.setUniqueIds()
92 return consensus
95def main():
96 parser = argparse.ArgumentParser(description="Align, link, and quantify across samples.")
97 parser.add_argument("inputs", nargs="+", help="mzML and/or featureXML files (2+)")
98 parser.add_argument("--out-prefix", default="consensus", help="Output file prefix")
99 parser.add_argument("--rt-tol", type=float, default=20.0, help="Max RT difference for linking (s)")
100 parser.add_argument("--mz-tol", type=float, default=10.0, help="Max m/z difference for linking")
101 parser.add_argument("--mz-unit", choices=["ppm", "Da"], default="ppm", help="m/z unit (default ppm)")
102 parser.add_argument("--ppm", type=float, default=10.0, help="ppm for feature detection (mzML inputs)")
103 parser.add_argument("--noise", type=float, default=1000.0, help="Noise threshold for detection")
104 parser.add_argument("--assume-centroided", action="store_true", help="Accept unknown MS1 type after verification")
105 parser.add_argument("--allow-unaligned", action="store_true", help="Explicitly continue if RT alignment fails")
106 args = parser.parse_args()
108 if len(args.inputs) < 2:
109 print("Error: provide at least 2 input files.")
110 sys.exit(1)
112 feature_maps = []
113 for path in args.inputs:
114 if not os.path.exists(path):
115 print(f"Error: file not found: {path}")
116 sys.exit(1)
117 print(f"Processing {path}...")
118 fm = load_or_detect(path, args.ppm, args.noise, args.assume_centroided)
119 fm.setUniqueIds()
120 print(f" {fm.size()} features")
121 feature_maps.append(fm)
123 print("Aligning retention times...")
124 ref_idx = align(feature_maps, args.allow_unaligned)
125 print(f" reference map: {args.inputs[ref_idx]}")
127 print("Linking features...")
128 consensus = link(feature_maps, args.inputs, args.rt_tol, args.mz_tol, args.mz_unit)
129 print(f" {consensus.size()} consensus features")
131 out_cxml = f"{args.out_prefix}.consensusXML"
132 ms.ConsensusXMLFile().store(out_cxml, consensus)
133 print(f"Wrote {out_cxml}")
135 # Quantification matrix: rows = consensus features, columns = samples.
136 # Concatenate positionally (index may be non-unique on degenerate maps).
137 import pandas as pd
138 intensity_df = consensus.get_intensity_df().reset_index(drop=True)
139 meta_df = consensus.get_metadata_df().reset_index(drop=True)
140 matrix = pd.concat([meta_df, intensity_df], axis=1)
141 out_csv = f"{args.out_prefix}_quant_matrix.csv"
142 matrix.to_csv(out_csv, index_label="consensus_id")
143 print(f"Wrote {out_csv} ({matrix.shape[0]} features x {intensity_df.shape[1]} samples)")
146if __name__ == "__main__":
147 main()