scverse / scverse/scanpy

`rank_genes_groups(method="t-test", mean_in_log_space=False)` runs the t-test on exponentiated values

Open
#4,336 0 comments 0 reactions 0 assignees View on GitHub

Nobody has claimed this yet.

Triage 🩺
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

Open the contributing guide

First steps

  1. Read the whole issue, then the project's contributing guide.
  2. Comment on the issue to say you are picking it up — it saves two people doing the same work.
  3. Fork the repository and make your change on a branch.
  4. 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

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.