bio-metabolomics-lipidomics

What this skill does

## Version Compatibility

Reference examples tested with: lipidr 2.16+, pygoslin 2.0+, MS-DIAL 5+

The achievable annotation level is fixed by the acquired evidence, not the software: sn-position and double-bond localization require EAD/OzID/PB/UVPD data that routine CID never produces, and class-resolved quantification requires one isotope-labeled internal standard per class. Verify both before trusting a name or a number.

Before using code patterns, verify installed versions match. If versions differ:
- R: `packageVersion('lipidr')` then `?function_name` to verify parameters
- Python: `pip show pygoslin` then `help(module.function)` to check signatures

If code throws ImportError, AttributeError, or TypeError, introspect the installed
package and adapt the example to match the actual API rather than retrying.

# Lipidomics Analysis

**"Analyze my lipidomics data"** -> Canonicalize names to the resolution level the evidence supports, quantify each class against its own standard, then run class/chain-aware differential and enrichment analysis.
- R: `lipidr::read_skyline()` / `as_lipidomics_experiment()`, `de_analysis()`, `lsea()`
- Nomenclature: `pygoslin` (Python) or `rgoslin` (R) for parsing/canonicalization
- Identification: MS-DIAL 5 (open) or LipidSearch (commercial) upstream

## The Single Most Important Insight -- A Lipid Name Is a Structural-Resolution Claim the Software Usually Overstates

The Liebisch/LIPID MAPS shorthand encodes, in its punctuation, exactly how much structure was measured: `PC 34:1` (space, sum composition) < `PC 16:0_18:1` (underscore, chains known) < `PC 16:0/18:1` (slash, sn-resolved) < `PC 16:0/18:1(9Z)` (double-bond position+geometry). The resolution level is a property of the evidence, not of the string. Tools manufacture overstatement three ways: a formatter that only knows `/`, an in-silico library entry authored at sn-level that a species-level match inherits, and "annotate to the nearest database structure" silently promoting a sum composition to a full structure. sn-position is almost never genuinely measured under CID, so treat every `/` as an unproven `_` until EAD/UVPD/derivatization evidence is in hand. The default rule is: when in doubt, drop a level.

## Structural-Resolution Hierarchy (Separator Semantics)

| Notation | Separator | What was measured | What may NOT be claimed |
|----------|-----------|-------------------|-------------------------|
| `PC 34:1` | space | class + total carbons:double-bonds (accurate mass + isotope + class diagnostic) | the two chains; sn; C=C position |
| `PC 16:0_18:1` | underscore `_` | the two acyl chains (MS/MS acyl losses, RT/ECN-consistent, not an in-source fragment) | which chain is sn-1 vs sn-2 |
| `PC 16:0/18:1` | slash `/` | sn-1/sn-2 assignment (EAD/UVPD/enzymatic - not a CID acyl-loss intensity guess) | C=C position/geometry |
| `PC 16:0/18:1(9Z)` | parentheses | exact double-bond position + cis/trans (OzID/PB/EAD/UVPD) | (full structure) |
| `PC O-34:1` / `PC P-34:1` | `O-` ether / `P-` plasmalogen | ether vs vinyl-ether linkage (diagnostic ion or acid-lability) | a sum composition alone cannot distinguish `P-34:1` from `O-34:2` (vinyl ether = ether + one C=C) |
| `Cer 18:1;O2/16:0` | `;O2` | sphingoid hydroxyl count (old `d18:1`) - measured, not assumed | backbone unsaturation if `d18:1` was a default rather than fragment-confirmed |

Canonicalize every name through Goslin before merging tables or querying LIPID MAPS; never string-match lipid names by hand. Goslin preserves a false `/` faithfully - it is necessary but not sufficient.

## Decision Tree by Question

