[auto-optimize]: Bugs found while testing PyDESeq2 2ca3b74

Aperta
  • Qualcosa non va?
#491 0 commenti 0 reazioni 0 assegnatari Vedi su GitHub

Nessuno ha ancora preso questa issue.

Valutazione

Difficoltà
5/5
Tempo stimato
Più di una settimana
Idoneità per principianti
20/100
Tipo di issue
Bug
Chiarezza
Specificata chiaramente
Stato di attività
Attiva
Stack tecnologico
python
Ambito
bioinformatics, data

Direzione di ricerca

Start with the reproduction for an individual bug and inspect the named entry point, such as _fit_parametric_dispersion_trend in dds.py or DeseqDataSet.cooks_outlier; compare the reported behavior with the DESeq2 result. Run the reproduction and verify the expected output after investigating that bug. The issue lists 15 separate bugs, and fixes are largely untested, so choose a specific item before estimating its scope.

Scritto dal modello di indicizzazione a partire dal testo della issue.

Descrizione

These bugs were found with auto-optimize, a tool developed by the Programmable Genomics Laboratory at UMass Chan Medical School: they turned up while we wrote a unit-test suite for PyDESeq2 and while we optimized its speed. Many of the bugs below are likely edge cases that a user may never experience in practice, and the table after this paragraph divides them by severity. Severity is our estimate of how likely a user is to meet the bug and what it costs them: high, likely to be met in ordinary use, with wrong output or a crash; medium, needs unusual input, options or an environment, or costs time rather than correctness; low, an edge case, malformed input, a rarely used option, or a documentation mismatch. This issue was written by Claude, so its wording may be more direct than necessary.

