Bio Applied Lc Ms Preprocessing

LC-MS Metabolomics Data Preprocessing

When to Use

• Converting a batch of raw or centroided mzML files into an aligned features x samples intensity matrix
• Picking chromatographic peaks (centWave) and correcting retention-time drift across injections (obiwarp)
• Normalizing untargeted metabolomics data with pooled QC injections (PQN, LOESS batch correction, CV filtering)
• Grouping adduct/isotope features ([M+H]+, [M+Na]+, [M+NH4]+) that belong to the same metabolite
• Preparing a feature table for downstream PCA, OPLS-DA, or differential abundance testing

Version Compatibility

• R >= 4.3, xcms >= 4.0 (Bioconductor >= 3.18), MSnbase >= 2.28, CAMERA >= 1.58
• Python >= 3.10, pyopenms >= 3.1, pandas >= 2.0, scikit-learn >= 1.4
• Raw data in mzML (centroided); vendor .raw/.d/.wiff must be converted first (ProteoWizard msconvert)

Prerequisites

• BiocManager::install(c("xcms", "CAMERA", "MSnbase")) in R; pip install pyopenms pandas scikit-learn in Python
• Concepts: m/z, retention time (RT), ESI polarity, MS1 vs MS/MS, pooled QC-sample design
• Related skills: bio-applied-metabolite-identification (spectral matching once you have a feature table), bio-applied-proteomics (shared MS instrumentation concepts)

Goal: Turn a directory of mzML files into an aligned, gap-filled feature table (rows = features, columns = samples). Approach: XCMS3 pipeline — readMSData → findChromPeaks (centWave) → adjustRtime (obiwarp) → groupChromPeaks (peak density correspondence) → fillChromPeaks.

library(xcms)
library(MSnbase)

#' Run the standard XCMS preprocessing pipeline on a set of mzML files
#'
#' @param mzml_files character vector of paths to centroided mzML files
#' @param sample_group character vector of group labels, same length/order as mzml_files
#' @return an XcmsExperiment with peaks picked, RT-aligned, and grouped across samples
run_xcms_pipeline <- function(mzml_files, sample_group) {
  pheno <- data.frame(
    sample_name  = sub("\\.mzML$", "", basename(mzml_files)),
    sample_group = sample_group
  )
  raw_data <- readMSData(mzml_files, pdata = new("NAnnotatedDataFrame", pheno), mode = "onDisk")

  # 1. Peak picking: centWave detects chromatographic peaks in profile-mode EIC traces
  cwp <- CentWaveParam(ppm = 15, peakwidth = c(5, 30), snthresh = 10, prefilter = c(3, 1000))
  xdata <- findChromPeaks(raw_data, param = cwp)

  # 2. Retention-time alignment: obiwarp warps each sample's RT axis onto a reference
  xdata <- adjustRtime(xdata, param = ObiwarpParam(binSize = 0.6))

  # 3. Correspondence: group peaks across samples into shared "features" by density in m/z-RT space
  pdp <- PeakDensityParam(sampleGroups = sample_group, bw = 5, minFraction = 0.5)
  xdata <- groupChromPeaks(xdata, param = pdp)

  # 4. Gap filling: integrate raw signal at feature coordinates where no peak was detected
  fillChromPeaks(xdata, param = ChromPeakAreaParam())
}

# xdata <- run_xcms_pipeline(list.files("data/mzML", "\\.mzML$", full.names = TRUE),
#                            rep(c("control", "treatment"), each = 6))
# feature_table <- featureValues(xdata, value = "into", method = "maxint")
# write.csv(cbind(featureDefinitions(xdata)[, c("mzmed", "rtmed")], feature_table), "feature_table.csv")

Goal: Correct systematic intensity drift and drop unreliable features before statistics. Approach: Probabilistic Quotient Normalization (PQN) removes sample-wise dilution effects; CV filtering on pooled QC injections removes features with unstable measurement.

import pandas as pd


