Saving/reading .h5ad file results in altered adata.uns['log1p']

[oligomyeggo]

• [X] I have checked that this issue has not already been reported.
• [X] 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.

---

Hello. I noticed an issue when trying to run tl.rank_genes_groups, which was the result of adata.uns['log1p'] being blank (i.e., an empty dictionary, {}). I have gone through my workflow and confirmed that adata.uns['log1p'] is the expected {'base': None} after each preprocessing step, so I don't think the issue is with any of preprocessing code. However, when I save my adata object to a .h5ad file using the .write() function and then read my .h5ad file using the sc.read() function, when I check adata.uns['log1p'] it is an empty dictionary - so maybe the issue is either in the writing or reading function? I am able to manually set adata.uns['log1p'] to {'base': None} after reading the file, and can then run downstream functions like tl.rank_genes_groups without issue. I have not had this problem previously when reading .h5ad files (into either the same Jupyter notebook or into a new Jupyter notebook). Since I can manually set adata.uns['log1p'] to {'base': None}, I don't think this issue is pressing. It's just a little strange to me. Thank you for any help/advice!

Note: Please read this guide detailing how to provide the necessary information for us to reproduce your bug.

Minimal code sample (that we can copy&paste without having any data)

This can all be run in one Jupyter notebook and should produce the issue (unless it's something exclusively on my end; I've been able to reproduce the error with my own data and one of the scanpy built-in test datasets). Sorry the code chunks are broken up/a little long; I am using the scran normalization approach outlined in the single cell tutorial.

adata = sc.datasets.pbmc3k()
sc.pp.filter_genes(adata, min_cells = 1)

# scran normalization
adata_pp = adata.copy()
sc.pp.normalize_per_cell(adata_pp, counts_per_cell_after = 1e6)
sc.pp.log1p(adata_pp)
sc.pp.pca(adata_pp, n_comps = 15)
sc.pp.neighbors(adata_pp)
sc.tl.leiden(adata_pp, key_added = 'groups', resolution = 0.5)
input_groups = adata_pp.obs['groups']
data_mat = adata.X.T

%%R -i data_mat -i input_groups -o size_factors
size_factors = sizeFactors(computeSumFactors(SingleCellExperiment(list(counts = data_mat)), 
                                             clusters = input_groups, 
                                             min.mean = 0.1))

del adata_pp
adata.obs['size_factors'] = size_factors
adata.layers['counts'] = adata.X.copy()
adata.X /= adata.obs['size_factors'].values[:,None]
sc.pp.log1p(adata)
adata.X = sp.sparse.csr_matrix(adata.X)
adata.raw = adata

sc.pp.highly_variable_genes(adata, flavor = 'cell_ranger', n_top_genes = 2000)
sc.pp.pca(adata, n_comps = 50, use_highly_variable = True, svd_solver = 'arpack')
sc.pp.neighbors(adata)
sc.tl.umap(adata)
adata.uns['log1p'] # this produces: {'base': None}
adata.write('test.h5ad')

adata = sc.read('test.h5ad')
adata.uns['log1p'] # the now produces: {}
sc.tl.leiden(adata, key_added='leiden_r1.0')
sc.tl.rank_genes_groups(adata, groupby='leiden_r1.0', key_added='rank_genes_all', method='wilcoxon') # this causes an error

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
Input In [13], in <cell line: 2>()
      1 # Calculate marker genes
----> 2 sc.tl.rank_genes_groups(adata, groupby='leiden_r1.0', key_added='rank_genes_all', method='wilcoxon', use_raw=False)

File /usr/local/lib/python3.8/dist-packages/scanpy/tools/_rank_genes_groups.py:590, in rank_genes_groups(adata, groupby, use_raw, groups, reference, n_genes, rankby_abs, pts, key_added, copy, method, corr_method, tie_correct, layer, **kwds)
    580 adata.uns[key_added] = {}
    581 adata.uns[key_added]['params'] = dict(
    582     groupby=groupby,
    583     reference=reference,
   (...)
    587     corr_method=corr_method,
    588 )
--> 590 test_obj = _RankGenes(adata, groups_order, groupby, reference, use_raw, layer, pts)
    592 if check_nonnegative_integers(test_obj.X) and method != 'logreg':
    593     logg.warning(
    594         "It seems you use rank_genes_groups on the raw count data. "
    595         "Please logarithmize your data before calling rank_genes_groups."
    596     )

File /usr/local/lib/python3.8/dist-packages/scanpy/tools/_rank_genes_groups.py:93, in _RankGenes.__init__(self, adata, groups, groupby, reference, use_raw, layer, comp_pts)
     82 def __init__(
     83     self,
     84     adata,
   (...)
     90     comp_pts=False,
     91 ):
---> 93     if 'log1p' in adata.uns_keys() and adata.uns['log1p']['base'] is not None:
     94         self.expm1_func = lambda x: np.expm1(x * np.log(adata.uns['log1p']['base']))
     95     else:

KeyError: 'base'

Versions

---

