Metabolomics Filtering Analysis

Summary

Filtered metabolomics data using a two-step process: blank subtraction followed by 80% cumulative signal threshold.

Original Data

Peaks
unique peaks before filtering

After Blank Filter

Peaks
peaks removed

Final (80% Filter)

Peaks
% of original kept
Two-Step Filtering: First, we remove peaks that appear in blanks (contamination). Then we apply 80% cumulative signal filtering to keep only the most significant peaks.

Blank Filtering Results

CategoryCountDescription
Sample OnlyPeaks only in samples, not in blanks (auto-kept)
Passed ValidationPassed both fold-change (≥3x) AND statistical test (p<0.05, FDR-corrected)
ContaminationFailed fold-change or statistical test (removed)
Blank OnlyPeaks only in blanks (not in samples)
No SignalNo detectable signal in any sample or blank

Statistical Validation Details

CriterionPassedFailedDescription
Fold-Change ≥3xSample mean ≥ 3× blank mean (pmp/Bioconductor)
Welch's t-testp < 0.05 (FDR-corrected, one-sided)
Insufficient DataNot enough replicates for t-test (used fold-change only)
Validation Method: Publication-quality blank subtraction using dual criteria: (1) fold-change ≥3x for biological significance (pmp/Bioconductor), and (2) Welch's t-test with Benjamini-Hochberg FDR correction for statistical significance (p < 0.05). For the univariate-testing approach used here, including the multiple-testing correction and the assumptions these tests rest on, see Vinaixa et al. 2012, Metabolites.

Source Data

Per-Tissue Concentration

Average number of peaks each tissue's samples need to reach 80% of total signal — fewer peaks means more concentrated signal.

80% Threshold Results - Per Sample

Each row shows one tissue sample and how many peaks were needed to account for 80% of its total signal (after blank filtering). Samples with fewer peaks needed have more concentrated signal.

How to Read This Table

Example: If a sample needed 500 peaks to reach 80% of its signal, those 500 peaks are the most abundant compounds. The smallest peak kept might contribute 0.02% - anything contributing less was filtered as noise.

Chemical Richness & Diversity

Comparison of metabolite richness (peak counts) and Shannon diversity index across treatments.

Data Source: Blank-filtered dataset ( peaks after 3x blank subtraction).
NOT the 80% cumulative-filtered dataset ( peaks).
Richness and diversity are calculated BEFORE the 80% filter to capture full metabolome complexity.

Chemical Richness (Peak Counts)

Number of detected metabolites per sample (mean ± SE) out of total peaks

Shannon Diversity Index (H)

H = -Σ(p × ln(p)) where p = relative abundance (see Petrén et al. 2024 on chemodiversity metrics)

Interpretation:
Chemical Richness: Higher values indicate more metabolites detected. Differences may reflect stress responses or metabolic shifts.
Shannon Diversity: Higher H indicates more even distribution of metabolite abundances. Lower H suggests dominance by fewer compounds.

Methodology (Exact Calculations)

Chemical Richness
Data source: df_blank_filtered (after 3x blank subtraction)

For each sample column (e.g., "BL - Drought"):
    richness = COUNT of rows WHERE value > 0.0

For each treatment:
    mean = SUM(sample_richness_values) / n_samples
    SE = STDEV(sample_richness_values) / SQRT(n_samples)
Shannon Diversity Index
Data source: df_blank_filtered (after 3x blank subtraction)

For each sample column:
    1. Get all peak abundances where value > 0.0
    2. total = SUM(all abundances)
    3. For each peak: p = abundance / total
    4. H = -SUM(p × ln(p))

For each treatment:
    mean = SUM(H_values) / n_samples
    SE = STDEV(H_values) / SQRT(n_samples)

Why blank-filtered, not 80%-filtered? The 80% cumulative filter removes low-abundance peaks to focus analysis. But richness and diversity metrics should capture the FULL metabolome complexity, so we use the larger blank-filtered dataset.

References:
• Richness, evenness and Shannon diversity as chemodiversity metrics, including their limitations: Petrén et al. 2024, Ecological Monographs — "Understanding the chemodiversity of plants: quantification, variation and ecological function"
• Shannon H and richness are standard formulas, written out above; the Petrén review covers what they do and do not capture (notably that neither accounts for structural disparity between compounds).

Compare Treatments

Toggle treatments to compare. Shows peaks that differ between selected treatments.

Select treatments:
m/z = mass-to-charge ratio. Molecular formulas assigned via MFAssignR (see Methods tab)
RT = retention time (minutes)
Abundance = normalized peak area (unitless, higher = more compound)

Treatment Overlap

Shows how peaks are distributed across treatment groups (Drought, Ambient, Watered). "Unique" means the peak is ONLY found in that treatment, not in the others.

Cross-Tissue Overlap

How peaks are shared between tissues ( unique peaks total).

% of peaks are found in all tissues
How to interpret: Peaks unique to a tissue (leaf-only or root-only) reflect tissue-specific metabolism, and that contrast is the one this dataset resolves: 835 features separate leaf from root after FDR correction.

