# Metabolomics

> Parse LC-MS mzML with pyOpenMS, PQN/LOESS-normalize feature tables, match m/z to HMDB/GNPS by ppm, run COBRApy FBA. Use when doing metabolomics preprocessing, metabolite ID, feature QC, MSEA enrichment, or flux modeling.

- Skill: `pavel-kravchenko/metabolomics` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/metabolomics`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/metabolomics/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Coding & Dev Tools
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/metabolomics

---


# metabolomics

## When to Use
- Preprocessing raw LC-MS/GC-MS data (mzML) into an aligned feature table: peak picking, RT alignment, gap filling.
- Normalizing and QC-filtering a metabolomics feature table (PQN, IS-based, QC-pool LOESS correction, CV filtering).
- Annotating features by exact mass (HMDB/KEGG) or MS/MS spectral match (GNPS/MassBank), and assigning MSI confidence levels.
- Testing for differential metabolite abundance between conditions and running pathway/MSEA enrichment.
- Modeling flux through a genome-scale metabolic network (FBA, gene essentiality, FVA) with COBRApy.

## Version Compatibility
- pyOpenMS ≥ 3.1, COBRApy ≥ 0.29 (COBRApy dropped `cobra.test`; use `cobra.io.load_model`), scikit-learn ≥ 1.3, RDKit ≥ 2023.09, Python ≥ 3.10.
- XCMS (R/Bioconductor) ≥ 4.0 for peak picking if not using pyOpenMS.

## Prerequisites
- `pip install pyopenms cobra rdkit pandas scikit-learn matplotlib`
- Familiarity with mass spectrometry basics (m/z, retention time, adducts) and basic statistics (t-test, multiple-testing correction).
- A downloaded reference database file for annotation (e.g., an HMDB metabolite CSV with `monoisotopic_mass`, `name`, `formula` columns).

## Quick Reference

| Step | Tool | Notes |
|------|------|-------|
| Peak picking | XCMS / MZmine / pyOpenMS | centWave algorithm for LC-MS |
| Normalization | PQN, IS-based, LOESS | Probabilistic Quotient Normalization |
| Exact mass match | HMDB / KEGG | ± 5 ppm tolerance |
| MS/MS matching | GNPS / MassBank | Cosine similarity threshold > 0.7 |
| Structure prediction | SIRIUS + CSI:FingerID | De novo from MS/MS |
| Statistics | t-test / limma | Log-transform before testing |
| Pathway enrichment | MetaboAnalyst MSEA | ORA or QEA |
| FBA modeling | COBRApy | `model.optimize()` |

**Goal:** Load raw LC-MS data and inspect MS1 spectra.
**Approach:** Use pyOpenMS to parse an mzML file, filter to MS1 scans, and pull the peak list from a spectrum.

```python
from pyopenms import MSExperiment, MzMLFile


def load_ms1_spectra(mzml_path):
    """Load an mzML file and return only MS1-level spectra."""
    exp = MSExperiment()
    MzMLFile().load(mzml_path, exp)
    ms1_spectra = [s for s in exp if s.getMSLevel() == 1]
    return ms1_spectra


spectra = load_ms1_spectra('sample.mzML')
print(f'MS1 spectra: {len(spectra)}')

spec = spectra[0]
mz, intensity = spec.get_peaks()
print(f'RT: {spec.getRT():.2f} s, Peaks: {len(mz)}')
```

**Goal:** Normalize a feature table (rows = samples, columns = features) and visualize sample separation.
**Approach:** Log-transform to stabilize variance, apply Probabilistic Quotient Normalization (PQN) to correct for dilution effects, then PCA to check for batch effects or group separation.

```python
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA


def pqn_normalize(df):
    """Probabilistic Quotient Normalization: scale each sample by its
    median quotient relative to a reference (median) sample."""
    reference = df.median(axis=0)
    quotients = df.div(reference)
    quotient_medians = quotients.median(axis=1)
    return df.div(quotient_medians, axis=0)


features = pd.read_csv('feature_table.csv', index_col=0)   # samples x features
metadata = pd.read_csv('metadata.csv', index_col=0)

features_log = np.log1p(features)
features_norm = pqn_normalize(features_log)

X_scaled = StandardScaler().fit_transform(features_norm)
pca = PCA(n_components=2)
scores = pca.fit_transform(X_scaled)

groups = metadata['condition']
for g in groups.unique():
    idx = (groups == g).values
    plt.scatter(scores[idx, 0], scores[idx, 1], label=g, alpha=0.7)
