pyopenms-mass-spectrometry

PyOpenMS — 质谱分析

概述

PyOpenMS 为 OpenMS C++ 库提供 Python 绑定，用于计算质谱分析。它支持蛋白质组学和代谢组学数据处理，包括 10 多种 MS 格式的文件读写、信号处理、特征检测、肽段/蛋白质鉴定，以及跨样本定量分析。

适用场景

• 处理用于蛋白质组学或代谢组学的原始 LC-MS/MS 数据（mzML、mzXML）
• 检测色谱特征，并将多个样本中的特征关联起来
• 使用 FDR 控制，根据 MS/MS 搜索引擎结果鉴定肽段和蛋白质
• 运行非靶向代谢组学工作流（峰提取 → 特征检测 → 对齐 → 注释）
• 在质谱文件格式之间转换（mzML、mzXML、featureXML、idXML）
• 对原始谱图数据进行平滑、过滤和质心化
• 对于简单的谱库匹配和代谢物鉴定，请改用 matchms
• 对于蛋白质序列分析（非质谱），请改用 biopython

前置条件

uv pip install pyopenms numpy pandas matplotlib

• Python 3.8+；峰数组操作需要 NumPy
• 输入数据：mzML 文件（标准 MS 格式）、FASTA 数据库（用于鉴定）
• 所有算法都遵循统一模式：algo = Algorithm(); params = algo.getParameters(); params.setValue(...); algo.setParameters(params)

预先访谈

编写任何分析代码前，先与用户确认以下事项。

decisions:
  - id: D1
    param: instrumentResolution
    kind: required
    source: data
    ask: "Is this high-resolution data, and are the spectra already centroided or still in profile mode?"
    default: null

  - id: D2
    param: signalToNoiseThreshold
    kind: required
    source: user
    depends_on: [D1]
    ask: "How far above noise must a peak rise to be kept during centroiding?"
    default: 1.0
    skip_if: "spectra already centroided by the acquisition software"

  - id: D3
    param: massTolerance
    kind: required
    source: data
    depends_on: [D1]
    ask: "What mass accuracy should feature detection and linking assume?"
    default: "10 ppm"

  - id: D4
    param: retentionTimeTolerance
    kind: required
    source: data
    ask: "How far apart in retention time may the same feature appear across runs?"
    default: "100 s"

  - id: D5
    param: chargeStateRange
    kind: required
    source: user
    ask: "Which charge states should features be searched for?"
    default: "1 to 3 - widen for intact protein or metabolite work"

  - id: D6
    param: intensityNormalization
    kind: optional
    source: user
    ask: "Should spectra be normalized before comparison, and to total ion current or to the base peak?"
    default: "not normalized"

  - id: D7
    param: smoothingFilter
    kind: optional_conditional
    source: data
    ask: "Does the signal need smoothing before peak picking, and with which filter width?"
    default: "no smoothing"

D1 会影响后续几乎所有步骤：跳过质心化的 profile 数据会导致每个扫描点都生成一个特征；将高分辨率容差用于低分辨率数据，则会关联并非同一化合物的特征。这两种情况都会生成完整的特征表。

快速入门

import pyopenms as ms

# Load mzML file
exp = ms.MSExperiment()
ms.MzMLFile().load("sample.mzML", exp)
print(f"Spectra: {exp.getNrSpectra()}, Chromatograms: {exp.getNrChromatograms()}")

# Examine first spectrum
spec = exp.getSpectrum(0)
mz, intensity = spec.get_peaks()
print(f"MS level: {spec.getMSLevel()}, RT: {spec.getRT():.2f}s, Peaks: {len(mz)}")

# Quick preprocessing: smooth + centroid
gauss = ms.GaussFilter()
p = gauss.getParameters(); p.setValue("gaussian_width", 0.1); gauss.setParameters(p)
gauss.filterExperiment(exp)

picker = ms.PeakPickerHiRes()
centroided = ms.MSExperiment()
picker.pickExperiment(exp, centroided)
print(f"Centroided spectra: {centroided.getNrSpectra()}")

核心 API

模块 1：文件读写与数据访问

以多种格式读取和写入质谱数据。

import pyopenms as ms

