GPU-accelerated single-cell analysis at scale with rapids-singlecell

Technical architecture

• CuPy sparse arrays preserve sparse data through analysis, avoiding matrix densification and repeated host-device transfers that would otherwise dominate runtime at scale.CuPy’s sparse implementations follow SciPy data structures and familiar NumPy and SciPy syntax.Source: Technical architecture · pp. 3–4
• RAPIDS libraries provide GPU tools for machine learning, graph computation, and vector search, while custom CUDA RawKernels streamline operations and reduce transfers and intermediate allocations.The custom kernels are JIT-compiled and support operations including gene selection and sparse gram matrix computation.Source: Technical architecture · p. 4
• RAPIDS Memory Manager can spill data to host memory when GPU memory is limited, allowing analyses to proceed with some speed reduction.The passage says this enables analyses of hundreds of thousands of cells on consumer-grade GPUs.Source: Technical architecture · p. 4
• Dask distributes chunked, out-of-core computation across GPUs or nodes, bypassing CuPy’s sparse-index ceiling and scaling analyses to millions of cells and beyond.The passage identifies aggregate memory and compute resources as the primary limits on this scaling.Source: Technical architecture · p. 4

Performance and reproducibility

• A 1-million-cell mouse-brain workflow took 26 seconds on a DGX B200 versus over 52 minutes on a 32-core workstation, a >120× speedup.The workflow included preprocessing, dimensionality reduction, neighborhood graph construction, and Leiden clustering.Source: Performance and reproducibility · p. 4
• RTX PRO 6000 achieved 120× and H200 110× speedups over CPU; Figure 2 identifies neighbor graph construction, clustering, and dimensionality reduction as the largest GPU gains.The benchmark used pooled memory allocation for H200, RTX PRO 6000, and B200.Source: Performance and reproducibility · p. 5, Figure 2 · p. 5
• For Harmony, tested GPUs achieved >250× speedup on 2.4-million-cell PCA embeddings; at 11.4 million cells they finished in under 25 seconds, while CPU execution was aborted after two hours.On datasets of approximately 200,000 cells, the reported speedup was roughly 100× across tested GPUs.Source: Performance and reproducibility · p. 5
• Dask-enabled out-of-core and multi-GPU execution analyzed the 100-million-cell Tahoe 100M dataset in less than 20 minutes.The integration uses chunk-based execution to handle datasets exceeding available GPU memory.Source: Performance and reproducibility · p. 5
• Stochastic methods need not be numerically identical across CPU and GPU, but downstream cluster assignments, marker-gene rankings, and cell-type annotations remain consistent.Exact numerical reproducibility is not guaranteed because of floating-point accumulation order, random number generation, and non-deterministic atomic operations.Source: Performance and reproducibility · p. 6
• GPU gains diminish below approximately 50,000 cells, for workflows dominated by small sequential operations, and on consumer GPUs with less than 8 GB VRAM.Unimplemented operations fall back to CPU and incur data-transfer overhead; mixed workflows should be profiled.Source: Performance and reproducibility · p. 6

Conclusion

• Reducing runtimes from hours to seconds makes parameter sweeps, threshold adjustments, and embedding comparisons practical for exploratory analysis.The passage also notes that acceleration can shorten the time from data generation to actionable results in clinical and translational settings.Source: Conclusion · p. 6
• The authors anticipate a growing role for GPU acceleration as single-cell atlases reach tens of millions of cells and spatial transcriptomics expands in resolution.Planned work includes spatial transcriptomics support, foundation-model integration, and broader hardware support.Source: Conclusion · pp. 6–7
• Scanpy API compatibility and a modular architecture support stepwise GPU adoption with minimal disruption to established workflows.The modular design is intended to accommodate new GPU-accelerated methods as the field evolves.Source: Conclusion · p. 7

Implementation of rapids-singlecell

• rapids-singlecell combines CuPy, RAPIDS libraries, RMM, and Dask, operating directly on AnnData within the scverse ecosystem.The stack includes cuML, cuGraph, and cuVS for machine learning, graph computation, and vector search, respectively.Source: Implementation of rapids-singlecell · p. 7
• CSR supports cell-wise operations, while CSC supports gene-wise operations; conversions are lazy and cached to reduce overhead.The passage gives normalization and QC metrics as cell-wise examples, and highly variable gene selection as a gene-wise example.Source: Implementation of rapids-singlecell · p. 7
• The implementations follow Scanpy algorithms to preserve equivalent biological results, with method differences described where GPU optimizations or library constraints require them.The passage specifies that the following methods sections focus primarily on those differing components.Source: Implementation of rapids-singlecell · p. 7

Preprocessing and normalization

• Preprocessing rescales each cell’s counts to a target sum, applies log(1 + x), and selects highly variable genes using per-gene means and variances from a sparse-matrix GPU reduction kernel.Some custom preprocessing kernels were reverse-ported to Scanpy via Numba to make those workloads multi-threaded on CPU.Source: Preprocessing and normalization · p. 7

Principal component analysis

