# Bio Applied Metabolite Identification

> Assign molecular formulas from accurate mass/adducts and match MS/MS spectra by cosine similarity to GNPS/MassBank/HMDB. Use for LC-MS peak annotation, MSI confidence scoring, or KEGG metabolite enrichment (MSEA).

- Skill: `pavel-kravchenko/bio-applied-metabolite-identification` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-applied-metabolite-identification`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-applied-metabolite-identification/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/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.

```python
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.

```python
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.

```python
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

| Level | Requirement | Typical method |
|-------|-------------|-----------------|
| 1 | Exact mass + MS/MS + retention time match to authentic standard | Reference standard run on the same instrument |
| 2 | Spectral library match (mass + MS/MS, no RT confirmation) | GNPS, MassBank, HMDB, NIST |
| 3 | Putative chemical class only | Mass + isotope pattern, no MS/MS match |
| 4 | Uncharacterized feature | Unknown, mass only |

## 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