# Read mzML (standard MS format)
exp = ms.MSExperiment()
ms.MzMLFile().load("data.mzML", exp)

# Indexed access for large files (memory-efficient)
loader = ms.IndexedMzMLFileLoader()
indexed_file = ms.OnDiscMSExperiment()
loader.load("large_data.mzML", indexed_file)
spec = indexed_file.getSpectrum(0)  # Load single spectrum on demand
print(f"Total spectra: {indexed_file.getNrSpectra()}")

# Read identification results (idXML)
protein_ids, peptide_ids = [], []
ms.IdXMLFile().load("results.idXML", protein_ids, peptide_ids)
print(f"Peptide IDs: {len(peptide_ids)}, Protein IDs: {len(protein_ids)}")

# Read feature map
fm = ms.FeatureMap()
ms.FeatureXMLFile().load("features.featureXML", fm)
print(f"Features: {fm.size()}")

# Write mzML with compression
exp_out = ms.MSExperiment()
# ... populate experiment ...
ms.MzMLFile().store("output.mzML", exp_out)

# Read FASTA database
entries = []
ms.FASTAFile().load("database.fasta", entries)
print(f"Proteins in DB: {len(entries)}")
for e in entries[:3]:
    print(f"  {e.identifier}: {e.sequence[:30]}...")

支持的格式：mzML、mzXML、mzData（谱图）；featureXML、consensusXML（特征）；idXML、mzIdentML、pepXML（鉴定结果）；TraML（转换）；mzTab（结果）；FASTA（序列）

模块 2：信号处理与峰提取

预处理原始谱图数据，以供后续分析使用。

import pyopenms as ms

exp = ms.MSExperiment()
ms.MzMLFile().load("raw.mzML", exp)

# Gaussian smoothing
gauss = ms.GaussFilter()
p = gauss.getParameters()
p.setValue("gaussian_width", 0.15)  # m/z width
gauss.setParameters(p)
gauss.filterExperiment(exp)

# Savitzky-Golay smoothing (alternative)
sg = ms.SavitzkyGolayFilter()
p = sg.getParameters()
p.setValue("frame_length", 15)  # Must be odd
sg.setParameters(p)
# sg.filterExperiment(exp)  # Use one smoother, not both

# Peak picking (centroiding) — required before feature detection
picker = ms.PeakPickerHiRes()
p = picker.getParameters()
p.setValue("signal_to_noise", 1.0)
picker.setParameters(p)
centroided = ms.MSExperiment()
picker.pickExperiment(exp, centroided)
print(f"Raw peaks in spec 0: {exp.getSpectrum(0).size()}")
print(f"Centroided peaks: {centroided.getSpectrum(0).size()}")

# Normalization
normalizer = ms.Normalizer()
p = normalizer.getParameters()
p.setValue("method", "to_one")  # "to_one" or "to_TIC"
normalizer.setParameters(p)
normalizer.filterPeakMap(centroided)

# Peak filtering — remove low-intensity noise
mower = ms.ThresholdMower()
p = mower.getParameters()
p.setValue("threshold", 100.0)  # Minimum intensity
mower.setParameters(p)
mower.filterPeakMap(centroided)

# Baseline reduction
morph = ms.MorphologicalFilter()
p = morph.getParameters()
p.setValue("struc_elem_length", 3.0)  # m/z window
morph.setParameters(p)
morph.filterExperiment(exp)

模块 3：特征检测与关联

检测色谱特征，并将多个样本中的特征关联起来。

import pyopenms as ms

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

# Feature detection (proteomics — centroided data)
ff = ms.FeatureFinder()
features = ms.FeatureMap()
seeds = ms.FeatureMap()
params = ms.FeatureFinder().getParameters("centroided")
ff.run("centroided", exp, features, params, seeds)
print(f"Detected {features.size()} features")

# Access feature properties
for f in features[:5]:
    print(f"  RT: {f.getRT():.1f}s, m/z: {f.getMZ():.4f}, "
          f"intensity: {f.getIntensity():.0f}, quality: {f.getOverallQuality():.3f}")

# Feature linking across samples — align retention times first
aligner = ms.MapAlignmentAlgorithmPoseClustering()
p = aligner.getParameters()
p.setValue("max_num_peaks_considered", 1000)
aligner.setParameters(p)