• Sparse PCA computes the mean-centered covariance implicitly and applies centering in chunks, avoiding explicit centering or densification of the sparse matrix.The workflow obtains top-k loadings through eigen-decomposition and transforms the data as Z = XV − 1(μ⊤V).Source: Principal component analysis · p. 7
• Uncentered truncated SVD is available when centering is unnecessary, and the GPU-engineered Dask sparse PCA algorithm was adapted for CPU use and became scanpy’s default large-scale sparse solver.Source: Principal component analysis · p. 8

Dimensionality reduction and graph analysis

• cuML accelerates UMAP and t-SNE, while cuGraph accelerates k-nearest neighbor graph construction and Leiden clustering; results are stored in AnnData .obsp slots for scanpy compatibility.Source: Dimensionality reduction and graph… · p. 8

Spatial autocorrelation metrics

• Moran’s I measures global spatial autocorrelation and Geary’s C measures local dissimilarity; both are computed with CuPy for dense and sparse matrices.Source: Spatial autocorrelation metrics · p. 8

Pathway and transcription factor activity inference

• GPU implementations of mlm, ulm, waggr, AUCell, and z-score enrichment methods use CuPy matrix operations directly on GPU-resident data, avoiding host-device transfers within rapids-singlecell workflows.Source: Pathway and transcription factor… · p. 8

Batch correction with Harmony

• Harmony operates on PCA embeddings and was re-engineered to eliminate dense one-hot encoding; it achieves several hundred-fold speedups on million-cell datasets while maintaining >95% Pearson correlation for corrected principal components versus the reference baseline.Source: Batch correction with Harmony · p. 8

• The E-distance class supports the E-distance metric for transcriptomic comparisons in perturbation experiments and accepts embeddings in adata.obsm or expression matrices in adata.layers.Its API mirrors pertpy and exposes pairwise, onesided_distances, and bootstrap methods.Source: Perturbation distance computation · p. 8
• E-distance is computed as twice the mean cross-group Euclidean distance minus the mean within-group distances, with O(n·m) cross-group and O(n² + m²) within-group pairwise distances.Source: Perturbation distance computation · p. 9
• A custom tiled CUDA kernel computes pairwise distance sums on the GPU, uses upper-triangle accumulation for within-group pairs, and normalizes sums by pair counts before assembling E[a,b] = 2·d[a,b] - d[a,a] - d[b,b].The kernel uses shared-memory data reuse and skips groups with fewer than two cells.Source: Perturbation distance computation · p. 9
• Multi-GPU execution distributes group pairs according to group sizes, computes assigned distance sums independently, and aggregates partial sum matrices on GPU 0.Source: Perturbation distance computation · p. 9
• Bootstrap confidence intervals resample cells and recompute distances each iteration, using O(n) rather than O(n²) memory and GPU-generated indices without host-device transfers.Bootstrap variance is approximated with the delta method.Source: Perturbation distance computation · p. 9
• onesided_distances compares one reference group with all other groups by constructing only required cross-group and diagonal pairs, reducing computation relative to a full pairwise matrix.Pairwise results are returned as pandas DataFrames or Series with pertpy-compatible naming and indexing.Source: Perturbation distance computation · p. 9

• CUDA RawKernels implement categorical scatter-add, centered Gram matrix construction, Euclidean and cosine pairwise distances, per-gene mean and variance reductions, and matrix symmetrization.The kernels use techniques including atomic operations, shared memory, float4 vectorized loads, and direct device-memory operation to avoid host-device transfers.Source: Custom GPU kernels · pp. 9–10

Memory management and scaling

• For datasets beyond single-GPU limits, RMM oversubscription spills to host memory, while Dask partitions AnnData across GPUs or nodes.The passage estimates single-GPU limits at about 1 million cells with int32 sparse indices; chunk sizes trade lower peak memory against throughput and scheduling overhead.Source: Memory management and scaling · p. 10

Benchmarking

• Benchmarks covered a 1.3-million-cell in-core workflow and Harmony batch correction on datasets of approximately 200,000, 2.4 million, and 11.4 million cells.Timing excluded I/O, each configuration was repeated three times, and mean runtimes were reported; CPU and GPU hardware configurations are specified.Source: Benchmarking · p. 10
• Under 20 minutes, the full analysis pipeline ran on the 100M-cell Tahoe dataset using eight NVIDIA DGX B300 GPUs.The pipeline included log-normalization, HVG selection, scaling, PCA, nearest-neighbor graph construction, UMAP, and Leiden clustering.Source: Benchmarking · pp. 10–11

• All rapids-singlecell implementations were validated against scanpy CPU baselines.Source: Numerical validation · p. 11
• Normalization, HVG selection, and PCA outputs matched within floating-point tolerance using numpy.allclose defaults.Source: Numerical validation · p. 11
• For UMAP, t-SNE, and Leiden, validation checked cluster-assignment adjusted Rand index and preservation of local embedding neighborhoods.Source: Numerical validation · p. 11

• Four authors are NVIDIA employees, and one is a LaminLabs employee.Source: Competing interests · p. 11
• FJT reports consulting for four companies and ownership interests in Dermagnostix GmbH and Cellarity.Source: Competing interests · p. 11