Pyopenms

PyOpenMS

Overview

PyOpenMS provides Python bindings to the OpenMS library for computational mass spectrometry, enabling analysis of proteomics and metabolomics data. Use it to read/write MS file formats, process raw spectra, detect and quantify features, identify peptides and proteins, and run end-to-end LC-MS/MS pipelines.

This skill ships ready-to-run scripts in `scripts/` covering the most common high-level workflows. Prefer running a script over writing new code—each is a parameterized CLI tool that handles loading, processing, and export. Drop into the Python API (and the references/) only when no script fits.

Installation

uv pip install pyopenms

Verify (note: __version__ works, but the bundled binary prints a one-line memory-status notice on import that is harmless):

import pyopenms as ms
print(ms.__version__)  # 3.5.0

Scripts (start here)

Run with python scripts/<name>.py --help for full options. All accept standard MS file formats and write featureXML/consensusXML/CSV/mzTab/PNG as appropriate.

Inspect & convert

Feature detection & quantification

Annotation

Identification

Chemistry

Targeted & visualization

Common script recipes

# Inspect a file
python scripts/inspect_ms_data.py sample.mzML --spectra-csv spectra.csv

# Untargeted metabolomics: features for one sample
python scripts/detect_features_metabo.py sample.mzML --out-csv features.csv

# Full multi-sample quantification study
python scripts/align_link_quantify.py s1.mzML s2.mzML s3.mzML --out-prefix study
python scripts/consensus_to_matrix.py study.consensusXML --out quant.csv --normalize median

# Peptide chemistry
python scripts/mass_calculator.py --peptide "PEPTIDEM(Oxidation)K" --charges 1 2 3 --isotopes 5
python scripts/digest_protein.py proteins.fasta --enzyme Trypsin --missed 2 --out peptides.csv

# Identification post-processing
python scripts/process_identifications.py search.idXML --fasta db.fasta --fdr 0.01 --out filtered.idXML --csv hits.csv

Key 3.5.0 API notes

These changed from older OpenMS releases—older tutorials and code will break:

• Feature finding: FeatureFinder("centroided") was removed. Use

FeatureFinderAlgorithmPicked (proteomics/centroided) or the MassTraceDetection → ElutionPeakDetection → FeatureFindingMetabo pipeline (metabolomics). See detect_features_*.py.

• idXML I/O: IdXMLFile().load/store require a ms.PeptideIdentificationList()

for peptide IDs (a plain Python list raises "can not handle type"). Protein IDs remain a plain list.

• Adduct decharging: the class is MetaboliteFeatureDeconvolution, and adducts

use Elements:Charge:Probability syntax (e.g. H:+:0.4, H-2O-1:0:0.05)—not bracket notation like [M+H]+.

• DataFrame columns: FeatureMap.get_df() uses lowercase rt/mz (not RT).

ConsensusMap provides get_intensity_df() and get_metadata_df().

• Bundled data caveat: the pip wheel ships HMDBMappingFile.tsv but not

HMDB2StructMapping.tsv; accurate_mass_search.py detects this and explains how to supply it.

Core data structures

• MSExperiment – collection of spectra and chromatograms
• MSSpectrum / MSChromatogram – a single spectrum / chromatographic trace
• Feature / FeatureMap – a detected LC-MS peak / collection of features
• ConsensusMap – features linked across samples (the quant table)
• PeptideIdentification / ProteinIdentification – search results
• AASequence / EmpiricalFormula – sequence and formula chemistry

For details: see references/data_structures.md.

Parameter management

Most algorithms expose an OpenMS Param object:

algo = ms.FeatureFindingMetabo()
p = algo.getDefaults()
for key in p.keys():
    print(key.decode(), "=", p.getValue(key), "|", p.getDescription(key))
p.setValue("charge_lower_bound", 1)
algo.setParameters(p)

Export to pandas

fm = ms.FeatureMap(); ms.FeatureXMLFile().load("features.featureXML", fm)
df = fm.get_df()             # columns include lowercase rt, mz, intensity, charge, quality

cm = ms.ConsensusMap(); ms.ConsensusXMLFile().load("study.consensusXML", cm)
intensities = cm.get_intensity_df()   # features x samples
metadata = cm.get_metadata_df()       # rt, mz, charge, quality, ...

Integration with other tools

Pandas (DataFrames), NumPy (peak arrays), scikit-learn (ML), Matplotlib/Seaborn (plots), and downstream tools via export: GNPS (FBMN), SIRIUS, and mzTab.

Resources

• Official docs (3.5.0): https://pyopenms.readthedocs.io/en/release-3.5.0/
• OpenMS: https://www.openms.org
• GitHub: https://github.com/OpenMS/OpenMS

References

• references/file_io.md – file format handling
• references/signal_processing.md – signal processing algorithms
• references/feature_detection.md – feature detection and linking
• references/identification.md – peptide and protein identification
• references/metabolomics.md – metabolomics-specific workflows
• references/data_structures.md – core objects and data structures

Core Data Structures

Overview

PyOpenMS uses C++ objects with Python bindings. Understanding these core data structures is essential for effective data manipulation.

Spectrum and Experiment Objects

MSExperiment

Container for complete LC-MS experiment data (spectra and chromatograms).

import pyopenms as ms

# Create experiment
exp = ms.MSExperiment()

# Load from file
ms.MzMLFile().load("data.mzML", exp)

# Access properties
print(f"Number of spectra: {exp.getNrSpectra()}")
print(f"Number of chromatograms: {exp.getNrChromatograms()}")

# Get RT range
rts = [spec.getRT() for spec in exp]
print(f"RT range: {min(rts):.1f} - {max(rts):.1f} seconds")

# Access individual spectrum
spec = exp.getSpectrum(0)

# Iterate through spectra
for spec in exp:
    if spec.getMSLevel() == 2:
        print(f"MS2 spectrum at RT {spec.getRT():.2f}")

# Get metadata
exp_settings = exp.getExperimentalSettings()
instrument = exp_settings.getInstrument()
print(f"Instrument: {instrument.getName()}")

MSSpectrum

Individual mass spectrum with m/z and intensity arrays.

# Create empty spectrum
spec = ms.MSSpectrum()

# Get from experiment
exp = ms.MSExperiment()
ms.MzMLFile().load("data.mzML", exp)
spec = exp.getSpectrum(0)

# Basic properties
print(f"MS level: {spec.getMSLevel()}")
print(f"Retention time: {spec.getRT():.2f} seconds")
print(f"Number of peaks: {spec.size()}")

# Get peak data as numpy arrays
mz, intensity = spec.get_peaks()
print(f"m/z range: {mz.min():.2f} - {mz.max():.2f}")
print(f"Max intensity: {intensity.max():.0f}")

# Access individual peaks
for i in range(min(5, spec.size())):  # First 5 peaks
    print(f"Peak {i}: m/z={mz[i]:.4f}, intensity={intensity[i]:.0f}")

# Precursor information (for MS2)
if spec.getMSLevel() == 2:
    precursors = spec.getPrecursors()
    if precursors:
        precursor = precursors[0]
        print(f"Precursor m/z: {precursor.getMZ():.4f}")
        print(f"Precursor charge: {precursor.getCharge()}")
        print(f"Precursor intensity: {precursor.getIntensity():.0f}")

# Set peak data
new_mz = [100.0, 200.0, 300.0]
new_intensity = [1000.0, 2000.0, 1500.0]
spec.set_peaks((new_mz, new_intensity))

MSChromatogram

Chromatographic trace (TIC, XIC, or SRM transition).

# Access chromatogram from experiment
for chrom in exp.getChromatograms():
    print(f"Chromatogram ID: {chrom.getNativeID()}")

    # Get data
    rt, intensity = chrom.get_peaks()

    print(f"  RT points: {len(rt)}")
    print(f"  Max intensity: {intensity.max():.0f}")

    # Precursor info (for XIC)
    precursor = chrom.getPrecursor()
    print(f"  Precursor m/z: {precursor.getMZ():.4f}")

Feature Objects

Feature

Detected chromatographic peak with 2D spatial extent (RT-m/z).

# Load features
feature_map = ms.FeatureMap()
ms.FeatureXMLFile().load("features.featureXML", feature_map)

# Access individual feature
feature = feature_map[0]

# Core properties
print(f"m/z: {feature.getMZ():.4f}")
print(f"RT: {feature.getRT():.2f} seconds")
print(f"Intensity: {feature.getIntensity():.0f}")
print(f"Charge: {feature.getCharge()}")

# Quality metrics
print(f"Overall quality: {feature.getOverallQuality():.3f}")
print(f"Width (RT): {feature.getWidth():.2f}")

# Convex hull (spatial extent)
hull = feature.getConvexHull()
print(f"Hull points: {hull.getHullPoints().size()}")

# Bounding box
bbox = hull.getBoundingBox()
print(f"RT range: {bbox.minPosition()[0]:.2f} - {bbox.maxPosition()[0]:.2f}")
print(f"m/z range: {bbox.minPosition()[1]:.4f} - {bbox.maxPosition()[1]:.4f}")

# Subordinate features (isotopes)
subordinates = feature.getSubordinates()
if subordinates:
    print(f"Isotopic features: {len(subordinates)}")
    for sub in subordinates:
        print(f"  m/z: {sub.getMZ():.4f}, intensity: {sub.getIntensity():.0f}")

# Metadata values
if feature.metaValueExists("label"):
    label = feature.getMetaValue("label")
    print(f"Label: {label}")

FeatureMap

Collection of features from a single LC-MS run.

# Create feature map
feature_map = ms.FeatureMap()

# Load from file
ms.FeatureXMLFile().load("features.featureXML", feature_map)

# Access properties
print(f"Number of features: {feature_map.size()}")

# Get unique features
print(f"Unique features: {feature_map.getUniqueId()}")

# Metadata
primary_path = feature_map.getPrimaryMSRunPath()
if primary_path:
    print(f"Source file: {primary_path[0].decode()}")

# Iterate through features
for feature in feature_map:
    print(f"Feature: m/z={feature.getMZ():.4f}, RT={feature.getRT():.2f}")

# Add new feature
new_feature = ms.Feature()
new_feature.setMZ(500.0)
new_feature.setRT(300.0)
new_feature.setIntensity(10000.0)
feature_map.push_back(new_feature)

# Sort features
feature_map.sortByRT()  # or sortByMZ(), sortByIntensity()

# Export to pandas
df = feature_map.get_df()
print(df.head())

