bio-workflows-metabolomics-pipeline

What this skill does

## Version Compatibility

Reference examples tested with: xcms 4.x+ (MsExperiment/XcmsExperiment), pmp 1.14+, ropls 1.34+, MetaboAnalystR 4.0+

Before using code patterns, verify installed versions match. If versions differ:
- R: `packageVersion('<pkg>')` then `?function_name` to verify parameters

This pipeline is only as honest as its weakest stage: a flawless feature table fed to a too-flexible drift model, or a clean OPLS-DA plot fed to background-free enrichment, produces confident wrong biology. Validate each stage against its own held-out check (QCs, permutation null, assay-coverage background), not against the next stage looking nice.

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

# Metabolomics Pipeline

**"Process my LC-MS metabolomics data end-to-end"** -> Chain xcms feature extraction, QC/normalization, confidence-stratified annotation, validated statistics, and background-aware pathway mapping, treating each stage's output as a hypothesis its component skill scrutinizes.
- R: `readMsExperiment()` -> `findChromPeaks()` -> `groupChromPeaks()` -> `fillChromPeaks()` -> `featureValues()` (xcms), then `QCRSC()`/`pqn_normalisation()` (pmp), `opls()` (ropls), `CalculateOraScore()`/`PerformPSEA()` (MetaboAnalystR)

## What Each Stage Decides and Where the Traps Are

This skill is an orchestrator: it sequences the five component skills and enforces the honest handoffs between them. It does not re-teach each stage's parameters -- those live in the component SKILLs cited per row.

| Stage | The decision it owns | The trap it must not paper over | Defers to |
|---|---|---|---|
| 1. Feature extraction | centWave/grouping/alignment parameters that set the detection floor | A feature table is a parameterized hypothesis; `fillChromPeaks` fabricates intensities; 1 compound = 5-15 features | metabolomics/xcms-preprocessing |
| 2. QC + normalization | drift correction, RSD/D-ratio filtering, dilution normalization, mechanism-aware imputation | Over-correction is invisible to QC RSD; half-min-impute-then-test inflates significance; confounded design is unrescuable | metabolomics/normalization-qc |
| 3. Annotation | the MSI/Schymanski confidence level of every name | A database hit is Level 4-5, not an identification; ambiguous m/z inflates downstream pathways | metabolomics/metabolite-annotation |
| 4. Statistics | univariate FDR + permutation-validated multivariate, reconciled | A clean PLS-DA score plot is the generic output of p>>n; R2 is no evidence; scaling changes conclusions | metabolomics/statistical-analysis |
| 5. Pathway mapping | ORA on IDs vs mummichog on m/z, with an explicit background | The background IS the null; enrichment launders annotation uncertainty into confident biology | metabolomics/pathway-mapping |

## Pipeline Flow

```
raw mzML (centroided)
   |  metabolomics/xcms-preprocessing
   v  readMsExperiment -> findChromPeaks -> adjustRtime -> (re)groupChromPeaks -> fillChromPeaks
features x samples table (+ is_filled flags, mzmed/rtmed)
   |  metabolomics/normalization-qc
   v  blank/detection filter -> within-batch drift (QCRSC) -> RSD/D-ratio filter -> PQN -> mechanism-aware impute
QC-clean, dilution-normalized matrix
   |  split: statistics  AND  annotation (independent axes)
   v
metabolomics/statistical-analysis            metabolomics/metabolite-annotation
permutation-validated hits + univariate FDR  confidence-stratified names (MSI level per feature)
   |                                                |
   +-------------------- join on feature_id --------+
   v  metabolomics/pathway-mapping
identified compounds -> ORA/MSEA   OR   raw m/z (no IDs) -> mummichog/PSEA (background = FULL table)
   v
pathways "consistent with perturbation", conditional on annotations + background
```

Stable-isotope tracing (flux) is a SEPARATE branch off labeled raw data, not a stage of this untargeted flow -- see metabolomics/isotope-tracing.

## Stage 1 -- Feature Extraction (modern xcms 4.x)

**Goal:** Turn centroided mzML into a features-by-samples table, carrying the parameters as part of the result.