Treatment-unique peaks are NOT drought biomarkers. With n=3–4 per group, zero features survive Benjamini–Hochberg correction for drought-vs-watered in either tissue, and the contrast does not beat a 300-iteration label-permutation test (leaf p ≈ 0.13, root p ≈ 0.74). A peak appearing in one treatment and not another at this sample size is expected from detection noise alone. Do not report treatment-unique peaks as markers of drought response.

Study limitations you should read first

Instrument class. These data were acquired on a Waters Xevo G2-XS Q-TOF with UPLC, not on an FT-ICR or Orbitrap instrument. The formula-assignment pipeline below (MFAssignR, Kendrick mass defect recalibration, elemental-ratio rules) was designed for ultrahigh-resolution FT-MS. Q-TOF resolution is roughly an order of magnitude lower, so more than one candidate formula can fit a given accurate mass, particularly above m/z 400. Treat formula assignments as provisional rather than confirmed identifications, and note that MS/MS-based annotation is the step that would raise identification confidence.
Positive-mode coverage gap. Data were acquired in positive ESI mode ([M+H]+) only. Sulfur-containing metabolites (glutathione, cysteine) and phosphorylated metabolites (nucleotide mono/di/tri-phosphates, sugar phosphates) ionize poorly in positive mode and are therefore under-represented here. Zero sulfur formulas were assigned; this reflects ionization chemistry and a conservative 34S confirmation rule, not the absence of those compounds from the plant. This is not recoverable by reprocessing. Capturing them requires a negative-mode acquisition. Everything positive mode does capture (the CHO/CHNO majority) is assigned reliably at a mean 1.27 ppm mass error.

Two-Step Filtering Process

We use a two-step filtering approach to ensure data quality: first removing contamination, then keeping only significant peaks.

Step 1: Blank Subtraction (Contamination Removal)

Blanks are samples run through the instrument with no plant material - they capture background contamination from solvents, plastics, and the instrument itself.

How it works:
For each peak found in BOTH samples and blanks:
1. Calculate fold-change: sample_mean / blank_mean
2. Perform Welch's t-test (one-sided: sample > blank)
3. Apply Benjamini-Hochberg FDR correction for multiple testing
4. KEEP only if BOTH criteria pass:
   - Fold-change ≥ 3x (biological significance, per pmp/Bioconductor)
   - FDR-adjusted p-value < 0.05 (statistical significance)

Peaks ONLY in samples (not in blanks) → auto-KEEP

Why dual criteria? The 3x fold-change threshold (pmp/Bioconductor standard) ensures biological relevance (a peak must be meaningfully higher in samples). The statistical test ensures the difference isn't due to random variation. FDR correction accounts for testing thousands of peaks simultaneously.

Statistical Details:
Welch's t-test: Compares sample vs blank means without assuming equal variance
One-sided test: Tests if sample > blank (not just different)
FDR correction: Benjamini-Hochberg method controls false discovery rate at 5% (Benjamini & Hochberg, 1995)

Step 2: 80% Cumulative Signal Threshold

After removing contamination, we filter to keep only the most abundant peaks.

How it works (for each sample separately):
1. Take all peaks and their area values for one sample
2. Sort peaks from LARGEST to SMALLEST
3. Add up the areas as you go down the list
4. Stop when you've added up 80% of the total
5. Everything above that line is kept

Important: A compound is kept if it makes the cut in ANY sample. This ensures we don't lose peaks that are important in specific tissues.

Why Two Steps?

StepPurposeWhat it removes
Blank SubtractionRemove contaminationPlasticizers, solvent impurities, instrument background
80% ThresholdRemove noiseLow-abundance peaks that contribute little to the biological profile
Note on the 80% cutoff: The algorithm keeps adding peaks until the cumulative sum crosses 80%. Two peaks with nearly identical contributions may get different treatment based on where the threshold falls. This is inherent to cumulative thresholds but ensures a consistent reduction in data complexity.

Step 3: Molecular Formula Assignment (MFAssignR)

After filtering, we assign molecular formulas to peaks using MFAssignR (Schum et al., 2020; GitHub), an R package for molecular formula assignment of ultrahigh resolution mass spectrometry data. The pipeline runs in five stages:

Stage 1: Noise Estimation (KMDNoise)

Estimates the noise level using Kendrick Mass Defect analysis. Peaks in a KMD region where no real analyte signals are expected provide a noise floor estimate (KMDN). All subsequent steps use a signal-to-noise threshold of 6 × KMDN.

Stage 2: Isotope Filtering (IsoFiltR)

Separates monoisotopic peaks from their polyisotopic partners (13C, 34S). Parameters: Carbrat = 60, Sulfrat = 30, mass accuracy = 5 ppm for both. This prevents isotope peaks from being assigned their own (incorrect) formulas.

Stage 3: Preliminary CHO Assignment (MFAssignCHO)

Assigns CHO-only formulas first, producing a set of high-confidence reference assignments used to select recalibration series.