ConsensusFeature

Feature linked across multiple samples.

# Load consensus map
consensus_map = ms.ConsensusMap()
ms.ConsensusXMLFile().load("consensus.consensusXML", consensus_map)

# Access consensus feature
cons_feature = consensus_map[0]

# Consensus properties
print(f"Consensus m/z: {cons_feature.getMZ():.4f}")
print(f"Consensus RT: {cons_feature.getRT():.2f}")
print(f"Consensus intensity: {cons_feature.getIntensity():.0f}")

# Get feature handles (individual map features)
feature_list = cons_feature.getFeatureList()
print(f"Present in {len(feature_list)} maps")

for handle in feature_list:
    map_idx = handle.getMapIndex()
    intensity = handle.getIntensity()
    mz = handle.getMZ()
    rt = handle.getRT()

    print(f"  Map {map_idx}: m/z={mz:.4f}, RT={rt:.2f}, intensity={intensity:.0f}")

# Get unique ID in originating map
for handle in feature_list:
    unique_id = handle.getUniqueId()
    print(f"Unique ID: {unique_id}")

ConsensusMap

Collection of consensus features across samples.

# Create consensus map
consensus_map = ms.ConsensusMap()

# Load from file
ms.ConsensusXMLFile().load("consensus.consensusXML", consensus_map)

# Access properties
print(f"Consensus features: {consensus_map.size()}")

# Column headers (file descriptions)
headers = consensus_map.getColumnHeaders()
print(f"Number of files: {len(headers)}")

for map_idx, description in headers.items():
    print(f"Map {map_idx}:")
    print(f"  Filename: {description.filename}")
    print(f"  Label: {description.label}")
    print(f"  Size: {description.size}")

# Iterate through consensus features
for cons_feature in consensus_map:
    print(f"Consensus feature: m/z={cons_feature.getMZ():.4f}")

# Export to DataFrame
df = consensus_map.get_df()

Identification Objects

PeptideIdentification

Identification results for a single spectrum.

# Load identifications
protein_ids = []                              # protein IDs: plain list
# pyOpenMS 3.5+: peptide IDs must be a PeptideIdentificationList, not a plain list
peptide_ids = ms.PeptideIdentificationList()
ms.IdXMLFile().load("identifications.idXML", protein_ids, peptide_ids)

# Access peptide identification
peptide_id = peptide_ids[0]

# Spectrum metadata
print(f"RT: {peptide_id.getRT():.2f}")
print(f"m/z: {peptide_id.getMZ():.4f}")

# Identification metadata
print(f"Identifier: {peptide_id.getIdentifier()}")
print(f"Score type: {peptide_id.getScoreType()}")
print(f"Higher score better: {peptide_id.isHigherScoreBetter()}")

# Get peptide hits
hits = peptide_id.getHits()
print(f"Number of hits: {len(hits)}")

for hit in hits:
    print(f"  Sequence: {hit.getSequence().toString()}")
    print(f"  Score: {hit.getScore()}")
    print(f"  Charge: {hit.getCharge()}")

PeptideHit

Individual peptide match to a spectrum.

# Access hit
hit = peptide_id.getHits()[0]

# Sequence information
sequence = hit.getSequence()
print(f"Sequence: {sequence.toString()}")
print(f"Mass: {sequence.getMonoWeight():.4f}")

# Score and rank
print(f"Score: {hit.getScore()}")
print(f"Rank: {hit.getRank()}")

# Charge state
print(f"Charge: {hit.getCharge()}")

# Protein accessions
accessions = hit.extractProteinAccessionsSet()
for acc in accessions:
    print(f"Protein: {acc.decode()}")

# Meta values (additional scores, errors)
if hit.metaValueExists("MS:1002252"):  # mass error
    mass_error = hit.getMetaValue("MS:1002252")
    print(f"Mass error: {mass_error:.4f} ppm")

ProteinIdentification

Protein-level identification information.

# Access protein identification
protein_id = protein_ids[0]

# Search engine info
print(f"Search engine: {protein_id.getSearchEngine()}")
print(f"Search engine version: {protein_id.getSearchEngineVersion()}")

# Search parameters
search_params = protein_id.getSearchParameters()
print(f"Database: {search_params.db}")
print(f"Enzyme: {search_params.digestion_enzyme.getName()}")
print(f"Missed cleavages: {search_params.missed_cleavages}")
print(f"Precursor tolerance: {search_params.precursor_mass_tolerance}")

# Protein hits
hits = protein_id.getHits()
for hit in hits:
    print(f"Accession: {hit.getAccession()}")
    print(f"Score: {hit.getScore()}")
    print(f"Coverage: {hit.getCoverage():.1f}%")

ProteinHit

Individual protein identification.

# Access protein hit
protein_hit = protein_id.getHits()[0]

# Protein information
print(f"Accession: {protein_hit.getAccession()}")
print(f"Description: {protein_hit.getDescription()}")
print(f"Sequence: {protein_hit.getSequence()}")

# Scoring
print(f"Score: {protein_hit.getScore()}")
print(f"Coverage: {protein_hit.getCoverage():.1f}%")

# Rank
print(f"Rank: {protein_hit.getRank()}")

Sequence Objects

AASequence

Amino acid sequence with modifications.

# Create sequence from string
seq = ms.AASequence.fromString("PEPTIDE")

# Basic properties
print(f"Sequence: {seq.toString()}")
print(f"Length: {seq.size()}")
print(f"Monoisotopic mass: {seq.getMonoWeight():.4f}")
print(f"Average mass: {seq.getAverageWeight():.4f}")

# Individual residues
for i in range(seq.size()):
    residue = seq.getResidue(i)
    print(f"Position {i}: {residue.getOneLetterCode()}")
    print(f"  Mass: {residue.getMonoWeight():.4f}")
    print(f"  Formula: {residue.getFormula().toString()}")

# Modified sequence
mod_seq = ms.AASequence.fromString("PEPTIDEM(Oxidation)K")
print(f"Modified: {mod_seq.isModified()}")

# Check modifications
for i in range(mod_seq.size()):
    residue = mod_seq.getResidue(i)
    if residue.isModified():
        print(f"Modification at {i}: {residue.getModificationName()}")

# N-terminal and C-terminal modifications
term_mod_seq = ms.AASequence.fromString("(Acetyl)PEPTIDE(Amidated)")

EmpiricalFormula

Molecular formula representation.

# Create formula
formula = ms.EmpiricalFormula("C6H12O6")  # Glucose

# Properties
print(f"Formula: {formula.toString()}")
print(f"Monoisotopic mass: {formula.getMonoWeight():.4f}")
print(f"Average mass: {formula.getAverageWeight():.4f}")

# Element composition
print(f"Carbon atoms: {formula.getNumberOf(b'C')}")
print(f"Hydrogen atoms: {formula.getNumberOf(b'H')}")
print(f"Oxygen atoms: {formula.getNumberOf(b'O')}")

# Arithmetic operations
formula2 = ms.EmpiricalFormula("H2O")
combined = formula + formula2  # Add water
print(f"Combined: {combined.toString()}")

Parameter Objects

Param

Generic parameter container used by algorithms.

# Get algorithm parameters
algo = ms.GaussFilter()
params = algo.getParameters()

# List all parameters
for key in params.keys():
    value = params.getValue(key)
    print(f"{key}: {value}")

# Get specific parameter
gaussian_width = params.getValue("gaussian_width")
print(f"Gaussian width: {gaussian_width}")

# Set parameter
params.setValue("gaussian_width", 0.2)

# Apply modified parameters
algo.setParameters(params)

# Copy parameters
params_copy = ms.Param(params)

Best Practices

Memory Management

# For large files, use indexed access instead of full loading
indexed_mzml = ms.IndexedMzMLFileLoader()
indexed_mzml.load("large_file.mzML")

# Access specific spectrum without loading entire file
spec = indexed_mzml.getSpectrumById(100)

Type Conversion

# Convert peak arrays to numpy
import numpy as np

mz, intensity = spec.get_peaks()
# These are already numpy arrays

# Can perform numpy operations
filtered_mz = mz[intensity > 1000]

Object Copying

# Create deep copy
exp_copy = ms.MSExperiment(exp)

# Modifications to copy don't affect original

Feature Detection and Linking

Overview

Feature detection identifies persistent signals (chromatographic peaks) in LC-MS data. Feature linking combines features across multiple samples for quantitative comparison.

Ready-to-run scripts: The skill ships CLIs that implement these workflows end to end: scripts/detect_features_metabo.py (metabolomics), scripts/detect_features_centroided.py (proteomics/centroided), scripts/align_link_quantify.py (alignment + linking + quant matrix), and scripts/detect_adducts.py (adduct grouping). Use them directly, or adapt the code below.

API note (pyOpenMS 3.5.0): The old FeatureFinder class and its run("centroided", ...) API were removed. Metabolomics now uses the MassTraceDetection -> ElutionPeakDetection -> FeatureFindingMetabo pipeline, and centroided/proteomics data uses FeatureFinderAlgorithmPicked. The patterns below reflect the current API.

Feature Detection Basics

A feature represents a chromatographic peak characterized by:

• m/z value (mass-to-charge ratio)
• Retention time (RT)
• Intensity
• Quality score
• Convex hull (spatial extent in RT-m/z space)

Feature Finding

Feature Finding for Metabolomics (FeatureFindingMetabo)

For small molecules, run the three-stage pipeline that replaced the removed FeatureFinder: detect mass traces, split them into elution peaks, then assemble isotope-grouped features.

import pyopenms as ms

# Load centroided data
exp = ms.MSExperiment()
ms.MzMLFile().load("centroided.mzML", exp)
exp.sortSpectra(True)

# Stage 1: mass trace detection
mtd = ms.MassTraceDetection()
p = mtd.getDefaults()
p.setValue("mass_error_ppm", 10.0)
p.setValue("noise_threshold_int", 1000.0)
mtd.setParameters(p)
mass_traces = []
mtd.run(exp, mass_traces, 0)

# Stage 2: elution peak detection
epd = ms.ElutionPeakDetection()
p = epd.getDefaults()
p.setValue("width_filtering", "fixed")
epd.setParameters(p)
mt_split = []
epd.detectPeaks(mass_traces, mt_split)

