Bio Applied Metabolite Identification

Metabolite Identification and Annotation

When to Use

• Assigning a molecular formula to an unknown LC-MS feature from its accurate m/z and adduct type
• Matching an experimental MS/MS spectrum against GNPS, MassBank, or HMDB reference libraries (cosine/modified-cosine scoring)
• Deciding what confidence level (MSI 1-4) an annotation deserves before reporting it
• Running differential abundance testing and volcano plots on an annotated metabolite table
• Running metabolite set enrichment analysis (MSEA/ORA) against KEGG pathways on a hit list

Version Compatibility

• Python >= 3.10, RDKit >= 2023.09, matchms >= 0.24, pyteomics >= 4.7
• statsmodels >= 0.14, scipy >= 1.11 (differential abundance, FDR)
• Reference libraries: GNPS (public spectral library, MSP/mgf export), MassBank of North America, HMDB 5.0

Prerequisites

• pip install rdkit matchms pyteomics statsmodels scipy numpy pandas
• Familiarity with centroided MS/MS spectra (m/z, intensity arrays) — see bio-applied-lc-ms-preprocessing for feature detection upstream of this skill
• A spectral library file in MSP or mzML/mgf format (e.g. downloaded from GNPS) for the matching step

Molecular Formula Assignment from Accurate Mass

Goal: given an observed m/z and an adduct, find candidate molecular formulas within a ppm tolerance and rank them by mass error and RDBE (rings + double bond equivalents).

Approach: enumerate candidate formulas with pyteomics' mass calculator, compute the adduct-corrected theoretical mass, filter by ppm error, then use RDKit only when you have a candidate SMILES/structure to confirm exact mass and RDBE.

from itertools import product
from pyteomics.mass import calculate_mass, isotopologues, most_probable_isotopologue

# Masses of the electron and common adduct ions (Da)
ADDUCTS = {
    "[M+H]+": 1.007276,
    "[M+Na]+": 22.989218,
    "[M+NH4]+": 18.033823,
    "[M-H]-": -1.007276,
}

def ppm_error(observed_mz: float, theoretical_mz: float) -> float:
    """Mass error in ppm between an observed and theoretical m/z."""
    return (observed_mz - theoretical_mz) / theoretical_mz * 1e6

def candidate_formulas(observed_mz: float, adduct: str, formulas: list[str],
                        tolerance_ppm: float = 5.0) -> list[dict]:
    """Score a list of candidate molecular formulas (e.g. ['C6H12O6', 'C9H8O2'])
    against an observed adduct m/z and keep only those within tolerance_ppm.
    """
    adduct_mass = ADDUCTS[adduct]
    hits = []
    for formula in formulas:
        neutral_mass = calculate_mass(formula=formula)
        theoretical_mz = neutral_mass + adduct_mass
        error = ppm_error(observed_mz, theoretical_mz)
        if abs(error) <= tolerance_ppm:
            hits.append({"formula": formula, "theoretical_mz": theoretical_mz, "ppm_error": error})
    return sorted(hits, key=lambda h: abs(h["ppm_error"]))

# Example: glucose [M+H]+ observed at m/z 181.0703
hits = candidate_formulas(181.0703, "[M+H]+", ["C6H12O6", "C7H12O5", "C5H8O7"])
print(hits)

MS/MS Spectral Library Matching (Cosine Similarity)

Goal: score an experimental MS/MS spectrum against a library of reference spectra to get putative identifications (MSI level 2).

Approach: use matchms to load spectra, apply standard preprocessing (normalize intensities, remove low-intensity noise), then compute cosine or modified-cosine similarity — the modified cosine also matches fragments shifted by the precursor mass difference, which recovers analogues/derivatives.

import numpy as np
from matchms import Spectrum
from matchms.similarity import CosineGreedy, ModifiedCosine
from matchms.filtering import normalize_intensities, select_by_relative_intensity

def build_spectrum(mz: np.ndarray, intensities: np.ndarray, precursor_mz: float) -> Spectrum:
    """Wrap raw MS/MS peak lists into a matchms Spectrum with basic cleanup."""
    spectrum = Spectrum(mz=mz, intensities=intensities,
                         metadata={"precursor_mz": precursor_mz})
    spectrum = normalize_intensities(spectrum)
    spectrum = select_by_relative_intensity(spectrum, intensity_from=0.01)  # drop <1% noise peaks
    return spectrum

