metabolomics/xcms-preprocessing/SKILL.md
Programmatic untargeted LC-MS feature extraction in R with the modern xcms 4.x MsExperiment/XcmsExperiment API, taking raw mzML to a feature table via CentWave peak detection, retention-time alignment, peak-density correspondence, gap-filling, CAMERA redundancy collapse, and built-in QC feature filtering. Use when converting centroided LC-MS runs into a features-by-samples matrix and deciding centWave/grouping/alignment parameters. For drift correction and QC/CV filtering execution see metabolomics/normalization-qc; for metabolite identification see metabolomics/metabolite-annotation; for the MS-DIAL GUI alternative with MS2Dec deconvolution see metabolomics/msdial-preprocessing; for downstream statistics see metabolomics/statistical-analysis.
npx skillsauth add GPTomics/bioSkills bio-metabolomics-xcms-preprocessingInstall this skill globally with one command. Works with Claude Code, Cursor, and Windsurf.
3 of 9 scanners reported clean
Some scanners were skipped, did not run, or reported a non-clean status. Review each row below.
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:
packageVersion('xcms') then ?CentWaveParam to verify parameter names and defaultsIf 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.
"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.
readMsExperiment() -> findChromPeaks() -> adjustRtime() -> groupChromPeaks() -> fillChromPeaks() (xcms)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.
| 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=...).
| 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 |
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.
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))
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.
# 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)
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.
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))
Goal: Integrate signal for features missing a detected peak in some samples.
Approach: fillChromPeaks with ChromPeakAreaParam; treat filled values as imputations, not measurements.
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
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.
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)
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.
qc <- sampleData(xdata)$sample_group == 'QC'
study <- sampleData(xdata)$sample_group %in% c('Control', 'Treatment')
xdata <- filterFeatures(xdata, filter = RsdFilter(threshold = 0.3, qcIndex = qc))
xdata <- filterFeatures(xdata, filter = DratioFilter(threshold = 0.5, qcIndex = qc, studyIndex = study))
ppm = 3 because the Orbitrap datasheet says 3 ppm.ppm is across-scan centroid scatter, which exceeds time-averaged mass accuracy; too tight fragments one ion into short ROIs that each fail prefilter.ppm to ~2-3x the empirical per-scan centroid scatter, not the datasheet number.c(20, 50) onto modern UHPLC.peakwidth ~ c(0.5x min, 2x max).snthresh while prefilter[I] already kills the trace.I=100 may both under-filter noise and kill trace metabolites.snthresh.prefilter[I] first for trace work; the lowest gate dominates.bw = 30 onto UHPLC, or choosing bw independently of alignment quality.bw=30 merges chromatographically resolved co-eluting compounds; with poor alignment a tight bw instead splits one compound across features.bw from residual post-alignment RT scatter (often 2-6 s on UHPLC); inspect EICs of merged/split features.is_filled; report per-feature filled fraction; for inference use unfilled values with MNAR-aware imputation (QRILC/GSimp), reserving the fill for dense exploratory PCA.| Threshold | Source | Rationale |
|-----------|--------|-----------|
| ppm Orbitrap/Q-Exactive 5-10, qTOF 15-30 | Tautenhahn 2008; instrument physics | ~2-3x measured across-scan centroid scatter, not spec accuracy |
| peakwidth UHPLC c(2,20), HPLC c(10,40), HILIC c(10,60) (s) | Smith 2006; chromatography | Must bracket measured EIC base-widths; default c(20,50) wrong for UHPLC |
| Points across peak >= ~6-7 | Zeng 2023 JASMS 34:1136 | Below this, peak-area precision degrades non-linearly; an acquisition limit no parameter recovers |
| prefilter = c(3, I) | Tautenhahn 2008 | Min 3 consecutive scans above intensity I; I set per instrument baseline |
| Grouping bw 2-6 s (UHPLC) | xcms vignette | Default 30 s merges resolved co-eluting compounds on fast chromatography |
| QC CV (RSD) < 0.20-0.30 | Broadhurst 2018 Metabolomics 14:72 | Features with high QC variance are unreliable |
| D-ratio < 0.5 | Broadhurst 2018 | Technical variance must sit well below biological |
| Blank flag k ~ 3-5 | Broadhurst 2018 | Test-sample mean must exceed k x blank mean |
| ~1 compound per 5-15 features | Mahieu 2017 Anal Chem 89:10397 | ~90% of detected features are adduct/isotope/fragment degeneracy |
| Error / symptom | Cause | Solution |
|-----------------|-------|----------|
| could not find function "readMSData" or legacy verbs missing | Using deprecated xcmsSet/readMSData API on xcms 4.x | Use readMsExperiment() + the findChromPeaks/groupChromPeaks verbs |
| unused argument (ppm = ...) in findChromPeaks | Passing loose args instead of a *Param object | Wrap in CentWaveParam(...) and pass via param = |
| Features defined on uncorrected RT | Skipped the regroup after adjustRtime | Call groupChromPeaks again after alignment |
| Garbage mass traces, almost no peaks | Profile (non-centroid) data fed to centWave | Centroid first (msconvert vendor peakPicking or Spectra::pickPeaks) |
| sampleGroups length/semantics error | Vector misaligned with sample order or missing | Pass sampleData(xdata)$group matching file order; it is mandatory |
| Three different binSize defaults confused | obiwarp (m/z, default 1) vs PeakDensity (m/z, 0.25) vs matchedFilter (m/z, 0.1) | Set each in its own *Param; they are not the same knob |
| as(xdata, 'xcmsSet') fails or warns | CAMERA expects the legacy container | Coerce the XcmsExperiment to xcmsSet only for CAMERA; keep modern objects upstream |
tools
End-to-end CLIP-seq pipeline from FASTQ to ENCODE-compliant binding sites, single-nucleotide crosslink maps, annotation, motifs, and (optionally) differential binding. Use when running the full Yeo lab eCLIP / iCLIP / iCLIP2 / iCLIP3 / irCLIP / PAR-CLIP analysis with SMInput control, protocol-specific UMI extraction, ENCODE STAR parameters, CLIPper or Skipper peak calling with stringent log2 FC and -log10 p thresholds, IDR rescue and self-consistency QC, and downstream motif registration with mCross or PEKA.
development
Detect, date, and contextualize whole-genome duplication (WGD / paleopolyploidy) events using wgd v2 (Chen et al 2024), KsRates (Sensalari 2022 substitution-rate-corrected Ks dating), DupGen_finder (Qiao 2019), MAPS (Li 2018 phylogenomic), POInT (Conant 2008 ordered-block), SLEDGe (2024 ML-based), Whale.jl (Bayesian DL+WGD), and synteny-anchored paranome construction. Use when identifying ancient polyploidy from Ks distributions and synteny block analysis, positioning WGD events relative to speciation, distinguishing tandem from segmental from WGD duplications, dating the 2R/3R vertebrate / fish / salmonid WGDs, building paranome and Ks-age mixture models, applying KsRates substitution-rate correction across lineages, or testing alternative biased-fractionation / dosage-balance models post-WGD.
tools
Build whole-genome alignments using Progressive Cactus (Armstrong 2020 reference-free clade-level WGA), Minigraph-Cactus (Hickey 2024 pangenome-aware), LASTZ chain/net (UCSC pipeline), MUMmer4 (Marçais 2018 pairwise), minimap2 -x asm5/10/20 (Li 2018 fast pairwise), AnchorWave (Song 2022 WGD-aware), and Mauve / progressiveMauve (bacterial). Operates the HAL toolkit (Hickey 2013) for downstream extraction including halSynteny, halLiftover, halBranchMutations, and hal2maf. Use when constructing multi-species alignments for comparative-annotation projection (TOGA), synteny detection, conservation analyses (phyloP / PhastCons), or pangenome graph construction; selecting between reference-free (Cactus) and reference-anchored (LASTZ chains/nets) approaches; tuning sensitivity for closely vs distantly related genomes; or producing HAL files for genome-wide downstream tools.
development
Detect syntenic blocks and structural rearrangements between genomes using MCScanX (Wang 2012), JCVI/MCScan (Tang 2008 Python), GENESPACE (Lovell 2022) for orthology-anchored riparian visualization, SyRI for structural variation, AnchorWave for sequence-level synteny, i-ADHoRe 3.0 for highly diverged species, SynNet for synteny networks, and ntSynt for multi-genome macrosynteny. Use when identifying collinear gene blocks across species, distinguishing macrosynteny from microsynteny, detecting inversions/translocations/duplications, anchoring orthology in WGD lineages, producing publication riparian plots, computing synteny block age via Ks (cross-references whole-genome-duplication), or running synteny-aware ortholog inference in polyploids.