# Stage 3: feature assembly with isotope grouping
ffm = ms.FeatureFindingMetabo()
p = ffm.getDefaults()
p.setValue("isotope_filtering_model", "metabolites (5% RMS)")  # or "none"
p.setValue("remove_single_traces", "true")
p.setValue("charge_lower_bound", 1)
p.setValue("charge_upper_bound", 3)
ffm.setParameters(p)
features = ms.FeatureMap()
chrom_out = []
ffm.run(mt_split, features, chrom_out)

print(f"Detected {features.size()} features")

# Save features
ms.FeatureXMLFile().store("features.featureXML", features)

Feature Finding for Proteomics (FeatureFinderAlgorithmPicked)

For centroided peptide data, use FeatureFinderAlgorithmPicked (replaces the removed FeatureFinder "centroided" workflow):

exp = ms.MSExperiment()
ms.MzMLFile().load("centroided.mzML", exp)
exp.sortSpectra(True)
exp.updateRanges()

ff = ms.FeatureFinderAlgorithmPicked()
params = ff.getDefaults()
params.setValue("isotopic_pattern:charge_low", 1)
params.setValue("isotopic_pattern:charge_high", 4)

features = ms.FeatureMap()
seeds = ms.FeatureMap()
# signature: run(input_map, output, param, seeds)
ff.run(exp, features, params, seeds)

print(f"Detected {features.size()} features")
ms.FeatureXMLFile().store("features.featureXML", features)

Accessing Feature Data

Iterate Through Features

# Load features
feature_map = ms.FeatureMap()
ms.FeatureXMLFile().load("features.featureXML", feature_map)

# Access individual features
for feature in feature_map:
    print(f"m/z: {feature.getMZ():.4f}")
    print(f"RT: {feature.getRT():.2f}")
    print(f"Intensity: {feature.getIntensity():.0f}")
    print(f"Charge: {feature.getCharge()}")
    print(f"Quality: {feature.getOverallQuality():.3f}")
    print(f"Width (RT): {feature.getWidth():.2f}")

    # Get convex hull
    hull = feature.getConvexHull()
    print(f"Hull points: {hull.getHullPoints().size()}")

Feature Subordinates (Isotope Pattern)

# Access isotopic pattern
for feature in feature_map:
    # Get subordinate features (isotopes)
    subordinates = feature.getSubordinates()

    if subordinates:
        print(f"Main feature m/z: {feature.getMZ():.4f}")
        for sub in subordinates:
            print(f"  Isotope m/z: {sub.getMZ():.4f}")
            print(f"  Isotope intensity: {sub.getIntensity():.0f}")

Export to Pandas

import pandas as pd

# Convert to DataFrame
df = feature_map.get_df()

print(df.columns)
# Columns are lowercase: rt, mz, intensity, charge, quality

# Analyze features
print(f"Mean intensity: {df['intensity'].mean()}")
print(f"RT range: {df['rt'].min():.1f} - {df['rt'].max():.1f}")

Feature Linking

Map Alignment

Align retention times before linking:

# Load multiple feature maps
fm1 = ms.FeatureMap()
fm2 = ms.FeatureMap()
ms.FeatureXMLFile().load("sample1.featureXML", fm1)
ms.FeatureXMLFile().load("sample2.featureXML", fm2)
feature_maps = [fm1, fm2]

# Pick the largest map as the alignment reference
aligner = ms.MapAlignmentAlgorithmPoseClustering()
ref_idx = max(range(len(feature_maps)), key=lambda i: feature_maps[i].size())
aligner.setReference(feature_maps[ref_idx])

# Align each non-reference map in place against the reference
transformer = ms.MapAlignmentTransformer()
for i, fm in enumerate(feature_maps):
    if i == ref_idx:
        continue
    trafo = ms.TransformationDescription()
    aligner.align(fm, trafo)
    transformer.transformRetentionTimes(fm, trafo, True)

Feature Linking Algorithm

Link features across samples:

# Create feature grouping algorithm
grouper = ms.FeatureGroupingAlgorithmQT()

# Configure parameters
params = grouper.getParameters()
params.setValue("distance_RT:max_difference", 30.0)  # Max RT difference (s)
params.setValue("distance_MZ:max_difference", 10.0)  # Max m/z difference (ppm)
params.setValue("distance_MZ:unit", "ppm")
grouper.setParameters(params)

# Prepare feature maps
feature_maps = [fm1, fm2, fm3]

# Create consensus map
consensus_map = ms.ConsensusMap()

# Link features (feature_maps is a list of FeatureMap)
grouper.group(feature_maps, consensus_map)

# Assign unique IDs before storing
consensus_map.setUniqueIds()

print(f"Created {consensus_map.size()} consensus features")

# Save consensus map
ms.ConsensusXMLFile().store("consensus.consensusXML", consensus_map)

Consensus Features

Access Consensus Data

# Load consensus map
consensus_map = ms.ConsensusMap()
ms.ConsensusXMLFile().load("consensus.consensusXML", consensus_map)

# Iterate through consensus features
for cons_feature in consensus_map:
    print(f"Consensus m/z: {cons_feature.getMZ():.4f}")
    print(f"Consensus RT: {cons_feature.getRT():.2f}")

    # Get features from individual maps
    for handle in cons_feature.getFeatureList():
        map_idx = handle.getMapIndex()
        intensity = handle.getIntensity()
        print(f"  Sample {map_idx}: intensity {intensity:.0f}")

Consensus Map Metadata

# Access file descriptions (map metadata)
file_descriptions = consensus_map.getColumnHeaders()

for map_idx, description in file_descriptions.items():
    print(f"Map {map_idx}:")
    print(f"  Filename: {description.filename}")
    print(f"  Label: {description.label}")
    print(f"  Size: {description.size}")

Building Quant Matrices from a ConsensusMap

ConsensusMap exposes two DataFrame helpers that make quantitative tables easy:

# Feature intensities, features (rows) x samples (columns)
intensity_df = consensus_map.get_intensity_df()

# Per-consensus-feature metadata: rt, mz, charge, quality
metadata_df = consensus_map.get_metadata_df()

# Join into a single annotated quant matrix
quant = metadata_df.join(intensity_df)

Adduct Detection

Identify different ionization forms of the same molecule. The class is MetaboliteFeatureDeconvolution (the old MetaboliteAdductDecharger does not exist in 3.5.0). Adducts are specified with Elements:Charge:Probability syntax, not bracket notation like [M+H]+:

# Create adduct deconvolution
mfd = ms.MetaboliteFeatureDeconvolution()

# Configure parameters
p = mfd.getDefaults()
p.setValue("potential_adducts", [b"H:+:0.4", b"Na:+:0.25", b"NH4:+:0.25", b"K:+:0.1", b"H-2O-1:0:0.05"])
p.setValue("charge_min", 1)
p.setValue("charge_max", 1)
mfd.setParameters(p)

# Detect adducts: compute(in, out, cons_groups, cons_edges)
fm_out = ms.FeatureMap()
groups = ms.ConsensusMap()
edges = ms.ConsensusMap()
mfd.compute(feature_map, fm_out, groups, edges)

Complete Feature Detection Workflow

End-to-End Example

import pyopenms as ms

def feature_detection_workflow(input_files, output_consensus):
    """
    Complete workflow: feature detection and linking across samples.

    Args:
        input_files: List of mzML file paths
        output_consensus: Output consensusXML file path
    """

    feature_maps = []

    # Step 1: Detect features in each file (metabolomics pipeline)
    for mzml_file in input_files:
        print(f"Processing {mzml_file}...")

        # Load experiment
        exp = ms.MSExperiment()
        ms.MzMLFile().load(mzml_file, exp)
        exp.sortSpectra(True)

        # Mass trace detection
        mtd = ms.MassTraceDetection()
        p = mtd.getDefaults()
        p.setValue("mass_error_ppm", 10.0)
        p.setValue("noise_threshold_int", 1000.0)
        mtd.setParameters(p)
        mass_traces = []
        mtd.run(exp, mass_traces, 0)

        # Elution peak detection
        epd = ms.ElutionPeakDetection()
        p = epd.getDefaults()
        p.setValue("width_filtering", "fixed")
        epd.setParameters(p)
        mt_split = []
        epd.detectPeaks(mass_traces, mt_split)

        # Feature assembly
        ffm = ms.FeatureFindingMetabo()
        p = ffm.getDefaults()
        p.setValue("isotope_filtering_model", "metabolites (5% RMS)")
        p.setValue("remove_single_traces", "true")
        p.setValue("charge_lower_bound", 1)
        p.setValue("charge_upper_bound", 3)
        ffm.setParameters(p)
        features = ms.FeatureMap()
        chrom_out = []
        ffm.run(mt_split, features, chrom_out)

        # Store filename in feature map
        features.setPrimaryMSRunPath([mzml_file.encode()])

        feature_maps.append(features)
        print(f"  Found {features.size()} features")

    # Step 2: Align retention times against the largest map
    print("Aligning retention times...")
    aligner = ms.MapAlignmentAlgorithmPoseClustering()
    ref_idx = max(range(len(feature_maps)), key=lambda i: feature_maps[i].size())
    aligner.setReference(feature_maps[ref_idx])
    transformer = ms.MapAlignmentTransformer()
    for i, fm in enumerate(feature_maps):
        if i == ref_idx:
            continue
        trafo = ms.TransformationDescription()
        aligner.align(fm, trafo)
        transformer.transformRetentionTimes(fm, trafo, True)

    # Step 3: Link features
    print("Linking features across samples...")
    grouper = ms.FeatureGroupingAlgorithmQT()
    params = grouper.getParameters()
    params.setValue("distance_RT:max_difference", 30.0)
    params.setValue("distance_MZ:max_difference", 10.0)
    params.setValue("distance_MZ:unit", "ppm")
    grouper.setParameters(params)

    consensus_map = ms.ConsensusMap()
    grouper.group(feature_maps, consensus_map)
    consensus_map.setUniqueIds()

    # Save results
    ms.ConsensusXMLFile().store(output_consensus, consensus_map)

    print(f"Created {consensus_map.size()} consensus features")
    print(f"Results saved to {output_consensus}")

    return consensus_map

# Run workflow
input_files = ["sample1.mzML", "sample2.mzML", "sample3.mzML"]
consensus = feature_detection_workflow(input_files, "consensus.consensusXML")

Feature Filtering

Filter by Quality

# Filter features by quality score
filtered_features = ms.FeatureMap()

for feature in feature_map:
    if feature.getOverallQuality() > 0.5:  # Quality threshold
        filtered_features.push_back(feature)

print(f"Kept {filtered_features.size()} high-quality features")