anndata 0.8.0rc1 scanpy 1.8.2 sinfo 0.3.4

PIL 9.0.1 anndata2ri 1.0.6 annoy NA asttokens NA backcall 0.2.0 backports NA beta_ufunc NA binom_ufunc NA cffi 1.15.0 cycler 0.10.0 cython_runtime NA dateutil 2.8.2 debugpy 1.5.1 decorator 5.1.1 defusedxml 0.7.1 dunamai 1.11.0 entrypoints 0.4 executing 0.8.3 get_version 3.5.4 h5py 3.6.0 hypergeom_ufunc NA igraph 0.9.9 ipykernel 6.9.1 ipython_genutils 0.2.0 ipywidgets 7.6.5 jedi 0.18.1 jinja2 3.0.3 joblib 1.1.0 jupyter_server 1.13.5 kiwisolver 1.4.0 leidenalg 0.8.9 llvmlite 0.38.0 markupsafe 2.1.1 matplotlib 3.5.1 matplotlib_inline NA mpl_toolkits NA natsort 8.1.0 nbinom_ufunc NA numba 0.55.1 numexpr 2.8.1 numpy 1.21.0 packaging 21.3 pandas 1.4.1 parso 0.8.3 pexpect 4.8.0 pickleshare 0.7.5 pkg_resources NA prompt_toolkit 3.0.28 ptyprocess 0.7.0 pure_eval 0.2.2 pycparser 2.21 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.11.2 pynndescent 0.5.6 pyparsing 3.0.7 pytz 2021.3 pytz_deprecation_shim NA rpy2 3.4.2 scipy 1.8.0 scrublet NA seaborn 0.11.2 sitecustomize NA six 1.14.0 skimage 0.19.2 sklearn 1.0.2 stack_data 0.2.0 statsmodels 0.13.2 tables 3.7.0 texttable 1.6.4 threadpoolctl 3.1.0 tornado 6.1 tqdm 4.63.0 traitlets 5.1.1 typing_extensions NA tzlocal NA umap 0.5.2 wcwidth 0.2.5 zmq 22.3.0

IPython 8.1.1 jupyter_client 7.1.2 jupyter_core 4.9.2 jupyterlab 3.3.1 notebook 6.4.8

Python 3.8.10 (default, Nov 26 2021, 20:14:08) [GCC 9.3.0] Linux-5.10.76-linuxkit-x86_64-with-glibc2.29 8 logical CPU cores, x86_64

Session information updated at 2022-03-17 17:22

[ghost]

Have the same, changing anndata to 0.7.8 version helps;

[bclopesrs]

Could someone please tell me how I can save a h5ad file from a Jupyter Notebook that I've made some modifications with another name and in the same path as before?

[Gbighe]

Hi, I found the answer to this question in another Q&A, please refer to (https://bytemeta.vip/repo/scverse/scanpy/issues/2239)

[maximilianh]

Got this on scanpy 1.9.1, did what #2239 suggested, as a temporary workaround

adata.uns['log1p']["base"] = None

[sbwiecko]

It looks like the adata.uns['log1p']["base"] = None is lost at the next reading.

Using an h5ad data set saved in result_file, the following happens:

adata = sc.read_h5ad(results_file)
adata.uns['log1p']['base'] is None
#---------------------------------------------------------------------------
#KeyError                                  Traceback (most recent call last)
#/home/sebastien/workshop_scRNA-Seq_UCLA/scanpy_tut3.py in line 2
#      [153](file:///home/sebastien/workshop_scRNA-Seq_UCLA/scanpy_tut3.py?line=152) adata = sc.read(results_file)
#----> [154](file:///home/sebastien/workshop_scRNA-Seq_UCLA/scanpy_tut3.py?line=153) adata.uns['log1p']['base'] is None

#
#KeyError: 'base'

adata.uns['log1p']["base"] = None # see https://github.com/scverse/scanpy/issues/2239
adata.uns['log1p']['base'] is None
# now returns True

adata.write(results_file)
adata.uns['log1p']['base'] is None
# still returns True

adata = sc.read_h5ad(results_file) # same with adata = sc.read(results_file)
adata.uns['log1p']['base'] is None
#
#KeyError: 'base'

This was with sc.__version__ 1.9.1

[hl-xue]

Hello, for information: I encountered the same problem with scanpy=1.9.1 either with anndata=0.8.0 or 0.7.8 or 0.7.6.

[Phil-Momoka]

I found it useful by calling scanpy.pp.log1p(adata) again before the function that returns the keyerror:base

For example, in the PBMC3K tutorial, calling this function again before step 43: Comparing to a single cluster. sc.pp.log1p(adata) sc.tl.rank_genes_groups(adata, 'leiden', groups=['0'], reference='1', method='wilcoxon') sc.pl.rank_genes_groups(adata, groups=['0'], n_genes=20)

I am new to scanpy so I do not know what exactly is going on here so all I can say is "for some reason" the code works this way

[flying-sheep]

The log1p→base thing was fixed in #2546