[jmschrei]

These tests were written with auto-optimize, a tool developed by the Programmable Genomics Laboratory at UMass Chan Medical School: this pull request adds a pytest suite for PyDESeq2's numerical core: size factors and normalization, the variance-stabilizing transform, dispersion estimation, GLM fitting, Cook's distances and refitting, Wald tests, p-value adjustment and apeglm shrinkage, through DeseqDataSet, DeseqStats, DefaultInference and the module-level functions behind them. It holds 300 tests in nine modules under tests/, with shared helpers in tests/helpers.py (a negative binomial count simulator that records its true parameters, hand-worked inputs, and reference implementations of the quantities PyDESeq2 computes), fixtures in tests/conftest.py, DESeq2 1.50.2 outputs with the R scripts that wrote them under tests/data/unit_a/ and tests/data/unit_b/, a tests/README_NEW_TESTS.md on running and extending the suite, and a tests/TARGETS.md that lists every public class, method and function with the tests that cover it. It changes no source, build or workflow file.

Tests and line coverage by part

Part	Tests before	Tests after	Line coverage before	Line coverage after
pydeseq2/datasets.py	46	46	88.9%	88.9%
pydeseq2/dds.py	55	209	87.4%	94.2%
pydeseq2/default_inference.py	48	213	96.1%	97.4%
pydeseq2/dispersions.py	41	157	96.8%	97.9%
pydeseq2/distributions.py	40	187	98.2%	98.2%
pydeseq2/ds.py	30	76	90.0%	94.5%
pydeseq2/glm.py	36	177	81.3%	100.0%
pydeseq2/grid_search.py	1	16	44.3%	98.6%
pydeseq2/inference.py	0	0	100.0%	100.0%
pydeseq2/__init__.py	0	0	100.0%	100.0%
pydeseq2/misc.py	55	234	86.7%	90.0%
pydeseq2/plotting.py	1	1	47.6%	47.6%
pydeseq2/preprocessing.py	44	141	100.0%	100.0%
pydeseq2/stats.py	36	131	90.2%	100.0%
pydeseq2/__version__.py	0	0	100.0%	100.0%
Total	65	327	86.0%	94.2%

"Before" is the existing tests, and "after" adds this suite. Tests per part count the tests that ran at least one line of the part, from per-test coverage; a test counts once per part, so the rows do not sum to the total.

Running it

pip install -e . pytest pytest-xdist
export OMP_NUM_THREADS=1 OPENBLAS_NUM_THREADS=1 MKL_NUM_THREADS=1 MPLBACKEND=agg
python -m pytest tests/test_unit_*.py                          # the new tests, serially
python -m pytest -n 4 tests/test_unit_*.py                     # with pytest-xdist
PYDESEQ2_TEST_CPUS=4 python -m pytest tests/test_unit_*.py     # four joblib workers per PyDESeq2 call
uvx hatch run hatch-test.py3.12-stable:run-cov -n auto         # every test, as CI runs it

• The new tests take about 20 s serially on one core, with PYDESEQ2_TEST_CPUS at its default of 1. The export line keeps BLAS and OpenMP to one thread per process.
• Python 3.12 or later, as the package requires. There is no new dependency: the tests use numpy, scipy, pandas, anndata and formulaic, which PyDESeq2 already depends on, and pytest from the test dependency group; -n needs pytest-xdist, which hatch's test environments install.
• PYDESEQ2_TEST_CPUS (default 1) is the n_cpus every new test passes to PyDESeq2. At 1 the per-gene solvers run in the test process; above 1 joblib starts that many workers for each call. No new test leaves n_cpus at its default of one worker per core.
• Every new module sets pytestmark = helpers.PANDAS_DEPRECATION_FILTERS, which ignores two pandas 3.1 deprecation messages that the repository's filterwarnings = error would otherwise turn into failures: anndata reading the future.infer_string option, and the inplace keyword of Series.drop in DeseqDataSet._fit_parametric_dispersion_trend. Nothing else is filtered, and the filters match by text, so they are also valid on pandas 2. The hatch environments resolve pandas 3.1.0 today, under which the existing tests fail on main (see Checked locally).
• The hatch test environments install pytest-randomly, which shuffles the test order and reseeds NumPy's global generator before each test. The tests do not depend on either: the simulator draws from its own seeded generators.
• In the first configuration under Checked locally, no test skips.

No workflow changes. testpaths is tests, so the Test workflow's hatch run hatch-test.<env>:run-cov -v --color=yes -n auto already collects the new modules in all three hatch environments.

How expected values are derived

Expected values are derived without PyDESeq2's code wherever practical:

• Inputs: counts from the simulator in tests/helpers.py, drawn from known size factors, coefficients and dispersions on a parametric trend, with all-zero genes, genes positive in a single sample, genes with means near 1e7 and planted outliers on request. It draws from NumPy's legacy RandomState streams, which NumPy keeps fixed across versions. Also small hand-worked matrices, and simulated datasets with DESeq2 outputs under tests/data/.
• Closed forms, worked by hand or written out in the tests: median-of-ratios and poscounts size factors, the negative binomial likelihood, the variance-stabilizing transform as the integral of the inverse standard deviation, trimmed means and variances, the hat diagonal, Cook's distances and Benjamini-Hochberg adjustment.
• Reference implementations in tests/helpers.py, written from DESeq2's definitions: the Cox-Reid adjusted dispersion objective minimized over a dense grid and then by bounded Brent minimization, the ridge-penalized GLM objective minimized by Newton's method, DESeq2's parametricDispersionFit, the Wald test with its sandwich covariance, apeglm's negative log posterior minimized by trust-region Newton from several starts, Cleveland's robust lowess, and DESeq2's independent filtering. Each has its own tests in tests/test_unit_harness.py.
• R: DESeq2 1.50.2 outputs for poscounts and iterative size factors, the variance-stabilizing transform, the dispersion trend, the dispersion prior variance and the dispersion outlier flags, with the R scripts that wrote them. DESeq2 1.50.2 reproduces the repository's existing single-factor fixtures, made with DESeq2 1.34.0, to within 6e-13 relative.
• Tolerances: each comparison is made at the precision of the method that produced the value, with the reason written next to the assertion: rtol=1e-3 on dispersions from L-BFGS-B, atol=1e-3 on natural-log LFCs from IRLS, atol=2e-2 on apeglm-shrunk LFCs, and 1e-8 relative on the Wald statistics, Cook's distances and hat diagonals recomputed from a fit's outputs. Integers, masks, NaN positions, shapes, index labels and dtypes are compared exactly. Tests that refit the same data with another gene order, n_cpus or batch_size compare at the solvers' precision, and leave out the genes whose MAP dispersion takes the grid fallback in one fit only. Whether L-BFGS-B reports failure on a gene-wise fit, and the gene takes the grid fallback, depends on rounding for a few genes and so on the platform; the VST tests compare at most two such genes at one step of that grid.
• Not pinned: private attributes, the order of internal calls, joblib usage, and intermediate arrays not exposed in obs, var, varm, obsm, uns, layers or a results table. A custom Inference that implements only the abstract methods runs deseq2, summary and lfc_shrink.
• Where the documentation or DESeq2 is silent: 1 of the 190 test functions or methods pin the current output instead, and each says so in its documentation or a comment.

What it covers

Module	Tests	What it checks
test_unit_stats.py	58	DeseqStats (Wald tests, summary filters and hypotheses, contrasts, lfc_shrink) and DeseqDataSet.cond, contrast and cooks_outlier, recomputed from public fields, unchanged by n_cpus or gene order.
test_unit_glm.py	53	nb_nll, the trimmed statistics, irls_solver and its fallbacks, grid_fit_beta, fit_LFC, Cook's distances, refit, and deseq2 end to end, unchanged by gene order or n_cpus.
test_unit_dispersions.py	40	fit_alpha_mle, rough and moments dispersions, dispersion_trend, the likelihoods and their gradient, grid_fit_alpha and the DefaultInference dispersion methods, against the reference minimizer.
test_unit_harness.py	36	The helpers themselves: the simulator's reproducibility and truth, the likelihood against scipy, and each reference implementation against hand-worked cases or its optimality conditions.
test_unit_data.py	33	fit_size_factors with every method and control genes, against hand-worked values, true size factors and DESeq2 1.50.2; deseq2_norm, fitted and applied; and the misc helpers.
test_unit_dispersions_dds.py	24	The DeseqDataSet dispersion steps: gene-wise, trend and MAP dispersions, the prior, outliers, all-zero and tiny datasets, and DESeq2 1.50.2's trend, prior and outlier flags.
test_unit_stats_functions.py	24	glm.wald_test against hand-worked values and the sandwich formula, nbinomFn and nbinomGLM against apeglm's posterior and its minimizer, grid_fit_shrink_beta, and stats.lowess against Cleveland's lowess.
test_unit_data_dataset.py	16	The DeseqDataSet constructor, its errors and deprecated arguments, variables, the design matrix against formulaic, to_picklable_anndata, n_cpus, and a custom Inference running the pipeline.
test_unit_data_vst.py	16	vst, vst_fit and vst_transform with both trends, blind and with the design, against a step-by-step transform and DESeq2 1.50.2, and new samples.

