GC-MS Data Processing

R Tools (XCMS for GC-MS)

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install(c("xcms", "MSnbase"))

Verify Installation

import matchms
import pyopenms as ms
print(f"matchms: {matchms.__version__}")
print(f"pyopenms: {ms.__version__}")

Core Capabilities

1. Loading GC-MS Data

Read raw GC-MS files (mzML, mzXML, CDF) and inspect TIC vs EIC.

import pyopenms as ms
import numpy as np

exp = ms.MSExperiment()
ms.MzMLFile().load("gcms_sample.mzML", exp)

# Build total ion chromatogram (TIC)
retention_times = []
tic_intensities = []
for spectrum in exp:
    if spectrum.getMSLevel() == 1:
        retention_times.append(spectrum.getRT())
        mz, intensity = spectrum.get_peaks()
        tic_intensities.append(np.sum(intensity))

print(f"Spectra: {len(retention_times)}")
print(f"RT range: {min(retention_times):.1f} - {max(retention_times):.1f} s")

Extract an extracted ion chromatogram (EIC) for a target m/z:

def extract_eic(experiment, target_mz, tolerance=0.5):
    rts = []
    intensities = []
    for spectrum in experiment:
        if spectrum.getMSLevel() == 1:
            rts.append(spectrum.getRT())
            mz_array, int_array = spectrum.get_peaks()
            mask = np.abs(mz_array - target_mz) <= tolerance
            intensities.append(np.sum(int_array[mask]) if np.any(mask) else 0.0)
    return np.array(rts), np.array(intensities)

# EIC for m/z 73 (common TMS fragment)
rt_arr, eic_arr = extract_eic(exp, target_mz=73.0, tolerance=0.5)

2. Peak Detection with XCMS (R, GC-MS Optimized)

CentWave parameters tuned for GC-MS: narrower peak widths, wider ppm tolerance for unit-mass-resolution instruments.

library(xcms)
library(MSnbase)

raw_files <- list.files("gcms_data", pattern = "\\.mzML$", full.names = TRUE)
raw_data <- readMSData(raw_files, mode = "onDisk")

# CentWave parameters optimized for GC-MS
# - peakwidth: GC peaks are typically 2-10 seconds (narrower than LC)
# - ppm: unit-mass GC-MS needs wider tolerance (~30-40 ppm);
#         for HR-GC-MS (e.g., GC-Orbitrap), use 5-10 ppm instead
# - noise: adjust to instrument baseline
cwp <- CentWaveParam(
    peakwidth = c(2, 10),
    ppm = 40,
    snthresh = 5,
    prefilter = c(3, 500),
    mzdiff = 0.5,
    noise = 500,
    integrate = 2
)

xdata <- findChromPeaks(raw_data, param = cwp)
cat("Peaks detected:", nrow(chromPeaks(xdata)), "\n")

3. Retention Index Calculation

Retention indices normalize retention times against an alkane ladder, enabling cross-instrument comparison.

Kovats Retention Index (isothermal GC, uses logarithmic interpolation):

import numpy as np

def kovats_ri(rt_compound, alkane_rts, alkane_carbons):
    """Calculate Kovats retention index for isothermal GC.

    Uses the original Kovats logarithmic formula:
    RI = 100 * (n + (log(tR(x)) - log(tR(n))) / (log(tR(n+1)) - log(tR(n))))

    Args:
        rt_compound: retention time of the target compound (seconds)
        alkane_rts: list of alkane retention times (seconds), sorted ascending
        alkane_carbons: list of carbon numbers matching alkane_rts
    """
    alkane_rts = np.array(alkane_rts, dtype=np.float64)
    alkane_carbons = np.array(alkane_carbons)

    idx = np.searchsorted(alkane_rts, rt_compound) - 1
    idx = max(0, min(idx, len(alkane_rts) - 2))

    n_z = alkane_carbons[idx]
    rt_z = alkane_rts[idx]
    rt_z1 = alkane_rts[idx + 1]

    ri = 100 * (n_z + (np.log(rt_compound) - np.log(rt_z)) / (np.log(rt_z1) - np.log(rt_z)))
    return ri

Linear Retention Index (LTPRI) (temperature-programmed GC, most common in metabolomics):