Stage 4: Internal Mass Recalibration (Recal)

Corrects systematic mass errors across the spectrum using homologous series from the CHO assignments. The spectrum is divided into 80 Da segments, and a polynomial correction is applied within each segment using up to 6 recalibrant series (auto-selected by quality score: lowest peak score + peak distance closest to 1, ≥4 observed members, >50 Da span). Expansion parameters: step_O = 3, step_H2 = 5, CalPeak = 150.

Stage 5: Full Formula Assignment (MFAssign)

Assigns molecular formulas to the recalibrated peaks. For each peak:

1. Calculate neutral mass: m/z - 1.007276 (remove proton from [M+H]+)
2. Find all formulas (C, H, O, N, S, P) matching within 3 ppm
3. Apply chemical validity rules:
   - H/C ratio: 0.3 to 3.0
   - O/C ratio: 0 to 1.2
   - Double bond equivalents (DBE-O): -13 to 13
   - Nitrogen rule (NMScut = on)
   - Sulfur isotope check (SulfCheck = on)
4. Confirm with isotope patterns (13C, 34S; iso_err = 3 ppm)
5. Select unambiguous best match

Elemental ratio constraints follow Koch et al. (2007) and the Seven Golden Rules of Kind & Fiehn (2007).

Parameters Used

ParameterValueMeaning
Ion ModePositive [M+H]+Compounds detected as protonated molecules
Mass Error3 ppmMax allowed difference between measured and theoretical mass
Isotope Error3 ppmMax error for isotope pattern confirmation
Mass Range100–1000 DaOnly assign formulas to peaks in this range
ElementsC, H, O, N≤5, S≤2, P≤3Allowed elements and maximum heteroatom counts
DeNovo300Above this mass, extend existing series rather than assign de novo
Recal mzRange80 DaWidth of each recalibration segment
S/N Threshold6 × KMDNSignal-to-noise cutoff for all stages

Doubly Charged Peaks (z=2)

Peaks with charge state z=2 are processed through the same pipeline separately. Due to a known limitation in MFAssignR (MFAssign_RMD does not support Zx=2), the final z=2 assignment uses MFAssignCHO, which assigns CHO-only formulas for doubly charged ions.

Understanding PPM Error

PPM (parts per million) measures how close the measured mass is to the theoretical formula mass:

ppm = (measured - theoretical) / theoretical × 1,000,000

Example: m/z 427.3778 vs C26H50O4+H theoretical 427.3782
         ppm = (427.3778 - 427.3782) / 427.3782 × 1,000,000 = -0.9 ppm

Lower ppm = higher confidence. Our assignments average ppm.

Formula Classes

ClassElementsTypical Compounds
CHOC, H, O onlySugars, fatty acids, terpenes
CHNO+ NitrogenAmino acids, alkaloids
CHNOS+ SulfurSulfur-containing amino acids
CHNOP+ PhosphorusPhospholipids, nucleotides
Important Limitation: A molecular formula (e.g., C26H50O4) tells you the atoms present, not the structure. Many different compounds (isomers) can share the same formula. Definitive identification requires MS/MS fragmentation data.
Assignment Results: peaks assigned formulas (% of filtered peaks). Mean mass error: ppm.
References:
• Schum, S.K., Brown, L.E., Mazzoleni, L.R. (2020). MFAssignR: Molecular formula assignment software for ultrahigh resolution mass spectrometry analysis of environmental complex mixtures. Environmental Research, 191, 110114. doi:10.1016/j.envres.2020.110114
• Koch, B.P., Dittmar, T., Witt, M., Kattner, G. (2007). Fundamentals of molecular formula assignment to ultrahigh resolution mass data of natural organic matter. Analytical Chemistry, 79(4), 1758–1763. doi:10.1021/ac061949s
• Kind, T., Fiehn, O. (2007). Seven Golden Rules for heuristic filtering of molecular formulas obtained by accurate mass spectrometry. BMC Bioinformatics, 8, 105. doi:10.1186/1471-2105-8-105

Filtered Metabolomics Data

Two-step filtered data: blank subtraction (3x threshold, pmp/Bioconductor) followed by 80% cumulative signal threshold.

Original Peaks

before any filtering

After Blank Filter

peaks removed

Final Peaks

% of original kept

What Do These Peak Names Mean?

Each compound is identified by a code like 3.90_564.1489n. This encodes two measurements:

PartExampleMeaning
First number3.90Retention time (minutes) - how long it took to pass through the column
Second number564.1489Mass (m/z) - the molecular weight detected
Suffixn or m/zJust notation style
Important: These are NOT identified compounds. We know something with mass 564.1489 eluted at 3.90 minutes, but we don't know what molecule it is yet.

How to Identify Peaks

  1. Database search - Look up the mass in METLIN, HMDB, or MassBank
  2. Run standards - Buy a pure compound and see if it matches
  3. MS/MS fragmentation - Break it apart and look at the pieces
  4. Literature - Check what others found in similar plants

What You Can Do Without Identification