msqrobpy

Python-native differential abundance analysis for quantitative LC-MS proteomics, inspired by the statistical design of the Bioconductor package msqrob2.

msqrobpy fits a linear model per feature (protein or peptide), moderates the residual variances with an empirical Bayes procedure, and tests user-defined linear contrasts. It works directly on polars DataFrames and does not reproduce the Bioconductor object model (QFeatures, SummarizedExperiment).

Features

Per-feature ordinary least squares (OLS) or robust M-estimation (Huber, IRLS) regression
Empirical Bayes moderation of residual variances, following the limma fitFDist / squeezeVar approach
Linear contrast testing with Benjamini-Hochberg multiple-testing correction
Peptide-to-protein aggregation utilities
A hurdle workflow that combines an abundance model with a detection (logistic) model
Cross-validation scripts that compare the output against msqrob2 in R

Installation

Requires Python 3.10 or later.

pip install git+https://github.com/CompOmics/MSqRobPy.git

For development:

git clone https://github.com/CompOmics/MSqRobPy.git
cd MSqRobPy
pip install -e .

Dependencies: numpy>=1.24, polars>=1.0, scipy>=1.10, statsmodels>=0.15.0, patsy>=0.5.6.

Dependency notes

Every public function takes and returns polars frames, and the numerical core is numpy. pandas is not imported anywhere in the package, but statsmodels declares it as a hard requirement, so pip still installs it. statsmodels receives plain numpy arrays, never a DataFrame.

The statsmodels lower bound is higher than the version it declares for itself: releases up to and including 0.14.4 import a private helper that scipy has removed, so import msqrobpy fails with ImportError: cannot import name '_lazywhere' from 'scipy._lib._util'. Verified broken with scipy 1.17.1, fixed by statsmodels>=0.15.0.

A stack verified end to end: numpy 2.4.6, scipy 1.17.1, polars 1.44.2, statsmodels 0.15.0, patsy 1.0.3. The declared floors polars==1.0.0 and patsy==0.5.6 were also checked against the full test suite.

On Windows terminals using a legacy code page, printing a polars frame can raise UnicodeEncodeError on the box-drawing characters. Either set PYTHONIOENCODING=utf-8 or switch the renderer:

import polars as pl
pl.Config.set_ascii_tables(True)

Data model

Polars has no row index, so both inputs carry their labels in an ordinary column:

intensity_df: one row per feature. One column holds the feature identifiers (feature_col, default "feature_id"); every other column is a sample. Values are assumed to be log-transformed and normalized. Missing values are null or NaN.
sample_metadata: one row per sample. One column holds the sample names (sample_col, default "sample"); the others hold the experimental factors used in the formula.

Sample alignment is an explicit left join of sample_metadata onto the sample columns of intensity_df, so the row order of sample_metadata does not matter. A sample column with no matching metadata row raises ValueError.

Tutorial

1. Input data

A synthetic dataset is included for testing:

from msqrobpy.example import simulate_protein_data

intensity_df, sample_metadata = simulate_protein_data(
    n_features=200,
    n_replicates_per_group=4,
)

print(intensity_df.select("feature_id", "S1", "S2", "S3", "S4").head(3))
# ┌────────────┬─────────┬─────────┬─────────┬─────────┐
# │ feature_id ┆ S1      ┆ S2      ┆ S3      ┆ S4      │
# │ str        ┆ f64     ┆ f64     ┆ f64     ┆ f64     │
# ╞════════════╪═════════╪═════════╪═════════╪═════════╡
# │ P1         ┆ 24.3779 ┆ 25.4321 ┆ 25.2482 ┆ 25.4380 │
# │ P2         ┆ 25.7072 ┆ 24.7912 ┆ 25.2925 ┆ 25.7186 │
# │ P3         ┆ 24.6030 ┆ 24.7451 ┆ 24.2956 ┆ 23.9691 │
# └────────────┴─────────┴─────────┴─────────┴─────────┘

