pymsstats

A pure-Python port of Bioconductor MSstats (Choi et al., Bioinformatics 2014) — vendor converters, pre-processing, normalization, peptide-to-protein summarization, per-protein LMM-based group comparison, sample-size design, and QC plotting for mass-spec proteomics.

• Full pipeline mirroring R: vendor → dataProcess (log2 + normalize + TMP/linear) → groupComparison → designSampleSize / plots
• No rpy2, no R install — every algorithm reimplemented directly in NumPy / SciPy / statsmodels / pandas
• 100% coverage of the R MSstats 4.14.2 public API — all 50 exported names (46 functions ported, 4 example datasets replaced by a synthetic generator)
• Bit-for-bit reproduction of TMP LogIntensities (5e-14), log2FC (8e-15), linear summarization (9e-14), designSampleSize numSample (exact), MSstatsContrastMatrix, checkRepeatedDesign and SDRFtoAnnotation against R MSstats 4.14.2
• 2.5×–3.9× faster than R MSstats on the 500-protein × 8-run benchmark (Pearson r = 1.0 on every output)
• Nine vendor converters (DIA-NN, Spectronaut, FragPipe, OpenMS, Skyline, Proteome Discoverer, Progenesis, OpenSWATH, DIA-Umpire) on top of MaxQuant
• Modular worker API (MSstatsPrepareFor*, MSstatsSummarize*, MSstatsGroupComparison*) + SDRF metadata helpers
• AnnData / pandas-friendly: accepts MSstats-format long DataFrame (10 canonical columns)

This is a standalone mirror of the canonical implementation that lives in omicverse. All algorithmic work is developed upstream in omicverse and synced here.

Install

pip install pymsstats

Depends only on numpy, scipy, pandas, and statsmodels.

Quick-start

import numpy as np
import pandas as pd
from pymsstats import data_process, group_comparison

# msstats_df is a long-format DataFrame with the 10 canonical MSstats columns:
#   ProteinName, PeptideSequence, PrecursorCharge, FragmentIon, ProductCharge,
#   IsotopeLabelType, Condition, BioReplicate, Run, Intensity
msstats_df = pd.read_csv("msstats_long.tsv", sep="\t")

# Pre-process: log2 → equalize-medians → drop low-obs features → TMP summarize.
processed = data_process(
    msstats_df,
    normalization="equalizeMedians",
    summary_method="TMP",
    log_base=2,
)
# processed columns: Protein, RUN, LogIntensities, GROUP, SUBJECT, n_features, n_obs

# Two-group test: group1 vs group0.
result = group_comparison(
    processed,
    contrast_matrix=np.array([[-1.0, 1.0]]),
    groups=["group0", "group1"],
)
# result columns: Protein, Label, log2FC, SE, df, t, pvalue, adj.pvalue, issue

Vendor converters

Six converters turn raw vendor output into MSstats long-format:

from pymsstats import (
    maxquant_to_msstats, diann_to_msstats, spectronaut_to_msstats,
    fragpipe_to_msstats, openms_to_msstats, skyline_to_msstats,
)

ann = pd.read_csv("annotation.csv")   # columns: Run, Condition, BioReplicate

# MaxQuant
pg = pd.read_csv("proteinGroups.txt", sep="\t")
ev = pd.read_csv("evidence.txt", sep="\t")
msstats_df = maxquant_to_msstats(pg, ev, ann)

# DIA-NN
report = pd.read_csv("report.tsv", sep="\t")
msstats_df = diann_to_msstats(report, ann, qvalue_cutoff=0.01)

# Spectronaut / FragPipe / OpenMS / Skyline
msstats_df = spectronaut_to_msstats(pd.read_csv("spectronaut.tsv", sep="\t"))
msstats_df = skyline_to_msstats(pd.read_csv("skyline.csv"), ann)

Sample-size design, quantification, plotting

from pymsstats import design_sample_size, quantification, group_comparison_plots

# how many replicates for 80% power at FC 1.25–1.5?
design = design_sample_size(processed, desired_fc=(1.25, 1.5), power=0.8)

# per-protein abundance matrix
quant = quantification(processed, type="Sample")

# volcano / heatmap on the comparison result
ax = group_comparison_plots(result, type="VolcanoPlot")

Low-level API

from pymsstats import equalize_medians, tmp_summarize, medpolish

# Each stage is exposed for users who want to slot one piece of the
# pipeline into a custom workflow.
norm_df = equalize_medians(feature_df, abundance_col="ABUNDANCE")
prot_df = tmp_summarize(norm_df)        # per (protein, run) LogIntensities
mp = medpolish(X)                       # raw R-stats::medpolish equivalent

R-parity status

The MSstats dataProcess → groupComparison pipeline is tested against R MSstats 4.14.2 on the canonical LFQ design (globally-unique BioReplicate IDs):

The OLS path (auto mode when SUBJECT is nested in GROUP) is bit-identical to R. The mixed path (paired/repeated-measures design where the same BioReplicate appears across groups) uses statsmodels.MixedLM; when its solver hits a singular Hessian (common for small balanced data) we fall back to a fixed-effects lm(ABUNDANCE ~ GROUP + SUBJECT), which gives slightly different SEs from R's lme4::lmer (the ranking — Pearson r ≥ 0.99 on pvalue — is unaffected).

Benchmark

500 proteins × 8 runs × 3 features, label-free DDA design:

Run it yourself:

python examples/benchmark.py --runs 3 --n-proteins 500

Reproducing R results exactly

tests/r_reference_driver.R runs MSstats::dataProcess + MSstats::groupComparison on the same TSV input the Python side dumps. tests/test_r_parity.py (skipped if no CMAP / MSstats env) checks every numerical field — see the R-parity table above.

# Run only the R-parity tests
pytest tests/test_r_parity.py -v

# Run everything (smoke + R-parity)
pytest tests/ -v

MSstats public API coverage (v0.3 — 100%)

All 50 exported names of R MSstats 4.14.2 are covered. The 46 functions are ported to Python; the 4 bundled example datasets are copyrighted Bioconductor data objects (not redistributed) and are replaced by a synthetic generator.

Core pipeline

Vendor converters

Statistical helpers

Modular worker API

SDRF helpers

Plotting

Bundled example datasets

The four R example datasets are copyrighted Bioconductor data objects; pymsstats provides make_example_dataset() / load_dda_example() / load_dia_example() / load_srm_example() which return a synthetic MSstats long-format DataFrame with the same 10 canonical columns so the pipeline can be exercised end-to-end without the original data.

Known limitations (within ported functions)

• The mixed-effects (lme4::lmer) path uses statsmodels.MixedLM with an OLS fallback on singular-Hessian cases (SE differs ≤ 5%; ranking unaffected — see R-parity table above).
• feature_subset='highQuality' in select_features falls back to 'all' with a warning (the leverage-based outlier rule is rarely exercised).
• TMT / labeled (heavy-light) workflows are out of scope (MSstatsTMT is a separate Bioconductor package).

Relationship to omicverse

Developed upstream in omicverse. This repo is a standalone mirror with the same code, same API, minus the omicverse packaging.

Citation

If you use this package, please cite the MSstats paper:

Choi, M., Chang, C.-Y., Clough, T., Broudy, D., Killeen, T., MacLean, B., & Vitek, O. MSstats: an R package for statistical analysis of quantitative mass spectrometry-based proteomic experiments. Bioinformatics 30 (17), 2524-2526 (2014). DOI: 10.1093/bioinformatics/btu305

…and acknowledge omicverse / this repo for the Python port.

License

Artistic-2.0 — matches the upstream Bioconductor package.