def linear_ri(rt_compound, alkane_rts, alkane_carbons):
    """Calculate linear (temperature-programmed) retention index.

    For temperature-programmed GC, retention times increase linearly with
    carbon number, so the simpler linear interpolation formula is used:
    RI = 100 * (n + (tR(x) - tR(n)) / (tR(n+1) - tR(n)))

    Args:
        rt_compound: retention time of the target compound (seconds)
        alkane_rts: list of alkane retention times (seconds), sorted ascending
        alkane_carbons: list of carbon numbers matching alkane_rts
    """
    alkane_rts = np.array(alkane_rts, dtype=np.float64)
    alkane_carbons = np.array(alkane_carbons)

    idx = np.searchsorted(alkane_rts, rt_compound) - 1
    idx = max(0, min(idx, len(alkane_rts) - 2))

    n_z = alkane_carbons[idx]
    rt_z = alkane_rts[idx]
    rt_z1 = alkane_rts[idx + 1]

    ri = 100 * (n_z + (rt_compound - rt_z) / (rt_z1 - rt_z))
    return ri

# Alkane ladder C8-C30
alkane_rts = [180.2, 240.5, 310.1, 385.7, 462.3, 540.8, 621.4, 703.9,
              788.1, 874.2, 962.0, 1051.5, 1142.8, 1235.6, 1330.1,
              1426.0, 1523.4, 1622.1, 1722.3, 1823.8, 1926.5, 2030.7, 2136.2]
alkane_carbons = list(range(8, 31))

ri = linear_ri(rt_compound=450.0, alkane_rts=alkane_rts, alkane_carbons=alkane_carbons)
print(f"Linear RI: {ri:.0f}")

4. Spectral Deconvolution and Library Matching with matchms

matchms supports GC-MS spectral comparison against reference libraries such as NIST and MassBank.

from matchms.importing import load_from_msp
from matchms.similarity import CosineGreedy
from matchms import calculate_scores
from matchms.filtering import (
    default_filters,
    normalize_intensities,
    select_by_mz,
    select_by_relative_intensity,
)

def process_spectrum(spectrum):
    spectrum = default_filters(spectrum)
    spectrum = normalize_intensities(spectrum)
    spectrum = select_by_mz(spectrum, mz_from=35, mz_to=500)
    spectrum = select_by_relative_intensity(spectrum, intensity_from=0.01)
    return spectrum

# Load reference library (NIST-format MSP export)
references = list(load_from_msp("nist_ei_library.msp"))
references = [process_spectrum(s) for s in references if s is not None]

# Load query spectra (deconvolved peaks exported as MSP)
queries = list(load_from_msp("sample_deconvolved.msp"))
queries = [process_spectrum(s) for s in queries if s is not None]

# Calculate cosine similarity scores
similarity = CosineGreedy(tolerance=0.5)
scores = calculate_scores(references, queries, similarity)

# Extract best matches per query
for i, query in enumerate(queries):
    best = scores.scores_by_query(query, "CosineGreedy_score", sort=True)
    # Each entry is (reference, score_object); access style depends on matchms version
    top_matches = []
    for ref, score_obj in best[:3]:
        sc = float(score_obj["score"] if isinstance(score_obj, dict) else score_obj[0])
        nm = int(score_obj["matches"] if isinstance(score_obj, dict) else score_obj[1])
        if sc > 0.7 and nm >= 6:
            top_matches.append((ref, sc, nm))
    if top_matches:
        ref, score, n_matches = top_matches[0]
        name = ref.get("compound_name", "Unknown")
        print(f"Query RT={query.get('retention_time')}: {name} "
              f"(score={score:.3f}, matched_peaks={n_matches})")

5. NIST Library Matching Workflow

For production-grade NIST matching, combine RI filtering with spectral similarity:

def nist_match_with_ri(query_spectrum, query_ri, references, ri_tolerance=20):
    """Filter candidates by RI window, then rank by spectral similarity."""
    similarity = CosineGreedy(tolerance=0.5)

    candidates = []
    for ref in references:
        ref_ri = ref.get("retention_index")
        if ref_ri is not None and abs(float(ref_ri) - query_ri) <= ri_tolerance:
            score_pair = similarity.pair(ref, query_spectrum)
            if score_pair is not None:
                # matchms >= 0.18 returns a dict-like object; older versions return a tuple
                score_val = float(score_pair["score"] if isinstance(score_pair, dict) else score_pair[0])
                n_matches = int(score_pair["matches"] if isinstance(score_pair, dict) else score_pair[1])
                if score_val > 0.6 and n_matches >= 5:
                    candidates.append({
                        "name": ref.get("compound_name", "Unknown"),
                        "score": score_val,
                        "matched_peaks": n_matches,
                        "ri_diff": abs(float(ref_ri) - query_ri),
                    })

    return sorted(candidates, key=lambda x: x["score"], reverse=True)

6. Peak Detection and Integration with PyOpenMS

import pyopenms as ms

exp = ms.MSExperiment()
ms.MzMLFile().load("gcms_sample.mzML", exp)