print(sample_metadata.head(3))
# ┌────────┬───────────┐
# │ sample ┆ condition │
# │ str    ┆ str       │
# ╞════════╪═══════════╡
# │ S1     ┆ control   │
# │ S2     ┆ control   │
# │ S3     ┆ control   │
# └────────┴───────────┘

2. Fit the models

fit_protein_model fits one linear model per row of intensity_df. The model is specified with a patsy formula referring to columns of sample_metadata.

from msqrobpy import fit_protein_model

fit = fit_protein_model(
    intensity_df,
    sample_metadata,
    formula="~ condition",
    feature_col="feature_id",   # default
    sample_col="sample",        # default
    robust=True,                # Huber M-estimation; set False for OLS
    empirical_bayes=True,       # moderate residual variances across features
)

print(fit.design_columns)
# ['Intercept', 'condition[T.treated]']

Features with fewer observed values than model parameters are skipped. Fitted coefficients are available as a feature x coefficient table:

print(fit.coefficients().head(3))
# ┌────────────┬───────────┬──────────────────────┐
# │ feature_id ┆ Intercept ┆ condition[T.treated] │
# │ str        ┆ f64       ┆ f64                  │
# ╞════════════╪═══════════╪══════════════════════╡
# │ P1         ┆ 25.1241   ┆ -0.0566              │
# │ P2         ┆ 25.3774   ┆ -0.5099              │
# │ P3         ┆ 24.4032   ┆ 0.0381               │
# └────────────┴───────────┴──────────────────────┘

3. Test a contrast

A contrast is a linear combination of design matrix columns. Pass the coefficient name directly, or an arithmetic expression over coefficient names.

res = fit.test_contrast("condition[T.treated]")

print(res.top(5).select("feature_id", "estimate", "std_error", "t_stat", "df", "p_value", "adj_p_value"))
# ┌────────────┬──────────┬───────────┬────────┬─────┬─────────┬─────────────┐
# │ feature_id ┆ estimate ┆ std_error ┆ t_stat ┆ df  ┆ p_value ┆ adj_p_value │
# │ str        ┆ f64      ┆ f64       ┆ f64    ┆ f64 ┆ f64     ┆ f64         │
# ╞════════════╪══════════╪═══════════╪════════╪═════╪═════════╪═════════════╡
# │ P74        ┆ 1.8037   ┆ 0.3179    ┆ 5.6746 ┆ inf ┆ 0.0000  ┆ 0.0000      │
# │ P66        ┆ 1.6577   ┆ 0.3257    ┆ 5.0895 ┆ inf ┆ 0.0000  ┆ 0.0000      │
# │ P56        ┆ 1.5946   ┆ 0.3258    ┆ 4.8941 ┆ inf ┆ 0.0000  ┆ 0.0001      │
# │ P151       ┆ 1.5915   ┆ 0.3359    ┆ 4.7376 ┆ inf ┆ 0.0000  ┆ 0.0001      │
# │ P41        ┆ 1.3777   ┆ 0.3169    ┆ 4.3473 ┆ inf ┆ 0.0000  ┆ 0.0006      │
# └────────────┴──────────┴───────────┴────────┴─────┴─────────┴─────────────┘

import polars as pl
print(res.table.filter(pl.col("adj_p_value") < 0.05).height)
# 14

res.table is a plain polars.DataFrame with these columns:

Column	Meaning
`feature_id`	Feature identifier, named after `feature_col`
`contrast`	Contrast name
`estimate`	Contrast estimate (log fold change on the input scale)
`std_error`	Standard error, using the moderated variance when available
`t_stat`	Moderated t-statistic
`df`	Posterior residual degrees of freedom (`inf` when the prior df is infinite)
`p_value`	Two-sided p-value
`adj_p_value`	Benjamini-Hochberg adjusted p-value
`sigma`, `sigma_posterior`	Raw and moderated residual standard deviation
`method`	`ols` or `rlm`