# Link features into consensus map
linker = ms.FeatureGroupingAlgorithmQT()
p = linker.getParameters()
p.setValue("distance_RT:max_difference", 60.0)  # seconds
p.setValue("distance_MZ:max_difference", 10.0)   # ppm
linker.setParameters(p)

consensus = ms.ConsensusMap()
linker.group([features_sample1, features_sample2, features_sample3], consensus)
print(f"Consensus features: {consensus.size()}")

# Export to pandas for downstream analysis
import pandas as pd
df = consensus.get_df()
print(f"Consensus table: {df.shape}")

模块 4：肽段与蛋白质鉴定

使用 FDR 控制处理搜索引擎结果并推断蛋白质。

import pyopenms as ms

# Load search engine results
protein_ids, peptide_ids = [], []
ms.IdXMLFile().load("search_results.idXML", protein_ids, peptide_ids)

# Examine peptide hits
for pep_id in peptide_ids[:3]:
    print(f"Spectrum: RT={pep_id.getRT():.1f}, MZ={pep_id.getMZ():.4f}")
    for hit in pep_id.getHits():
        seq = hit.getSequence()
        print(f"  {seq} score={hit.getScore():.4f} charge={hit.getCharge()}")

# FDR filtering (target-decoy approach)
fdr = ms.FalseDiscoveryRate()
fdr.apply(peptide_ids)

# Filter at 1% FDR
filtered = []
for pep_id in peptide_ids:
    hits = [h for h in pep_id.getHits() if h.getScore() <= 0.01]
    if hits:
        pep_id.setHits(hits)
        filtered.append(pep_id)
print(f"Peptide IDs at 1% FDR: {len(filtered)}")

# Protein inference
inference = ms.BasicProteinInferenceAlgorithm()
inference.run(peptide_ids, protein_ids)

for prot_id in protein_ids:
    for hit in prot_id.getHits()[:5]:
        print(f"Protein: {hit.getAccession()}, score: {hit.getScore():.4f}")

# Peptide sequence handling
seq = ms.AASequence.fromString("PEPTIDER")
print(f"Molecular weight: {seq.getMonoWeight():.4f}")
print(f"Formula: {seq.getFormula()}")

# Modified sequence
mod_seq = ms.AASequence.fromString("PEPTM(Oxidation)DER")
print(f"Modified weight: {mod_seq.getMonoWeight():.4f}")

# Enzymatic digestion
digestor = ms.ProteaseDigestion()
digestor.setEnzyme("Trypsin")
digest = []
digestor.digest(ms.AASequence.fromString("MKWVTFISLLLLFSSAYSRGVFRR"), digest)
print(f"Tryptic peptides: {len(digest)}")

模块 5：代谢组学流程

完成从原始数据到特征表的非靶向代谢组学工作流。

import pyopenms as ms

# Step 1: Load and centroid raw data
exp = ms.MSExperiment()
ms.MzMLFile().load("metabolomics_sample.mzML", exp)

picker = ms.PeakPickerHiRes()
centroided = ms.MSExperiment()
picker.pickExperiment(exp, centroided)

# Step 2: Feature detection for metabolomics (small molecules)
ff = ms.FeatureFinder()
features = ms.FeatureMap()
seeds = ms.FeatureMap()
params = ms.FeatureFinder().getParameters("centroided")
params.setValue("isotopic_pattern:charge_low", 1)
params.setValue("isotopic_pattern:charge_high", 3)
ff.run("centroided", centroided, features, params, seeds)
print(f"Detected features: {features.size()}")

# Step 3: Adduct detection (group related adducts)
decharger = ms.MetaboliteAdductDecharger()
p = decharger.getParameters()
p.setValue("potential_adducts", "H:+:0.6;Na:+:0.3;K:+:0.1")  # Positive mode
decharger.setParameters(p)
# decharger.compute(features, feature_map_out, consensus_map_out)

# Step 4: RT alignment across samples
import pyopenms as ms
import pandas as pd

# Assuming feature maps from multiple samples
sample_files = ["sample1.featureXML", "sample2.featureXML", "sample3.featureXML"]
feature_maps = []
for f in sample_files:
    fm = ms.FeatureMap()
    ms.FeatureXMLFile().load(f, fm)
    feature_maps.append(fm)

