bio-metabolomics-xcms-preprocessing

What this skill does

## Version Compatibility

Reference examples tested with: xcms 4.x+ (MsExperiment/XcmsExperiment containers), Spectra 1.12+, CAMERA 1.58+

Before using code patterns, verify installed versions match. If versions differ:
- R: `packageVersion('xcms')` then `?CentWaveParam` to verify parameter names and defaults

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

A feature table is only meaningful alongside its full processing specification: which xcms version, every `*Param` value, and the fill/filter ordering. The table is a parameterized hypothesis about which molecules exist, not the data.

# XCMS Untargeted LC-MS Preprocessing

**"Turn my raw LC-MS files into a feature table"** -> Detect chromatographic peaks per file, align retention times across runs, group corresponding peaks into features, fill gaps, then collapse adduct/isotope redundancy.
- R: `readMsExperiment()` -> `findChromPeaks()` -> `adjustRtime()` -> `groupChromPeaks()` -> `fillChromPeaks()` (xcms)

## The Single Most Important Insight -- The Feature Table Is a Model-Dependent Artifact, Not Ground Truth

Every cell in the table is the output of a detection + grouping + filling model with chosen parameters. Two analysts with different centWave/grouping settings produce materially different tables from identical raw files, so "not detected" is a statement about the parameters, not the sample. Three consequences reorganize the whole workflow: (1) preprocessing parameters silently set the detection floor - a compound absent from results may be present in the raw data but excluded by `noise`/`prefilter`/`peakwidth`/`snthresh`; (2) `fillChromPeaks` integrates whatever signal sits in a feature window even when no peak exists, fabricating a positive number where the honest answer is "below detection"; (3) one compound yields 5-15 features (adducts, isotopologues, in-source fragments, multimers), so a 10,000-feature table is plausibly ~1,000 compounds (Mahieu 2017). Report parameters as part of the result, inspect EICs and alignment of every hit, and collapse redundancy before annotation.

## API Generations -- Use Modern, Not Legacy

| Path | Containers | Verbs | Status |
|------|-----------|-------|--------|
| Modern (xcms 4.x) | `MsExperiment` (raw, Spectra backend) / `XcmsExperiment` (result) | `findChromPeaks` / `adjustRtime` / `groupChromPeaks` / `fillChromPeaks` driven by `*Param` objects | Preferred |
| Legacy (xcms <3) | `xcmsSet` / `xcmsRaw` | `findPeaks` / `group` / `retcor` / `fillPeaks`; `readMSData(mode='onDisk')` | Deprecated - do not use in new code |

Parameters are objects, not loose args: `findChromPeaks(data, param = CentWaveParam(...))`, never `findChromPeaks(data, ppm=..., peakwidth=...)`.

## Decision Tree by Scenario

| Situation | Do | Why |
|-----------|----|----|
| High-res centroid (Orbitrap, Q-Exactive, qTOF) | `CentWaveParam` | Wavelet on real mass traces, no fixed binning |
| Low-res / quadrupole / profile-only | `MatchedFilterParam` | Model-peak on binned EICs tolerates poor resolution |
| Profile data of any kind | Centroid first (msconvert vendor peakPicking, or `Spectra::pickPeaks`) | centWave requires centroids; profile input yields garbage mass traces |
| Many shared, well-behaved peaks across samples | `PeakGroupsParam` (after an initial `groupChromPeaks`) | Loess on universal anchor peaks; gentle and fast |
| Few shared peaks / sparse / strong nonlinear drift | `ObiwarpParam` | Full-profile warping needs no prior peaks |
| Cohort with large case/control compositional differences | `ObiwarpParam`, or `PeakGroupsParam` with `subset =` QC indices | Few universal anchors mis-register the condition-specific metabolome |
| New instrument, no parameter priors | AutoTuner / IPO for a starting neighborhood, then verify against EIC FWHM | Optimizers maximize a surrogate, not biology (McLean 2020) |
| GC-EI data | Deconvolution tools, not xcms peak picking -> metabolomics/msdial-preprocessing | Co-elution + universal fragmentation require component separation first |

## Peak Detection

**Goal:** Detect chromatographic peaks in each centroided file.