4. Multi-factor designs and contrast arithmetic

Add covariates to the formula and combine coefficients in the contrast expression:

fit = fit_protein_model(intensity_df, sample_metadata, "~ condition + batch", robust=True)
print(fit.design_columns)
# ['Intercept', 'condition[T.B]', 'condition[T.C]', 'batch[T.b2]']

res = fit.test_contrast("condition[T.C] - condition[T.B]", name="C_vs_B")
print(res.top(3).select("feature_id", "contrast", "estimate", "std_error", "t_stat", "adj_p_value"))
# ┌────────────┬──────────┬──────────┬───────────┬─────────┬─────────────┐
# │ feature_id ┆ contrast ┆ estimate ┆ std_error ┆ t_stat  ┆ adj_p_value │
# │ str        ┆ str      ┆ f64      ┆ f64       ┆ f64     ┆ f64         │
# ╞════════════╪══════════╪══════════╪═══════════╪═════════╪═════════════╡
# │ P17        ┆ C_vs_B   ┆ 1.0588   ┆ 0.4720    ┆ 2.2431  ┆ 0.3319      │
# │ P37        ┆ C_vs_B   ┆ -1.0731  ┆ 0.4720    ┆ -2.2735 ┆ 0.3319      │
# │ P39        ┆ C_vs_B   ┆ 1.2421   ┆ 0.4720    ┆ 2.6314  ┆ 0.3319      │
# └────────────┴──────────┴──────────┴───────────┴─────────┴─────────────┘

A contrast may also be given as a sequence ordered like design_columns, or as a mapping from coefficient name to weight. Both are convenient for programmatically generated contrasts:

fit.test_contrast([0.0, -1.0, 1.0, 0.0], name="C_vs_B")
fit.test_contrast({"condition[T.C]": 1.0, "condition[T.B]": -1.0}, name="C_vs_B")

lme4-style random effect terms such as (1 | run) are stripped from the formula before the design matrix is built. Mixed models are not fitted; only the fixed-effect part is used.

5. Peptide-to-protein aggregation

If the input is a long peptide table, aggregate it to a feature x sample frame first:

import polars as pl
from msqrobpy import aggregate_peptides

long_df = pl.DataFrame({
    "protein":   ["P1", "P1", "P1", "P1", "P2", "P2"],
    "peptide":   ["pepA", "pepB", "pepA", "pepB", "pepC", "pepC"],
    "sample":    ["S1", "S1", "S2", "S2", "S1", "S2"],
    "intensity": [20.1, 19.8, 21.3, 21.0, 18.4, 18.9],
})

protein_df = aggregate_peptides(
    long_df,
    protein_col="protein",
    peptide_col="peptide",
    sample_col="sample",
    intensity_col="intensity",
    min_peptides=2,          # drop proteins with fewer distinct peptides
)
print(protein_df)
# ┌─────────┬─────────┬─────────┐
# │ protein ┆ S1      ┆ S2      │
# │ str     ┆ f64     ┆ f64     │
# ╞═════════╪═════════╪═════════╡
# │ P1      ┆ 19.9500 ┆ 21.1500 │
# └─────────┴─────────┴─────────┘

The result feeds straight into the model. Name the identifier column through feature_col, and make sure sample_metadata covers the sample columns the aggregation produced:

protein_meta = pl.DataFrame({"sample": ["S1", "S2"], "condition": ["control", "treated"]})
fit = fit_protein_model(protein_df, protein_meta, "~ condition", feature_col="protein")

aggregate_features does the same for any grouping column (peptide ID, PTM site), using min_observations in place of min_peptides.

The default summary is the median (robust_summary), which runs natively in polars. Pass any callable through summary_func to replace it; it receives the group values as a numpy array and must return a scalar.

6. Hurdle model

fit_hurdle_model fits the abundance model on observed intensities and, in parallel, a logistic regression on the detection pattern (observed versus missing). Evidence from both components is combined with Stouffer's method.