def pqn_normalize(feature_table: pd.DataFrame) -> pd.DataFrame:
    """Probabilistic Quotient Normalization for an LC-MS feature table.

    feature_table: samples (rows) x features (columns) of raw intensities.
    Each sample is first total-intensity-normalized, then scaled by the median
    ratio ("quotient") of its features to a reference spectrum (the per-feature
    median across samples) -- corrects dilution without being skewed by a few
    high-intensity outlier features.
    """
    row_sums = feature_table.sum(axis=1)
    integral_norm = feature_table.div(row_sums, axis=0)
    reference_spectrum = integral_norm.median(axis=0)
    quotients = integral_norm.div(reference_spectrum, axis=1)
    median_quotient = quotients.median(axis=1)
    return feature_table.div(median_quotient, axis=0)


def filter_by_qc_cv(feature_table: pd.DataFrame, qc_samples: list, cv_threshold: float = 0.30) -> pd.DataFrame:
    """Drop features whose coefficient of variation in pooled QC samples exceeds cv_threshold.

    Standard metabolomics QC: CV > 30% in repeated QC injections marks a feature
    as too noisy/unstable to trust for biological comparisons.
    """
    qc_data = feature_table.loc[qc_samples]
    cv = qc_data.std(axis=0) / qc_data.mean(axis=0)
    keep_features = cv[cv <= cv_threshold].index
    return feature_table[keep_features]

Goal: Collapse redundant adduct/isotope features into one entry per metabolite before annotation. Approach: Features that co-elute (same RT) and differ by a known adduct mass shift almost always come from the same underlying molecule.

import pandas as pd

# Mass shift of each adduct relative to [M+H]+, in Da
ADDUCT_MASS_SHIFT = {"[M+Na]+": 21.9819, "[M+K]+": 37.9559, "[M+NH4]+": 17.0265}


def group_adducts(features: pd.DataFrame, mz_col: str = "mz", rt_col: str = "rt",
                   rt_tolerance: float = 0.05, ppm_tolerance: float = 10) -> pd.DataFrame:
    """Assign an adduct_group id to co-eluting, mass-shift-related features.

    features: DataFrame indexed by feature id with numeric mz_col/rt_col columns.
    rt_tolerance is in the same units as rt_col (minutes for XCMS 'rtmed').
    Returns a copy of features with an added 'adduct_group' integer column.
    """
    features = features.copy()
    features["adduct_group"] = -1
    assigned, group_id = set(), 0
    for i, row in features.iterrows():
        if i in assigned:
            continue
        features.loc[i, "adduct_group"] = group_id
        assigned.add(i)
        for shift in ADDUCT_MASS_SHIFT.values():
            expected_mz = row[mz_col] + shift
            ppm_window = expected_mz * ppm_tolerance / 1e6
            match = features[
                (features[rt_col].sub(row[rt_col]).abs() <= rt_tolerance)
                & (features[mz_col].sub(expected_mz).abs() <= ppm_window)
            ]
            for j in match.index:
                if j not in assigned:
                    features.loc[j, "adduct_group"] = group_id
                    assigned.add(j)
        group_id += 1
    return features

Pitfalls

• Mass accuracy drift: calibrate/recalibrate m/z before peak picking (ppm in CentWaveParam too tight/loose otherwise causes split or merged peaks and false IDs)
• Adduct confusion: the same metabolite yields [M+H]+, [M+Na]+, [M+K]+ ions — group adducts (see above) before any statistical test, or you will double-count one metabolite as several
• Missing value handling: zeros/NA in a feature table usually mean below-detection-limit, not truly absent — run fillChromPeaks first, then impute remaining gaps with kNN or min/2, never zero-fill
• CentWave parameters are instrument-specific: ppm, peakwidth, and snthresh tuned for Orbitrap data will misfire on QTOF data — check EICs of known standards before running the full batch
• Skipping QC injections: without pooled QC samples interspersed through the run you cannot distinguish batch drift from biology; always include them and use their CV/LOESS trend for correction

See Also

• bio-applied-metabolite-identification — spectral matching, formula assignment, MSEA once you have a feature table
• bio-applied-proteomics — shared LC-MS instrumentation and quantification concepts
• bio-applied-statistics-for-bioinformatics — t-test/limma/FDR for differential abundance
• bio-applied-dimensionality-reduction — PCA/OPLS-DA on the normalized feature table