Filter by Intensity

# Keep only intense features
min_intensity = 10000

filtered_features = ms.FeatureMap()
for feature in feature_map:
    if feature.getIntensity() >= min_intensity:
        filtered_features.push_back(feature)

Filter by m/z Range

# Extract features in specific m/z range
mz_min = 200.0
mz_max = 800.0

filtered_features = ms.FeatureMap()
for feature in feature_map:
    mz = feature.getMZ()
    if mz_min <= mz <= mz_max:
        filtered_features.push_back(feature)

Feature Annotation

Add Identification Information

# Annotate features with peptide identifications
# Load identifications
# pyOpenMS 3.5+: peptide IDs must be a PeptideIdentificationList, not a plain list
protein_ids = []
peptide_ids = ms.PeptideIdentificationList()
ms.IdXMLFile().load("identifications.idXML", protein_ids, peptide_ids)

# Create ID mapper
mapper = ms.IDMapper()

# Map IDs to features
mapper.annotate(feature_map, peptide_ids, protein_ids)

# Check annotations
for feature in feature_map:
    peptide_ids_for_feature = feature.getPeptideIdentifications()
    if peptide_ids_for_feature:
        print(f"Feature at {feature.getMZ():.4f} m/z identified")

Best Practices

Parameter Optimization

Optimize parameters for your data type:

# Test different mass-trace tolerance values (metabolomics pipeline)
mz_tolerances = [5.0, 10.0, 20.0]  # ppm

for tol in mz_tolerances:
    mtd = ms.MassTraceDetection()
    p = mtd.getDefaults()
    p.setValue("mass_error_ppm", tol)
    p.setValue("noise_threshold_int", 1000.0)
    mtd.setParameters(p)
    mass_traces = []
    mtd.run(exp, mass_traces, 0)

    epd = ms.ElutionPeakDetection()
    mt_split = []
    epd.detectPeaks(mass_traces, mt_split)

    ffm = ms.FeatureFindingMetabo()
    features = ms.FeatureMap()
    chrom_out = []
    ffm.run(mt_split, features, chrom_out)

    print(f"Tolerance {tol} ppm: {features.size()} features")

Visual Inspection

Export features for visualization:

# Convert to DataFrame for plotting
df = feature_map.get_df()

import matplotlib.pyplot as plt

plt.figure(figsize=(10, 6))
plt.scatter(df['rt'], df['mz'], s=df['intensity']/1000, alpha=0.5)
plt.xlabel('Retention Time (s)')
plt.ylabel('m/z')
plt.title('Feature Map')
plt.colorbar(label='Intensity (scaled)')
plt.show()

File I/O and Data Formats

Overview

PyOpenMS supports multiple mass spectrometry file formats for reading and writing. This guide covers file handling strategies and format-specific operations.

Supported Formats

Spectrum Data Formats

• mzML: Standard XML-based format for mass spectrometry data
• mzXML: Earlier XML-based format
• mzData: XML format (deprecated but supported)

Identification Formats

• idXML: OpenMS native identification format
• mzIdentML: Standard XML format for identification data
• pepXML: X! Tandem format
• protXML: Protein identification format

Feature and Quantitation Formats

• featureXML: OpenMS format for detected features
• consensusXML: Format for consensus features across samples
• mzTab: Tab-delimited format for reporting

Sequence and Library Formats

• FASTA: Protein/peptide sequences
• TraML: Transition lists for targeted experiments

Reading mzML Files

In-Memory Loading

Load entire file into memory (suitable for smaller files):

import pyopenms as ms

# Create experiment container
exp = ms.MSExperiment()

# Load file
ms.MzMLFile().load("sample.mzML", exp)

# Access data
print(f"Spectra: {exp.getNrSpectra()}")
print(f"Chromatograms: {exp.getNrChromatograms()}")

Indexed Access

Efficient random access for large files:

# Create indexed access
indexed_mzml = ms.IndexedMzMLFileLoader()
indexed_mzml.load("large_file.mzML")

# Get specific spectrum by index
spec = indexed_mzml.getSpectrumById(100)

# Access by native ID
spec = indexed_mzml.getSpectrumByNativeId("scan=5000")

Streaming / On-Disc Access

Memory-efficient processing for very large files uses OnDiscMSExperiment, which parses the index of an indexed mzML and loads spectra on demand instead of holding the whole run in memory. (The old MSExperimentConsumer subclassing pattern is not available in pyOpenMS 3.5.)

# Requires an indexed mzML. Write one with the write-index option set:
exp = ms.MSExperiment()
ms.MzMLFile().load("large.mzML", exp)
f = ms.MzMLFile()
opt = f.getOptions(); opt.setWriteIndex(True); f.setOptions(opt)
f.store("large_indexed.mzML", exp)

# Now access spectra lazily, one at a time
od = ms.OnDiscMSExperiment()
if od.openFile("large_indexed.mzML"):
    count = 0
    for i in range(od.getNrSpectra()):
        spec = od.getSpectrum(i)   # loaded from disk on demand
        if spec.getMSLevel() == 2:
            count += 1
    print(f"Processed {count} MS2 spectra")

Cached Access

A cached binary representation trades a little disk space for faster repeated reads. Use the static CachedmzML.store/load methods:

# Write a cached representation
exp = ms.MSExperiment()
ms.MzMLFile().load("sample.mzML", exp)
ms.CachedmzML().store("sample.cachedmzML", exp)

# Load it back for on-demand spectrum access
cached = ms.CachedmzML()
ms.CachedmzML().load("sample.cachedmzML", cached)
print(f"{cached.getNrSpectra()} spectra")
spec = cached.getSpectrum(0)

Writing mzML Files

Basic Writing

# Create or modify experiment
exp = ms.MSExperiment()
# ... add spectra ...

# Write to file
ms.MzMLFile().store("output.mzML", exp)

Compression Options

# Configure compression
file_handler = ms.MzMLFile()

options = ms.PeakFileOptions()
options.setCompression(True)  # Enable compression
file_handler.setOptions(options)

file_handler.store("compressed.mzML", exp)

Reading Identification Data

idXML Format

# Load identification results
protein_ids = []                              # protein IDs: plain list
# pyOpenMS 3.5+: peptide IDs must be a PeptideIdentificationList, not a plain list
peptide_ids = ms.PeptideIdentificationList()

ms.IdXMLFile().load("identifications.idXML", protein_ids, peptide_ids)

# Access peptide identifications
for peptide_id in peptide_ids:
    print(f"RT: {peptide_id.getRT()}")
    print(f"MZ: {peptide_id.getMZ()}")

    # Get peptide hits
    for hit in peptide_id.getHits():
        print(f"  Sequence: {hit.getSequence().toString()}")
        print(f"  Score: {hit.getScore()}")
        print(f"  Charge: {hit.getCharge()}")

mzIdentML Format

# Read mzIdentML
protein_ids = []                              # protein IDs: plain list
peptide_ids = ms.PeptideIdentificationList()  # pyOpenMS 3.5+: not a plain list

ms.MzIdentMLFile().load("results.mzid", protein_ids, peptide_ids)

pepXML Format

# Load pepXML
protein_ids = []                              # protein IDs: plain list
peptide_ids = ms.PeptideIdentificationList()  # pyOpenMS 3.5+: not a plain list

ms.PepXMLFile().load("results.pep.xml", protein_ids, peptide_ids)

Reading Feature Data

featureXML

# Load features
feature_map = ms.FeatureMap()
ms.FeatureXMLFile().load("features.featureXML", feature_map)

# Access features
for feature in feature_map:
    print(f"RT: {feature.getRT()}")
    print(f"MZ: {feature.getMZ()}")
    print(f"Intensity: {feature.getIntensity()}")
    print(f"Quality: {feature.getOverallQuality()}")

consensusXML

# Load consensus features
consensus_map = ms.ConsensusMap()
ms.ConsensusXMLFile().load("consensus.consensusXML", consensus_map)

# Access consensus features
for consensus_feature in consensus_map:
    print(f"RT: {consensus_feature.getRT()}")
    print(f"MZ: {consensus_feature.getMZ()}")

    # Get feature handles (sub-features from different maps)
    for handle in consensus_feature.getFeatureList():
        map_index = handle.getMapIndex()
        intensity = handle.getIntensity()
        print(f"  Map {map_index}: {intensity}")

Reading FASTA Files

# Load protein sequences
fasta_entries = []
ms.FASTAFile().load("database.fasta", fasta_entries)

for entry in fasta_entries:
    print(f"Identifier: {entry.identifier}")
    print(f"Description: {entry.description}")
    print(f"Sequence: {entry.sequence}")

Reading TraML Files

# Load transition lists for targeted experiments
targeted_exp = ms.TargetedExperiment()
ms.TraMLFile().load("transitions.TraML", targeted_exp)

# Access transitions
for transition in targeted_exp.getTransitions():
    print(f"Precursor MZ: {transition.getPrecursorMZ()}")
    print(f"Product MZ: {transition.getProductMZ()}")

Writing mzTab Files

# Create mzTab for reporting
mztab = ms.MzTab()

# Add metadata
metadata = mztab.getMetaData()
metadata.mz_tab_version.set("1.0.0")
metadata.title.set("Proteomics Analysis Results")

# Add protein data
protein_section = mztab.getProteinSectionRows()
# ... populate protein data ...

# Write to file
ms.MzTabFile().store("report.mzTab", mztab)

Format Conversion

mzXML to mzML

# Read mzXML
exp = ms.MSExperiment()
ms.MzXMLFile().load("data.mzXML", exp)

# Write as mzML
ms.MzMLFile().store("data.mzML", exp)

Extract Chromatograms from mzML

# Load experiment
exp = ms.MSExperiment()
ms.MzMLFile().load("data.mzML", exp)

# Extract specific chromatogram
for chrom in exp.getChromatograms():
    if chrom.getNativeID() == "TIC":
        rt, intensity = chrom.get_peaks()
        print(f"TIC has {len(rt)} data points")

File Metadata

Access mzML Metadata

# Load file
exp = ms.MSExperiment()
ms.MzMLFile().load("sample.mzML", exp)

# Get experimental settings
exp_settings = exp.getExperimentalSettings()

# Instrument info
instrument = exp_settings.getInstrument()
print(f"Instrument: {instrument.getName()}")
print(f"Model: {instrument.getModel()}")

# Sample info
sample = exp_settings.getSample()
print(f"Sample name: {sample.getName()}")