| Question / situation | Approach | Why |
|----------------------|----------|-----|
| Accurate class-level quantification, high throughput | Shotgun (direct infusion) or HILIC-LC-MS | constant concentration / class bands -> clean ratio to a co-eluting class IS |
| Resolve isobars/isomers, deep low-abundance coverage | RP-LC-MS (± ion mobility) | RT axis adds an identity coordinate; co-elution flags in-source fragments |
| Double-bond position, sn-position, ether/plasmalogen | LC-MS + EAD/OzID/PB/UVPD (± IM) | only these break C=C / glycerol backbone; CID is blind to them |
| Spatial localization | MS-imaging (MS-DIAL 5 spatial mode) | tissue context with predicted-CCS database |
| Need PC acyl chains | negative-mode formate/acetate adduct -> `[M-CH3]-` | `[M+H]+` gives only the m/z 184 head-group ion (class, no chains) |
| Neutral lipids (TG/DG) chains | `[M+NH4]+` adduct | drives neutral-loss-of-fatty-acid fragmentation |
| Suspicious elevated LPC / DG / FA pool | RT co-elution test vs the parent class | an LPC eluting at a PC's RT is an in-source fragment, not biology |
| An apparent odd-chain species (`PC 33:1`) | require MS/MS chain confirmation | usually an in-source fragment or 13C-isotope artifact of an even neighbor |
| Merge names across tools / before a DB lookup | Goslin canonicalization first | abbreviations and separators are tool-specific; hand string-matching corrupts merges |
| Untargeted oxidized-lipid claim | escalate to a targeted, standard-anchored oxylipin panel | untargeted oxidized-lipid IDs are hypotheses; auto-oxidation in the tube fabricates them |

## Load, Normalize, and Run Differential Analysis (lipidr)

**Goal:** Import a quantified lipid table, normalize within class, and find lipids that differ between groups with class/chain-aware output.

**Approach:** Read a Skyline/matrix export into a `LipidomicsExperiment`, attach sample groups, normalize (PQN or class internal standard), then `de_analysis` with an explicit contrast; visualize as a class-faceted volcano.

```r
library(lipidr)

# data_normalized ships with lipidr (PQN-normalized, log2); substitute a real import:
#   d <- read_skyline(list.files(datadir, 'data.csv', full.names = TRUE))
#   d <- add_sample_annotation(d, 'clinical.csv')
#   d <- normalize_pqn(d, measure = 'Area', exclude = 'blank', log = TRUE)
data(data_normalized)

# Contrast references sample-group labels directly; group_col defaults to the first annotation
de_results <- de_analysis(data_normalized, HighFat_water - NormalDiet_water, measure = 'Area')

# logFC.cutoff is on the log2 scale used by limma's topTable inside de_analysis
sig <- significant_molecules(de_results, p.cutoff = 0.05, logFC.cutoff = 1)

plot_results_volcano(de_results, show.labels = FALSE)
```

## Class-Based Internal-Standard Quantification (the non-negotiable)

**Goal:** Convert per-class signal to comparable abundances without baking in class-dependent ionization error.

**Approach:** Ratio each species to a stable-isotope-labeled standard of its OWN class, spiked before extraction so it shares the class's recovery loss; never quantify one class with another class's standard.

```r
# normalize_istd divides each lipid by the internal standard of its matched class.
# Requires one labeled IS per class present in the data (e.g. SPLASH/EquiSPLASH covers ~13 classes).
d_istd <- normalize_istd(data_normalized, measure = 'Area', exclude = 'blank', log = TRUE)

# Class-level summary is only valid WITHIN a class unless per-class response factors were calibrated:
# cross-class molar ratios (e.g. 'PE is 3x PC') carry head-group response bias and are not licensed here.
plot_lipidclass(d_istd, 'sd')
```

## Honest Annotation-Level Assignment (Goslin)

**Goal:** Downgrade any name to the level the evidence supports and verify the claimed level is internally consistent.

**Approach:** Parse with Goslin, read the perceived level, and re-emit at SPECIES (or MOLECULAR_SPECIES) unless sn/C=C evidence exists.

```python
from pygoslin.parser.Parser import LipidParser
from pygoslin.domain.LipidLevel import LipidLevel

parser = LipidParser()
lipid = parser.parse('PC 16:0/18:1')      # a slash-claimed name from a tool export

claimed_level = lipid.lipid.info.level    # LipidLevel enum the string asserts
# Without EAD/UVPD evidence, re-emit at the honest molecular-species level (drops the unproven s