Line coverage

How it was measured:

• Tool: coverage.py 7.17.0 through pytest-cov 7.1.0, with --cov=pydeseq2 --cov-context=test and the repository's [tool.coverage] settings, on Python 3.12.13 with pandas 3.0.6, numpy 2.5.4 and scipy 1.18.1, and -n 4.
• Covered: a statement of src/pydeseq2 that ran at least once; branch coverage is off.
• n_cpus 1 in both runs: the existing tests leave n_cpus at its default of one worker per core, so they ran with the CPU count Python reports set to 1, and the new tests with PYDESEQ2_TEST_CPUS=1. The per-gene solvers then run in the pytest process, where coverage.py credits them to a test.
• The two suites: the existing tests (tests/test_pydeseq2.py, tests/test_edge_cases.py and tests/test_stats.py) ran in a checkout of 2ca3b74; "this suite" is the new modules alone, on this branch.
• Tests per part: from coverage.py's per-test contexts. A test counts for a file when it runs at least one of its lines; lines that run when a module is imported count for no test, which is why inference.py, whose covered lines are its class and method definitions, shows none.
• Lines counted: each source line counts once: 1,276 executable lines in 15 files.

Tests	Lines covered
Existing tests	1,097 (86.0%)
This suite alone	1,133 (88.8%)
Existing tests and this suite	1,202 (94.2%)

What the suite does not cover

