highly variable genes + batch_key --> reciprocal condition number error

[iamsalil]

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 master branch of scanpy.

What happened?

Note: This bug seems to have been mentioned on 2/13/2022 in this discourse.scverse.org thread: https://discourse.scverse.org/t/error-in-highly-variable-gene-selection/276. However, I can't seem to access scverse.org so I can't see what the cause/resolution was.

Issue: I am running highly_variable_genes with flavor='seurat_v3'. When I do not include a batch_key, the function runs fine. When I add a batch_key, I get a numerical error.

Minimal code sample

sc.pp.highly_variable_genes(adata, flavor='seurat_v3', n_top_genes=1000, batch_key="Covariate")

Error output

----> 1 sc.pp.highly_variable_genes(adata, flavor='seurat_v3', n_top_genes=1000, batch_key="Covariate")

File ~/miniconda3/envs/singlecell/lib/python3.8/site-packages/scanpy/preprocessing/_highly_variable_genes.py:428, in highly_variable_genes(adata, layer, n_top_genes, min_disp, max_disp, min_mean, max_mean, span, n_bins, flavor, subset, inplace, batch_key, check_values)
    422     raise ValueError(
    423         '`pp.highly_variable_genes` expects an `AnnData` argument, '
    424         'pass `inplace=False` if you want to return a `pd.DataFrame`.'
    425     )
    427 if flavor == 'seurat_v3':
--> 428     return _highly_variable_genes_seurat_v3(
    429         adata,
    430         layer=layer,
    431         n_top_genes=n_top_genes,
    432         batch_key=batch_key,
    433         check_values=check_values,
    434         span=span,
    435         subset=subset,
    436         inplace=inplace,
    437     )
    439 if batch_key is None:
    440     df = _highly_variable_genes_single_batch(
    441         adata,
    442         layer=layer,
   (...)
    449         flavor=flavor,
    450     )

File ~/miniconda3/envs/singlecell/lib/python3.8/site-packages/scanpy/preprocessing/_highly_variable_genes.py:85, in _highly_variable_genes_seurat_v3(adata, layer, n_top_genes, batch_key, check_values, span, subset, inplace)
     83 x = np.log10(mean[not_const])
     84 model = loess(x, y, span=span, degree=2)
---> 85 model.fit()
     86 estimat_var[not_const] = model.outputs.fitted_values
     87 reg_std = np.sqrt(10**estimat_var)

File _loess.pyx:899, in _loess.loess.fit()

ValueError: b'reciprocal condition number  6.4708e-16'

Versions

Details

-----
anndata     0.9.1
scanpy      1.9.3
-----
PIL                         9.5.0
asttokens                   NA
backcall                    0.2.0
cffi                        1.15.1
comm                        0.1.2
cycler                      0.10.0
cython_runtime              NA
dateutil                    2.8.2
debugpy                     1.5.1
decorator                   5.1.1
defusedxml                  0.7.1
executing                   0.8.3
h5py                        3.8.0
igraph                      0.10.4
importlib_resources         NA
ipykernel                   6.19.2
ipython_genutils            0.2.0
jedi                        0.18.1
joblib                      1.2.0
jupyter_server              1.23.4
kiwisolver                  1.4.4
leidenalg                   0.9.1
llvmlite                    0.40.0
matplotlib                  3.7.1
mpl_toolkits                NA
natsort                     8.3.1
numba                       0.57.0
numpy                       1.24.3
packaging                   23.0
pandas                      2.0.1
parso                       0.8.3
pexpect                     4.8.0
pickleshare                 0.7.5
pkg_resources               NA
platformdirs                2.5.2
prompt_toolkit              3.0.36
psutil                      5.9.0
ptyprocess                  0.7.0
pure_eval                   0.2.2
pydev_ipython               NA
pydevconsole                NA
pydevd                      2.6.0
pydevd_concurrency_analyser NA
pydevd_file_utils           NA
pydevd_plugins              NA
pydevd_tracing              NA
pygments                    2.15.1
pyparsing                   3.0.9
pytz                        2022.7
scipy                       1.10.1
session_info                1.0.0
setuptools                  66.0.0
six                         1.16.0
sklearn                     1.2.2
stack_data                  0.2.0
texttable                   1.6.7
threadpoolctl               3.1.0
tornado                     6.2
traitlets                   5.7.1
typing_extensions           NA
wcwidth                     0.2.5
yaml                        6.0.1
zipp                        NA
zmq                         25.0.2
-----
IPython             8.12.0
jupyter_client      8.1.0
jupyter_core        5.3.0
jupyterlab          3.5.3
notebook            6.5.4
-----
Python 3.8.16 (default, Mar  2 2023, 03:21:46) [GCC 11.2.0]
Linux-4.15.0-212-generic-x86_64-with-glibc2.17
-----
Session information updated at 2023-09-19 02:18

[DawnChou]

I am encountering the same error. Have you fixed it?

[Aeget1000]

I encountered the same problem

[VladimirShitov]

Hi, everyone, @DawnChou , @Aeget1000 , @iamsalil . I faced the same problem as well. TLDR: choose a higher span value in sc.pp.highly_variable_genes. The default is 0.3, which caused an error for me as well. 0.5 worked fine in my case. The information below might be interesting for developers or anyone who wants to understand this error more deeply.

I got the error when using HLCA data. If scanpy developers are interested, I can point to the dataset to reproduce this problem. It is quite big, but I don't know any other example yet.

The error is caused by this line. I got it when selecting HVGs by "dataset" batch key in HLCA. Batch "Sims_2019" caused the problem. Surprisingly, relationships between mean and var as well as between x and y seemed ok:

However, something was still causing the problem. I tried to locate the error in the loess calucation in the original package but did not succeed. Anyway, this is a bit out of the scope of scanpy. Setting span to a higher value (0.5) solved the problem for me. If there is no strong argument against it, I suggest changing the default value from 0.3 to 0.5.

By the way, there is another potential bug in this function. If all the values are constant and not_const only consists of False, kernel dies when trying to run model.fit(). Maybe it is prevented previously, but in case it isn't, you might want to check that.