# Smooth spectra
gaussian = ms.GaussFilter()
params = gaussian.getParameters()
params.setValue("gaussian_width", 0.15)
gaussian.setParameters(params)
gaussian.filterExperiment(exp)

# Peak picking (centroiding)
picker = ms.PeakPickerHiRes()
pick_params = picker.getParameters()
pick_params.setValue("signal_to_noise", 3.0)
picker.setParameters(pick_params)

centroided = ms.MSExperiment()
picker.pickExperiment(exp, centroided)
print(f"Centroided spectra: {centroided.getNrSpectra()}")

Common GC-MS Issues and Solutions

TIC vs EIC Selection

• TIC (Total Ion Chromatogram): use for overall sample overview and retention time alignment. Prone to masking low-abundance compounds.
• EIC (Extracted Ion Chromatogram): use for targeted quantification. Select characteristic fragment ions (e.g., m/z 73, 147 for TMS derivatives; m/z 117 for amino acid TMS).

Derivatization Artifacts

TMS (trimethylsilyl) derivatization is the most common approach for GC-MS metabolomics. Common reagents include MSTFA (N-Methyl-N-(trimethylsilyl)trifluoroacetamide) and BSTFA (N,O-Bis(trimethylsilyl)trifluoroacetamide), often with 1% TMCS as catalyst. TMS derivatization produces multiple derivatives per metabolite:

# Common TMS artifacts to flag and exclude
TMS_ARTIFACTS = {
    "hexamethylcyclotrisiloxane": 207,   # D3 siloxane, column bleed
    "octamethylcyclotetrasiloxane": 281,  # D4 siloxane, column bleed
    "trimethylsilanol": 75,               # TMS reagent byproduct
}

def flag_artifacts(peak_list, artifact_mz_values, tolerance=0.5):
    """Flag peaks matching known derivatization artifacts."""
    flagged = []
    for peak in peak_list:
        base_peak_mz = peak.get("base_peak_mz", 0)
        is_artifact = any(
            abs(base_peak_mz - mz) <= tolerance
            for mz in artifact_mz_values
        )
        flagged.append({**peak, "is_artifact": is_artifact})
    return flagged

Column Bleed Detection

Column bleed produces characteristic siloxane fragments that increase with temperature:

COLUMN_BLEED_MZ = [207.0, 281.0, 355.0, 429.0]

def detect_column_bleed(spectrum_mz, spectrum_intensity, threshold=0.05):
    """Check if column bleed fragments dominate the spectrum."""
    total_intensity = np.sum(spectrum_intensity)
    if total_intensity == 0:
        return False
    bleed_intensity = sum(
        spectrum_intensity[np.argmin(np.abs(spectrum_mz - mz))]
        for mz in COLUMN_BLEED_MZ
        if np.min(np.abs(spectrum_mz - mz)) < 0.5
    )
    return (bleed_intensity / total_intensity) > threshold

Workflow: Complete GC-MS Metabolomics Pipeline

# 1. Load raw data
exp = ms.MSExperiment()
ms.MzMLFile().load("gcms_sample.mzML", exp)

# 2. Smooth and centroid
gaussian = ms.GaussFilter()
g_params = gaussian.getParameters()
g_params.setValue("gaussian_width", 0.15)
gaussian.setParameters(g_params)
gaussian.filterExperiment(exp)

picker = ms.PeakPickerHiRes()
centroided = ms.MSExperiment()
picker.pickExperiment(exp, centroided)

# 3. Export spectra for library matching
ms.MzMLFile().store("gcms_centroided.mzML", centroided)

# 4. Calculate RIs for detected peaks using alkane ladder
# 5. Match against NIST/MassBank using matchms with RI filtering
# 6. Export annotated feature table

Best Practices

• Always run alkane standards in the same batch for reliable RI calculation
• Set mass range to 35-500 m/z for EI spectra (below 35 is noise)
• Require minimum 6 matched peaks and cosine score > 0.7 for confident library matches
• Check for multiple TMS derivatives of the same metabolite (mono-, di-, tri-TMS)
• Use EIC for quantification, TIC for qualitative screening
• Blank subtraction is critical: run solvent blanks and derivatization blanks
• Monitor column bleed at high-temperature segments of the temperature program
• RI tolerance of 10-20 units for intra-lab matching, 30-50 for cross-lab comparison

Resources

• matchms documentation: https://matchms.readthedocs.io
• OpenMS documentation: https://www.openms.org
• NIST Chemistry WebBook: https://webbook.nist.gov
• Golm Metabolome Database (GMD): http://gmd.mpimp-golm.mpg.de
• MassBank spectral database: https://massbank.eu