# Align retention times
aligner = ms.MapAlignmentAlgorithmPoseClustering()
aligner.setReference(0)  # Use first sample as reference

# Step 5: Link features to consensus map
linker = ms.FeatureGroupingAlgorithmQT()
p = linker.getParameters()
p.setValue("distance_RT:max_difference", 30.0)
p.setValue("distance_MZ:max_difference", 10.0)
linker.setParameters(p)

consensus = ms.ConsensusMap()
linker.group(feature_maps, consensus)

# Step 6: Export metabolite table
df = consensus.get_df()
print(f"Feature table: {df.shape[0]} features × {df.shape[1]} columns")
print(df[["RT", "mz", "intensity_0", "intensity_1", "intensity_2"]].head())

核心概念

算法参数模式

所有 PyOpenMS 算法都遵循相同的 3 步模式：

algo = ms.AlgorithmClass()           # 1. Instantiate
params = algo.getParameters()         # 2. Get parameters
params.setValue("param_name", value)   # 3. Set values
algo.setParameters(params)            # 4. Apply
# algo.process(input, output)         # 5. Execute

查看可用参数：

for key in params.keys():
    print(f"{key}: {params.getValue(key)} (type: {params.getDescription(key)})")

核心数据结构

对象	描述	主要方法
MSExperiment	谱图与色谱图的集合	getNrSpectra()、getSpectrum(i)、addSpectrum()
MSSpectrum	单个质谱（m/z、强度）	get_peaks()、getMSLevel()、getRT()、size()
MSChromatogram	色谱迹线	get_peaks()、getNativeID()
Feature	检测到的色谱峰	getRT()、getMZ()、getIntensity()、getOverallQuality()
FeatureMap	特征集合	size()、get_df()、迭代
ConsensusMap	跨样本关联的特征	size()、get_df()
PeptideIdentification	单个谱图的搜索结果	getHits()、getRT()、getMZ()
AASequence	带修饰的氨基酸序列	fromString()、getMonoWeight()、getFormula()

支持的文件格式

格式	读取器类	内容
mzML	MzMLFile	原始谱图（标准）
mzXML	MzXMLFile	原始谱图（旧版）
featureXML	FeatureXMLFile	检测到的特征
consensusXML	ConsensusXMLFile	关联的特征
idXML	IdXMLFile	肽段/蛋白质鉴定结果
mzIdentML	MzIdentMLFile	肽段/蛋白质鉴定结果（标准）
FASTA	FASTAFile	蛋白质序列
TraML	TraMLFile	MRM 转换
mzTab	MzTabFile	定量结果

常见工作流

工作流 1：蛋白质组学特征检测流程

import pyopenms as ms

# Load → Smooth → Centroid → Detect → Export
exp = ms.MSExperiment()
ms.MzMLFile().load("proteomics.mzML", exp)

# Preprocessing
gauss = ms.GaussFilter()
p = gauss.getParameters(); p.setValue("gaussian_width", 0.1); gauss.setParameters(p)
gauss.filterExperiment(exp)

picker = ms.PeakPickerHiRes()
centroided = ms.MSExperiment()
picker.pickExperiment(exp, centroided)

# Feature detection
ff = ms.FeatureFinder()
features = ms.FeatureMap()
params = ff.getParameters("centroided")
ff.run("centroided", centroided, features, params, ms.FeatureMap())

# Export to pandas
df = features.get_df()
print(f"Features: {df.shape}")
df.to_csv("proteomics_features.csv", index=False)

工作流 2：多样本代谢组学比较

1. 加载所有样本的 mzML 文件（核心 API 模块 1）
2. 对每个样本进行质心化处理（核心 API 模块 2 — PeakPickerHiRes）
3. 检测每个样本的特征（核心 API 模块 3 — FeatureFinder）
4. 对齐保留时间（核心 API 模块 3 — MapAlignmentAlgorithmPoseClustering）
5. 将特征关联到共识图（核心 API 模块 3 — FeatureGroupingAlgorithmQT）
6. 将共识表导出到 pandas，并筛选重复样本间 CV 符合条件的条目

工作流 3：鉴定结果分析

import pyopenms as ms
import pandas as pd

