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

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.

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.

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.

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.

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):

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

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.