**Approach:** Build a `CentWaveParam` with `ppm` and `peakwidth` set from the actual instrument and chromatography (see Quantitative Thresholds), then call `findChromPeaks`.

```r
library(xcms)
# spectraFiles: centroided mzML paths; pd: data.frame with one row per file
raw <- readMsExperiment(spectraFiles = mzml_files, sampleData = pd)

# ppm is across-scan centroid scatter (~2-3x measured error), NOT the spec mass accuracy.
# peakwidth is c(min, max) in SECONDS, measured from EIC base-widths of known peaks.
cwp <- CentWaveParam(ppm = 10, peakwidth = c(2, 20), snthresh = 10,
                     prefilter = c(3, 1000), noise = 1000, mzdiff = -0.001,
                     integrate = 1L, mzCenterFun = 'wMean')
xdata <- findChromPeaks(raw, param = cwp)
nrow(chromPeaks(xdata))
```

## Retention-Time Alignment

**Goal:** Remove cross-run RT drift so the same compound lands at the same RT in every sample.

**Approach:** Choose obiwarp (no prior peaks) or peakGroups (anchor-based); align to a pooled QC, never to file #1. Regroup afterward because RTs changed.

```r
# obiwarp: full-profile warping. binSize here is the m/z profile bin (default 1),
# distinct from PeakDensityParam$binSize and MatchedFilterParam$binSize.
xdata <- adjustRtime(xdata, param = ObiwarpParam(binSize = 0.6))

# peakGroups alternative needs an initial correspondence and good universal anchors:
# xdata <- groupChromPeaks(xdata, param = pdp_anchor)
# xdata <- adjustRtime(xdata, param = PeakGroupsParam(minFraction = 0.85, span = 0.4,
#     subset = which(sampleData(xdata)$sample_type == 'QC'), subsetAdjust = 'average'))
plotAdjustedRtime(xdata)
```

## Correspondence (Grouping)

**Goal:** Match peaks across samples into consensus features.

**Approach:** Peak-density grouping in m/z slices; `bw` is the dominant knob and must reflect residual post-alignment RT scatter, not raw peak width.

```r
pdp <- PeakDensityParam(sampleGroups = sampleData(xdata)$sample_group,
                        bw = 5, minFraction = 0.5, minSamples = 1, binSize = 0.025)
xdata <- groupChromPeaks(xdata, param = pdp)
nrow(featureDefinitions(xdata))
```

## Gap-Filling

**Goal:** Integrate signal for features missing a detected peak in some samples.

**Approach:** `fillChromPeaks` with `ChromPeakAreaParam`; treat filled values as imputations, not measurements.

```r
xdata <- fillChromPeaks(xdata, param = ChromPeakAreaParam())
filled <- chromPeakData(xdata)$is_filled   # logical flag; lives in chromPeakData, not chromPeaks
feat <- featureValues(xdata, value = 'into')        # features x samples matrix
defs <- featureDefinitions(xdata)                   # mzmed / rtmed / npeaks per feature
```

## Redundancy Collapse

**Goal:** Group the same compound's adducts/isotopes/fragments back toward compound spectra before annotation.

**Approach:** CAMERA in order groupFWHM -> groupCorr -> findIsotopes -> findAdducts (isotopes before adducts). Correlation grouping needs enough samples to be meaningful and can over- or under-merge - verify against the table size.

```r
library(CAMERA)
xsa <- xsAnnotate(as(xdata, 'xcmsSet'))
xsa <- groupFWHM(xsa, perfwhm = 0.6)
xsa <- groupCorr(xsa)
xsa <- findIsotopes(xsa, mzabs = 0.01, ppm = 10)
xsa <- findAdducts(xsa, polarity = 'positive')
peaklist <- getPeaklist(xsa)
```

## QC Feature Filtering (Preprocessing/QC Bridge)

**Goal:** Drop features that fail conventional QC, operationalizing Broadhurst 2018 inside the xcms object.

**Approach:** `filterFeatures` with `RsdFilter` (CV in QCs), `DratioFilter` (sd_QC/sd_sample), `PercentMissingFilter`, `BlankFlag`. Drift correction and the full QC pipeline live in metabolomics/normalization-qc.

```r
qc <- sampleData(xdata)$sample_group == 'QC'
study <- sampleData(xdata)$sample_group %in% c('Control', 'Trea