MaldiAMRKit - Peak Sets

This notebook covers the peak sets: compact (m/z, intensity) containers (PeakSet / PeakList) that represent each spectrum as a variable-length set of peaks for downstream peak-set models. It shows how to extract peak sets, rank and pad peaks, persist and cache them, and apply fit-free peak alignment.

Import Libraries

import numpy as np

from maldiamrkit import MaldiSet
from maldiamrkit.alignment import align_peaks
from maldiamrkit.detection import MaldiPeakDetector, PeakList, create_peakset_input

Load Dataset

data = MaldiSet.from_directory(
    "../data/",
    "../data/metadata/metadata.csv",
    aggregate_by=dict(antibiotics="Drug"),
)
X = data.X
print(f"Binned spectra shape: {X.shape}")

Binned spectra shape: (29, 6000)

Extract a PeakList from binned spectra

MaldiPeakDetector.transform_peaklist detects peaks per spectrum and returns a PeakList - one PeakSet per spectrum, each a set of (m/z, intensity) peaks. Extraction is stateless (a pure per-spectrum function), so the result is identical whether computed over the whole dataset or inside a single cross-validation fold. The m/z axis is recovered from the binned column labels.

detector = MaldiPeakDetector(method="local", prominence=1e-4)
peaks = detector.transform_peaklist(X, top_k=50)

print(peaks)
print(f"Spectra: {len(peaks)}")
print(f"Per-spectrum peak counts: {peaks.n_peaks[:8]} ...")
print(f"warped (pre-alignment): {peaks.meta['warped']}")

PeakList(n=29, method=local, top_k=50)
Spectra: 29
Per-spectrum peak counts: [50 50 50 50 50 50 50 50] ...
warped (pre-alignment): False

Inspect a single PeakSet

Peaks are stored sorted by ascending m/z. Index the PeakList by position to get a PeakSet; as_array() returns an (n_peaks, 2) array of (m/z, intensity).

ps = peaks[0]
print(ps)
print(f"n_peaks: {ps.n_peaks}")
print(f"first 5 (m/z, intensity):\n{ps.as_array()[:5]}")

PeakSet(n_peaks=50, mz=[2012.0..9530.0])
n_peaks: 50
first 5 (m/z, intensity):
[[2.01200000e+03 6.87894941e-04]
 [2.03000000e+03 3.54902016e-03]
 [2.04500000e+03 9.10884239e-04]
 [2.05100000e+03 1.97589893e-03]
 [2.06000000e+03 9.63189503e-04]]

Rank and keep the strongest peaks

top_k keeps the highest-ranked peaks. The ranking is selected at detection time via rank_by: "intensity" (default), "prominence", or "persistence". For a non-intensity ranking, each PeakSet carries the per-peak score, so any later truncation (e.g. padding) stays consistent with how peaks were detected.

peaks_prom = detector.transform_peaklist(X, top_k=20, rank_by="prominence")
ps_prom = peaks_prom[0]

print(f"rank_by: {peaks_prom.meta['rank_by']}")
print(f"per-peak score carried: {ps_prom.score is not None}")
print(f"5 strongest m/z in spectrum 0: {ps_prom.top_k(5).mz}")

rank_by: prominence
per-peak score carried: True
5 strongest m/z in spectrum 0: [2030. 2201. 4361. 5375. 6284.]

Persistent-homology ranking

Persistent homology (method="ph") ranks peaks by topological persistence, which is robust to noise. Persistence ranking requires the PH detector.

ph_detector = MaldiPeakDetector(method="ph", persistence_threshold=5e-4)
peaks_ph = ph_detector.transform_peaklist(X, top_k=20, rank_by="persistence")

print(peaks_ph)
print(f"persistence scores (spectrum 0, first 5): {peaks_ph[0].score[:5]}")

PeakList(n=29, method=ph, top_k=20)
persistence scores (spectrum 0, first 5): [0.00354396 0.00167008 0.00120384 0.00132015 0.0027002 ]

