GC-MS peak tables lie: align before you model

GC-MS data has a property that most instrument formats do not: the x-axis is not a physical constant. A chromatographic peak at 12.41 minutes in Monday's run can be at 12.38 minutes on Wednesday, and 12.52 after a column trim. Every step downstream — library search, peak-table comparison, multivariate modelling — inherits that instability.

The temptation is to take the vendor peak table and start modelling. Here is the order that actually works, and why each step has to come before the next.

First: get out of the vendor container

Vendors trap the useful outputs — TIC, extracted-ion chromatograms, scan spectra — inside proprietary containers: Agilent .D folders (the binary .ms file plus method metadata), Thermo .raw, Shimadzu .qgd. Two practical routes:

• Vendor export to netCDF or mzML (open formats), then parse in Python. Agilent ships an AIA/netCDF exporter; ProteoWizard's msconvert converts several vendor formats if installed with the right readers.
• Direct read of the .D folder from Python is possible (pymsfilereader on Windows uses the vendor libraries), but it ties your pipeline to a machine with the vendor stack installed.

The key point about .D folders: copy the whole folder, not just the .ms file. The method (temperature programme, scan rate, solvents delay) lives in the sibling files, and without it a scan index is not a time.

Once you have mzML, pyteomics gets you the TIC without ceremony:

from pyteomics import mzml
import numpy as np

times, tics = [], []
with mzml.read("run.mzML") as reader:
    for spec in reader:
        scan = spec.get("scanList", {}).get("scan") or [{}]
        rt = scan[0].get("scan start time")
        tic = spec.get("total ion current")
        if rt is not None and tic is not None:
            times.append(float(rt))
            tics.append(float(tic))

times, tics = np.asarray(times), np.asarray(tics)
# The unit lives in scanList.scan[0]['unitName'] — minutes or seconds.
# Assert it once, or your axis is silently off by 60x.
print(f"{len(times)} scans, {times[0]:.2f} -> {times[-1]:.2f}")

Step 1 — subtract the blank before anything else

Background ions, column bleed and carryover are shared across runs; they are the one component that *does* align by construction. If you detect peaks before blank subtraction, your peak table is partly a table of the column's history. Do it on the trace, not on the peak table: a peak that only exists in the blank should never become a feature.

# blank and sample as aligned arrays on the same time grid
net = sample_tics - blank_tics
net[net < 0] = 0

Record that you did this. "Reported as blank-subtracted" is a method statement, and reviewers will ask.

Step 2 — align retention times against anchors

Alignment means: find the time shift that best overlays shared peaks, then apply it. For simple drift on a single column, a robust anchor approach beats a full warping algorithm:

def estimate_shift(reference, run, times, window_idx):
    """Median shift (in samples) between an anchor peak in ref and a run."""
    def apex(y, t):
        lo = int(np.searchsorted(t, window_idx[0]))
        hi = int(np.searchsorted(t, window_idx[1]))
        return t[lo + np.argmax(y[lo:hi])]
    return apex(run, times) - apex(reference, times)

anchor = (11.8, 13.2)          # minutes around a peak you trust in every run
shift_samples = estimate_shift(ref_tics, net, times, anchor)
aligned = np.interp(times - shift_samples * np.median(np.diff(times)), times, net)

That handles constant drift. For non-linear drift (temperature-programme variability, column ageing), use a piecewise or correlation-optimized warping method — but the discipline is the same: align on peaks you can identify, verify the residual shift after alignment, and report it. If the median absolute residual is comparable to your peak width at half height, the alignment failed and every downstream comparison is suspect.

Where the chemistry allows, avoid the problem entirely: add a retention-index standard series and express peaks as Kovats indices (RI = 100·n + 100·(t − t_n)/(t_{n+1} − t_n) between bracketing alkanes). RI values are stable across runs, columns and instruments in a way raw minutes never are.

Step 3 — detect peaks on the aligned signal

Only now detect. On the aligned, blank-subtracted trace:

from scipy.signal import find_peaks

peak_idx, props = find_peaks(net, prominence=0.05 * net.max(),
                             width=(3, None))  # min 3 scans wide
peak_times = times[peak_idx]

Two rules of thumb worth writing into the documentation:

• Width gate. Anything narrower than a few scans is noise or a spike; anything that spans most of the run is column bleed.
• Co-elution is the default. A clean Gaussian peak is the exception. Extraction for quantitation should use extracted-ion chromatograms at qualifying masses, not the TIC apex — otherwise the "peak area" includes whatever else eluted at the same time.

Step 4 — the library match is a hypothesis

The NIST match factor (the familiar 0–999 number) is a normalized dot product between your spectrum and a library entry. It says "these two spectra are similar", not "this compound is present". Under co-elution, a 900+ match against the wrong compound is routine — library search optimizes similarity, and the library contains everything.

Minimum discipline for a call you will act on:

1. Qualifier ions. Confirm the proposed compound with 2–3 characteristic fragments, not just the base peak.
2. Standards or retention. The same spectrum at the right retention index is a much stronger claim than the same spectrum alone.
3. A written uncertainty. "Tentatively identified, match 912, co-elutes with an unidentified isomer" is an honest sentence. "Identified as X" from the same data is not.

Step 5 — the feature matrix, structured for modelling

If the goal is a model — batch classification, process monitoring, QC — build the matrix from aligned, blank-subtracted data:

# Bin the mass axis once, then aggregate every scan into (time bin, m/z bin).
# m/z binning wider than the mass calibration error is safer than matching exact masses.

Because vendor mass-axis calibration differs, exact-mass matching across instruments produces spurious features; a tolerance or a binned axis is the honest version. Then the usual small-data rules apply — see /tutorials/small-data for the evaluation discipline, and the conformal post for why every prediction should carry an interval.

The short version

Order matters: blank → baseline → align → detect → identify → model. Every skipped step reappears as unexplained variance in the model, where it is much harder to diagnose. The vendor peak table is an output of an unknown pipeline — useful for a sanity check, not for a training set.

The free parser at /tools/gcms reads supported GC-MS files into a parsed TIC with peaks and derived analytics (in memory, nothing stored) — useful for a one-file look before you commit to a pipeline. /docs/instruments/gcms covers the format specifics; if you mostly have IR or Raman data, the FTIR baseline post covers the spectral side of the same discipline.

Honest limits

• Alignment targets drift; aggressive warping can also distort peak shapes and areas. If quantitation matters, correct with internal standards rather than warping.
• Kovats indices require the alkane series and a stable temperature programme — they fix comparability, not peak assignment.
• Vendor-format direct readers depend on vendor DLLs and OS; the portable path is a conversion to mzML/netCDF recorded in the pipeline log.
• A correct analysis of a badly sampled run (overloaded column, saturated detector) is still a bad measurement. That is upstream of every step here.