Metadata-Version: 2.1
Name: maldiassist
Version: 1.1.0
Summary: Mathematical utilities for MALDI-TOF mass spectrometry (Python port of the MALDIassist R package)
Author: Wonseok Oh
License: MIT
Project-URL: Homepage, https://github.com/hiows/MALDIassist-py
Project-URL: Repository, https://github.com/hiows/MALDIassist-py
Project-URL: Documentation, https://github.com/hiows/MALDIassist-py#readme
Project-URL: Changelog, https://github.com/hiows/MALDIassist-py/releases
Project-URL: Bug Tracker, https://github.com/hiows/MALDIassist-py/issues
Requires-Python: >=3.9
Requires-Dist: numpy>=1.23
Requires-Dist: scipy>=1.9
Requires-Dist: pandas>=1.5
Provides-Extra: viz
Requires-Dist: matplotlib>=3.6; extra == "viz"
Provides-Extra: test
Requires-Dist: pytest>=7; extra == "test"
Description-Content-Type: text/markdown

# maldiassist (Python)

[![PyPI](https://img.shields.io/pypi/v/maldiassist.svg)](https://pypi.org/project/maldiassist/)

> Python package **v1.1.0** · Original R package: [hiows/MALDIassist](https://github.com/hiows/MALDIassist) (v1.0.2)

**maldiassist** is a Python port of the [MALDIassist](https://github.com/hiows/MALDIassist) R package
(v1.0.2 API surface; R + Rcpp/C++). It reproduces the **same algorithms and numerical results** within
floating-point tolerance. The workflow covers the full path from raw Bruker spectra to a
cohort-level peak matrix:

- Loading Bruker MALDI-TOF spectra
- Smoothing and baseline correction
- Gaussian KDE-based peak detection (including shoulder peaks)
- Peak-quality metrics and filtering (intensity / prominence / strength)
- Spectrum alignment to internal standards (linear / lowess)
- Cohort feature analysis (frequent m/z discovery, overlap filtering, matched peak matrix, two-group significance testing)
- Visualization (matplotlib)

Performance-critical kernel-regression routines are implemented in C++ (via
[nanobind](https://github.com/wjakob/nanobind), a direct port of the original Rcpp core).

## Example

The figure below shows the core workflow on real Bruker MALDI-TOF data: a raw spectrum (gray)
is smoothed and baseline-corrected (red), then peaks are detected and filtered (blue).

![MALDIassist workflow example](docs/figures/workflow_example.png)

---

## Installation

Requires Python 3.9+.

The package ships a small C++ extension (a nanobind port of the original Rcpp
kernel-regression core) that accelerates peak detection. Pre-built **wheels**
bundle this extension, so no compiler is needed:

```bash
pip install maldiassist
pip install "maldiassist[viz]"   # with visualization (matplotlib)
```

Alternatively, install a pre-built wheel from the
[GitHub Releases](https://github.com/hiows/MALDIassist-py/releases) page, or
build from source (requires a C++17 compiler and CMake, handled automatically
by `scikit-build-core`):

```bash
pip install "git+https://github.com/hiows/MALDIassist-py.git"
```

> The C++ extension is optional at runtime: if it cannot be imported the package
> transparently falls back to the pure-Python implementation, producing identical
> results (just slower).

For development (clone first, then editable install):

```bash
git clone https://github.com/hiows/MALDIassist-py.git
cd MALDIassist-py
pip install -e .              # core (numpy, scipy, pandas) + C++ extension
pip install -e ".[viz]"       # with visualization (matplotlib)
pip install -e ".[viz,test]"  # visualization + test tooling
```

The original R package can be installed from GitHub (CRAN submission pending):

```r
# install.packages("remotes")
remotes::install_github("hiows/MALDIassist")
```

---

## Quick start

The pipeline follows six steps: **load → preprocess → KDE → peak picking → peak filtering → visualization**.

```python
import maldiassist as ma
```

### 1. Load

Point `load_maldi_spectra()` at a directory of Bruker flex data files. It returns a
name-keyed dict of raw spectra.

```python
raw_spectra = ma.load_maldi_spectra("data/")
```

### 2. Preprocess

Apply Savitzky-Golay smoothing and baseline subtraction.

```python
preprocessed_spectra = ma.preprocess_maldi_spectra(
    raw_spectra,
    hws_sg=10,           # half-window size for Savitzky-Golay
    pno_sg=3,            # polynomial order
    baseline_type="snip",  # baseline algorithm
    iter_snip=50,        # SNIP iterations
)
```

Build Gaussian KDE spectra for peak filtering:

```python
kde_spectra = ma.build_kde_spectra(preprocessed_spectra, bw=1)
kde_spectrum_only = {k: v["spectrum"] for k, v in kde_spectra.items()}
```

### 3. Peak picking

`find_peaks_spectra()` detects peaks (including shoulder peaks) using a Gaussian KDE approach.

```python
peaks_list = ma.find_peaks_spectra(
    preprocessed_spectra,
    bw=1,                          # KDE bandwidth
    hws_peaks=10,
    weight_type="raw",
    cutoff_kappa_peak_strength=0.3,
    peak_retention_fraction=0.25,
)
```

### 4. Peak filtering

`filter_peaks_spectra()` removes low-quality peaks by intensity, prominence, and strength cutoffs.
Pass the **KDE spectra** (not the preprocessed spectra) as the `spectra` argument, matching the R API.

```python
filtered_peaks_list = ma.filter_peaks_spectra(
    kde_spectrum_only,
    peaks_list,
    cutoff_peak_intensity=100,
    cutoff_peak_prominence=50,
    cutoff_peak_strength=0.5,
    normalization_type="raw",
)
```

### 5. Visualization

Overlay a raw spectrum (gray), its preprocessed version (red), and the filtered peaks (blue)
in the same style as the Example figure (requires the `[viz]` extra). The exact spectrum
depends on your `data/` directory and which sample is selected.

```python
import matplotlib.pyplot as plt

example_range = (12000, 15000)
sample = next(iter(raw_spectra))
spec = raw_spectra[sample]
pp = preprocessed_spectra[sample]
fp = filtered_peaks_list[sample]

lo, hi = example_range
x_raw = spec.loc[(spec["mz"] >= lo) & (spec["mz"] <= hi), "mz"]
y_raw = spec.loc[(spec["mz"] >= lo) & (spec["mz"] <= hi), "intensity"]
x_pp = pp.loc[(pp["mz"] >= lo) & (pp["mz"] <= hi), "mz"]
y_pp = pp.loc[(pp["mz"] >= lo) & (pp["mz"] <= hi), "intensity"]

fig, ax = plt.subplots()
ax.plot(x_raw, y_raw, color="gray", lw=1.5, label="raw")
ax.plot(x_pp, y_pp, color="red", lw=2, label="preprocessed")
fp_r = fp[(fp["mz"] >= lo) & (fp["mz"] <= hi)]
ax.vlines(fp_r["mz"], 0, fp_r["intensity"], color="blue", lw=1.5, label="peaks")
ax.set_xlabel("m/z")
ax.set_ylabel("Intensity")
ax.set_ylim(0, y_raw.max())
plt.show()
```

You can also overlay all spectra in a single plot:

```python
ma.visualize_spectra(preprocessed_spectra)
```

---

## Cohort analysis

The example below follows a two-species MALDI-TOF cohort from
[PRIDE PXD058284](https://www.ebi.ac.uk/pride/archive/projects/PXD058284): load Bruker spectra
and sample metadata, preprocess and pick peaks, align across samples, build a matched-peak
matrix, and test for group-discriminating m/z features.

### 1. Load spectra and metadata

```python
import pandas as pd

raw_spectra = ma.load_maldi_spectra("data/raw/")
metadata = pd.read_excel("data/sample_metadata.xlsx")

# two-level species grouping (one label per sample)
sample_group = metadata.set_index("SampleID").reindex(raw_spectra.keys())["Species"].to_numpy()
```

### 2. Preprocess, KDE, and peak picking

```python
preprocessed_spectra = ma.preprocess_maldi_spectra(
    raw_spectra,
    hws_sg=10,
    pno_sg=3,
    baseline_type="snip",
    iter_snip=50,
)

kde_spectra = ma.build_kde_spectra(preprocessed_spectra, bw=1)
kde_spectrum_only = {k: v["spectrum"] for k, v in kde_spectra.items()}

peaks_list = ma.find_peaks_spectra(
    preprocessed_spectra,
    bw=1,
    hws_peaks=10,
    weight_type="raw",
    cutoff_kappa_peak_strength=0.3,
    peak_retention_fraction=0.25,
)

filtered_peaks_list = ma.filter_peaks_spectra(
    kde_spectrum_only,
    peaks_list,
    cutoff_peak_intensity=100,
    cutoff_peak_prominence=50,
    cutoff_peak_strength=0.5,
    normalization_type="raw",
)
```

### 3. Align spectra

`align_spectra()` corrects m/z drift using internal standards selected from frequent,
high-intensity peaks. Choose `"linear"` (two-point) or `"lowess"` (multi-point) alignment.
It returns one aligned `spectrum` / `peaks` pair per sample (in `alignment_results`) plus
the reference `standard_mz` values used as anchors.

```python
aligned = ma.align_spectra(
    kde_spectrum_only,
    filtered_peaks_list,
    bin_width=20,
    alignment_mode="linear",  # or "lowess"
    hws_alignment=50,
)

aligned_peaks = {k: v["peaks"] for k, v in aligned["alignment_results"].items()}
exclude_mz = list(aligned["standard_mz"].values())  # alignment anchors, excluded below
```

### 4. Find frequent m/z values

`find_frequent_mz()` scans pooled peak m/z values across the aligned samples and refines each
bin location with Gaussian KDE. Pass the alignment anchors to `exclude_mz` so the internal
standards are dropped from the feature set.

```python
freq_mz = ma.find_frequent_mz(
    aligned_peaks,
    bin_width=20,
    exclude_mz=exclude_mz,
)

# optional: drop nearby redundant features (R 1.0.2+)
freq_mz = ma.filter_overlap_mz(freq_mz, hws_selection=5)
```

### 5. Build a matched peak matrix

`build_matched_matrix()` matches each sample's peaks to the frequent m/z references and
returns a detection matrix (`detected_matrix`) and a signed m/z-difference matrix
(`delta_mz_matrix`), both sample-by-marker.

```python
matched = ma.build_matched_matrix(
    aligned_peaks,
    reference_mz=freq_mz["mz"].to_numpy(),
    hws_match=10,
)

mat = matched["detected_matrix"]
```

Visualize the matrix with `heatmap_matched_matrix()`, optionally annotated by a per-sample
grouping (requires `[viz]`):

```python
ma.heatmap_matched_matrix(
    mat,
    group=sample_group,   # one entry per sample (row)
    hide_rownames=True,
    hide_colnames=True,
)
```

### 6. Test for significant m/z features

`estimate_significance()` runs a per-feature two-group comparison (t-test or Wilcoxon) on a
sample-by-marker matrix and returns raw and adjusted p-values. Subset the matrix to the
significant markers to highlight the discriminating features.

```python
sig = ma.estimate_significance(
    mat,
    group=sample_group,      # two-level grouping, one entry per sample
    stat_method="t.test",    # or "wilcox"
    adj_method="BH",         # "none", "BH", or "bonferroni"
)

sig_cols = sig.loc[sig["adj_pvalue"] < 0.01, "feat_names"]
ma.heatmap_matched_matrix(
    mat[sig_cols],
    group=sample_group,
    hide_rownames=True,
    hide_colnames=True,
)
```

Applied to the PXD058284 two-species cohort (*E. coli* vs *K. pneumoniae*), the significant
markers cleanly separate the samples by species:

![Cohort matched-peak heatmap](results/heatmap_significant.png)

---

## R/Python parity

The core numerical pipeline was validated step-by-step against **MALDIassist** R v1.0.0 on the
same Bruker cohort and parameter set (Quick start / cohort analysis above). The public API
surface tracks R **v1.0.2**, with two intentional SemVer-friendly differences: heatmap
clustering still defaults to `"complete"` (R 1.0.2 defaults to `"ward.D2"`; pass
`clustering_method="ward.D2"` to match), and `estimate_peak_boundaries` remains available
but deprecated (removed in R 1.0.2). Summary of the v1.0.0 validation:

| Step | Match |
| --- | --- |
| Loading / preprocessing / peak detection / filtering / alignment | m/z and intensity within ~5×10⁻¹¹ |
| `find_frequent_mz` | identical row count; count and freq_ratio exact |
| `build_matched_matrix` | detection matrix exact (difference 0) |
| `estimate_significance` | p-value difference ~10⁻¹⁶; R and Python agree on adj. p < 0.01 calls |

Core algorithmic details (Nadaraya–Watson kernel regression with 1st–3rd derivatives,
SNIP/TopHat baselines, Savitzky–Golay boundary coefficients, R `hist`/`pretty` binning,
R `lowess`, `p.adjust`, Wilcoxon continuity correction, etc.) are reproduced identically.

> Note: kernel summation is computed in the **same sequential accumulation order** as R's
> Rcpp loop. NumPy's default pairwise summation can introduce tiny floating-point
> differences over symmetric windows that flip tie-breaking in extremum selection.

### Performance

The kernel-regression hot paths (Gaussian KDE grid evaluation with 1st–3rd derivatives and
bisection root finding) are implemented in C++ via nanobind (`src/spectrum_math_cpp.cpp`, a
direct port of the original Rcpp `spectrum_math.cpp`). The compiled backend uses the identical
sequential summation order as the pure-Python reference, so results match to floating-point
tolerance while running substantially faster. When the extension is unavailable,
`maldiassist.spectrum_math` falls back to the pure-Python path automatically. Set
`MALDIASSIST_DISABLE_CPP=1` to force the pure-Python backend (used by the parity tests in
`tests/test_kde_parity.py`).

---

## Main functions

| Function | Purpose |
| --- | --- |
| `load_maldi_spectra()` | Load Bruker raw spectra from a directory |
| `preprocess_maldi_spectra()` | Smooth and baseline-correct spectra |
| `find_peaks()` / `find_peaks_spectra()` | Detect ordinary and shoulder peaks (single / dict) |
| `find_peaks_fast()` / `find_peaks_spectra_fast()` | Fast local peak detection (single / dict) |
| `filter_peaks()` / `filter_peaks_spectra()` | Filter peaks by intensity, prominence, and strength |
| `build_kde_spectrum()` / `build_kde_spectra()` | Build Gaussian KDE spectra (single / dict) |
| `find_frequent_mz()` | Find frequent m/z values across a cohort |
| `filter_overlap_mz()` | Remove overlapping nearby frequent m/z features |
| `align_spectra()` | Align spectra to internal standards (linear / lowess) |
| `build_matched_matrix()` | Assemble a cohort peak intensity matrix |
| `estimate_significance()` | Two-group significance testing per m/z feature |
| `visualize_spectrum()` / `visualize_spectra()` | Visualize spectra with matplotlib (requires `[viz]`) |
| `heatmap_matched_matrix()` | Heatmap of a matched-peak matrix (requires `[viz]`) |

> `estimate_peak_boundaries()` is still exported but **deprecated** (removed from R 1.0.2)
> and will be removed in a future major release.

### Correspondence with the R package

| Step | R function | Python function |
| --- | --- | --- |
| Loading | `load_maldi_spectra` | `load_maldi_spectra` |
| Preprocessing | `preprocess_maldi_spectra` | `preprocess_maldi_spectra` |
| KDE | `build_kde_spectrum` / `build_kde_spectra` | `build_kde_spectrum` / `build_kde_spectra` |
| Peak detection | `find_peaks` / `find_peaks_spectra` | `find_peaks` / `find_peaks_spectra` |
| Peak filtering | `filter_peaks` / `filter_peaks_spectra` | `filter_peaks` / `filter_peaks_spectra` |
| Alignment | `align_spectra` | `align_spectra` |
| Frequent m/z | `find_frequent_mz` | `find_frequent_mz` |
| Overlap filter | `filter_overlap_mz` | `filter_overlap_mz` |
| Matched matrix | `build_matched_matrix` | `build_matched_matrix` |
| Significance test | `estimate_significance` | `estimate_significance` |
| Visualization | `visualize_spectrum/spectra`, `heatmap_matched_matrix` | `visualize_spectrum/spectra`, `heatmap_matched_matrix` |

---

## Author

**Wonseok Oh** ([ORCID: 0009-0002-0687-8466](https://orcid.org/0009-0002-0687-8466))

## How to cite

If you use **maldiassist** (this Python package) in your research, please cite:

```bibtex
@software{maldiassist_py,
  title   = {maldiassist: Mathematical Utilities for MALDI-TOF Mass Spectrometry (Python)},
  author  = {Wonseok Oh},
  year    = {2026},
  version = {1.1.0},
  url     = {https://github.com/hiows/MALDIassist-py}
}
```

The underlying R package **MALDIassist** can be cited separately (from R:
`citation("MALDIassist")`). Its v1.0.0 archive is
[10.5281/zenodo.21376279](https://doi.org/10.5281/zenodo.21376279); the R concept DOI is
[10.5281/zenodo.21219451](https://doi.org/10.5281/zenodo.21219451).

## License

MIT © 2026 Wonseok Oh. See [LICENSE](LICENSE) for details.

This project is a Python port of the [MALDIassist](https://github.com/hiows/MALDIassist)
R package (v1.0.2), which is also released under the MIT License (© 2026 Wonseok Oh).
See [CHANGELOG.md](CHANGELOG.md) for release notes.