plt.xlabel(f'PC1 ({pca.explained_variance_ratio_[0]*100:.1f}%)')
plt.ylabel(f'PC2 ({pca.explained_variance_ratio_[1]*100:.1f}%)')
plt.legend(); plt.title('PCA of metabolomics data')
```

**Goal:** Annotate a feature's m/z with candidate HMDB metabolites.
**Approach:** Subtract the adduct mass to get the neutral monoisotopic mass, then filter a reference database within a ppm tolerance window.

```python
import pandas as pd

ADDUCT_MASS = {
    '[M+H]+': 1.00728, '[M+Na]+': 22.98922,
    '[M-H]-': -1.00728, '[M+NH4]+': 18.03437,
}


def match_exact_mass(observed_mz, hmdb_df, adduct='[M+H]+', ppm_tol=5):
    """Match an observed m/z to reference metabolites within ppm tolerance.

    hmdb_df must have a 'monoisotopic_mass' column (neutral mass).
    """
    neutral_mass = observed_mz - ADDUCT_MASS.get(adduct, 0)
    ppm_error = abs((hmdb_df['monoisotopic_mass'] - neutral_mass) / neutral_mass * 1e6)
    matches = hmdb_df[ppm_error <= ppm_tol].copy()
    matches['ppm_error'] = ppm_error[ppm_error <= ppm_tol]
    return matches.sort_values('ppm_error')


hmdb = pd.read_csv('hmdb_metabolites.csv')  # monoisotopic_mass, name, formula
hits = match_exact_mass(180.0634, hmdb, adduct='[M+H]+', ppm_tol=5)
print(hits.head())
```

**Goal:** Run flux balance analysis (FBA) on a genome-scale metabolic model and screen gene essentiality.
**Approach:** Load a curated model with COBRApy, optimize the biomass objective, then knock out each gene one at a time and flag those that collapse growth.

```python
import cobra
from cobra.flux_analysis import flux_variability_analysis


def essential_genes(model, growth_cutoff=0.01):
    """Return gene IDs whose single knockout drops growth below cutoff."""
    essential = []
    for gene in model.genes:
        with model:
            gene.knock_out()
            sol = model.optimize()
            if sol.objective_value < growth_cutoff:
                essential.append(gene.id)
    return essential


model = cobra.io.load_model('e_coli_core')  # bundled BiGG test model

solution = model.optimize()
print(f'Growth rate: {solution.objective_value:.4f} h^-1')
print(solution.fluxes[solution.fluxes.abs() > 0.1].sort_values())

print('Essential genes:', essential_genes(model))

fva = flux_variability_analysis(model, fraction_of_optimum=0.9)
print(fva.sort_values('maximum', ascending=False).head(10))
```

## Differential Abundance in R (limma)

For multi-factor designs or moderated variance estimates on log-transformed feature tables, `limma` is the standard choice (same engine used for microarray/proteomics):

```r
library(limma)

# feature_matrix: features x samples, already log-transformed and normalized
design <- model.matrix(~ condition, data = metadata)
fit <- lmFit(feature_matrix, design)
fit <- eBayes(fit)

# top differential metabolites, BH-corrected
results <- topTable(fit, coef = 2, number = Inf, adjust.method = "BH")
sig <- results[results$adj.P.Val < 0.05, ]
```

## MSI Identification Levels

| Level | Evidence | Example |
|-------|----------|---------|
| 1 | Exact mass + MS/MS + RT vs authentic standard | Glucose confirmed by standard |
| 2 | Spectral match to library (no standard) | GNPS library hit (cosine > 0.7) |
| 3 | Putative: chemical class from MS/MS fragments | "Sphingolipid" from fragment pattern |
| 4 | Uncharacterized | "Feature 1234" |

## Pitfalls
- **In-source fragmentation**: fragments appear as extra MS1 features; check MS2 spectra before annotating as a distinct metabolite.
- **Adduct confusion**: `[M+Na]+` vs `[M+H]+` differ by 21.98 Da; annotate all plausible adducts, don't assume `[M+H]+`.
- **Batch effects**: inject QC pooled samples every ~10 injections; LOESS-correct feature intensities against QC drift.
- **Missing values**: features missing in > 50% of samples are usually noise, not MCAR; use min/2 imputation only for low-abundance MCAR features.
- **Log transform before t-test**: raw MS intensities are approximately log-normal; testing on raw values inflates variance and skews p-values.
- **`cobra.test` is gone**: modern COBRApy loads bundled models via `cobra.io.load_model('e_coli_core')`, not `cobra.test.create_test_model`.

## See Also
- `bio-applied-lc-ms-preprocessing` — deeper XCMS/MZmine peak picking and alignment workflows.
- `bio-applied-metabolite-identification` — expanded annotation and MSEA pathway workflows.
- `bio-applied-metabolic-flux` — extended FBA, gene knockout, and 13C-MFA workflows.
- `bio-core-pathways` — KEGG/Reactome pathway enrichment shared across omics types.