`rank_genes_groups(method="t-test", mean_in_log_space=False)` runs the t-test on exponentiated values
Nobody has claimed this yet.
- Dominant language
- Python
- Stars
- 2.6k
- Forks
- 779
- Avg merge
- 1d 4h
- Merged PRs (30d)
- 27
Description
Please make sure these conditions are met
- I have checked that this issue has not already been reported.
- I have confirmed this bug exists on the latest version of scanpy.
- (optional) I have confirmed this bug exists on the main branch of scanpy.
What happened?
Expected: sc.tl.rank_genes_groups(..., method="t-test") returns the same scores and pvals whether mean_in_log_space is True or False. The parameter is documented as choosing how logfoldchanges is computed ("Whether to do log(mean(e^x)) (False) or log(e^mean(x)) (True)"), and the 1.13.0a1 release note says the same.
Got on main (a656a33b) and on 1.13.0a2: every t statistic differs between the two settings (script below; the assertion fails). On 1.12.4, the latest release, the same script prints a difference of 0.0 and passes (1.12.x only accepts mean_in_log_space through **kwds).
Why: compute_statistics calls _basic_stats(exponentiate_values=not mean_in_log_space, need_var=True) for the t-test methods, so the per-group means and variances that t_test passes to scipy.stats.ttest_ind_from_stats are computed on expm1(X). With mean_in_log_space=False the test is therefore Welch's t on linear-scale values, although the function "expects logarithmized data" (I checked: the False scores match scipy.stats.ttest_ind(expm1(x_a), expm1(x_b), equal_var=False) to 4e-5 on pbmc68k_reduced). Shrinking the example showed that the data do not matter: any log1p matrix and any two groups reproduce it, and method="wilcoxon" is unaffected because ranks are always taken on X.
ScanpyV1 defaults mean_in_log_space=True, so default calls are fine, but the docstring recommends False as the "accurate" option and ScanpyV2Preview defaults to False, so method="t-test" under the v2 preset silently runs a different test. tests/test_rank_genes_groups.py::test_mean_in_log_space asserts only logfoldchanges, which is why this was not caught.
Proposed fix: compute the test statistics on X and, when mean_in_log_space=False, the exponentiated means for the fold change separately, as the wilcoxon branch already does, plus a test that scores/pvals are identical for both settings. A branch with fix, test and release-note fragment is ready; PR follows.
Found in a source-level correctness audit of research software (methods and harnesses: https://github.com/cindykrafft/research-software-audit/tree/claude/software-package-audit-ablwee/audits/scanpy)
Generated by Claude Code
Minimal code sample
# /// script
# requires-python = ">=3.12"
# dependencies = [
# "scanpy@git+https://github.com/scverse/scanpy.git@main",
# ]
# ///
import numpy as np
import scanpy as sc
from anndata import AnnData
rng = np.random.default_rng(0)
adata = AnnData(np.log1p(rng.poisson(1.0, (200, 50)).astype(np.float32))) # log1p of random counts
adata.obs["group"] = ["a"] * 50 + ["b"] * 150
t = {}
for mean_in_log_space in (True, False):
sc.tl.rank_genes_groups(
adata, "group", groups=["a"], method="t-test", mean_in_log_space=mean_in_log_space
)
df = sc.get.rank_genes_groups_df(adata, "a").set_index("names")
t[mean_in_log_space] = df["scores"].sort_index()
print("max |t(True) - t(False)| =", float(np.abs(t[True] - t[False]).max()))
np.testing.assert_allclose(t[True], t[False], atol=1e-4) # expected to pass: same test, same data
Error output
max |t(True) - t(False)| = 0.644777774810791
Traceback (most recent call last):
File "sc1.py", line 24, in <module>
np.testing.assert_allclose(t[True], t[False], atol=1e-4) # expected to pass: same test, same data
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File ".../numpy/testing/_private/utils.py", line 1780, in assert_allclose
assert_array_compare(compare, actual, desired, err_msg=str(err_msg),
File ".../numpy/testing/_private/utils.py", line 988, in assert_array_compare
raise AssertionError(msg)
AssertionError:
Not equal to tolerance rtol=1e-07, atol=0.0001
Mismatched elements: 50 / 50 (100%)
First 5 mismatches are at indices:
[0]: -0.806300699710846 (ACTUAL), -0.8120160698890686 (DESIRED)
[1]: -0.890629768371582 (ACTUAL), -0.8614875674247742 (DESIRED)
[2]: 0.34413713216781616 (ACTUAL), -2.9856526491845326e-15 (DESIRED)
[3]: -2.134220838546753 (ACTUAL), -2.4044182300567627 (DESIRED)
[4]: 0.5714088082313538 (ACTUAL), 0.4728967547416687 (DESIRED)
Max absolute difference among violations: 0.6447778
Max relative difference among violations: 2.526569e+14
ACTUAL: array([-0.806301, -0.89063 , 0.344137, -2.134221, 0.571409, -0.639508,
0.616735, 2.074773, -0.64225 , 0.657275, -3.797236, 1.079398,
-0.315629, 1.587186, -1.364316, -0.596969, -1.449556, -0.230574,...
DESIRED: array([-8.120161e-01, -8.614876e-01, -2.985653e-15, -2.404418e+00,
4.728968e-01, -3.775894e-01, 6.650800e-01, 1.948812e+00,
-8.193644e-01, 4.609343e-01, -4.056521e+00, 1.219001e+00,...
Versions
Python 3.12.3
scanpy 1.14.0.dev1+ga656a33b0 (main @ a656a33b)
anndata 0.13.3.post0
numpy 2.5.2
scipy 1.18.1
pandas 3.0.5
numba 0.67.0
scikit-learn 1.9.0
statsmodels 0.15.0
Contributor guide
First steps
- Read the whole issue, then the project's contributing guide.
- Comment on the issue to say you are picking it up — it saves two people doing the same work.
- Fork the repository and make your change on a branch.
- Open a pull request that references the issue number.
Research direction
Start with the compute_statistics and _basic_stats entry points, then compare the t-test and wilcoxon branches described in the issue. Extend tests/test_rank_genes_groups.py::test_mean_in_log_space with coverage for scores and pvals under both settings. Done means the t-test statistics and p-values are identical while logfoldchanges retain their intended behavior.
Written by the indexing model from the issue text.
Assessment
- Tech stack
- python
- Domain
- bioinformatics
- Issue type
- Bug
- Difficulty
- 3/5
- Estimated time
- 1-2 days
- Activity status
- Active
- Clarity
- Clearly specified
- Newbie friendliness
- 72/100