# Load and filter identifications
protein_ids, peptide_ids = [], []
ms.IdXMLFile().load("search.idXML", protein_ids, peptide_ids)

# FDR filtering
fdr = ms.FalseDiscoveryRate()
fdr.apply(peptide_ids)

# Build results table
rows = []
for pep_id in peptide_ids:
    for hit in pep_id.getHits():
        if hit.getScore() <= 0.01:  # 1% FDR
            rows.append({
                "sequence": str(hit.getSequence()),
                "score": hit.getScore(),
                "charge": hit.getCharge(),
                "rt": pep_id.getRT(),
                "mz": pep_id.getMZ()
            })

df = pd.DataFrame(rows)
print(f"Identified peptides at 1% FDR: {len(df)}")
print(f"Unique sequences: {df['sequence'].nunique()}")
df.to_csv("peptide_identifications.csv", index=False)

关键参数

参数	模块	默认值	范围/选项	作用
gaussian_width	GaussFilter	0.2	0.05–1.0	m/z 平滑窗口宽度
frame_length	SavitzkyGolayFilter	11	5–31（奇数）	平滑窗口点数
signal_to_noise	PeakPickerHiRes	1.0	0.5–10.0	质心化的最小 S/N
isotopic_pattern:charge_low	FeatureFinder	1	1–6	最小电荷态
isotopic_pattern:charge_high	FeatureFinder	3	1–10	最大电荷态
distance_RT:max_difference	FeatureGroupingAlgorithmQT	100	10–300 s	特征关联的 RT 容差
distance_MZ:max_difference	FeatureGroupingAlgorithmQT	10	1–50 ppm	特征关联的 m/z 容差
threshold	ThresholdMower	0.0	0–10000	最小峰强度
method	Normalizer	"to_one"	"to_one"/"to_TIC"	归一化方法

最佳实践

1. 特征检测前始终进行质心化处理 — FeatureFinder 需要质心化数据。首先使用 PeakPickerHiRes
2. 对大型文件使用索引加载 — OnDiscMSExperiment 按需加载谱图，避免处理大于 1 GB 的文件时出现内存问题
3. 保留原始数据 — 在处理前存储原始实验：orig = ms.MSExperiment(exp)。处理过程具有破坏性
4. 注意 Profile 与质心数据的区别 — 检查数据类型：spec.getType() 返回 1（profile）或 2（质心）。某些算法需要特定的数据类型
5. 发现参数 — 调整参数前先打印参数：for k in params.keys(): print(k, params.getValue(k))
6. 尽早导出到 pandas — 使用 features.get_df() 或 consensus.get_df()，借助 pandas/numpy 进行统计分析

常用示例

示例 1：谱图统计与质量检查

import pyopenms as ms
import numpy as np

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

ms1_count = sum(1 for s in exp if s.getMSLevel() == 1)
ms2_count = sum(1 for s in exp if s.getMSLevel() == 2)
rt_range = (exp.getSpectrum(0).getRT(), exp.getSpectrum(exp.getNrSpectra()-1).getRT())
peak_counts = [s.size() for s in exp]
print(f"MS1: {ms1_count}, MS2: {ms2_count}")
print(f"RT range: {rt_range[0]:.1f}–{rt_range[1]:.1f}s")
print(f"Peaks per spectrum: mean={np.mean(peak_counts):.0f}, "
      f"median={np.median(peak_counts):.0f}")

示例 2：理论谱图生成

import pyopenms as ms

# Generate theoretical fragment spectrum for a peptide
tsg = ms.TheoreticalSpectrumGenerator()
p = tsg.getParameters()
p.setValue("add_b_ions", "true")
p.setValue("add_y_ions", "true")
p.setValue("add_metainfo", "true")
tsg.setParameters(p)

spec = ms.MSSpectrum()
seq = ms.AASequence.fromString("PEPTIDER")
tsg.getSpectrum(spec, seq, 1, 2)  # charge 1 to 2
mz, intensity = spec.get_peaks()
print(f"Theoretical fragments for PEPTIDER: {len(mz)} ions")

示例 3：格式转换（mzXML → mzML）

import pyopenms as ms

