vsn 0.1.0

Pure Python implementation of Variance Stabilization and Normalization

VSN in Python

A NumPy/SciPy implementation of variance stabilization and normalization (VSN), based on the Bioconductor vsn package by Wolfgang Huber and contributors.

The implementation fits an additive-multiplicative error model, estimates sample-specific calibration parameters, applies the generalized logarithm, and uses the same robust least-trimmed-squares workflow as the R implementation.

The corrected implementation reproduces the supplied R vsn::justvsn() reference values to floating-point precision on the validated sparse seven-sample dataset.

Project components

File	Purpose
vsn2.py	Reusable Python VSN implementation
vsn2mo.py	Interactive Marimo normalization dashboard
run.py	Run VSN on selected columns of a CSV or TSV matrix
compare.py	Compare saved Python VSN output with R output
run_and_compare.py	Select matching samples, run Python VSN, save output, and compare with R
run.qmd	Quarto Shiny dashboard using the R vsn package

Requirements

For vsn2.py:

numpy
scipy

Install with:

python -m pip install numpy scipy

The matrix runner and comparison utilities also require pandas:

python -m pip install numpy scipy pandas

The Marimo dashboard declares its dependencies inline and can be run with uv.

Quick start

import numpy as np
from vsn2 import vsn_matrix

x = np.array(
    [
        [1200.0, 1500.0, np.nan],
        [2500.0, 2300.0, 2700.0],
        [5000.0, np.nan, 4800.0],
        [9000.0, 8700.0, 9200.0],
    ]
)

result = vsn_matrix(
    x,
    min_data_points_per_stratum=0,
)

transformed = result.hx

For ordinary datasets, retain the default min_data_points_per_stratum=42. The lower value above is only required for very small examples.

Input and missing values

The input must be a floating-point matrix with:

• rows representing features, proteins, or probes
• columns representing samples
• missing observations represented by np.nan

The core vsn_matrix() function does not automatically interpret zero as missing. Convert zeros before fitting when zero represents an undetected intensity:

x = np.asarray(x, dtype=float)
x[x == 0] = np.nan

Rows that are entirely missing across the selected samples are excluded from model fitting. Their positions are retained, and their transformed output remains entirely np.nan.

Rows that are only partially missing are retained. Every finite input value receives a VSN-transformed value, while missing positions remain missing.

Transformation

For sample j and feature i, the fitted natural-scale transformation is:

h_ij = asinh(exp(log_b_j) * x_ij + a_j)

The returned standard VSN output is:

hx_ij = h_ij / log(2) - hoffset

where:

hoffset = log2(2 * exp(mean(log_b)))

For affine calibration, each sample has its own offset a_j and log-scale log_b_j.

Model fitting

The implementation uses:

• profile maximum likelihood
• analytical gradients matching the R/C likelihood
• L-BFGS-B optimization through scipy.optimize.minimize
• robust least-trimmed-squares iterations
• five intensity slices per LTS iteration
• R-compatible quantile and missing-value behavior
• R-compatible starting parameters and final hoffset

Default optimizer parameters

DEFAULT_OPTIMPAR = {
    "factr": 5e7,
    "pgtol": 2e-4,
    "maxit": 60000,
    "trace": 0,
    "cvg_niter": 7,
    "cvg_eps": 0.0,
}

SciPy receives:

maxcor = 5
ftol   = factr * machine_epsilon
gtol   = pgtol
maxiter = maxit

Offsets are unbounded. Log-scale parameters are bounded to [-100, 100].

Correct R-compatible LTS partitioning

The important compatibility correction is the implementation of:

cut(rank(hmean, na.last = TRUE), breaks = 5)

When R receives a scalar number of breaks, R first creates equally spaced internal boundaries over the original rank range. R then expands only the first and last endpoints by 0.1% of that range.

The matching Python implementation is:

cut_breaks = np.linspace(rank_min, rank_max, n_slices + 1)
cut_breaks[0] -= rank_span * 0.001
cut_breaks[-1] += rank_span * 0.001

slice_labels = (
    np.searchsorted(
        cut_breaks,
        rank_hmean,
        side="left",
    )
    - 1
)
slice_labels = np.clip(slice_labels, 0, n_slices - 1)

The earlier implementation expanded the complete range before generating all boundaries. That shifted every internal boundary and changed which sparse rows entered the LTS fit.

In the validated sparse dataset:

Rows in input:                     7,816
Samples:                               7
Finite values compared:           18,023
Old Python RMSE versus R:       0.002289762308143096
Corrected Python RMSE versus R: 1.5963303695428504e-15
Corrected maximum error:        1.0658141036401503e-14

The corrected differences are at floating-point precision.

Missing values during LTS selection

Row means are calculated from available transformed values.

Residual variance is calculated with missing values propagated, matching R:

squared_residuals = (hy - hmean[:, np.newaxis]) ** 2
rvar = np.sum(squared_residuals, axis=1)

This intentionally does not use np.nansum().

A partially missing row therefore has:

finite transformed values
finite row mean
finite rank and intensity slice
NaN residual variance

Such a row is normally omitted from later trimmed fitting iterations. Rows in the lowest-intensity slice are retained by the explicit R-compatible selection rule. Regardless of LTS selection, the final fitted transformation is applied to every finite input value.

API

vsn_matrix()

vsn_matrix(
    x,
    reference=None,
    strata=None,
    lts_quantile=0.9,
    subsample=0,
    verbose=False,
    return_data=True,
    calib="affine",
    pstart=None,
    min_data_points_per_stratum=42,
    optimpar=None,
    defaultpar=None,
)

Parameters:

• x: two-dimensional NumPy array, with features in rows and samples in columns
• reference: optional fitted VsnResult used for reference normalization
• strata: optional one-based integer labels for row strata
• lts_quantile: retained LTS fraction, between 0.5 and 1.0
• subsample: number of rows sampled per stratum; 0 uses all rows
• verbose: print fitting diagnostics
• return_data: calculate and store the transformed matrix in result.hx
• calib: "affine" or "none"
• pstart: optional custom starting coefficients
• min_data_points_per_stratum: minimum number of rows required per stratum
• optimpar: overrides selected optimizer settings
• defaultpar: overrides the default optimizer dictionary

VsnResult

The result contains:

• hx: transformed matrix when return_data=True
• coefficients: fitted offsets and log-scales
• mu: transformed row means
• sigsq: estimated residual variance
• hoffset: final stratum-specific VSN offset
• strata: row stratum labels
• lbfgsb: optimizer status code; zero indicates success
• calib: calibration mode

Coefficient layout:

offsets = result.coefficients[:, :, 0]
log_scales = result.coefficients[:, :, 1]
scales = np.exp(log_scales)

Direct transformation

import numpy as np
from vsn2 import vsn2_trsf

transformed = vsn2_trsf(
    x=x,
    p=result.coefficients,
    strata=np.ones(x.shape[0], dtype=int),
    hoffset=result.hoffset,
    calib="affine",
)

Optimizer overrides

result = vsn_matrix(
    x,
    optimpar={
        "factr": 5e7,
        "pgtol": 2e-4,
        "maxit": 60000,
        "trace": 0,
        "cvg_niter": 7,
        "cvg_eps": 0.0,
    },
)

Python uses underscores in cvg_niter and cvg_eps, corresponding to R's cvg.niter and cvg.eps.

Matrix runner

Run VSN on explicitly named columns:

python run.py input.tsv output.tsv \
  --id-columns "Protein.Group,Protein.Names,Genes" \
  --intensity-columns "Sample1,Sample2,Sample3,Sample4"

Select intensity columns with regular expressions:

python run.py input.tsv output.tsv \
  --id-columns "Protein.Group,Protein.Names,Genes" \
  --intensity-regex "^F:"

The runner:

• reads CSV or TSV input
• preserves the requested identifier columns
• converts nonnumeric and nonfinite intensity entries to missing values
• treats zero as missing by default
• preserves finite values in partially missing rows
• writes VSN-transformed values
• writes fitted parameters to a separate CSV file

Use --keep-zero only when zero is a genuine measured intensity.

Comparing with R

Compare a saved Python result with R output:

python compare.py python_vsn.tsv r_vsn.csv

Optional outputs:

python compare.py python_vsn.tsv r_vsn.csv \
  --details-output comparison_by_sample.csv \
  --differences-output differences_by_row.csv

The comparison reports:

• matched intensity columns
• number of finite paired values
• RMSE
• mean absolute error
• mean Python-minus-R difference
• maximum absolute error
• per-sample Pearson correlation
• per-sample R-squared
• missing-value counts

Both fits must use the same input samples and rows for a meaningful numerical comparison. Fitting Python on all samples and comparing only a subset with an R fit made on that subset is not an equivalent test.

Generate and compare in one command

python run_and_compare.py raw_matrix.tsv r_vsn.csv \
  --python-output python_vsn_matched_samples.tsv

This workflow:

• identifies samples represented in the R output
• extracts the matching raw columns
• fits Python VSN on exactly those samples
• saves the Python-transformed matrix
• saves fitted parameters
• compares Python and R values

Marimo dashboard

Run locally:

uv run marimo run vsn2mo.py

The dashboard supports:

• CSV and TSV uploads
• intensity-column prefix selection
• editable VSN and optimizer parameters
• raw and transformed diagnostics
• timestamped result columns when output names already exist
• parameterized download filenames

The Marimo core contains the same corrected R-compatible rank-slice partitioning as vsn2.py.

Quarto Shiny dashboard

Run locally:

quarto preview run.qmd

The dashboard uses the R vsn package directly and exposes:

• calibration mode
• LTS quantile
• subsampling
• minimum rows per stratum
• factr, pgtol, and maximum iterations
• optimizer trace level
• maximum LTS iterations
• LTS convergence tolerance

It also provides raw and transformed mean-SD plots, distributions, sample boxplots, regression diagnostics, timestamp-safe result columns, and parameterized downloads.

R reference workflow

The equivalent R normalization is:

library(vsn)

x <- as.matrix(data[, intensity_columns, drop = FALSE])
x[x == 0] <- NA_real_

normalized <- vsn::justvsn(
  x,
  minDataPointsPerStratum = 0
)

For ordinary datasets, retain the package default minimum unless a lower threshold is intentional.

References

Huber W, von Heydebreck A, Sültmann H, Poustka A, Vingron M. Variance stabilization applied to microarray data calibration and to the quantification of differential expression. Bioinformatics. 2002;18(Suppl 1):S96-S104.

Bioconductor package: vsn, by Wolfgang Huber and contributors.

Current status

The following components have been independently checked in the current implementation:

• profile negative log-likelihood
• analytical gradient
• affine parameterization
• R-compatible starting parameters
• R-compatible hoffset
• L-BFGS-B memory setting
• missing-value propagation during LTS residual selection
• R-compatible quantile interpolation
• R-compatible five-slice rank partitioning
• restoration of all-missing rows in the final output
• application of the final transformation to partially observed rows

The formerly documented residual RMSE values of approximately 0.000432 and 0.00229 are not irreducible differences. The sparse-data discrepancy was caused by incorrect internal cut() boundaries and is fixed in the current vsn2.py and vsn2mo.py.