from msqrobpy import fit_hurdle_model

hurdle = fit_hurdle_model(intensity_df, sample_metadata, "~ condition", robust=True)
res = hurdle["test_contrast"]("condition[T.treated]")

print(res.top(3).select(
    "feature_id", "abundance_estimate", "abundance_p_value",
    "detection_estimate", "detection_p_value", "combined_p_value", "adj_combined_p_value",
))
# ┌────────────┬──────────────┬──────────────┬─────────────┬─────────────┬─────────────┬─────────────┐
# │ feature_id ┆ abundance_es ┆ abundance_p_ ┆ detection_e ┆ detection_p ┆ combined_p_ ┆ adj_combine │
# │            ┆ timate       ┆ value        ┆ stimate     ┆ _value      ┆ value       ┆ d_p_value   │
# │ str        ┆ f64          ┆ f64          ┆ f64         ┆ f64         ┆ f64         ┆ f64         │
# ╞════════════╪══════════════╪══════════════╪═════════════╪═════════════╪═════════════╪═════════════╡
# │ P74        ┆ 1.8037       ┆ 0.0000       ┆ NaN         ┆ NaN         ┆ 0.0000      ┆ 0.0002      │
# │ P141       ┆ 1.2951       ┆ 0.0003       ┆ NaN         ┆ NaN         ┆ 0.0003      ┆ 0.0294      │
# │ P21        ┆ -1.0330      ┆ 0.0026       ┆ NaN         ┆ NaN         ┆ 0.0026      ┆ 0.1319      │
# └────────────┴──────────────┴──────────────┴─────────────┴─────────────┴─────────────┴─────────────┘

The returned dictionary holds abundance_fit, detection_models, and the test_contrast callable. The result table reports abundance_estimate, detection_estimate, their p-values, and combined_p_value with a BH-adjusted counterpart.

ContrastResult.top() ranks the hurdle table on adj_combined_p_value, because the abundance-model column adj_p_value is absent here. Pass sort_by to rank on any other column. Features whose detection pattern is constant across samples are dropped from the detection component and reported on the abundance evidence alone, which is why detection_estimate is NaN in the rows above. A perfectly separated detection pattern makes the logistic fit diverge instead, which produces large estimates with uninformative p-values.

API reference

Object	Purpose
`fit_protein_model(intensity_df, sample_metadata, formula, feature_col="feature_id", sample_col="sample", robust=True, empirical_bayes=True, maxiter=5)`	Fit per-feature models; returns a `ProteinModelFit`
`ProteinModelFit.test_contrast(contrast, name=None)`	Test a contrast given as a name, expression, sequence or mapping; returns a `ContrastResult`
`ProteinModelFit.coefficients()`	Feature x coefficient frame
`ContrastResult.table` / `.top(n, sort_by)`	Full result table / top-ranked features
`fit_hurdle_model(...)`	Combined abundance and detection workflow
`aggregate_peptides(...)`, `aggregate_features(...)`, `robust_summary(...)`	Peptide-to-protein and generic feature summarization
`FeatureModelResult`	Per-feature coefficients, unscaled covariance, sigma, df, weights

FeatureModelResult.coef and .vcov_unscaled are numpy arrays ordered by .design_columns. Use .coef_dict() for a name-keyed view.

Migrating from the pandas API

Earlier versions took pandas DataFrames with the labels in the index. The changes:

Before	Now
Feature labels in `intensity_df.index`	Column named by `feature_col`, default `"feature_id"`
Sample labels in `sample_metadata.index`	Column named by `sample_col`, default `"sample"`
`feature_axis=1` for features in columns	Removed. Transpose with `intensity_df.transpose(...)` before the call
`FeatureModelResult.coef` is a `pd.Series`	numpy array plus `.design_columns`, or `.coef_dict()`
`FeatureModelResult.vcov_unscaled` is a `pd.DataFrame`	numpy array indexed by position in `.design_columns`
`ContrastResult.table` is a `pd.DataFrame`	`polars.DataFrame`
`summary_func` receives a `pd.Series`	receives a numpy array