# Source files
for source_file in exp_settings.getSourceFiles():
    print(f"Source: {source_file.getNameOfFile()}")

Best Practices

Memory Management

For large files: 1. Use indexed or streaming access instead of full in-memory loading 2. Process data in chunks 3. Clear data structures when no longer needed

# Good for large files
indexed_mzml = ms.IndexedMzMLFileLoader()
indexed_mzml.load("huge_file.mzML")

# Process spectra one at a time
for i in range(indexed_mzml.getNrSpectra()):
    spec = indexed_mzml.getSpectrumById(i)
    # Process spectrum
    # Spectrum automatically cleaned up after processing

Error Handling

try:
    exp = ms.MSExperiment()
    ms.MzMLFile().load("data.mzML", exp)
except Exception as e:
    print(f"Failed to load file: {e}")

File Validation

# Check if file exists and is readable
import os

if os.path.exists("data.mzML") and os.path.isfile("data.mzML"):
    exp = ms.MSExperiment()
    ms.MzMLFile().load("data.mzML", exp)
else:
    print("File not found")

Peptide and Protein Identification

Overview

PyOpenMS supports peptide/protein identification through integration with search engines and provides tools for post-processing identification results including FDR control, protein inference, and annotation.

The skill ships ready-to-use CLI scripts: scripts/process_identifications.py performs FDR filtering and export on idXML files, and scripts/inspect_ms_data.py summarizes idXML (and other MS) files.

pyOpenMS 3.5+ API note: IdXMLFile().load()/store() require the peptide-IDs argument to be a ms.PeptideIdentificationList(), not a plain Python list. Passing a plain list raises Exception: can not handle type of (...). The protein-IDs argument is still a plain Python list.

Supported Search Engines

PyOpenMS integrates with these search engines:

• Comet: Fast tandem MS search
• Mascot: Commercial search engine
• MSGFPlus: Spectral probability-based search
• XTandem: Open-source search tool
• OMSSA: NCBI search engine
• Myrimatch: High-throughput search
• MSFragger: Ultra-fast search

Reading Identification Data

idXML Format

import pyopenms as ms

# Load identification results
protein_ids = []                              # protein IDs: plain list
# pyOpenMS 3.5+: peptide IDs must be a PeptideIdentificationList, not a plain list
peptide_ids = ms.PeptideIdentificationList()

ms.IdXMLFile().load("identifications.idXML", protein_ids, peptide_ids)

print(f"Protein identifications: {len(protein_ids)}")
print(f"Peptide identifications: {len(peptide_ids)}")

Access Peptide Identifications

# Iterate through peptide IDs
for peptide_id in peptide_ids:
    # Spectrum metadata
    print(f"RT: {peptide_id.getRT():.2f}")
    print(f"m/z: {peptide_id.getMZ():.4f}")

    # Get peptide hits (ranked by score)
    hits = peptide_id.getHits()
    print(f"Number of hits: {len(hits)}")

    for hit in hits:
        sequence = hit.getSequence()
        print(f"  Sequence: {sequence.toString()}")
        print(f"  Score: {hit.getScore()}")
        print(f"  Charge: {hit.getCharge()}")
        print(f"  Mass error (ppm): {hit.getMetaValue('mass_error_ppm')}")

        # Get modifications
        if sequence.isModified():
            for i in range(sequence.size()):
                residue = sequence.getResidue(i)
                if residue.isModified():
                    print(f"    Modification at position {i}: {residue.getModificationName()}")

Access Protein Identifications

# Access protein-level information
for protein_id in protein_ids:
    # Search parameters
    search_params = protein_id.getSearchParameters()
    print(f"Search engine: {protein_id.getSearchEngine()}")
    print(f"Database: {search_params.db}")

    # Protein hits
    hits = protein_id.getHits()
    for hit in hits:
        print(f"  Accession: {hit.getAccession()}")
        print(f"  Score: {hit.getScore()}")
        print(f"  Coverage: {hit.getCoverage()}")
        print(f"  Sequence: {hit.getSequence()}")

False Discovery Rate (FDR)

FDR Filtering

Apply FDR filtering to control false positives:

# Create FDR object
fdr = ms.FalseDiscoveryRate()

# Apply FDR at PSM level
fdr.apply(peptide_ids)

# Filter by FDR threshold
fdr_threshold = 0.01  # 1% FDR
filtered_peptide_ids = ms.PeptideIdentificationList()  # use push_back to add IDs

for peptide_id in peptide_ids:
    # Keep hits below FDR threshold
    filtered_hits = []
    for hit in peptide_id.getHits():
        if hit.getScore() <= fdr_threshold:  # Lower score = better
            filtered_hits.append(hit)

    if filtered_hits:
        peptide_id.setHits(filtered_hits)
        filtered_peptide_ids.push_back(peptide_id)

print(f"Peptides passing FDR: {len(filtered_peptide_ids)}")

Score Transformation

Convert scores to q-values:

# Apply score transformation
fdr.apply(peptide_ids)

# Access q-values
for peptide_id in peptide_ids:
    for hit in peptide_id.getHits():
        q_value = hit.getMetaValue("q-value")
        print(f"Sequence: {hit.getSequence().toString()}, q-value: {q_value}")

Protein Inference

ID Mapper

Map peptide identifications to proteins:

# Create mapper
mapper = ms.IDMapper()

# Map to features
feature_map = ms.FeatureMap()
ms.FeatureXMLFile().load("features.featureXML", feature_map)

# Annotate features with IDs
mapper.annotate(feature_map, peptide_ids, protein_ids)

# Check annotated features
for feature in feature_map:
    pep_ids = feature.getPeptideIdentifications()
    if pep_ids:
        for pep_id in pep_ids:
            for hit in pep_id.getHits():
                print(f"Feature {feature.getMZ():.4f}: {hit.getSequence().toString()}")

Protein Grouping

Group proteins by shared peptides:

# Create protein inference algorithm
inference = ms.BasicProteinInferenceAlgorithm()

# Run inference
inference.run(peptide_ids, protein_ids)

# Access protein groups
for protein_id in protein_ids:
    hits = protein_id.getHits()
    if len(hits) > 1:
        print("Protein group:")
        for hit in hits:
            print(f"  {hit.getAccession()}")

Peptide Sequence Handling

AASequence Object

Work with peptide sequences:

# Create peptide sequence
seq = ms.AASequence.fromString("PEPTIDE")

print(f"Sequence: {seq.toString()}")
print(f"Monoisotopic mass: {seq.getMonoWeight():.4f}")
print(f"Average mass: {seq.getAverageWeight():.4f}")
print(f"Length: {seq.size()}")

# Access individual amino acids
for i in range(seq.size()):
    residue = seq.getResidue(i)
    print(f"Position {i}: {residue.getOneLetterCode()}, mass: {residue.getMonoWeight():.4f}")

Modified Sequences

Handle post-translational modifications:

# Sequence with modifications
mod_seq = ms.AASequence.fromString("PEPTIDEM(Oxidation)K")

print(f"Modified sequence: {mod_seq.toString()}")
print(f"Mass with mods: {mod_seq.getMonoWeight():.4f}")

# Check if modified
print(f"Is modified: {mod_seq.isModified()}")

# Get modification info
for i in range(mod_seq.size()):
    residue = mod_seq.getResidue(i)
    if residue.isModified():
        print(f"Residue {residue.getOneLetterCode()} at position {i}")
        print(f"  Modification: {residue.getModificationName()}")

Peptide Digestion

Simulate enzymatic digestion:

# Create digestion enzyme
enzyme = ms.ProteaseDigestion()
enzyme.setEnzyme("Trypsin")

# Set missed cleavages
enzyme.setMissedCleavages(2)

# Digest protein sequence
protein_seq = "MKTAYIAKQRQISFVKSHFSRQLEERLGLIEVQAPILSRVGDGTQDNLSGAEKAVQVKVKALPDAQFEVVHSLAKWKRQTLGQHDFSAGEGLYTHMKALRPDEDRLSPLHSVYVDQWDWERVMGDGERQFSTLKSTVEAIWAGIKATEAAVSEEFGLAPFLPDQIHFVHSQELLSRYPDLDAKGRERAIAKDLGAVFLVGIGGKLSDGHRHDVRAPDYDDWSTPSELGHAGLNGDILVWNPVLEDAFELSSMGIRVDADTLKHQLALTGDEDRLELEWHQALLRGEMPQTIGGGIGQSRLTMLLLQLPHIGQVQAGVWPAAVRESVPSLL"

# Get peptides
peptides = []
enzyme.digest(ms.AASequence.fromString(protein_seq), peptides)

print(f"Generated {len(peptides)} peptides")
for peptide in peptides[:5]:  # Show first 5
    print(f"  {peptide.toString()}, mass: {peptide.getMonoWeight():.2f}")

Theoretical Spectrum Generation

Fragment Ion Calculation

Generate theoretical fragment ions:

# Create peptide
peptide = ms.AASequence.fromString("PEPTIDE")

# Generate b and y ions
fragments = []
ms.TheoreticalSpectrumGenerator().getSpectrum(fragments, peptide, 1, 1)

print(f"Generated {fragments.size()} fragment ions")

# Access fragments
mz, intensity = fragments.get_peaks()
for m, i in zip(mz[:10], intensity[:10]):  # Show first 10
    print(f"m/z: {m:.4f}, intensity: {i}")

Complete Identification Workflow

End-to-End Example

import pyopenms as ms