**Approach:** Use the `MsExperiment`/`XcmsExperiment` containers with `*Param` objects; align to pooled QC, regroup after alignment, and treat filled values as imputations. Full parameter rationale (ppm, peakwidth, bw, prefilter) lives in metabolomics/xcms-preprocessing.

```r
library(xcms)
# pd: data.frame, one row per file, with a sample_group column ('QC'/'Control'/'Treatment')
raw <- readMsExperiment(spectraFiles = mzml_files, sampleData = pd)

cwp <- CentWaveParam(ppm = 10, peakwidth = c(2, 20), snthresh = 10,
                     prefilter = c(3, 1000), noise = 1000)   # set from instrument; see xcms-preprocessing
xdata <- findChromPeaks(raw, param = cwp)
xdata <- adjustRtime(xdata, param = ObiwarpParam(binSize = 0.6))
pdp <- PeakDensityParam(sampleGroups = sampleData(xdata)$sample_group,
                        bw = 5, minFraction = 0.5, binSize = 0.025)
xdata <- groupChromPeaks(xdata, param = pdp)        # regroup on corrected RT
xdata <- fillChromPeaks(xdata, param = ChromPeakAreaParam())

feat <- featureValues(xdata, value = 'into')        # features x samples; filled cells are imputations
defs <- featureDefinitions(xdata)                   # mzmed / rtmed per feature, for annotation + mummichog
```

## Stage 2 -- QC, Drift, Normalization (not naive median + half-min)

**Goal:** Filter junk features, flatten injection-order drift, normalize per-sample dilution, and impute by mechanism -- before any test sees the data.

**Approach:** Follow the normalization-qc pipeline order: blank/detection filter -> within-batch drift correction (QCRSC) -> RSD/D-ratio filter -> PQN -> mechanism-aware imputation. Do NOT silently half-min-impute and feed limma; validate drift correction on held-out QCs, not on QC clustering.

```r
library(pmp)
# feature_matrix: features in ROWS, samples in COLUMNS (pmp convention); transpose featureValues output
fm <- t(feat)

filtered <- filter_peaks_by_fraction(fm, classes = sample_class, min_frac = 0.5, qc_label = 'QC')
corrected <- QCRSC(df = filtered, order = injection_order, batch = batch_id,
                   classes = sample_class, spar = 0, minQC = 5, qc_label = 'QC')  # CV-selected spline
rsd_filtered <- filter_peaks_by_rsd(corrected, max_rsd = 30, classes = sample_class, qc_label = 'QC')
normalized <- pqn_normalisation(rsd_filtered, classes = sample_class, qc_label = 'QC')
# Impute only the sparse residual holes, by mechanism (QRILC for MNAR / left-censored); see normalization-qc.
```

Drift correction should lower QC RSD AND leave biological-sample RSD unchanged; if biological RSD rises, the spline absorbed signal. Mechanism-aware imputation (QRILC/GSimp for left-censored zeros) replaces the old half-min step, which collapses imputed-subset variance and inflates false significance.

## Stage 3 -- Annotation Before Claiming IDs

**Goal:** Attach an MSI/Schymanski confidence level to each feature so the pathway stage knows what it is allowed to claim.

**Approach:** Match MS/MS to a library (Level 2a) or run SIRIUS/CSI:FingerID (formula Level 4, structure Level 2b/3); a bare m/z is Level 5. Collapse ion families (CAMERA) first so adducts of one compound are not counted as separate metabolites. Mechanics and thresholds live in metabolomics/metabolite-annotation. Annotation and statistics are independent axes -- run them in parallel and join on feature_id.

## Stage 4 -- Statistics (univariate FDR + validated multivariate)

**Goal:** Decide which metabolites genuinely differ, with neither a score plot nor an unadjusted p-value standing alone.

**Approach:** Transform (if heteroscedastic), pick a scaling explicitly (run >=1 alternative and check the conclusion is not scaling-fragile), run a Welch/Mann-Whitney univariate test with BH FDR, AND a permutation-validated OPLS-DA (`permI >= 1000`), then r