Behaviour that contradicts the documentation or DESeq2 is left untested, so that fixing it will not break the suite. tests/README_NEW_TESTS.md describes every such case under "Known issues in the current code (not tested)", together with one case the tests do check: with three or fewer residual degrees of freedom, fit_dispersion_prior warns and keeps the trigamma formula where DESeq2 1.50.2 simulates. Docstrings or comments in the tests describe 13 of them where a test avoids them; the other four (zero counts in deseq2_norm_transform, the trend's first fit, min_beta in irls_solver, and the circular import) are described only in the README. There are 17 such behaviours of main. #491 reports 13 of them with reproductions, and describes one more, gene-wise dispersions held at min_disp when the starting value sits there, as a limitation DESeq2 shares. It cites the other three, which are already reported: lfc_null and alt_hypothesis applied per coefficient (#483), the prior curvature added to every row of the shrinkage Hessian (#470), and the MAP dispersion fallback dropping the prior (#471).

Public names no test file names (3 of 74): DeseqDataSet.plot_dispersions, make_MA_plot and make_scatter.

Lines no test reaches:

• Every file runs in at least one test.
• The largest gaps in files that are reached are in pydeseq2/dds.py (29 lines not run, 94.2% covered), pydeseq2/plotting.py (22, 47.6% covered), pydeseq2/ds.py (11, 94.5% covered), pydeseq2/datasets.py (3, 88.9% covered), pydeseq2/misc.py (3, 90.0% covered) and pydeseq2/default_inference.py (2, 97.4% covered).

Platforms: Linux on x86_64 here, and on GitHub Actions in a fork with ubuntu-latest, windows-latest and macos-latest (arm64) in all three hatch test environments. With pandas below 3.1 every test passed on all three, and under pandas 3.1.0 the only failures, on every OS and in every environment, were the 49 existing tests that fail and the 2 that error on main (see Checked locally). Other architectures (Linux on aarch64, macOS on Intel) were not tried. Expected values are compared at the precision of each solver rather than bitwise, because other numpy and scipy builds round differently; the suite passes with the lowest scipy the package allows (table below).

Checked locally

The suite was run on Ubuntu 24.04.5 LTS (x86_64, AMD EPYC 9575F 64-Core Processor, 256 cores), each configuration limited to eight cores with one BLAS and OpenMP thread per process. The CI configurations use the hatch test environments CI builds (uvx hatch env create hatch-test.<env>, resolved on 2026-10-28: numpy 2.5.4, scipy 1.18.1, anndata 0.13.4, pytest 9.2.0), and the others copies of them with the versions their names give:

Configuration	Result
Python 3.12, pandas 3.0.6, serial (default)	300 passed, each of three runs (19-20 s)
Python 3.12, pandas 3.0.6, -n 4	300 passed (22 s)
Python 3.12, pandas 3.0.6, serial, PYDESEQ2_TEST_CPUS=4	300 passed (47 s)
Python 3.12, pandas 3.0.6, serial, file order (-p no:randomly)	300 passed (24 s)
CI hatch-test.py3.12-stable with pandas 3.0.6, whole tests directory	365 passed (46 s)
CI hatch-test.py3.14-stable with pandas 3.0.6, whole tests directory	365 passed (36 s)
CI hatch-test.py3.12-stable as resolved today (pandas 3.1.0), new tests	300 passed (26 s)
CI hatch-test.py3.14-stable as resolved today (pandas 3.1.0), new tests	300 passed (19 s)
CI hatch-test.py3.14-pre as resolved today (pandas 3.1.0), new tests	300 passed (20 s)
CI hatch-test.py3.12-stable as resolved today (pandas 3.1.0), whole tests directory, two pandas deprecation messages ignored	365 passed (92 s)
Python 3.12, pandas 2.3.3, whole tests directory, -n 4	365 passed (40 s)
Python 3.12, formulaic 1.0.2 (lowest allowed) and pandas 2.3.3, whole tests directory, -n 4	365 passed (30 s)
Python 3.12, scipy 1.15.0 (lowest allowed), numpy 2.4.6 and pandas 3.1.0, new tests, -n 4	300 passed (14 s)
Python 3.12, pandas 3.0.6, -n 4, merged with #479	300 passed (14 s)
Python 3.12, pandas 3.0.6, -n 4, merged with #477	300 passed (14 s)
Python 3.12, pandas 3.0.6, -n 4, merged with #472	300 passed (15 s)
Python 3.12, pandas 3.0.6, -n 4, merged with #471	300 passed (16 s)
Python 3.12, pandas 3.0.6, -n 4, merged with #470	300 passed (14 s)
Python 3.12, pandas 3.0.6, -n 4, merged with #451	300 passed (14 s)

• pandas 3.1.0, released on 2026-10-22 after the last passing scheduled CI run (2026-10-15), is what all three CI environments resolve today, and under it the existing tests fail on main (2ca3b74): 49 fail and 2 error. The repository's filterwarnings = error turns two Pandas4Warnings into errors: anndata 0.13.4 reads the deprecated future.infer_string option whenever an AnnData is built, and DeseqDataSet._fit_parametric_dispersion_trend calls Series.drop(..., inplace=True). The Test workflow will therefore fail on this pull request, in the existing tests, until main handles these two warnings. The configurations that run the whole tests/ directory in the stable environments use pandas 3.0.6, under which the existing tests pass; the new tests pass under pandas 3.1.0 in all three CI environments, and with the two messages ignored on the command line the whole directory passes there too.
• In the configurations that run the whole tests/ directory, the CPU count Python reports was capped at 4, because the existing tests leave n_cpus at its default of one worker per core and the machine has 256 cores.
• Open pull requests, as of 2026-10-31: six change code these tests reach and merge cleanly with this branch, and merged with any one of Speed up iterative size factor fitting by not slicing AnnData in the objective #479, Store dispersion trend coefficients as a dict so the dataset can be written to h5ad #477, Honor the requested coefficient in shrinkage grid fallback #472, Preserve dispersion regularization in optimizer fallback #471, Fix prior curvature in shrinkage Hessian #470 or Add tximport-style transcript-length normalization and pytximport support #451 the suite passes (table above). Store dispersion trend coefficients as a dict so the dataset can be written to h5ad #477 stores uns["trend_coeffs"] and uns["vst_trend_coeffs"] as a dict instead of a pandas Series; the tests read the two coefficients by key, through helpers.trend_coeffs, which works with either. With Store dispersion trend coefficients as a dict so the dataset can be written to h5ad #477 merged, the whole tests/ directory also passes with Python 3.12 and pandas 3.0.6 (CPU count capped at 4, as in the other runs of the whole directory): the 300 new tests and 66 existing tests, the same 66 that pass on main merged with Store dispersion trend coefficients as a dict so the dataset can be written to h5ad #477 alone. Eleven other open pull requests that change PyDESeq2's code (feat: add likelihood ratio test (LRT) to DeseqStats #452, add cupy backend #447, Feature/gpu inference #445, jax deseq2 #441, Normalization factors and DESeqDataSetFromTximport equivalent #412, fix: improve the numerical stability of lfc_shrink #370, Add Relative Log Expression (RLE) Plots #348, ENH volcano plots #317, ENH implement local fit for trend curve #234, Adding an example with real data #180 and FEAT: Implement Likelihood-Ratio test #178) conflict with main and were not tried; Chore(deps): bump the actions-deps group across 1 directory with 2 updates #473 and [pre-commit.ci] pre-commit autoupdate #449 change no Python code and no test settings.
• The repository's pre-commit hooks pass on every changed file (ruff 0.15.21 check and format, and the pre-commit-hooks checks); mypy and zizmor check files this pull request does not change and were not run. The Check Build workflow's uv build and twine check --strict pass.

🤖 Generated with Claude Code