`score_genes` does not work on a dataset loaded with `backed="r+"`

[mfiers]

Hi - excellent software, thanks! -

but I do have a problem. If i load a disk backed dataset, I cannot run sc.tl.score_genes.

Given these two sets:

ad = sc.read_h5ad('scdataset.h5ad', backed='r+')
ad2 = sc.read_h5ad('scdataset.h5ad')

and

random_genes = list(ad.var_names.to_series().sample(100))

this works perfectly:

sc.tl.score_genes(ad2, random, score_name="random100", random_state=42)

but, this:

sc.tl.score_genes(ad, random, score_name="random100", random_state=42)

yields the following error:

-----------------------------------------------------------------
ValueError                      Traceback (most recent call last)
<ipython-input-113-9cb28e089b25> in <module>
----> 1 sc.tl.score_genes(ad, random, score_name="random100", random_state=42)

~/.pyenv/versions/mfpy372/lib/python3.7/site-packages/scanpy/tools/_score_genes.py in score_genes(adata, gene_list, ctrl_size, gene_pool, n_bins, score_name, random_state, copy, use_raw)
     90     else:
     91         obs_avg = pd.Series(
---> 92             np.nanmean(_adata[:, gene_pool].X, axis=0), index=gene_pool)  # average expression of genes
     93 
     94     obs_avg = obs_avg[np.isfinite(obs_avg)] # Sometimes (and I don't know how) missing data may be there, with nansfor

<__array_function__ internals> in nanmean(*args, **kwargs)

~/.pyenv/versions/mfpy372/lib/python3.7/site-packages/numpy/lib/nanfunctions.py in nanmean(a, axis, dtype, out, keepdims)
    949     cnt = np.sum(~mask, axis=axis, dtype=np.intp, keepdims=keepdims)
    950     tot = np.sum(arr, axis=axis, dtype=dtype, out=out, keepdims=keepdims)
--> 951     avg = _divide_by_count(tot, cnt, out=out)
    952 
    953     isbad = (cnt == 0)

~/.pyenv/versions/mfpy372/lib/python3.7/site-packages/numpy/lib/nanfunctions.py in _divide_by_count(a, b, out)
    216         else:
    217             if out is None:
--> 218                 return a.dtype.type(a / b)
    219             else:
    220                 # This is questionable, but currently a numpy scalar can

ValueError: setting an array element with a sequence.

thanks
Mark

[flying-sheep]

Weird, this sure looks like a bug! Can you give use code to reproduce it with the builtin datasets?

[bkmartinjr]

previously filed on anndata: scverse/anndata#234

[mfiers]

Not sure if still required (the code in the anndata repo is close to this), but here is a minimal example:

import scanpy as sc
adata =  sc.datasets.pbmc3k()
h5adfile = 'pbmc3k.h5ad'
adata.write(h5adfile)
a1 = sc.read_h5ad(h5adfile)
a2 = sc.read_h5ad(h5adfile, backed='r+')
sc.tl.score_genes(a1, ['KIR3DL2-1', 'AL590523.1', 'CT476828.1'])
# OK
sc.tl.score_genes(a2, ['KIR3DL2-1', 'AL590523.1', 'CT476828.1'])
# ERROR 

thanks..

[Liripo]

Has anyone solved it?

[flying-sheep]

This is now handled explicitly and doesn’t fail with a confusing error anymore:

---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)
Cell In[4], [line 1](vscode-notebook-cell:?execution_count=4&line=1)
----> [1](vscode-notebook-cell:?execution_count=4&line=1) sc.tl.score_genes(a2, ['KIR3DL2-1', 'AL590523.1', 'CT476828.1'])

    [... skipping hidden 1 frame]

File ~/Dev/Python/scverse/scanpy/src/scanpy/tools/_score_genes.py:126, in score_genes(adata, gene_list, ctrl_as_ref, ctrl_size, gene_pool, n_bins, score_name, random_state, copy, use_raw, layer)
    124 if is_backed_type(adata.X) and not use_raw:
    125     msg = f"score_genes is not implemented for matrices of type {type(adata.X)}"
--> [126](https://file+.vscode-resource.vscode-cdn.net/home/phil/Dev/Python/scverse/scanpy/data/~/Dev/Python/scverse/scanpy/src/scanpy/tools/_score_genes.py:126)     raise NotImplementedError(msg)
    128 if random_state is not None:
    129     np.random.seed(random_state)

NotImplementedError: score_genes is not implemented for matrices of type <class 'anndata._core.sparse_dataset._CSRDataset'>