[QST] Question about how to save the files

[liangyf24]

What is your question?
Hi there,

First of all, thank you very much for maintaining this repository. I really appreciate the effort that has gone into making GPU-accelerated single-cell analysis possible.

I am writing to ask for advice regarding an issue I encountered when saving results.

Dataset context
• Dataset size: ~880,000 cells
• Workflow involves GPU-based preprocessing (via rapids_singlecell multi GPU workflow)
• One of the layers (e.g. adata.layers["scaled"]) is very large
• According to Dask diagnostics, this layer corresponds to ~94 GB of data

The issue:
When trying to save the AnnData object, both h5ad and zarr formats fail. For example, when writing to h5ad, I encounter errors like:
File ~/miniconda3/envs/rapids_singlecell/lib/python3.13/site-packages/dask/local.py:191, in start_state_from_dask(dsk, cache, sortkey, keys)
189 if task is None:
190 if dependents[key] and not cache.get(key):
--> 191 raise ValueError(
192 "Missing dependency {} for dependents {}".format(
193 key, dependents[key]
194 )
195 )
196 continue
197 elif isinstance(task, DataNode):

ValueError: Missing dependency ('scale_kernel_center-0940b9c96e20102f885e498c537f15dd', 9, 0) for dependents {('store-map-cc0b083db7973a40ae52a6e4cdb83699', 9, 0)}
Error raised while writing key 'scaled' of <class 'h5py._hl.group.Group'> to /layers

I also tried:
• Writing to Zarr instead of H5AD
• adata.layers['scaled'].compute() is not feasible in my case due to memory constraints (the layer alone is ~94 GB).

My question:
What would be the recommended way to store final results in this situation?

Specifically:
• Is it expected that large intermediate layers (e.g. scaled data) should be dropped before saving?
• Is there a recommended pattern for persisting results when working with very large GPU/Dask-backed layers?
• Are there best practices for saving AnnData objects at this scale that I may be missing?

Any guidance on how you would handle this in a production-scale workflow would be greatly appreciated.

Thank you again for your time and for maintaining this project.

[liangyf24]

I tried subsetting to highly variable genes (HVGs). After that, everything worked fine and I was able to save the .h5ad file.
However, when I tried to run the same workflow using only a single GPU, I couldn’t even complete the scaling step, as it resulted in a CUDA out-of-memory error. I would like to know how this problem can be solved.

I also have a question regarding the behavior of layers in a multi-GPU setting. I stored the raw counts in adata.layers["counts"] by assigning adata.layers["counts"] = adata.X.copy(). However, after running normalization and log1p, I noticed that this layer was modified as well. Under normal circumstances, I would expect this layer to be isolated and remain unchanged. Could you explain why this happens and how to properly preserve raw counts in this case?

[Intron7]

It would really help if you could share the exact code you’re running (or a minimal reproducible snippet). Right now I don’t have enough to reproduce the failure on my side or to pinpoint where the graph is getting into a bad state.

One suspicion is that the issue is scale/chunk related: if you’re densifying (explicitly or implicitly) and the resulting chunks are very large or irregular, the write will generate a huge store-map-* task graph and can hit edge cases. If you can share the chunking, we can test whether rechunking to a smaller target (or computing a better chunk size before writing) resolves it.

If you can post that snippet + those prints, I can try to reproduce and suggest the smallest change that makes the write robust.

[liangyf24]

Of course, here is my code and I will highlight the part where the error occurred.

In a multi-GPU setup using Dask + CuPy + RAPIDS (rapids_singlecell), I observed two unexpected behaviors:
1. A layer intended to store raw counts (adata.layers["counts"]) is modified after running normalize_total and log1p, even though no layer argument is provided.
2. Saving the .h5ad file fails unless I subset to HVGs beforehand. If HVG subsetting is applied, saving works normally.

Below is a minimal reproducible example illustrating both issues.

Multi-GPU environment

import rmm
from rmm.allocators.cupy import rmm_cupy_allocator
import cupy as cp

rmm.reinitialize(
managed_memory=True,
pool_allocator=False,
devices=[0, 1, 2],
)
cp.cuda.set_allocator(rmm_cupy_allocator)

from dask_cuda import LocalCUDACluster
from dask.distributed import Client

cluster = LocalCUDACluster(
CUDA_VISIBLE_DEVICES="0,1,2",
n_workers=3,
threads_per_worker=8,
device_memory_limit="77GB",
memory_limit="77GB",
rmm_pool_size="5GB",
rmm_managed_memory=True,
jit_unspill=True,
protocol="tcp",
)
client = Client(cluster)

Data construction

import dask.array as da
import anndata as ad
import rapids_singlecell as rsc
import numpy as np

X = da.from_array(
cp.random.poisson(1.0, size=(100_000, 5_000)).astype(cp.float32),
chunks=(33_334, 5_000),
)

adata = ad.AnnData(X)

in my dataset the data is lookalike

Store raw counts

adata.X = adata.X.persist()
adata.layers["counts"] = adata.X.copy()

Counts layer diagnostics (before normalization)

Shape: (863085, 29718)
Type: scipy.sparse.csr_matrix
Dtype: float32
Sparsity: 92.62%

Min: 1.0
Max: 15145.0
Mean: 2.9397
Median: 1.0000

Negative: 0
NaN: 0
Inf: 0

Is integer dtype: No
Are all values integers: Yes

Preprocessing

rsc.pp.normalize_total(adata, target_sum=1e4)
rsc.pp.log1p(adata)

Counts layer diagnostics (after normalization)

Shape: (863085, 29359)
Type: scipy.sparse.csr_matrix
Dtype: float32
Sparsity: 92.53%

Min: 0.2503129
Max: 6046.8872
Mean: 4.5580
Median: 2.3015

Negative: 0
NaN: 0
Inf: 0

Is integer dtype: No
Are all values integers: No
adata.layers['counts'].data
array([2.9770765, 2.9770765, 2.9770765, ..., 2.295684 , 2.295684 ,
2.295684 ], shape=(1893551058,), dtype=float32)

adata.X.data
array([1.380547 , 1.380547 , 1.380547 , ..., 1.1926137, 1.1926137,
1.1926137], shape=(1893551058,), dtype=float32)

HVG step (optional)

rsc.pp.highly_variable_genes(adata)

Uncommenting the following line avoids save error

####### adata = adata[:, adata.var.highly_variable].copy()

Rechunk before scaling

n_workers = 3
rows_per_worker = (adata.n_obs + n_workers - 1) // n_workers
adata.X = adata.X.rechunk((rows_per_worker, adata.n_vars)).persist()

Scale (dense materialization)

rsc.pp.scale(adata, max_value=10, layer="None")

persist scale chunk

Saving

Saving fails if HVGs are NOT subset

adata.write("test.h5ad")

Here is some packages info from my environment:
=== Python session info ===
Python: 3.13.11 | packaged by Anaconda, Inc. | (main, Dec 10 2025, 21:28:48) [GCC 14.3.0]
Platform: Linux-3.10.0-1160.el7.x86_64-x86_64-with-glibc2.35

=== Core packages ===
anndata: 0.12.7
scanpy: 1.11.5
rapids_singlecell: 0.13.5
dask: 2025.9.1
distributed: 2025.9.1
dask_cuda: 25.12.000
cupy: 13.6.0
rmm: 25.12.00
numpy: 2.2.6
scipy: 1.16.3
zarr: 3.1.5
h5py: 3.15.1
cuda: 12.9.4

[liangyf24]

I confirm that only normalize step will influence the layers in my environment.