Severity # Bug Category Effect
high 1 Parametric trend fitted on genes DESeq2 leaves out wrong output wrong output (trend coefficients, and through them every MAP dispersion)
medium 2 Cook's three-count exception applied to every design wrong output wrong output (p-values kept that DESeq2 sets to NA)
3 Genes all zero after replacement get p-value NaN wrong output wrong output (pvalue)
4 One-sided and lessAbs p-values stop at 0.5 wrong output wrong output (pvalue, padj)
6 poscounts leaves out genes whose positive counts are 1 wrong output wrong output (size factors)
7 Zero counts enter deseq2_norm_transform's median wrong output wrong output (size factors of 0; NaN and inf normalized counts)
8 Shrinkage stops early at near-Poisson dispersions wrong output wrong output (shrunk LFC)
9 Refitted genes keep first-pass MAP_dispersions wrong output wrong output (var["MAP_dispersions"])
13 Duplicated gene names make deseq2() raise crash crash
14 Series.drop(inplace=True) fails the tests on pandas 3.1 crash crash (with warnings as errors, as in the repository's tests)
low 5 prior_LFC_var is squared in the Wald test's ridge wrong output wrong output (lfcSE, stat, pvalue)
10 lowess returns 0 on an exact fit wrong output wrong output
11 irls_solver ignores min_beta when IRLS converges wrong output wrong output (non-default bounds); documentation mismatch
12 trimmed_variance(..., axis=1) centres along the wrong axis wrong output wrong output (square input); crash (other 2-D input)
15 import pydeseq2.grid_search fails in a fresh interpreter crash crash

All 15 were reproduced on PyDESeq2 main at commit 2ca3b74 (source lines refer to it), installed from source with Python 3.12.13, pandas 3.0.6, numpy 2.5.4, scipy 1.18.1, anndata 0.13.4 and formulaic 1.2.2 on Linux x86-64; bug 14 needs pandas 3.1.0 and was run with it in place of 3.0.6. DESeq2 values come from DESeq2 1.50.2 (apeglm 1.32.0, R 4.5.3), run with the R code shown under a bug on the counts its reproduction writes. Bug 13 turned up while we optimized PyDESeq2 and the others while we wrote the tests. Three more bugs we found are already reported and are not repeated here: lfc_null and alt_hypothesis applied to each coefficient instead of the contrast (#483), the prior curvature broadcast across the rows of the shrinkage Hessian (#470), and the MAP dispersion fallback dropping the prior (#471).

A limitation shared with DESeq2. One more behaviour is not in the list because DESeq2 does the same. When a gene's starting dispersion, the smaller of its rough and moments estimates, is clipped to min_disp (1e-8), fit_alpha_mle returns 1e-8 with converged=True even when the Cox-Reid adjusted negative log-likelihood is lower well above it: _fit_MoM_dispersions clips the start (dds.py:1116-1121) and fit_alpha_mle starts L-BFGS-B there (dispersions.py:211-219), where the objective is flat in log dispersion and the projected gradient is already below tolerance. In the reproduction below the fit from 1e-8 returns 1e-08 (objective 27.1724) and the fit from 0.1 returns 0.03234 (objective 26.9205); DESeq2 1.50.2's estimateDispersionsGeneEst also returns 1e-08 for this gene with the same size factors. DESeq2's documentation describes the gene-wise dispersion as the maximum likelihood estimate. Starting above the floor when the clipped start sits on it, or comparing with the grid fit, would return that estimate, but PyDESeq2 would then differ from DESeq2 on these genes.

Reproduction

It needs only PyDESeq2, not the setup below.

python - <<'EOF'
import numpy as np
from pydeseq2.dispersions import fit_alpha_mle, fit_moments_dispersions, fit_rough_dispersions
from pydeseq2.distributions import nb_nll
counts = np.array([18, 2, 10, 4, 14, 6, 19, 7, 5, 5])
size_factors = np.array([0.87, 0.82, 0.99, 1.02, 1.54, 1.30, 1.18, 1.36, 0.52, 1.07])
design = np.column_stack([np.ones(10), np.tile([0.0, 1.0], 5)])
normed = counts / size_factors
group = design[:, 1] == 1
mu = size_factors * np.where(group, normed[group].mean(), normed[~group].mean())


def objective(alpha):  # Cox-Reid adjusted negative log-likelihood
	W = mu / (1 + mu * alpha)
	return nb_nll(counts, mu, alpha) + 0.5 * np.linalg.slogdet((design.T * W) @ design)[1]


rough = fit_rough_dispersions(normed[:, None], design)[0]
moments = fit_moments_dispersions(normed[:, None], size_factors)[0]
print("starting value min(rough, moments):", min(rough, moments))
for start in [1e-8, 0.1]:
	alpha, converged = fit_alpha_mle(counts, design, mu, start, 1e-8, 10.0)
	print(f"start {start:g}: dispersion {alpha:.4g}, converged {converged}, objective {objective(alpha):.4f}")
EOF

Output:

starting value min(rough, moments): 0.0
start 1e-08: dispersion 1e-08, converged True, objective 27.1724
start 0.1: dispersion 0.03234, converged True, objective 26.9205

DESeq2 1.50.2 with the same counts and size factors:

suppressMessages(library(DESeq2))
y <- c(18, 2, 10, 4, 14, 6, 19, 7, 5, 5)
sf <- c(0.87, 0.82, 0.99, 1.02, 1.54, 1.30, 1.18, 1.36, 0.52, 1.07)
coldata <- data.frame(condition = factor(rep(c("A", "B"), 5)), row.names = paste0("s", 1:10))
dds <- DESeqDataSetFromMatrix(matrix(y, nrow = 1, dimnames = list("gene", rownames(coldata))), coldata, ~condition)
sizeFactors(dds) <- sf
dds <- estimateDispersionsGeneEst(dds, quiet = TRUE)
cat("DESeq2", as.character(packageVersion("DESeq2")), "dispGeneEst:", mcols(dds)$dispGeneEst, "iter:", mcols(dds)$dispGeneIter, "\n")
# DESeq2 1.50.2 dispGeneEst: 1e-08 iter: 10

The fixes we suggest are untested, except those for bugs 1, 6 and 8, which we checked on their reproductions only. Each section below says what happens, what should happen and which code we believe is responsible, and its collapsed block holds a reproduction with the expected and actual output; the reproductions run in a directory prepared by the setup block after the table, unless their section says otherwise.

Setup used by the reproductions

In an empty directory, with python running PyDESeq2 at 2ca3b74 (for example pip install git+https://github.com/scverse/PyDESeq2@2ca3b74c9b2b994b44c071f17ec76ad1cf1595e9). sim.py simulates negative binomial counts from a seed and builds a DeseqDataSet with n_cpus=1; load.R reads the counts a reproduction writes into R, for the DESeq2 comparisons shown under some bugs.

cat > sim.py <<'EOF'
import numpy as np
import pandas as pd
from pydeseq2.dds import DeseqDataSet


def simulate(n_genes=200, groups=(4, 4), seed=0):
	"""Counts (samples x genes) and metadata; dispersion 0.05 + 1 / mean."""
	rng = np.random.default_rng(seed)
	group = np.repeat(np.arange(len(groups)), groups)
	base = rng.lognormal(4, 1.5, n_genes)
	lfc = rng.normal(0, 0.5, (len(groups), n_genes))
	lfc[0] = 0
	mu = base * np.exp(lfc[group])
	disp = 0.05 + 1 / base
	counts = rng.negative_binomial(1 / disp, 1 / (1 + disp * mu))
	samples = [f"s{i}" for i in range(len(group))]
	counts = pd.DataFrame(counts, index=samples, columns=[f"g{j}" for j in range(n_genes)])
	metadata = pd.DataFrame({"condition": np.array(list("ABCDEF"))[group]}, index=samples)
	return counts, metadata


def dataset(counts, metadata, **kwargs):
	return DeseqDataSet(counts=counts, metadata=metadata, n_cpus=1, quiet=True, **kwargs)
EOF
cat > load.R <<'EOF'
suppressMessages(library(DESeq2))
counts <- t(as.matrix(read.csv("counts.csv", row.names = 1, check.names = FALSE)))
coldata <- read.csv("metadata.csv", row.names = 1, stringsAsFactors = TRUE)
EOF

Wrong output

1. The parametric dispersion trend is fitted on genes that DESeq2 leaves out, and a gene dropped from it never returns

One gene with a single positive count changes the parametric dispersion trend (the default fit_type), and with it the prior centre of every gene's MAP dispersion: 100 simulated genes in two groups of six give coefficients (a0, a1) of (0.0507, 1.7973), and adding a gene whose only positive count is 30 gives (0.0389, 3.0042); DESeq2 1.50.2 gives (0.0505, 1.8002) with or without that gene, because estimateDispersionsFit keeps only genes with dispGeneEst > 100 * minDisp and parametricDispersionFit starts from the curve (0.1, 1) and, before every fit, selects from all of them the genes whose ratio to the current curve lies between 1e-4 and 15; the single-sample gene is 23.7 times the starting curve and never enters. _fit_parametric_dispersion_trend (dds.py:1170-1219, which also fits vst_trend_coeffs) fits first on every non-zero gene and then drops genes outside the window without reconsidering them, so a gene that bent the first fit toward itself can stay (9.5 times the final curve here). Selecting from the genes above 100 * min_disp before each fit, from the starting curve (0.1, 1), would fix it; that loop with PyDESeq2's own dispersion_trend_gamma_glm gives (0.0507, 1.7973) on these 101 genes. Related: #482, which reported the same difference on public data and was closed by its author before review.

Reproduction 
python - <<'EOF'
from sim import simulate, dataset
counts, metadata = simulate(n_genes=100, groups=(6, 6), seed=1)
counts["single"] = [0] * 11 + [30]
counts.to_csv("counts.csv")
metadata.to_csv("metadata.csv")
for genes in [counts.columns[:-1], counts.columns]:
	dds = dataset(counts[genes], metadata)
	dds.fit_size_factors()
	dds.fit_genewise_dispersions()
	dds.fit_dispersion_trend()
	print(len(genes), "genes, trend a0, a1:", dds.uns["trend_coeffs"].round(4).tolist())
EOF

Expected: the same coefficients with and without the single-sample gene, as in DESeq2.

Actual:

100 genes, trend a0, a1: [0.0507, 1.7973]
101 genes, trend a0, a1: [0.0389, 3.0042]

DESeq2 1.50.2 on the same counts:

source("load.R")
for (genes in list(setdiff(rownames(counts), "single"), rownames(counts))) {
	dds <- estimateSizeFactors(DESeqDataSetFromMatrix(counts[genes, ], coldata, ~condition))
	dds <- estimateDispersionsFit(estimateDispersionsGeneEst(dds, quiet = TRUE), quiet = TRUE)
	print(round(attr(dispersionFunction(dds), "coefficients"), 4))
}
# asymptDisp  extraPois
#     0.0505     1.8002     (both times)
2. The three-larger-counts exception to Cook's filtering is applied to every design

In a design whose only variable has three levels, a gene above the Cook's cutoff keeps its p-value when three or more samples have a larger count than the sample with its largest Cook's distance: g_out (largest distance 38.73 in s0, cutoff 5.417, six larger counts) keeps p-value 0.000126 for B vs A; DESeq2 1.50.2 reports NA, because results() applies this exception only when the design has a single variable that is a two-level factor (length(designVars) == 1 and nlevels(var) == 2). DeseqDataSet.cooks_outlier (dds.py:1057-1061) applies it for every design; restricting it the same way would fix it.

Reproduction 
python - <<'EOF'
from scipy.stats import f
from pydeseq2.ds import DeseqStats
from sim import simulate, dataset
counts, metadata = simulate(n_genes=100, groups=(6, 6, 6), seed=5)
counts["g_out"] = [3000, 50, 52, 48, 51, 49] + [50, 47, 53, 49, 52, 51] + [3500, 3400, 3600, 3550, 3450, 3500]
counts.to_csv("counts.csv")
metadata.to_csv("metadata.csv")
dds = dataset(counts, metadata)
dds.deseq2()
cooks, x = dds.layers["cooks"][:, -1], counts["g_out"].to_numpy()
print(f"largest Cook's distance {cooks.max():.2f} in {counts.index[cooks.argmax()]}, cutoff {f.ppf(0.99, 3, 15):.3f},",
	f"samples with a larger count: {(x > x[cooks.argmax()]).sum()}")
ds = DeseqStats(dds, contrast=["condition", "B", "A"], n_cpus=1, quiet=True)
ds.summary()
print("g_out pvalue:", ds.results_df.loc["g_out", "pvalue"].round(6))
EOF

Expected: g_out pvalue: nan, as in DESeq2.

Actual:

largest Cook's distance 38.73 in s0, cutoff 5.417, samples with a larger count: 6
g_out pvalue: 0.000126

DESeq2 1.50.2 on the same counts:

source("load.R")
dds <- DESeq(DESeqDataSetFromMatrix(counts, coldata, ~condition), quiet = TRUE)
results(dds, contrast = c("condition", "B", "A"))["g_out", "pvalue"]
# [1] NA
3. Genes that are all zero after outlier replacement get p-value NaN instead of 1

A gene whose only positive count is replaced as a Cook's outlier becomes all zero. run_wald_test gives it p-value 1 (ds.py:322-325), and Cook's filtering in summary() then sets that p-value to NaN; DESeq2 1.50.2 reports 1 for the same gene, since results() sets the p-value of genes that are zero after replacement (nowZero) to 1 after its Cook's filtering. cooks_outlier reads replace_cooks (dds.py:1043-1050), where only refitted genes have their distances zeroed (dds.py:1413-1414), so a gene that is replaced but not refitted stays flagged and _cooks_filtering (ds.py:515) overwrites its p-value. Leaving new_all_zeroes_genes out of cooks_outlier, or setting their p-values to 1 after Cook's filtering, would fix it. Related: #239, which set these genes' p-values to 1 in run_wald_test.

Reproduction 
python - <<'EOF'
from pydeseq2.ds import DeseqStats
from sim import simulate, dataset
counts, metadata = simulate(n_genes=60, groups=(8, 8), seed=4)
counts["single"] = [0] * 15 + [100]
counts.to_csv("counts.csv")
metadata.to_csv("metadata.csv")
dds = dataset(counts, metadata)
dds.deseq2()
print("new all-zero genes:", list(dds.new_all_zeroes_genes))
ds = DeseqStats(dds, contrast=["condition", "B", "A"], n_cpus=1, quiet=True)
ds.run_wald_test()
print("pvalue after run_wald_test:", ds.p_values["single"])
ds.summary()
print("pvalue after summary:", ds.results_df.loc["single", "pvalue"])
EOF

Expected: pvalue after summary: 1.0, as in DESeq2.

Actual:

new all-zero genes: ['single']
pvalue after run_wald_test: 1.0
pvalue after summary: nan

DESeq2 1.50.2 on the same counts:

source("load.R")
dds <- DESeq(DESeqDataSetFromMatrix(counts, coldata, ~condition), quiet = TRUE)
results(dds, contrast = c("condition", "B", "A"))["single", "pvalue"]
# [1] 1
4. One-sided and lessAbs p-values are 0.5 for every gene on the far side of the threshold

With alt_hypothesis set to "greater", "less" or "lessAbs", every gene whose statistic is clipped to 0 gets p-value 0.5: in the reproduction 74, 72 and 54 genes do, where DESeq2's formulas on the same log2 fold changes and standard errors give 0.5498 to 1, 0.5453 to 1 and 0.5015 to 1; DESeq2's results() clips the statistic but takes the p-value from the unclipped upper tail (pfunc((LFC - T)/SE) for greater, pfunc((-T - LFC)/SE) for less, the larger of pfunc((T - LFC)/SE) and pfunc((LFC + T)/SE) for lessAbs). The repository's fixture tests/data/single_factor/r_test_res_greater.csv has p-value 0.999938 for gene3, whose statistic is 0, and test_alt_hypothesis compares p-values only where the statistic is not 0. greater, less and less_abs in glm.wald_test (glm.py:249-267) take the p-value of the clipped statistic; computing it from the unclipped (lfc - lfc_null) / wald_se would fix it (less_abs also reports a signed statistic where DESeq2's is never negative). Related: #483 (lfc_null and alt_hypothesis applied per coefficient), a separate bug in the same function.

Reproduction 
python - <<'EOF'
import numpy as np
from scipy.stats import norm
from pydeseq2.ds import DeseqStats
from sim import simulate, dataset
dds = dataset(*simulate(n_genes=100, groups=(4, 4), seed=3))
dds.deseq2()
for alternative, lfc_null in [("greater", 0.5), ("less", -0.5), ("lessAbs", 0.5)]:
	ds = DeseqStats(dds, contrast=["condition", "B", "A"], lfc_null=lfc_null, alt_hypothesis=alternative, n_cpus=1, quiet=True)
	ds.summary()
	res = ds.results_df
	lfc, se, t = res["log2FoldChange"], res["lfcSE"], abs(lfc_null)
	deseq2 = {"greater": norm.sf((lfc - t) / se), "less": norm.sf((-t - lfc) / se),
		"lessAbs": np.maximum(norm.sf((t - lfc) / se), norm.sf((lfc + t) / se))}[alternative]
	zero = (res["stat"] == 0).to_numpy()
	print(f"{alternative}: {zero.sum()} genes with stat 0, pvalue {sorted(set(res['pvalue'][zero].round(4)))},",
		f"DESeq2's formula {deseq2[zero].min():.4f} to {deseq2[zero].max():.4f}")
EOF

Expected: p-values equal to DESeq2's formula, between 0.5 and 1, for the genes with statistic 0.

Actual:

greater: 74 genes with stat 0, pvalue [0.5], DESeq2's formula 0.5498 to 1.0000
less: 72 genes with stat 0, pvalue [0.5], DESeq2's formula 0.5453 to 1.0000
lessAbs: 54 genes with stat 0, pvalue [0.5], DESeq2's formula 0.5015 to 1.0000
5. prior_LFC_var is squared in the Wald test's ridge

DeseqStats(prior_LFC_var=...) uses the reciprocal of the square of each value as the Wald test's ridge: with prior_LFC_var=[4, 4] the lowest-count gene's standard error is 1.098596, the one with ridge I / 16, not 0.841471, the one with ridge I / 4; the docstring calls it the "Prior variance for LFCs, used for ridge regularization", and a normal prior with variance v has penalty 1 / v, as in DESeq2's fitGLMsWithPrior (lambda <- 1/betaPriorVar). ds.py:293 sets ridge_factor = np.diag(1 / self.prior_LFC_var**2); np.diag(1 / self.prior_LFC_var) would match the docstring, which could also give the scale (DESeq2's prior variance is on the log2 scale and fitNbinomGLMs divides lambda by log(2)^2, while this ridge applies to natural-log coefficients).

Reproduction 
python - <<'EOF'
import numpy as np
from pydeseq2.ds import DeseqStats
from pydeseq2.glm import wald_test
from sim import simulate, dataset
dds = dataset(*simulate(n_genes=50, groups=(4, 4), seed=6))
dds.deseq2()
ds = DeseqStats(dds, contrast=["condition", "B", "A"], prior_LFC_var=np.array([4.0, 4.0]), n_cpus=1, quiet=True)
ds.run_wald_test()
gene = dds.var["_normed_means"].idxmin()
X = dds.obsm["design_matrix"].to_numpy()
lfc = dds.varm["LFC"].loc[gene].to_numpy()
mu = dds.obs["size_factors"].to_numpy() * np.exp(X @ lfc)
for v in [4.0, 16.0]:
	se = wald_test(X, dds.var.loc[gene, "dispersions"], lfc, mu, np.eye(2) / v, np.array([0.0, 1.0]), 0.0, None)[2]
	print(f"Wald SE with ridge I / {v:g}: {se:.6f}")
print(f"DeseqStats.SE with prior_LFC_var=[4, 4]: {ds.SE[gene]:.6f}")
EOF

Expected: DeseqStats.SE equal to the standard error with ridge I / 4.

Actual:

Wald SE with ridge I / 4: 0.841471
Wald SE with ridge I / 16: 1.098596
DeseqStats.SE with prior_LFC_var=[4, 4]: 1.098596
6. poscounts size factors leave out genes whose positive counts are all 1

fit_size_factors("poscounts") leaves out genes whose positive counts are all 1: the size factors are 0.4781, 1.2134, 1.3907 and 1.2395 with and without the gene ones (counts 1, 0, 1, 1); DESeq2 1.50.2 gives 0.6721, 1.3354, 0.9174 and 1.2146 with it, since its geoMeanNZ leaves out only all-zero genes, and the docstring says the method "will only exclude genes which have all-0 values". dds.py:623 keeps genes with (~np.isinf(logmeans)) & (logmeans > 0), and logmeans, which counts zeros as log 0, is 0 for such a gene. Keeping the genes with at least one positive count, (self.X > 0).any(axis=0), would fix it; it gives DESeq2's values on this matrix.

Reproduction 
python - <<'EOF'
import pandas as pd
from sim import dataset
counts = pd.DataFrame(
	{"g1": [10, 20, 30, 40], "g2": [0, 15, 25, 5], "g3": [7, 0, 9, 30],
	 "g4": [100, 80, 0, 60], "ones": [1, 0, 1, 1], "g6": [3, 9, 27, 0]},
	index=["s1", "s2", "s3", "s4"])
metadata = pd.DataFrame({"condition": ["A", "A", "B", "B"]}, index=counts.index)
counts.to_csv("counts.csv")
metadata.to_csv("metadata.csv")
for genes in [counts.columns, counts.columns.drop("ones")]:
	dds = dataset(counts[genes], metadata)
	dds.fit_size_factors("poscounts")
	print(len(genes), "genes:", dds.obs["size_factors"].round(4).tolist())
EOF

Expected: with 6 genes, DESeq2's size factors.

Actual:

6 genes: [0.4781, 1.2134, 1.3907, 1.2395]
5 genes: [0.4781, 1.2134, 1.3907, 1.2395]

DESeq2 1.50.2 on the same counts:

source("load.R")
round(sizeFactors(estimateSizeFactors(DESeqDataSetFromMatrix(counts, coldata, ~condition), type = "poscounts")), 4)
#     s1     s2     s3     s4
# 0.6721 1.3354 0.9174 1.2146
7. deseq2_norm_transform takes zero counts of new samples into the median as -inf

deseq2_norm_transform takes the median of log(count) - logmean over the fitted genes with zero counts included as -inf, so a new sample with zeros in half the fitted genes gets size factor 0 and NaN or inf normalized counts, and one with fewer zeros gets a smaller size factor (0.9913 against 0.9957 over its positive counts); the fitted genes have no zeros in the data they were fitted on, but new data, such as the counts passed to DeseqDataSet.vst_transform(counts), can, and DESeq2's estimateSizeFactorsForMatrix(counts, geoMeans = ...) takes the median over is.finite(loggeomeans) & cnts > 0. preprocessing.py:87-94 takes the median over every fitted gene; taking each sample's median over the fitted genes where its count is positive would fix it.

Reproduction 
python - <<'EOF' 2>/dev/null
import numpy as np
from pydeseq2.preprocessing import deseq2_norm_fit, deseq2_norm_transform
train = np.array([[10, 20, 30, 40, 50], [12, 18, 33, 41, 45], [9, 25, 28, 39, 55]])
new = np.array([[11, 21, 29, 42, 48], [0, 0, 0, 40, 50], [0, 20, 30, 40, 50]])
logmeans, filtered = deseq2_norm_fit(train)
normed, size_factors = deseq2_norm_transform(new, logmeans, filtered)
print("size factors:", size_factors.round(4))
print("normalized counts of sample 2:", normed[1])
with np.errstate(divide="ignore"):
	ratios = np.log(new[:, filtered]) - logmeans[filtered]
print("median over positive counts:", np.array([np.exp(np.median(r[np.isfinite(r)])) for r in ratios]).round(4))
EOF

Expected: size factors equal to the medians over positive counts, and finite normalized counts.

Actual:

size factors: [1.0096 0.     0.9913]
normalized counts of sample 2: [nan nan nan inf inf]
median over positive counts: [1.0096 1.0018 0.9957]
8. lfc_shrink returns a wrong MAP, reported as converged, at dispersions of 1e-6 and below

nbinomGLM, which lfc_shrink runs for each gene, stops away from the minimum and reports convergence at dispersions of 1e-6 and below (min_disp is 1e-8): it returns (2.1157, 1.3388) at 1e-6 and (0.8963, 0.5049) at 1e-8, where nbinomFn is lower by 1.44 and 168.6 at (2.3746, 1.0349); apeglm 1.32.0's nbinomCR method, which lfcShrink(type = "apeglm") uses, gives (2.3747, 1.0348) and (2.3747, 1.0349) for the same gene and priors, and nbinomGLM itself finds the minimum at dispersion 0.1. nbinomGLM scales the objective by its value at beta = 0 (distributions.py:216-226), which grows with size through the beta-free part of (counts + size) * logaddexp(xbeta + offset, log(size)) (1.1e8 at dispersion 1e-6), so L-BFGS-B's relative ftol of 1e-8 stops it once the unscaled objective changes by about 1. Writing that term as (counts + size) * np.log1p(np.exp(xbeta + offset) / size) gave (2.3746, 1.0349) at both dispersions in our check with the same optimizer settings. Related: #347 and #349, which tightened the tolerances for genes with large counts; this case remains at 2ca3b74.

Reproduction 
python - <<'EOF'
import numpy as np
from scipy.optimize import minimize
from pydeseq2.distributions import nbinomGLM, nbinomFn
design = np.column_stack([np.ones(8), np.repeat([0.0, 1.0], 4)])
counts = np.array([10, 12, 9, 11, 30, 28, 33, 31.0])
offset = np.zeros(8)
for disp in [1e-1, 1e-6, 1e-8]:
	size = 1 / disp
	beta, _, converged = nbinomGLM(design, counts, size, offset, 15.0, 1.0)
	objective = lambda b: nbinomFn(b, design, counts, size, offset, 15.0, 1.0, 1)
	best = minimize(objective, [2.4, 1.0], method="Nelder-Mead", options={"xatol": 1e-10, "fatol": 1e-12}).x
	print(f"dispersion {disp:.0e}: nbinomGLM {beta.round(4)} converged={converged}, objective {objective(beta):.4f};",
		f"minimum {best.round(4)}, objective {objective(best):.4f}")
EOF

Expected: nbinomGLM at the minimum for every dispersion, or converged=False.

Actual:

dispersion 1e-01: nbinomGLM [2.3997 0.9851] converged=True, objective 332.2967; minimum [2.3997 0.9851], objective 332.2967
dispersion 1e-06: nbinomGLM [2.1157 1.3388] converged=True, objective 110526000.6893; minimum [2.3746 1.0349], objective 110525999.2452
dispersion 1e-08: nbinomGLM [0.8963 0.5049] converged=True, objective 14736547433.8335; minimum [2.3745 1.035 ], objective 14736547265.1892
9. MAP_dispersions of refitted genes keep the value from before outlier replacement

After deseq2() with refit_cooks=True, a refitted gene's MAP_dispersions is still the value from before outlier replacement, while dispersions comes from the refit: g1 has MAP_dispersions 0.59995 with and without the refit, and dispersions 0.11337 after it; DESeq2 1.50.2 reports dispMAP and dispersion of 0.11336 for g1, because refitWithoutOutliers copies every column of the refit back (mcols(object)[refitReplace, idx] <- mcols(objectSub)). _refit_without_outliers (dds.py:1398-1409) copies _normed_means, LFC, genewise_dispersions, fitted_dispersions and dispersions, but not MAP_dispersions or the _genewise_converged, _MAP_converged, _LFC_converged and _outlier_genes columns; copying those too would fix it.

Reproduction 
python - <<'EOF'
from sim import simulate, dataset
counts, metadata = simulate(n_genes=100, groups=(7, 7), seed=2)
counts.loc["s0", ["g1", "g2"]] *= 20
counts.to_csv("counts.csv")
metadata.to_csv("metadata.csv")
for refit in [True, False]:
	dds = dataset(counts, metadata, refit_cooks=refit)
	dds.deseq2()
	columns = ["MAP_dispersions", "dispersions"] if refit else ["MAP_dispersions"]
	print(f"refit_cooks={refit}:", dds.var.loc[["g1", "g2", "g4"], columns].round(5).to_dict())
EOF

Expected: after the refit, MAP_dispersions equal to dispersions for g1 and g2, as in DESeq2.

Actual:

refit_cooks=True: {'MAP_dispersions': {'g1': 0.59995, 'g2': 0.41028, 'g4': 0.06673}, 'dispersions': {'g1': 0.11337, 'g2': 0.07503, 'g4': 0.06673}}
refit_cooks=False: {'MAP_dispersions': {'g1': 0.59995, 'g2': 0.41028, 'g4': 0.06673}}

DESeq2 1.50.2 on the same counts:

source("load.R")
dds <- DESeq(DESeqDataSetFromMatrix(counts, coldata, ~condition), quiet = TRUE)
round(as.data.frame(mcols(dds)[c("g1", "g2", "g4"), c("dispMAP", "dispersion")]), 5)
#    dispMAP dispersion
# g1 0.11336    0.11336
# g2 0.07500    0.07500
# g4 0.06669    0.06669
10. lowess returns 0 for points of an exact fit after robustness iterations

After a robustness iteration, stats.lowess weights points by their residual over six times the median absolute residual; on an exact fit the residuals are rounding noise, so some points get weight 0, and where a whole neighbourhood does, lstsq returns beta = 0 and the fitted value is 0, as at points 46 to 49 of the line y = 3 - 7.5x with frac=0.1, iter=5; a local linear regression of points on a line should return the line, as the same call does with iter=1. stats.py:116-121 treats only a median of exactly 0 as an exact fit; treating a median below a small multiple of the targets' scale as zero would fix it. lowess smooths the rejection counts of independent filtering (ds.py:480), whose plateaus are locally exact fits; we have not looked for a dataset where this changes padj.

Reproduction

The points hit depend on the numpy build; this output is from numpy 2.5.4.

python - <<'EOF'
import numpy as np
from pydeseq2.stats import lowess
x = np.linspace(0.1, 0.95, 50)
y = 3 - 7.5 * x
for iterations in [1, 5]:
	fit = lowess(x, y, frac=0.1, iter=iterations)
	off = np.flatnonzero(np.abs(fit - y) > 1e-6)
	print(f"iter={iterations}: points off the line {off.tolist()}")
	if len(off):
		print("fitted:", fit[off].round(4).tolist())
		print("line:", y[off].round(4).tolist())
EOF

Expected: no points off the line for either number of iterations.

Actual:

iter=1: points off the line []
iter=5: points off the line [44, 45, 46, 47, 48, 49]
fitted: [-3.373, -3.4016, 0.0, 0.0, 0.0, 0.0]
line: [-3.4745, -3.6046, -3.7347, -3.8648, -3.9949, -4.125]
11. irls_solver does not enforce min_beta when IRLS converges, and documents max_beta's default as -30

irls_solver hands over to the bounded optimizer only when a coefficient exceeds max_beta in absolute value, so with bounds that are not symmetric IRLS can converge below min_beta: min_beta=-2 returns an LFC of -3.1146 with converged=True; the docstrings of irls_solver, Inference.irls and DefaultInference.irls call min_beta the "Lower-bound on LFC" (and give max_beta's default as -30, where it is 30). The default bounds, -30 and 30, are symmetric and not affected. glm.py:100 tests np.abs(beta_hat) > max_beta; testing (beta_hat < min_beta) | (beta_hat > max_beta) would apply both bounds, and the default is misprinted at glm.py:49 and inference.py:75.

Reproduction 
python - <<'EOF'
import numpy as np
from pydeseq2.glm import irls_solver
size_factors = np.array([0.8, 1.2, 1.0, 0.9, 1.1, 1.05])
design = np.column_stack([np.ones(6), [0, 0, 0, 1, 1, 1]])
counts = np.array([48, 55, 50, 2, 3, 2])
beta, mu, hat, converged = irls_solver(counts, size_factors, design, 0.05, min_beta=-2, max_beta=30)
print("beta", beta.round(4), "converged", converged)
print([line.strip() for line in irls_solver.__doc__.splitlines() if "Upper-bound" in line])
EOF

Expected: an LFC no lower than -2, and a default of 30 in the docstring.

Actual:

beta [ 3.9448 -3.1146] converged True
['Upper-bound on LFC. (default: ``-30``).']
12. trimmed_variance(..., axis=1) centres along the wrong axis

dispersions.trimmed_variance(x, trim, axis=1) subtracts the row means without restoring the reduced axis, so a square input gives wrong values (55.37, 290.93 and 1957.97 for a 3 x 3 input with trim=0, against 272.81, 1.0067 and 1269.41 from 1.51 * np.var(x, axis=1)) and any other 2-D input raises a broadcasting ValueError; the docstring calls axis the "Dimension along which to compute variance". PyDESeq2 itself uses only axis=0, which is right. dispersions.py:111-112 computes (x - rm) ** 2 with rm = trimmed_mean(x, trim=trim, axis=axis); subtracting np.expand_dims(rm, axis) would fix it.

Reproduction 
python - <<'EOF'
import numpy as np
from pydeseq2.dispersions import trimmed_variance
x = np.array([[1.0, 2, 30], [4, 5, 6], [70, 8, 9]])
print("axis=1:", trimmed_variance(x, trim=0, axis=1))
print("transposed, axis=0:", trimmed_variance(x.T, trim=0, axis=0))
print("1.51 * np.var(x, axis=1):", 1.51 * np.var(x, axis=1))
try:
	trimmed_variance(np.arange(16.0).reshape(2, 8), axis=1)
except ValueError as error:
	print("2 x 8, axis=1:", error)
EOF

Expected: axis=1 equal to the transposed axis=0 result, and a result for the 2 x 8 input.

Actual:

axis=1: [  55.36666667  290.92666667 1957.96666667]
transposed, axis=0: [2.72806667e+02 1.00666667e+00 1.26940667e+03]
1.51 * np.var(x, axis=1): [2.72806667e+02 1.00666667e+00 1.26940667e+03]
2 x 8, axis=1: operands could not be broadcast together with shapes (2,8) (2,)

Crashes

13. deseq2() raises a pandas reindexing error when gene names are duplicated

Counts with a duplicated gene name build a DeseqDataSet with only anndata's warning that variable names are not unique, and deseq2() then raises ValueError: cannot reindex on an axis with duplicate labels in the parametric trend fit (dds.py:1171), or InvalidIndexError in fit_dispersion_prior (dds.py:830) with fit_type="mean"; DESeq2 1.50.2 renames duplicated row names with a warning ("1 duplicate rownames were renamed by adding numbers") and runs. PyDESeq2 should do the same, or refuse the names when the DeseqDataSet is built. Several steps index var by gene name (dds.py:830, dds.py:1171-1178); checking var_names.is_unique in DeseqDataSet.__init__, or calling var_names_make_unique() there with a warning, would fix it.

Reproduction 
python - <<'EOF' 2>/dev/null
from sim import simulate, dataset
counts, metadata = simulate(n_genes=50, groups=(4, 4), seed=0)
counts.columns = list(counts.columns[:-1]) + ["g0"]
for fit_type in ["parametric", "mean"]:
	dds = dataset(counts, metadata, fit_type=fit_type)
	try:
		dds.deseq2()
		print(fit_type, "ok")
	except Exception as error:
		print(f"{fit_type}: {type(error).__name__}: {error}")
EOF

Expected: deseq2() runs, or the DeseqDataSet constructor raises a clear error about the duplicated name.

Actual:

parametric: ValueError: cannot reindex on an axis with duplicate labels
mean: InvalidIndexError: Reindexing only valid with uniquely valued Index objects
14. The parametric trend fit calls Series.drop(inplace=True), deprecated in pandas 3.1, which fails the test suite

With pandas 3.1.0 (released 2026-10-22) every parametric trend fit warns Pandas4Warning: The inplace keyword in Series.drop is deprecated and will be removed in a future version (PDEP-8) from dds.py:1212 and dds.py:1216, and the repository's pytest setting filterwarnings = ["error", ...] turns it into an error, which fails 33 tests of the existing suite and the scverse integration-testing runs of #484 and #485; PyDESeq2's own code should not use an API that pandas versions it supports (pandas>=2.3) deprecate, since these calls will fail once the keyword is removed. dds.py:1185-1186, 1212-1215 and 1216-1219 call Series.drop(..., inplace=True); assigning the result instead would fix it. Two other Pandas4Warnings the suite meets under pandas 3.1.0 are not PyDESeq2's: future.infer_string, read by anndata 0.13.4 (anndata.compat.pandas_as_str, from DeseqDataSet.__init__ at dds.py:208; fixed on anndata's main branch), and DataFrame.rename(..., inplace=True) in test_design_matrix_init itself. Related: #484 and #485, whose integration-testing logs show this warning.

Reproduction

Not rerun before filing, unlike the other reproductions: it needs pandas 3.1.0; the output shown is from an earlier run.

With pandas 3.0.6 the script prints only its first line.

python - <<'EOF'
import warnings
import pandas as pd
from sim import simulate, dataset
print("pandas", pd.__version__)
dds = dataset(*simulate(n_genes=100, groups=(4, 4), seed=1))
with warnings.catch_warnings(record=True) as caught:
	warnings.simplefilter("always")
	dds.deseq2()
for where, message in sorted({(f"{w.filename.rsplit('/', 1)[-1]}:{w.lineno}", str(w.message).splitlines()[0]) for w in caught if "pydeseq2" in w.filename}):
	print(where, message)
dds = dataset(*simulate(n_genes=100, groups=(4, 4), seed=1))
warnings.simplefilter("error", DeprecationWarning)
try:
	dds.deseq2()
except DeprecationWarning as error:
	print("with warnings as errors:", type(error).__name__, str(error).splitlines()[0])
EOF

Expected: no warning from PyDESeq2's code.

Actual:

pandas 3.1.0
dds.py:1212 The inplace keyword in Series.drop is deprecated and will be removed in a future version (PDEP-8).
dds.py:1216 The inplace keyword in Series.drop is deprecated and will be removed in a future version (PDEP-8).
with warnings as errors: Pandas4Warning The inplace keyword in Series.drop is deprecated and will be removed in a future version (PDEP-8).
15. import pydeseq2.grid_search fails in a fresh interpreter

import pydeseq2.grid_search as the first PyDESeq2 import raises ImportError: cannot import name 'grid_fit_shrink_beta' from partially initialized module 'pydeseq2.grid_search' (most likely due to a circular import), while importing pydeseq2.dds first avoids it; grid_search is listed in the API reference (docs/api.md), so it should import on its own. grid_search.py:4 imports pydeseq2.distributions (for nbinomFn at grid_search.py:250) and distributions.py:8 imports grid_fit_shrink_beta from pydeseq2.grid_search; moving either import into the function that uses it would break the cycle.

Reproduction 
python -c "import pydeseq2.grid_search" 2>&1 | tail -n 1
python -c "import pydeseq2.dds, pydeseq2.grid_search; print('import after pydeseq2.dds: ok')"

Expected: both imports succeed.

Actual:

ImportError: cannot import name 'grid_fit_shrink_beta' from partially initialized module 'pydeseq2.grid_search' (most likely due to a circular import) (/path/to/pydeseq2/grid_search.py)
import after pydeseq2.dds: ok

🤖 Generated with Claude Code

Lingua principale
Python
Stelle
768
Fork
94
Merge medio
16m
PR unite (30g)
1

Preparare l'ambiente

Come iniziare

  1. Leggi tutta la issue e poi la guida ai contributi del progetto.
  2. Commenta sulla issue per dire che te ne occupi tu — evita che due persone facciano lo stesso lavoro.
  3. Fai un fork del repository e lavora su un branch.
  4. Apri una pull request che faccia riferimento al numero della issue.

Altre issue di scverse/PyDESeq2

Tutte le issue di scverse/PyDESeq2

Issue simili

Altre issue su Python

Ricevi le nuove issue nella tua casella

Un breve riepilogo di issue GitHub adatte ai principianti.