A pandas input converts with pl.from_pandas(df.reset_index(names="feature_id")).

Implementation notes

Robust fitting uses Huber M-estimation with tuning constant k = 1.345 and a MAD-based scale re-estimated at each IRLS iteration. The default of 5 iterations matches the msqrobLm default in msqrob2.
The unscaled covariance matrix is (X'WX)^-1, matching .vcovUnscaled() in msqrob2. For robust fits, the residual degrees of freedom are sum(weights) - rank.
Variance moderation follows limma: the centred log-variances are used to isolate the prior trigamma term, which is inverted to recover the prior degrees of freedom d0 and prior variance s0^2. When the estimated between-feature variance is non-positive, d0 is set to infinity and the moderated t-statistic is referred to a normal distribution.
patsy builds the design matrix from a dict of numpy arrays taken from the metadata frame, so no pandas object is constructed.

Validating against msqrob2

The tests/cross_validation/ directory contains scripts that run both implementations on the same synthetic dataset and compare the output:

cd tests/cross_validation
python run_all.py

This requires R with msqrob2 installed (BiocManager::install("msqrob2")). The script enforces a 5% relative tolerance on the regression quantities and 20% on the empirical Bayes quantities and the statistics derived from them.

Both bounds are far looser than the measured agreement. On the bundled dataset (100 features, 8 samples, 5% missing), run against R 4.3.3 with msqrob2 1.10.0 and limma 3.58.1:

Quantity	Measured max relative difference
Coefficients, sigma, residual df, unscaled vcov, logFC	0% (OLS), at most 0.0062% (RLM)
Posterior variance and df	0.22% (OLS), 0.0001% (RLM)
Standard error and t-statistic	0.11% (OLS), at most 0.0064% (RLM)
p-value and adjusted p-value	1.5% worst case, 0.099% median (OLS); 0.0019% (RLM)

The largest single gap is the OLS prior variance, 0.153619 in R against 0.153959 here, which propagates into the OLS standard errors and p-values. The RLM differences come from the IRLS implementation rather than the moderation: MASS::rlm and _rlm_fit land about 6e-5 apart in relative terms. These margins hold for this dataset; the enforced tolerances stay loose because a different variance structure can widen the gap between the two moderation strategies. See tests/cross_validation/README.md for the full table.

generate_test_data.py now names the identifier column of each exported matrix (feature_id, sample) instead of leaving the header blank. run_msqrob2_r.R reads those files with row.names = 1, which takes the first column whatever its header, so the R side is unaffected.

Tests

pip install pytest
pytest tests/test_basic.py

License

Apache License 2.0. See LICENSE.

Navigation Menu

Repository files navigation

msqrobpy

Features

Installation

Dependency notes

Data model

Tutorial

1. Input data

2. Fit the models

3. Test a contrast

4. Multi-factor designs and contrast arithmetic

5. Peptide-to-protein aggregation

6. Hurdle model

API reference

Migrating from the pandas API

Implementation notes

Validating against msqrob2

Tests

License

About

Releases

Packages

Contributors

Languages

Name		Name	Last commit message	Last commit date
dist		dist
src		src
tests		tests
LICENSE		LICENSE
README.md		README.md
demo.py		demo.py
pyproject.toml		pyproject.toml

Latest commit

History

Folders and files

Repository files navigation

msqrobpy

Features

Installation

Dependency notes

Data model

Tutorial

1. Input data

2. Fit the models

3. Test a contrast

4. Multi-factor designs and contrast arithmetic

5. Peptide-to-protein aggregation

6. Hurdle model

API reference

Migrating from the pandas API

Implementation notes

Validating against msqrob2

Tests

License

About

Resources

Stars

Watchers

Forks

Releases

Packages

Contributors

Languages