Metadata-Version: 2.4
Name: vsn
Version: 0.1.0
Summary: Pure Python implementation of Variance Stabilization and Normalization
Project-URL: Repository, https://github.com/animesh/vsn
Project-URL: Issues, https://github.com/animesh/vsn/issues
Requires-Python: >=3.9
Description-Content-Type: text/markdown
Requires-Dist: numpy>=1.23
Requires-Dist: scipy>=1.9

# 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`:

```text
numpy
scipy
```

Install with:

```bash
python -m pip install numpy scipy
```

The matrix runner and comparison utilities also require pandas:

```bash
python -m pip install numpy scipy pandas
```

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

## Quick start

```python
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:

```python
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:

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

The returned standard VSN output is:

```text
hx_ij = h_ij / log(2) - hoffset
```

where:

```text
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

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

SciPy receives:

```text
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:

```r
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:

```python
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:

```text
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:

```python
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:

```text
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()`

```python
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:

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

### Direct transformation

```python
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

```python
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:

```bash
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:

```bash
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:

```bash
python compare.py python_vsn.tsv r_vsn.csv
```

Optional outputs:

```bash
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

```bash
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:

```bash
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:

```bash
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:

```r
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`.