exp = ms.MSExperiment()
ms.MzXMLFile().load("old_data.mzXML", exp)
ms.MzMLFile().store("converted.mzML", exp)
print(f"Converted {exp.getNrSpectra()} spectra to mzML format")

故障排除

问题	原因	解决方案
加载时出现 FileNotFound	路径或格式错误	确认文件存在；使用正确的读取器类（例如，读取 .mzML 时使用 MzMLFile）
检测后特征图为空	数据是轮廓数据，而非质心数据	在运行 FeatureFinder 前，先运行 PeakPickerHiRes
检测到的特征很少	S/N 阈值过于严格	将 signal_to_noise 降至 0.5–1.0；调整 isotopic_pattern 的电荷范围
处理大型文件时出现内存错误	一次性加载整个文件	使用 OnDiscMSExperiment 进行索引访问
setValue 类型错误	参数的值类型错误	检查预期类型：params.getDescription(key)。数值使用浮点数，枚举使用字符串
RT 对齐失败	公共特征太少	增大 max_num_peaks_considered；确认样本来自同一实验
FDR 值全为 1.0	搜索结果中没有诱饵命中	确保搜索使用了目标-诱饵数据库；检查分数方向
特征链接过于激进	RT/m/z 容差过大	减小 distance_RT:max_difference 和 distance_MZ:max_difference
处理速度慢	处理了所有谱图	先按 MS 级别过滤：[s for s in exp if s.getMSLevel() == 1]

捆绑资源

• references/data_structures_reference.md — PyOpenMS 核心对象的全面文档（MSExperiment、MSSpectrum、Feature、FeatureMap、ConsensusMap、PeptideIdentification、AASequence、Param）。涵盖：对象创建、属性访问、迭代模式、内存管理、类型转换。

  • 涵盖：全部 13 个以上核心数据结构类，以及创建、访问和迭代示例
  • 已内联：Key Concepts 中的核心数据结构摘要表、Core API Module 4 中的 AASequence 用法
  • 已省略：与 PyOpenMS API 文档重复的详尽属性列表

• references/signal_feature_processing.md — 整合自两个原始参考文档的信号处理算法和特征检测/链接工作流。

  • 涵盖：平滑（Gaussian、Savitzky-Golay）、峰拾取（HiRes、CWT）、归一化、峰过滤（阈值、窗口、N-largest）、基线降低、反卷积、谱图合并、RT 对齐、质量校准、特征查找（质心数据、代谢组学）、特征链接、加合物检测、质量控制
  • 已内联：Core API Module 2 中的 GaussFilter、PeakPickerHiRes、Normalizer、ThresholdMower；Core API Module 3 中的 FeatureFinder、MapAlignment、FeatureGroupingAlgorithmQT
  • 已省略：CWT 峰拾取器的详细参数（相较于 HiRes，使用频率较低）；同位素小波变换（专用场景）

• references/identification_metabolomics.md — 整合自两个原始参考文档的肽/蛋白质鉴定和非靶向代谢组学流水线。

  • 涵盖：搜索引擎集成、FDR 计算、蛋白质推断、酶解、谱图库搜索、理论谱图生成、非靶向代谢组学流水线（峰拾取 → 特征检测 → 加合物检测 → 对齐 → 链接 → 缺失值填充 → 注释 → 导出）、化合物鉴定（基于质量、基于 MS/MS）、归一化（TIC）、质量控制（CV 过滤、空白过滤）、MetaboAnalyst 导出
  • 已内联：Core API Module 4 中的 FDR 过滤、蛋白质推断、AASequence 处理；Core API Module 5 中的代谢组学流水线
  • 已省略：Comet/Mascot 的详细参数配置（因搜索引擎而异）；MetaboAnalyst 导出格式详情（外部工具）

相关技能

• matchms-spectral-matching — 基于 MS/MS 的谱图库匹配和代谢物鉴定
• biopython-molecular-biology — 蛋白质序列分析（FASTA 解析、BLAST）

参考资料

• PyOpenMS 文档：https://pyopenms.readthedocs.io
• OpenMS 项目：https://www.openms.org
• Röst et al. (2016) OpenMS：一种灵活的质谱数据分析开源软件平台。Nature Methods，DOI: 10.1038/nmeth.3959
• GitHub：https://github.com/OpenMS/OpenMS