Single-cell RNA-seq pipelines rarely fail because PCA or Leiden clustering are conceptually difficult. They fail because data volume grows faster than workstation memory, intermediate objects become too large, and exploratory scripts are promoted into production analysis without any execution discipline.

At small scale, both Seurat and Scanpy feel interactive and forgiving. At large scale, the question changes from “how do I cluster cells?” to “how do I process hundreds of thousands of cells without crashing the node or losing reproducibility?”

Scanpy logo from the official documentation

Image source: Scanpy documentation.

What Changes at Scale

The core analysis stages remain familiar:

  • quality control
  • normalization
  • highly variable gene selection
  • dimensionality reduction
  • graph construction and clustering
  • marker analysis and annotation

What changes is the compute model. Sparse matrices, chunked I/O, and batch-aware scheduling become more important than the exact plotting function you use.

For a rough resource model, memory pressure grows approximately with cell count $n$, feature count $p$, and the number of stored representations:

$$ M \approx k \cdot n \cdot p_{effective} $$

where $k$ reflects data type, sparsity, and how many layers or embeddings you keep resident in memory.

Seurat and Scanpy Solve Similar Problems Differently

Seurat is an R package built for QC, analysis, and exploration of single-cell RNA-seq data, and Seurat v5 adds infrastructure aimed at multimodal and million-cell scale workflows. Scanpy provides a Python toolkit built with anndata and is explicitly designed for scalable gene expression analysis, including workflows that can extend beyond one million cells.

In practice:

  • use Seurat when your group already standardizes on R-based analysis, reference mapping, and existing lab workflows
  • use Scanpy when you want tighter Python ecosystem integration, easier pipeline scripting, and stronger interoperability with broader ML tooling

The engineering principle is the same in both cases: keep the count matrix sparse for as long as possible and avoid copying full objects unnecessarily.

A Cluster-Friendly Workflow Layout

Do not run the entire analysis as one notebook cell sequence. Split the workflow into coarse stages that can be retried independently.

project/
|-- data/
|-- results/
|-- scripts/
|   |-- 01_qc.py
|   |-- 02_normalize.py
|   `-- 03_cluster.py
`-- envs/
    |-- scanpy.yml
    `-- seurat.yml

That structure matters because cluster execution is mostly about restartability. If neighbor graph construction fails after two hours, you should not need to recompute raw QC metrics from scratch.

Example: Scanpy on a Shared Cluster

For large datasets, write intermediate .h5ad files between stages and submit each stage as a separate batch job. A compact preprocessing script can look like this:

import scanpy as sc

adata = sc.read_10x_mtx("data/pbmc/", var_names="gene_symbols")
adata.var_names_make_unique()

sc.pp.filter_cells(adata, min_genes=300)
sc.pp.filter_genes(adata, min_cells=10)
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, n_top_genes=3000)

adata = adata[:, adata.var.highly_variable].copy()
sc.pp.pca(adata)
sc.pp.neighbors(adata)
sc.tl.leiden(adata, resolution=0.8)

adata.write("results/pbmc_clustered.h5ad")

And the corresponding Slurm job can stay minimal:

#!/bin/bash
#SBATCH --job-name=scrna-scanpy
#SBATCH --cpus-per-task=16
#SBATCH --mem=64G
#SBATCH --time=08:00:00

module load miniconda
conda activate scanpy
python scripts/03_cluster.py

This is usually enough for mid-sized analyses. For much larger datasets, you start optimizing around on-disk representations, chunk sizes, and batchwise processing instead of full in-memory transforms.

Example: Seurat for High-Cell-Count Workflows

Seurat v5 explicitly documents scalable workflows, including sketch-based analysis and support for backends that help with million-cell-scale data. A lean clustering setup still follows the same operational pattern: create the object, normalize, find variable features, reduce dimensions, cluster, save.

library(Seurat)

counts <- Read10X(data.dir = "data/pbmc/")
obj <- CreateSeuratObject(counts = counts, min.features = 300)

obj <- NormalizeData(obj)
obj <- FindVariableFeatures(obj, nfeatures = 3000)
obj <- ScaleData(obj)
obj <- RunPCA(obj)
obj <- FindNeighbors(obj, dims = 1:30)
obj <- FindClusters(obj, resolution = 0.8)

saveRDS(obj, file = "results/pbmc_seurat.rds")

The main scaling discipline is not the function sequence. It is controlling object size, writing checkpoints, and choosing methods that do not densify the matrix behind your back.

Cluster Sizing Heuristic

For operational planning, think in terms of throughput per stage rather than one monolithic runtime. If a stage processes $c$ cells per CPU-hour, total runtime is roughly:

$$ T \approx \frac{N_{cells}}{c \cdot N_{cpu}} $$

This is crude, but useful when requesting resources on shared infrastructure. Over-requested memory leaves jobs pending forever; under-requested memory gets them killed.

Practical Rules That Prevent Pain

The highest-value rules are not exotic:

  • keep raw data immutable and write new outputs per stage
  • store sparse formats, not CSV exports of matrices
  • checkpoint after QC, normalization, and clustering
  • separate exploratory plotting from production preprocessing
  • pin package versions for Seurat, Scanpy, and their storage backends

If you follow those rules, scaling from 20,000 to 500,000 cells becomes mostly an infrastructure problem, not a methodological crisis.

References