def identification_workflow(spectrum_file, fasta_file, output_file):
    """
    Complete identification workflow with FDR control.

    Args:
        spectrum_file: Input mzML file
        fasta_file: Protein database (FASTA)
        output_file: Output idXML file
    """

    # Step 1: Load spectra
    exp = ms.MSExperiment()
    ms.MzMLFile().load(spectrum_file, exp)
    print(f"Loaded {exp.getNrSpectra()} spectra")

    # Step 2: Configure search parameters
    search_params = ms.SearchParameters()
    search_params.db = fasta_file
    search_params.precursor_mass_tolerance = 10.0  # ppm
    search_params.fragment_mass_tolerance = 0.5  # Da
    search_params.enzyme = "Trypsin"
    search_params.missed_cleavages = 2
    search_params.modifications = ["Oxidation (M)", "Carbamidomethyl (C)"]

    # Step 3: Run the database search.
    # NOTE: OpenMS search engines (Comet, MSGFPlus, XTandem, ...) are exposed as
    # command-line TOPP tools, NOT pyOpenMS Python classes (there is no
    # ms.CometAdapter). Run them as a subprocess and read the idXML they emit, e.g.:
    #   subprocess.run(["CometAdapter", "-in", "spectra.mzML",
    #                   "-database", fasta_file, "-out", "raw_identifications.idXML"])
    # These adapter executables ship with an OpenMS (not pyOpenMS) installation.

    # Here we load the pre-computed search results
    protein_ids = []                              # protein IDs: plain list
    # pyOpenMS 3.5+: peptide IDs must be a PeptideIdentificationList, not a plain list
    peptide_ids = ms.PeptideIdentificationList()
    ms.IdXMLFile().load("raw_identifications.idXML", protein_ids, peptide_ids)

    print(f"Initial peptide IDs: {len(peptide_ids)}")

    # Step 4: Apply FDR filtering
    fdr = ms.FalseDiscoveryRate()
    fdr.apply(peptide_ids)

    # Filter by 1% FDR
    filtered_peptide_ids = ms.PeptideIdentificationList()  # use push_back to add IDs
    for peptide_id in peptide_ids:
        filtered_hits = []
        for hit in peptide_id.getHits():
            q_value = hit.getMetaValue("q-value")
            if q_value <= 0.01:
                filtered_hits.append(hit)

        if filtered_hits:
            peptide_id.setHits(filtered_hits)
            filtered_peptide_ids.push_back(peptide_id)

    print(f"Peptides after FDR (1%): {len(filtered_peptide_ids)}")

    # Step 5: Protein inference
    inference = ms.BasicProteinInferenceAlgorithm()
    inference.run(filtered_peptide_ids, protein_ids)

    print(f"Identified proteins: {len(protein_ids)}")

    # Step 6: Save results
    ms.IdXMLFile().store(output_file, protein_ids, filtered_peptide_ids)
    print(f"Results saved to {output_file}")

    return protein_ids, filtered_peptide_ids

# Run workflow
protein_ids, peptide_ids = identification_workflow(
    "spectra.mzML",
    "database.fasta",
    "identifications_fdr.idXML"
)

Spectral Library Search

Library Matching

# Load spectral library
library = ms.MSPFile()
library_spectra = []
library.load("spectral_library.msp", library_spectra)

# Load experimental spectra
exp = ms.MSExperiment()
ms.MzMLFile().load("data.mzML", exp)

# Compare spectra
spectra_compare = ms.SpectraSTSimilarityScore()

for exp_spec in exp:
    if exp_spec.getMSLevel() == 2:
        best_match_score = 0
        best_match_lib = None

        for lib_spec in library_spectra:
            score = spectra_compare.operator()(exp_spec, lib_spec)
            if score > best_match_score:
                best_match_score = score
                best_match_lib = lib_spec

        if best_match_score > 0.7:  # Threshold
            print(f"Match found: score {best_match_score:.3f}")

Best Practices

Decoy Database

Use target-decoy approach for FDR calculation:

# Generate decoy database
decoy_generator = ms.DecoyGenerator()

# Load target database
fasta_entries = []
ms.FASTAFile().load("target.fasta", fasta_entries)

# Generate decoys
decoy_entries = []
for entry in fasta_entries:
    decoy_entry = decoy_generator.reverseProtein(entry)
    decoy_entries.append(decoy_entry)

# Save combined database
all_entries = fasta_entries + decoy_entries
ms.FASTAFile().store("target_decoy.fasta", all_entries)

Score Interpretation

Understand score types from different engines:

# Interpret scores based on search engine
for peptide_id in peptide_ids:
    search_engine = peptide_id.getIdentifier()

    for hit in peptide_id.getHits():
        score = hit.getScore()

        # Score interpretation varies by engine
        if "Comet" in search_engine:
            # Comet: higher E-value = worse
            print(f"E-value: {score}")
        elif "Mascot" in search_engine:
            # Mascot: higher score = better
            print(f"Ion score: {score}")

Metabolomics Workflows

Overview

PyOpenMS provides specialized tools for untargeted metabolomics analysis including feature detection optimized for small molecules, adduct grouping, compound identification, and integration with metabolomics databases.

Code examples target pyOpenMS 3.5.0. APIs removed in 3.5.0 (e.g. FeatureFinder().run("centroided", ...) and MetaboliteAdductDecharger) are not used here. The skill also ships ready-to-run scripts implementing these workflows end-to-end: scripts/detect_features_metabo.py, scripts/align_link_quantify.py, scripts/consensus_to_matrix.py, scripts/detect_adducts.py, scripts/accurate_mass_search.py, and scripts/export_gnps_sirius.py.

Untargeted Metabolomics Pipeline

Complete Workflow

import pyopenms as ms

def metabolomics_pipeline(input_files, output_dir):
    """
    Complete untargeted metabolomics workflow.

    Args:
        input_files: List of mzML file paths (one per sample)
        output_dir: Directory for output files
    """

    # Step 1: Feature detection (MassTraceDetection ->
    # ElutionPeakDetection -> FeatureFindingMetabo)
    feature_maps = []

    for mzml_file in input_files:
        print(f"Processing {mzml_file}...")

        # Load data
        exp = ms.MSExperiment()
        ms.MzMLFile().load(mzml_file, exp)
        exp.sortSpectra(True)

        # Detect mass traces
        mtd = ms.MassTraceDetection()
        p = mtd.getDefaults()
        p.setValue("mass_error_ppm", 10.0)
        p.setValue("noise_threshold_int", 1000.0)
        mtd.setParameters(p)
        mass_traces = []
        mtd.run(exp, mass_traces, 0)

        # Split traces into elution peaks
        epd = ms.ElutionPeakDetection()
        p = epd.getDefaults()
        p.setValue("width_filtering", "fixed")
        epd.setParameters(p)
        mt_split = []
        epd.detectPeaks(mass_traces, mt_split)

        # Assemble metabolite features
        ffm = ms.FeatureFindingMetabo()
        p = ffm.getDefaults()
        p.setValue("isotope_filtering_model", "metabolites (5% RMS)")
        p.setValue("remove_single_traces", "true")
        p.setValue("charge_lower_bound", 1)
        p.setValue("charge_upper_bound", 3)
        ffm.setParameters(p)
        features = ms.FeatureMap()
        chrom_out = []
        ffm.run(mt_split, features, chrom_out)
        features.setUniqueIds()

        features.setPrimaryMSRunPath([mzml_file.encode()])
        feature_maps.append(features)

        print(f"  Detected {features.size()} features")

    # Step 2: Adduct detection and grouping
    print("Detecting adducts...")
    adduct_grouped_maps = []

    for fm in feature_maps:
        mfd = ms.MetaboliteFeatureDeconvolution()
        p = mfd.getDefaults()
        # potential_adducts uses Elements:Charge:Probability syntax
        p.setValue("potential_adducts",
                   [b"H:+:0.4", b"Na:+:0.25", b"NH4:+:0.25",
                    b"K:+:0.1", b"H-2O-1:0:0.05"])
        p.setValue("charge_min", 1)
        p.setValue("charge_max", 1)
        mfd.setParameters(p)

        fm_out = ms.FeatureMap()
        groups = ms.ConsensusMap()
        edges = ms.ConsensusMap()
        mfd.compute(fm, fm_out, groups, edges)  # 4 args
        adduct_grouped_maps.append(fm_out)

    # Step 3: RT alignment (PoseClustering: pick a reference,
    # then align each map in place against it)
    print("Aligning retention times...")
    aligner = ms.MapAlignmentAlgorithmPoseClustering()
    reference = adduct_grouped_maps[0]
    aligner.setReference(reference)

    transformer = ms.MapAlignmentTransformer()
    for fm in adduct_grouped_maps:
        trafo = ms.TransformationDescription()
        aligner.align(fm, trafo)
        transformer.transformRetentionTimes(fm, trafo, True)

    aligned_maps = adduct_grouped_maps

    # Step 4: Feature linking
    print("Linking features...")
    grouper = ms.FeatureGroupingAlgorithmQT()

    params = grouper.getParameters()
    params.setValue("distance_RT:max_difference", 60.0)  # seconds
    params.setValue("distance_MZ:max_difference", 5.0)  # ppm
    params.setValue("distance_MZ:unit", "ppm")
    grouper.setParameters(params)

    consensus_map = ms.ConsensusMap()
    grouper.group(aligned_maps, consensus_map)  # list of FeatureMap
    consensus_map.setUniqueIds()

    print(f"Created {consensus_map.size()} consensus features")

    # Step 5: Export results
    consensus_file = f"{output_dir}/consensus.consensusXML"
    ms.ConsensusXMLFile().store(consensus_file, consensus_map)

    # Export quant matrix for downstream analysis
    # get_intensity_df() -> features x samples; get_metadata_df() -> rt/mz/charge/quality
    intensities = consensus_map.get_intensity_df()
    metadata = consensus_map.get_metadata_df()
    csv_file = f"{output_dir}/metabolite_table.csv"
    metadata.join(intensities).to_csv(csv_file)

    print(f"Results saved to {output_dir}")

    return consensus_map

# Run pipeline
input_files = ["sample1.mzML", "sample2.mzML", "sample3.mzML"]
consensus = metabolomics_pipeline(input_files, "output")

Adduct Detection

Configure Adduct Types

Adducts are configured with the Elements:Charge:Probability syntax (e.g. b"Na:+:0.25"), not bracket notation like [M+Na]+. Probabilities across the list should sum to ~1.0.

# Create adduct detector
mfd = ms.MetaboliteFeatureDeconvolution()

# Configure common adducts
params = mfd.getDefaults()

# Positive mode adducts (Elements:Charge:Probability)
positive_adducts = [
    b"H:+:0.4",
    b"Na:+:0.25",
    b"NH4:+:0.25",
    b"K:+:0.1",
    b"H-2O-1:0:0.05",  # neutral water loss
]

# Negative mode adducts (set negative_mode="true" when using these)
negative_adducts = [
    b"H-1:-:0.6",
    b"Cl:-:0.2",
]

# Set for positive mode
params.setValue("potential_adducts", positive_adducts)
params.setValue("charge_min", 1)
params.setValue("charge_max", 1)
mfd.setParameters(params)

# Apply adduct detection (4 args: in, out, groups, edges)
feature_map_out = ms.FeatureMap()
groups = ms.ConsensusMap()
edges = ms.ConsensusMap()
mfd.compute(feature_map, feature_map_out, groups, edges)