def match_against_library(query: Spectrum, library: list[Spectrum],
                           tolerance: float = 0.005, min_matches: int = 3) -> list[tuple]:
    """Return (library_index, score, n_matched_peaks) sorted by score, using
    modified-cosine so precursor mass shifts (analogues) are still matched.
    """
    similarity = ModifiedCosine(tolerance=tolerance)
    results = []
    for i, ref in enumerate(library):
        score, n_matched = similarity.pair(query, ref)[()]
        if n_matched >= min_matches:
            results.append((i, float(score), int(n_matched)))
    return sorted(results, key=lambda r: r[1], reverse=True)

# query = build_spectrum(mz_array, intensity_array, precursor_mz=181.0703)
# best_hits = match_against_library(query, reference_library_spectra)

Note: CosineGreedy is the plain cosine (peak m/z must match directly); ModifiedCosine additionally allows a constant offset equal to the precursor mass difference — prefer it when searching for structural analogues, not exact standards.

Differential Abundance and MSEA

Goal: test which annotated metabolites differ between two groups, then check whether the significant hits are enriched in KEGG pathways.

Approach: Mann-Whitney U per metabolite with Benjamini-Hochberg FDR, then a hypergeometric over-representation test (ORA) against a KEGG pathway-to-compound mapping.

import numpy as np
import pandas as pd
from scipy.stats import mannwhitneyu, hypergeom
from statsmodels.stats.multitest import multipletests

def differential_abundance(intensity_df: pd.DataFrame, group_a: list[str],
                            group_b: list[str]) -> pd.DataFrame:
    """Mann-Whitney U test per metabolite (rows) between two sample groups (columns),
    with log2 fold change and BH-FDR-corrected p-values.
    """
    records = []
    for metabolite, row in intensity_df.iterrows():
        a, b = row[group_a].values, row[group_b].values
        stat, pval = mannwhitneyu(a, b, alternative="two-sided")
        log2fc = np.log2((a.mean() + 1e-9) / (b.mean() + 1e-9))
        records.append({"metabolite": metabolite, "log2fc": log2fc, "pvalue": pval})
    result = pd.DataFrame(records)
    result["padj"] = multipletests(result["pvalue"], method="fdr_bh")[1]
    return result.sort_values("padj")

def msea_ora(hit_metabolites: set[str], pathway_to_compounds: dict[str, set[str]],
             background_size: int) -> pd.DataFrame:
    """Hypergeometric over-representation test: for each KEGG pathway, is the
    overlap with hit_metabolites bigger than expected by chance?
    """
    rows = []
    for pathway, compounds in pathway_to_compounds.items():
        overlap = hit_metabolites & compounds
        if not overlap:
            continue
        pval = hypergeom.sf(len(overlap) - 1, background_size, len(compounds), len(hit_metabolites))
        rows.append({"pathway": pathway, "n_hits": len(overlap), "n_pathway": len(compounds), "pvalue": pval})
    result = pd.DataFrame(rows)
    result["padj"] = multipletests(result["pvalue"], method="fdr_bh")[1]
    return result.sort_values("padj")

MSI Identification Confidence Levels

Pitfalls

• Mass accuracy thresholds: always compare in ppm, not absolute Da — 5 ppm at m/z 500 is only 0.0025 Da, so a fixed Da cutoff over-tolerates high-mass ions and under-tolerates low-mass ones.
• Schymanski/MSI confidence levels: level 1 requires a co-injected authentic standard; most annotations from library search alone are level 2-3 at best — do not report level-1 confidence from a database hit.
• Ion mode bias: ESI+ and ESI- ionize different metabolite classes (e.g. amines favor positive mode, carboxylic acids favor negative) — run and annotate both modes rather than assuming one is sufficient.
• Adduct ambiguity: the same neutral mass can fit several adduct/formula combinations at a given m/z — always check multiple adducts, not just [M+H]+, before accepting a formula.
• Cosine score inflation with few peaks: a high cosine score from only 2-3 matched fragments is not reliable — require a minimum matched-peak count (e.g. >=3) alongside the score threshold.

See Also

• bio-applied-lc-ms-preprocessing — feature detection and alignment upstream of identification
• bio-applied-metabolic-flux — using identified/quantified metabolites in flux models
• bio-applied-rdkit-basics — molecular structure handling for candidate confirmation
• bio-applied-statistics-for-bioinformatics — general hypothesis testing and FDR background