Pad to a dense array for modelling

Peak-set models consume a fixed-width tensor. to_padded packs the variable-length sets into a dense (n_spectra, P, 2) array plus a boolean mask flagging real (non-padded) peaks. Spectra with more than P peaks keep their top-ranked ones.

values, mask = peaks.to_padded(max_peaks=50)

print(f"values shape: {values.shape}  (n_spectra, P, 2 = m/z, intensity)")
print(f"mask shape:   {mask.shape}")
print(f"real peaks in spectrum 0: {int(mask[0].sum())} / {mask.shape[1]}")

values shape: (29, 50, 2)  (n_spectra, P, 2 = m/z, intensity)
mask shape:   (29, 50)
real peaks in spectrum 0: 50 / 50

Index and slice a PeakList

Integer access returns a PeakSet; slicing returns a sub-PeakList that keeps the matching sample index and metadata.

subset = peaks[:5]
print(subset)
print(f"index: {list(subset.index)}")

PeakList(n=5, method=local, top_k=50)
index: ['10s', '11s', '12s', '13s', '14s']

Persist and cache

save / load round-trips a PeakList to .npz + .json. Because extraction is stateless and leak-safe, a PeakList can also be cached: pass cache_dir= and a repeat call with the same data and configuration reuses the stored result instead of recomputing.

import tempfile
from pathlib import Path

tmp = Path(tempfile.mkdtemp())

# save / load round-trip
peaks.save(tmp / "peaks")
reloaded = PeakList.load(tmp / "peaks")
print(f"reloaded {len(reloaded)} spectra (method={reloaded.meta['method']})")

# content + config keyed cache: the second call is a cache hit
cache = tmp / "cache"
detector.transform_peaklist(X, top_k=50, cache_dir=cache)
detector.transform_peaklist(X, top_k=50, cache_dir=cache)
print(f"cached files: {sorted(p.name for p in cache.glob('peaklist_*'))}")

reloaded 29 spectra (method=local)
cached files: ['peaklist_9593a4e7dea1d102_bdbc80b1d9a8bd5b.json', 'peaklist_9593a4e7dea1d102_bdbc80b1d9a8bd5b.npz']

Binning-free extraction from raw files

create_peakset_input is the faithful, binning-free path: it preprocesses each raw spectrum and detects peaks at full m/z resolution (no binning), returning an un-warped PeakList. It mirrors create_raw_input from the alignment module.

peaks_raw = create_peakset_input("../data/", top_k=50)

print(peaks_raw)
print(f"source: {peaks_raw.meta['source']}")
print(f"warped: {peaks_raw.meta['warped']}")

PeakList(n=29, method=local, top_k=50)
source: create_peakset_input
warped: False

Fit-free peak alignment

align_peaks warps a peak set’s m/z onto a reference peak list using the shared alignment strategies. It is fit-free: the caller supplies the reference (estimated on the training split, in-fold), so it never leaks. Intensities are carried through unchanged. Here we shift a spectrum onto a reference deliberately offset by 5 Da.

sample = peaks_raw[0]
ref_mz = sample.mz + 5.0  # generate an offset of 5 Da

aligned = align_peaks(sample, ref_mz, method="shift", max_shift_da=20.0)

print(f"median m/z shift applied: {np.median(aligned.mz - sample.mz):.2f} Da")
print(f"n_peaks unchanged: {aligned.n_peaks == sample.n_peaks}")
print(f"intensities unchanged: {np.array_equal(sample.intensity, aligned.intensity)}")

median m/z shift applied: 5.00 Da
n_peaks unchanged: True
intensities unchanged: True

Avoiding data leakage

Peak extraction is a pure per-spectrum function (meta["warped"] = False), so precomputing or caching a PeakList over the whole dataset cannot leak. Any fitted step (an align_peaks reference, standardization) must still be estimated on the training split, inside the fold.