Access Adduct Information

# Check adduct annotations
for feature in feature_map_out:
    # Get adduct type if annotated
    if feature.metaValueExists("adduct"):
        adduct = feature.getMetaValue("adduct")
        neutral_mass = feature.getMetaValue("neutral_mass")
        print(f"m/z: {feature.getMZ():.4f}")
        print(f"  Adduct: {adduct}")
        print(f"  Neutral mass: {neutral_mass:.4f}")

Compound Identification

Accurate Mass Search (HMDB)

The built-in AccurateMassSearchEngine annotates a FeatureMap against HMDB and writes results as mzTab.

Caveat: the pip wheel ships HMDBMappingFile.tsv but not

HMDB2StructMapping.tsv, so engine.init() fails unless you supply the

struct file yourself. Download it from

https://github.com/OpenMS/OpenMS/blob/develop/share/OpenMS/CHEMISTRY/HMDB2StructMapping.tsv

and point db:struct at it.

engine = ms.AccurateMassSearchEngine()
p = engine.getDefaults()
p.setValue("mass_error_value", 5.0)
p.setValue("mass_error_unit", "ppm")
p.setValue("ionization_mode", "positive")
# If the struct file is missing from the wheel, supply it:
# p.setValue("db:struct", b"/path/to/HMDB2StructMapping.tsv")
engine.setParameters(p)
engine.init()

mztab = ms.MzTab()
engine.run(feature_map, mztab)
ms.MzTabFile().store("out.mzTab", mztab)

Mass-Based Annotation (custom database)

# Annotate features against your own compound list.
# (Plain Python, no pyOpenMS-specific DB required.)

# Load compound database (example structure)
# In practice, use external database like HMDB, METLIN

compound_db = [
    {"name": "Glucose", "formula": "C6H12O6", "mass": 180.0634},
    {"name": "Citric acid", "formula": "C6H8O7", "mass": 192.0270},
    # ... more compounds
]

# Annotate features
mass_tolerance = 5.0  # ppm

for feature in feature_map:
    observed_mz = feature.getMZ()

    # Calculate neutral mass (assuming [M+H]+)
    neutral_mass = observed_mz - 1.007276  # Proton mass

    # Search database
    for compound in compound_db:
        mass_error_ppm = abs(neutral_mass - compound["mass"]) / compound["mass"] * 1e6

        if mass_error_ppm <= mass_tolerance:
            print(f"Potential match: {compound['name']}")
            print(f"  Observed m/z: {observed_mz:.4f}")
            print(f"  Expected mass: {compound['mass']:.4f}")
            print(f"  Error: {mass_error_ppm:.2f} ppm")

MS/MS-Based Identification

# Load MS2 data
exp = ms.MSExperiment()
ms.MzMLFile().load("data_with_ms2.mzML", exp)

# Extract MS2 spectra
ms2_spectra = []
for spec in exp:
    if spec.getMSLevel() == 2:
        ms2_spectra.append(spec)

print(f"Found {len(ms2_spectra)} MS2 spectra")

# Match to spectral library
# (Requires external tool or custom implementation)

Data Normalization

Total Ion Current (TIC) Normalization

import numpy as np

# Load consensus map
consensus_map = ms.ConsensusMap()
ms.ConsensusXMLFile().load("consensus.consensusXML", consensus_map)

# Calculate TIC per sample
n_samples = len(consensus_map.getColumnHeaders())
tic_per_sample = np.zeros(n_samples)

for cons_feature in consensus_map:
    for handle in cons_feature.getFeatureList():
        map_idx = handle.getMapIndex()
        tic_per_sample[map_idx] += handle.getIntensity()

print("TIC per sample:", tic_per_sample)

# Normalize to median TIC
median_tic = np.median(tic_per_sample)
normalization_factors = median_tic / tic_per_sample

print("Normalization factors:", normalization_factors)

# Apply normalization
consensus_map_normalized = ms.ConsensusMap(consensus_map)
for cons_feature in consensus_map_normalized:
    feature_list = cons_feature.getFeatureList()
    for handle in feature_list:
        map_idx = handle.getMapIndex()
        normalized_intensity = handle.getIntensity() * normalization_factors[map_idx]
        handle.setIntensity(normalized_intensity)

Quality Control

Coefficient of Variation (CV) Filtering

import pandas as pd
import numpy as np

# Export intensities to pandas (features x samples)
df = consensus_map.get_intensity_df()

# Assume QC samples are columns with 'QC' in name
qc_cols = [col for col in df.columns if 'QC' in col]

if qc_cols:
    # Calculate CV for each feature in QC samples
    qc_data = df[qc_cols]
    cv = (qc_data.std(axis=1) / qc_data.mean(axis=1)) * 100

    # Filter features with CV < 30% in QC samples
    good_features = df[cv < 30]

    print(f"Features before CV filter: {len(df)}")
    print(f"Features after CV filter: {len(good_features)}")

Blank Filtering

# Remove features present in blank samples
blank_cols = [col for col in df.columns if 'Blank' in col]
sample_cols = [col for col in df.columns if 'Sample' in col]

if blank_cols and sample_cols:
    # Calculate mean intensity in blanks and samples
    blank_mean = df[blank_cols].mean(axis=1)
    sample_mean = df[sample_cols].mean(axis=1)

    # Keep features with 3x higher intensity in samples than blanks
    ratio = sample_mean / (blank_mean + 1)  # Add 1 to avoid division by zero
    filtered_df = df[ratio > 3]

    print(f"Features before blank filtering: {len(df)}")
    print(f"Features after blank filtering: {len(filtered_df)}")

Missing Value Imputation

import pandas as pd
import numpy as np

# Load intensities (features x samples)
df = consensus_map.get_intensity_df()

# Replace zeros with NaN
df = df.replace(0, np.nan)

# Count missing values
missing_per_feature = df.isnull().sum(axis=1)
print(f"Features with >50% missing: {sum(missing_per_feature > len(df.columns)/2)}")

# Simple imputation: replace with minimum value
for col in df.columns:
    if df[col].dtype in [np.float64, np.int64]:
        min_val = df[col].min() / 2  # Half minimum
        df[col].fillna(min_val, inplace=True)

Metabolite Table Export

Create Analysis-Ready Table

import pandas as pd

def create_metabolite_table(consensus_map, output_file):
    """
    Create metabolite quantification table for statistical analysis.
    """

    # Get column headers (file descriptions)
    headers = consensus_map.getColumnHeaders()

    # Initialize data structure
    data = {
        'mz': [],
        'rt': [],
        'feature_id': []
    }

    # Add sample columns
    for map_idx, header in headers.items():
        sample_name = header.label or f"Sample_{map_idx}"
        data[sample_name] = []

    # Extract feature data
    for idx, cons_feature in enumerate(consensus_map):
        data['mz'].append(cons_feature.getMZ())
        data['rt'].append(cons_feature.getRT())
        data['feature_id'].append(f"F{idx:06d}")

        # Initialize intensities
        intensities = {map_idx: 0.0 for map_idx in headers.keys()}

        # Fill in measured intensities
        for handle in cons_feature.getFeatureList():
            map_idx = handle.getMapIndex()
            intensities[map_idx] = handle.getIntensity()

        # Add to data structure
        for map_idx, header in headers.items():
            sample_name = header.label or f"Sample_{map_idx}"
            data[sample_name].append(intensities[map_idx])

    # Create DataFrame
    df = pd.DataFrame(data)

    # Sort by RT
    df = df.sort_values('rt')

    # Save to CSV
    df.to_csv(output_file, index=False)

    print(f"Metabolite table with {len(df)} features saved to {output_file}")

    return df

# Create table
df = create_metabolite_table(consensus_map, "metabolite_table.csv")

Integration with External Tools

Export for MetaboAnalyst

def export_for_metaboanalyst(df, output_file):
    """
    Format data for MetaboAnalyst input.

    Requires sample names as columns, features as rows.
    """

    # Transpose DataFrame
    # Remove metadata columns
    sample_cols = [col for col in df.columns if col not in ['mz', 'rt', 'feature_id']]

    # Extract sample data
    sample_data = df[sample_cols]

    # Transpose (samples as rows, features as columns)
    df_transposed = sample_data.T

    # Add feature identifiers as column names
    df_transposed.columns = df['feature_id']

    # Save
    df_transposed.to_csv(output_file)

    print(f"MetaboAnalyst format saved to {output_file}")

# Export
export_for_metaboanalyst(df, "for_metaboanalyst.csv")

Best Practices

Sample Size and Replicates

• Include QC samples (pooled sample) every 5-10 injections
• Run blank samples to identify contamination
• Use at least 3 biological replicates per group
• Randomize sample injection order

Parameter Optimization

Test parameters on pooled QC sample:

# Test different mass trace detection parameters
mass_errors_ppm = [3.0, 5.0, 10.0]
noise_thresholds = [500.0, 1000.0, 2000.0]

exp.sortSpectra(True)

for mass_err in mass_errors_ppm:
    for noise in noise_thresholds:
        mtd = ms.MassTraceDetection()
        p = mtd.getDefaults()
        p.setValue("mass_error_ppm", mass_err)
        p.setValue("noise_threshold_int", noise)
        mtd.setParameters(p)
        mass_traces = []
        mtd.run(exp, mass_traces, 0)

        epd = ms.ElutionPeakDetection()
        p = epd.getDefaults()
        p.setValue("width_filtering", "fixed")
        epd.setParameters(p)
        mt_split = []
        epd.detectPeaks(mass_traces, mt_split)

        ffm = ms.FeatureFindingMetabo()
        p = ffm.getDefaults()
        p.setValue("isotope_filtering_model", "metabolites (5% RMS)")
        p.setValue("remove_single_traces", "true")
        ffm.setParameters(p)
        features = ms.FeatureMap()
        chrom_out = []
        ffm.run(mt_split, features, chrom_out)

        print(f"mass_error_ppm={mass_err}, noise={noise}: "
              f"{features.size()} features")

Retention Time Windows

Adjust based on chromatographic method:

# For 10-minute LC gradient
params.setValue("distance_RT:max_difference", 30.0)  # 30 seconds

# For 60-minute LC gradient
params.setValue("distance_RT:max_difference", 90.0)  # 90 seconds

Signal Processing

Overview

PyOpenMS provides algorithms for processing raw mass spectrometry data including smoothing, filtering, peak picking, centroiding, normalization, and deconvolution.

Algorithm Pattern

Most signal processing algorithms follow a standard pattern:

import pyopenms as ms

# 1. Create algorithm instance (GaussFilter shown as a concrete example)
algo = ms.GaussFilter()

# 2. Get and modify parameters
params = algo.getParameters()
params.setValue("gaussian_width", 0.2)
algo.setParameters(params)

# 3. Apply to data
algo.filterExperiment(exp)  # or filterSpectrum(spec)

Tip: scripts/process_spectra.py runs a configurable smoothing →

centroiding → normalization → thresholding chain from the command line, so you

rarely need to wire these steps up by hand.

Smoothing

Gaussian Filter

Apply Gaussian smoothing to reduce noise:

# Create Gaussian filter
gaussian = ms.GaussFilter()

# Configure parameters
params = gaussian.getParameters()
params.setValue("gaussian_width", 0.2)  # Width in m/z or RT units
params.setValue("ppm_tolerance", 10.0)  # For m/z dimension
params.setValue("use_ppm_tolerance", "true")
gaussian.setParameters(params)

# Apply to experiment
gaussian.filterExperiment(exp)

# Or apply to single spectrum
spec = exp.getSpectrum(0)
gaussian.filterSpectrum(spec)

Savitzky-Golay Filter

Polynomial smoothing that preserves peak shapes:

# Create Savitzky-Golay filter
sg_filter = ms.SavitzkyGolayFilter()

# Configure parameters
params = sg_filter.getParameters()
params.setValue("frame_length", 11)  # Window size (must be odd)
params.setValue("polynomial_order", 4)  # Polynomial degree
sg_filter.setParameters(params)

# Apply smoothing
sg_filter.filterExperiment(exp)

Peak Picking and Centroiding

Peak Picker High Resolution

Detect peaks in high-resolution data:

# Create peak picker
peak_picker = ms.PeakPickerHiRes()

# Configure parameters
params = peak_picker.getParameters()
params.setValue("signal_to_noise", 3.0)  # S/N threshold
params.setValue("spacing_difference", 1.5)  # Minimum peak spacing
peak_picker.setParameters(params)

# Pick peaks
exp_picked = ms.MSExperiment()
peak_picker.pickExperiment(exp, exp_picked)

Iterative Peak Picker

The CWT-based PeakPickerCWT was removed in modern OpenMS. For data where PeakPickerHiRes struggles (e.g. broader or low-resolution peaks), use PeakPickerIterative, which refits peak widths over several iterations:

# Create iterative peak picker
it_picker = ms.PeakPickerIterative()

# Configure parameters
params = it_picker.getParameters()
params.setValue("signal_to_noise_", 1.0)
params.setValue("peak_width", 0.15)         # expected peak width
params.setValue("nr_iterations_", 5)
it_picker.setParameters(params)

# Pick peaks
exp_picked = ms.MSExperiment()
it_picker.pickExperiment(exp, exp_picked)

Normalization

Normalizer

Normalize peak intensities within spectra:

# Create normalizer
normalizer = ms.Normalizer()

# Configure normalization method
params = normalizer.getParameters()
params.setValue("method", "to_one")  # Options: "to_one", "to_TIC"
normalizer.setParameters(params)

# Apply normalization
normalizer.filterExperiment(exp)

Peak Filtering

Threshold Mower

Remove peaks below intensity threshold:

# Create threshold filter
mower = ms.ThresholdMower()

# Configure threshold
params = mower.getParameters()
params.setValue("threshold", 1000.0)  # Absolute intensity threshold
mower.setParameters(params)

# Apply filter
mower.filterExperiment(exp)

Window Mower

Keep only highest peaks in sliding windows:

# Create window mower
window_mower = ms.WindowMower()

# Configure parameters
params = window_mower.getParameters()
params.setValue("windowsize", 50.0)  # Window size in m/z
params.setValue("peakcount", 2)  # Keep top N peaks per window
window_mower.setParameters(params)

# Apply filter
window_mower.filterExperiment(exp)

N Largest Peaks

Keep only the N most intense peaks:

# Create N largest filter
n_largest = ms.NLargest()

# Configure parameters
params = n_largest.getParameters()
params.setValue("n", 200)  # Keep 200 most intense peaks
n_largest.setParameters(params)

# Apply filter
n_largest.filterExperiment(exp)

Baseline Reduction

Morphological Filter

Remove baseline using morphological operations:

# Create morphological filter
morph_filter = ms.MorphologicalFilter()

# Configure parameters
params = morph_filter.getParameters()
params.setValue("struc_elem_length", 3.0)  # Structuring element size
params.setValue("method", "tophat")  # Method: "tophat", "bothat", "erosion", "dilation"
morph_filter.setParameters(params)

# Apply filter
morph_filter.filterExperiment(exp)

Spectrum Merging

Spectra Merger

Combine multiple spectra into one:

# Create merger
merger = ms.SpectraMerger()

# Configure parameters
params = merger.getParameters()
params.setValue("average_gaussian:spectrum_type", "profile")
params.setValue("average_gaussian:rt_FWHM", 5.0)  # RT window
merger.setParameters(params)

# Merge spectra
merger.mergeSpectraBlockWise(exp)

Deconvolution

Charge Deconvolution

Determine charge states and convert to neutral masses:

# Create feature deconvoluter
deconvoluter = ms.FeatureDeconvolution()

# Configure parameters
params = deconvoluter.getParameters()
params.setValue("charge_min", 1)
params.setValue("charge_max", 4)
params.setValue("potential_charge_states", "1,2,3,4")
deconvoluter.setParameters(params)

# Apply deconvolution. Input is a FeatureMap (not an MSExperiment); the two
# ConsensusMaps receive the charge groups and the connecting edges.
feature_map_out = ms.FeatureMap()
groups = ms.ConsensusMap()
edges = ms.ConsensusMap()
deconvoluter.compute(feature_map, feature_map_out, groups, edges)

Deisotoping a Spectrum

The IsotopeWaveletTransform algorithm was removed. To collapse isotope envelopes in a centroided spectrum to monoisotopic peaks, use the static Deisotoper.deisotopeAndSingleCharge:

spec = exp.getSpectrum(0)
spec.sortByPosition()
# Positional args: spectrum, fragment_tolerance, fragment_unit_ppm, min_charge,
# max_charge, keep_only_deisotoped, min_isopeaks, max_isopeaks,
# make_single_charged, annotate_charge, annotate_iso_peak_count,
# use_decreasing_model, start_intensity_check, add_up_intensity, annotate_features
ms.Deisotoper.deisotopeAndSingleCharge(
    spec, 10.0, True, 1, 3, True, 2, 10, True, True, False, True, 3, False, False
)

Retention Time Alignment

Map Alignment

Align retention times across multiple runs:

# Create map aligner
aligner = ms.MapAlignmentAlgorithmPoseClustering()

# Load multiple experiments
exp1 = ms.MSExperiment()
exp2 = ms.MSExperiment()
ms.MzMLFile().load("run1.mzML", exp1)
ms.MzMLFile().load("run2.mzML", exp2)

# Create reference
reference = ms.MSExperiment()

# Align experiments
transformations = []
aligner.align(exp1, exp2, transformations)

# Apply transformation
transformer = ms.MapAlignmentTransformer()
transformer.transformRetentionTimes(exp2, transformations[0])

Mass Calibration

Internal Calibration

Calibrate mass axis using known reference masses:

# Create internal calibration
calibration = ms.InternalCalibration()

# Set reference masses
reference_masses = [500.0, 1000.0, 1500.0]  # Known m/z values

# Calibrate
calibration.calibrate(exp, reference_masses)

Quality Control

Spectrum Statistics

Calculate quality metrics:

# Get spectrum
spec = exp.getSpectrum(0)

# Calculate statistics
mz, intensity = spec.get_peaks()

# Total ion current
tic = sum(intensity)

# Base peak
base_peak_intensity = max(intensity)
base_peak_mz = mz[intensity.argmax()]

print(f"TIC: {tic}")
print(f"Base peak: {base_peak_mz} m/z at {base_peak_intensity}")

Spectrum Preprocessing Pipeline

Complete Preprocessing Example

import pyopenms as ms

def preprocess_experiment(input_file, output_file):
    """Complete preprocessing pipeline."""

    # Load data
    exp = ms.MSExperiment()
    ms.MzMLFile().load(input_file, exp)

    # 1. Smooth with Gaussian filter
    gaussian = ms.GaussFilter()
    gaussian.filterExperiment(exp)

    # 2. Pick peaks
    picker = ms.PeakPickerHiRes()
    exp_picked = ms.MSExperiment()
    picker.pickExperiment(exp, exp_picked)

    # 3. Normalize intensities
    normalizer = ms.Normalizer()
    params = normalizer.getParameters()
    params.setValue("method", "to_TIC")
    normalizer.setParameters(params)
    normalizer.filterExperiment(exp_picked)

    # 4. Filter low-intensity peaks
    mower = ms.ThresholdMower()
    params = mower.getParameters()
    params.setValue("threshold", 10.0)
    mower.setParameters(params)
    mower.filterExperiment(exp_picked)

    # Save processed data
    ms.MzMLFile().store(output_file, exp_picked)

    return exp_picked

# Run pipeline
exp_processed = preprocess_experiment("raw_data.mzML", "processed_data.mzML")

Best Practices

Parameter Optimization

Test parameters on representative data:

# Try different Gaussian widths
widths = [0.1, 0.2, 0.5]

for width in widths:
    exp_test = ms.MSExperiment()
    ms.MzMLFile().load("test_data.mzML", exp_test)

    gaussian = ms.GaussFilter()
    params = gaussian.getParameters()
    params.setValue("gaussian_width", width)
    gaussian.setParameters(params)
    gaussian.filterExperiment(exp_test)

    # Evaluate quality
    # ... add evaluation code ...

Preserve Original Data

Keep original data for comparison:

# Load original
exp_original = ms.MSExperiment()
ms.MzMLFile().load("data.mzML", exp_original)

# Create copy for processing
exp_processed = ms.MSExperiment(exp_original)

# Process copy
gaussian = ms.GaussFilter()
gaussian.filterExperiment(exp_processed)

# Original remains unchanged

Profile vs Centroid Data

Check data type before processing:

# Check if spectrum is centroided
spec = exp.getSpectrum(0)

if spec.isSorted():
    # Likely centroided
    print("Centroid data")
else:
    # Likely profile
    print("Profile data - apply peak picking")