Scanpy vs Seurat: Choosing a Single-Cell Analysis Stack

Choose by ecosystem, not by accuracy. Scanpy (Wolf, Angerer & Theis, Genome Biology 2018) and Seurat (Hao et al., Cell 2021) implement the same canonical single-cell RNA-seq pipeline — quality control, library-size normalisation, highly variable gene selection, PCA, a nearest-neighbour graph, graph-based clustering, UMAP, differential expression — and neither is the more correct implementation of it. What differs is the language, the object you carry your analysis in, the packages that plug into that object, and, most consequentially for anyone comparing two analyses, the default value of nearly every parameter in the pipeline. Run both at their defaults on one dataset and you should expect different cluster counts and different marker lists, because you did not run the same analysis twice. Pick the stack whose downstream packages you need, fix the parameters explicitly, and record them.

The honest framing for this page: the .bloom* platform runs Python. Scanpy, scvi-tools and the other scverse packages we have verified run here as jobs.

Seurat does not run here, because the platform does not run R. Everything below about Seurat is about Seurat as an external package you would run in your own R session — which is a legitimate choice, and often the right one.

At a Glance#

Every row in the "default" half of that table is a parameter, not a property of the method. That distinction is the whole article.

The Data Structures#

AnnData: one matrix, annotated on both axes#

Scanpy operates on AnnData, an annotated-matrix container whose slots map one-to-one onto the things a single-cell analysis accumulates: the expression matrix in X, per-cell metadata in obs, per-gene metadata in var, cell embeddings in obsm, cell–cell graphs in obsp, alternative matrices in layers, and everything unstructured in uns. A single .h5ad file therefore holds the entire analysis state, which is why "send me the h5ad" is a complete handoff in Python single-cell work. The slot-by-slot breakdown is in the platform's Scanpy docs.

The practical consequence of the AnnData design is that the ecosystem is

composable by convention. scvi-tools writes its latent space into obsm['X_scVI']; Scanpy's pp.neighbors reads whatever representation you point it at; a clustering function reads obsp['connectivities']. Packages do not need to know about each other, only about the slot contract.

The Seurat object: assays, reductions and graphs#

Seurat carries a Seurat object composed of one or more assays (RNA, ADT, SCT, integrated), each with its own counts, normalised data and scaled data; plus DimReduc objects for PCA and UMAP, Graph objects for the nearest-neighbour and shared-nearest-neighbour graphs, and meta.data for per-cell annotations. The multi-assay structure is the object's real strength: multimodal data — CITE-seq protein alongside RNA, or ATAC alongside RNA — is a first-class citizen rather than an extension, which is the architecture the Seurat v4 paper was written around (Hao et al., Cell 2021).

Interoperability, and why it is not free#

Conversion between the two is routine but lossy at the edges. SeuratDisk (mojaveazure/seurat-disk) provides SaveH5Seurat, Convert and LoadH5Seurat for "rapid and on-disk conversion between h5Seurat and AnnData objects"; note that the package carries an experimental lifecycle badge and is not on CRAN. On the other side,

zellkonverter (Bioconductor) and sceasy (cellgeni/sceasy) move between .h5ad and SingleCellExperiment/Seurat objects, and anndata2ri (theislab/anndata2ri) bridges AnnData and R within one session.

What survives a round trip reliably: the count matrix, cell and gene metadata, and embeddings. What needs checking every time: graphs, nested assay layers, command history, factor levels, and anything a package stashed in uns or misc. Treat the conversion as a step that needs its own assertion — count cells, count genes, and check that the cluster labels you expected came across — rather than as a format change.

The Standard Pipeline, Step by Step#

Both packages walk the same path. Naming it in both vocabularies is the fastest way to read someone else's script in the other language.

Two structural differences hide inside that symmetry.

Scanpy mutates one object in place; Seurat accumulates slots. Scanpy functions write into the AnnData you passed, so order matters and a re-run on an already-processed object is a real hazard — log1p applied twice is not an error, just wrong. Seurat returns a modified object with each step added as a named slot plus a command log, which makes "what did I actually run?" answerable from the object itself.

Scanpy's neighbour graph is a UMAP kernel; Seurat's is a shared-NN graph. sc.pp.neighbors performs a k-NN search and then derives connectivities with a kernel selected by method, whose default is 'umap' — the fuzzy simplicial set construction from the UMAP paper (McInnes, Healy & Melville, arXiv:1802.03426). FindNeighbors builds a k-NN graph with annoy and then, by default, computes a shared-nearest-neighbour graph with prune.SNN = 1/15. Leiden and Louvain are then partitioning different graphs built from the same PCA space. This is the single most under-appreciated source of divergence between the two pipelines, and it is invisible in a side-by-side reading of the two scripts.

Same Data, Different Clusters#

This is the part worth internalising. The disagreement between a Scanpy analysis and a Seurat analysis of one dataset is almost never a bug in either package. It is the accumulated effect of a dozen defaults, each individually defensible, compounding through a pipeline where every step feeds the next.

The defaults below were read from the package sources (scverse/scanpy, satijalab/seurat) in October 2026. Defaults move between releases — in both packages — so treat this as a list of

knobs to pin, and read them off the version you are actually running.

Where each divergence actually bites#

The normalisation target. Seurat's NormalizeData divides each cell by its total counts, multiplies by scale.factor = 1e4, and natural-log transforms with log1p — a fixed target for every dataset. Scanpy's normalize_total takes target_sum = None by default, which normalises to the median counts per cell of the dataset in front of it. Pass target_sum=1e4 and the two agree; leave it alone and your normalised matrix is on a dataset-dependent scale. The log base, for what it is worth, is the one place these two agree without being asked: both apply log1p with a natural log, and Scanpy's base argument lets you change that silently for everything downstream, so it is worth stating rather than assuming.

Highly variable genes. Scanpy's flavor is not cosmetic: 'seurat' implements the dispersion-based selection of Satija et al., Nature Biotechnology 2015, 'cell_ranger' the Cell Ranger variant, and 'seurat_v3' the normalised-variance method of Stuart et al., Cell 2019. Scanpy's own documentation is explicit about the correspondence: with no batch key, flavor='seurat' mimics Seurat's FindVariableFeatures(method = 'mean.var.plot'), while flavor='seurat_v3' mimics FindVariableFeatures(method = 'vst'). Since vst is Seurat's default, the Scanpy flavor that matches a default Seurat run is 'seurat_v3' — and the two flavors expect different input: the dispersion flavors want logarithmised data, 'seurat_v3' wants raw counts. Choosing the wrong one does not error; it silently selects a different gene set, and every downstream step inherits it.

The number of PCs. FindNeighbors defaults to dims = 1:10. Ten components versus thirty or fifty is a different geometry, and the k-NN graph built in it is a different graph. There is no universally right answer here, which is exactly why it should never be left implicit.

Leiden versus Louvain. This one has a published magnitude. The Leiden paper's central claim is that Louvain can return badly connected or even disconnected communities: in the authors' experimental analysis, "up to 25% of the communities are badly connected and up to 16% are disconnected", and they prove Leiden's communities are guaranteed connected while running faster (Traag, Waltman & van Eck, Scientific Reports 2019; the Louvain algorithm itself is Blondel et al., J. Stat. Mech. 2008). Those figures are measured on the paper's benchmark and real-world networks, not on a single-cell k-NN graph, so do not transfer them to your dataset as a number. Transfer the direction: Leiden is the better-behaved algorithm, Scanpy defaults to it, and Seurat's default algorithm = 1 is the original Louvain. Seurat can run Leiden — it is algorithm = 4, and it requires the leidenalg Python module, which is a small irony worth noting when the comparison is framed as Python versus R.

Resolution. Scanpy's Leiden defaults to resolution = 1, Seurat's FindClusters to 0.8. Resolution is the most direct dial on cluster count there is. Two analyses that differ only here will report different numbers of cell types, and neither is wrong.

Differential expression. Three separate things differ. The test: Scanpy documents 't-test' as its default while Seurat's is a Wilcoxon rank-sum test, and the choice matters most for small clusters and zero-inflated genes (see Soneson & Robinson, Nature Methods 2018 for a systematic treatment of test behaviour on scRNA-seq data). The correction: Benjamini–Hochberg in Scanpy versus Bonferroni in Seurat, which is strictly more conservative, so the same p-values yield shorter Seurat marker lists. And the pre-filter: FindMarkers discards genes below logfc.threshold and min.pct before testing, whereas rank_genes_groups ranks every gene — so the two tools' outputs are not even over the same gene universe before you compare them.

What to do about it#

• Pin every parameter in the table above explicitly, in both stacks, and keep them in the script rather than the notebook.
• Compare partitions, not labels. Cluster IDs are arbitrary. Use an adjusted Rand index or a contingency table between the two labellings.
• Match one pipeline to the other deliberately if you need agreement: flavor='seurat_v3', the same PC count, the same k, the same resolution, and the same clustering algorithm. Agreement is achievable; it is just not the default.
• Record the versions. A default that changed between minor releases is indistinguishable from a different analysis.

Scale and Memory#

Historically this was Scanpy's clearest advantage, and it has narrowed.

Scanpy / AnnData. Sparse matrices throughout, plus a backed mode that reads .h5ad from disk rather than loading X into memory, which is what makes atlas-scale objects tractable in a single process. The cost is that not every function works on a backed object, so a large analysis becomes a sequence of load-compute-write steps rather than one interactive session.

Seurat. Seurat v5 addressed this directly with on-disk layers via

BPCells (bnprks/BPCells) and a sketch-based strategy: subsample a representative set of cells, run the expensive integration or clustering on the sketch, then project the result back onto the full dataset. That is the contribution of Hao et al., Nature Biotechnology 2023, the dictionary-learning paper behind v5. If your objection to Seurat was "it cannot hold my data", check the current version before repeating it.

A fair summary: for a few thousand to a few hundred thousand cells on a laptop, both are fine and the bottleneck is your patience with plotting. Above that, both need a deliberate strategy, and the strategies are different — Scanpy's is backed/chunked access, Seurat's is on-disk layers plus sketching. Neither is free, and the choice is better made on which one your collaborators can actually operate.

The orthogonal lever, in both stacks, is to stop doing PCA-and-k-NN at all and learn the representation instead. scvi-tools (Gayoso et al., Nature Biotechnology 2022; scverse/scvi-tools) trains the variational model of Lopez et al., Nature Methods 2018 on raw counts, producing a batch-corrected latent space that drops into obsm and is clustered by the same Leiden call. It is GPU work and it scales with batch size rather than with available RAM, which makes it the usual answer for multi-sample datasets where the batch effect, not the cell count, is the problem. See the platform's scvi-tools docs for the model variants.

The Surrounding Ecosystems#

Neither package is used alone, and the packages around them are the real decision.

scverse, around AnnData. scvi-tools for deep generative modelling, batch correction and label transfer; squidpy (scverse/squidpy) for spatial transcriptomics neighbourhood and image analysis; scVelo (Bergen et al., Nature Biotechnology 2020; theislab/scvelo) for RNA velocity;

muon for multimodal AnnData. They share the slot contract, so chaining them is a matter of pointing the next tool at the right obsm key.

Bioconductor and the Seurat wrappers, around the Seurat object.

Signac (stuart-lab/signac) for scATAC-seq, Azimuth (satijalab/azimuth) for reference-based annotation, plus the much larger Bioconductor universe reachable through SingleCellExperiment. If your downstream analysis is statistical genomics — pseudobulk differential expression with established count models, annotation packages, curated reference data — R has the deeper bench, and it is not close.

The honest asymmetry. New single-cell methods ship in Python first more often than not, because that is where the deep-learning tooling is. Mature statistical methods have R implementations that have been stress-tested for a decade. That is the trade, and it is why the strongest single-cell groups run both and convert.

For integration method choice specifically, the benchmark to read before picking is Luecken et al., Nature Methods 2021, which evaluates atlas-level integration methods on both batch correction and biological-signal conservation rather than on either alone.

How to Choose#

Choose Scanpy when#

• Your downstream work is deep learning. scvi-tools, scANVI label transfer, totalVI, and essentially every new model ships as a Python package that expects an AnnData.
• You are working at atlas scale and want backed, chunked access rather than a sketching strategy.
• Your pipeline is already Python — an ML stack, a workflow engine, a compute platform that runs Python jobs. Keeping one language removes a conversion step that has to be asserted every time.
• You need spatial neighbourhood analysis through squidpy, or RNA velocity through scVelo.

Choose Seurat when#

• Your collaborators write R. This is the single best reason on either list. A pipeline nobody else on the team can read is a liability that no benchmark fixes.
• Your data is multimodal — CITE-seq, ATAC + RNA — and you want the assay structure and the v4/v5 integration machinery rather than assembling it.
• Your downstream analysis is statistical genomics that lives in Bioconductor: pseudobulk count models, curated references, annotation packages.
• You want publication-ready plots with less work. Seurat's ggplot2-based defaults need less tuning than Matplotlib's. This is a real cost difference, not a matter of taste.

Run both when#

You are reproducing a published analysis, or the cluster structure is the claim. Two pipelines agreeing on a partition is far stronger evidence than one pipeline at its defaults — and when they disagree, the table above tells you where to look first.

Limitations to Keep in View#

• Neither package validates your clusters. A graph-based clustering at any resolution returns clusters. Whether they are cell types, states, cell-cycle phase or an uncorrected batch is a biological judgement, and it needs marker evidence or an orthogonal measurement — surface protein, flow cytometry, imaging — not a different resolution.
• UMAP distances are not distances. The UMAP paper argues the method preserves more global structure than t-SNE and runs faster, which is not a claim that inter-cluster gaps on the plot are quantitative (McInnes et al., arXiv:1802.03426). Do not read cluster separation off the embedding; read it off the graph or the differential expression.
• Default QC thresholds are not biology. A mitochondrial-percent cutoff tuned for PBMCs will delete real cells from cardiac or hepatic tissue, where high mitochondrial content is expected.
• Cluster-level differential expression inflates significance. Testing cells as independent replicates after using the same cells to define the clusters is a known circularity in both packages' default marker workflows. For condition-level claims across samples, pseudobulk per sample.
• Conversion is a step that can fail quietly. See the interoperability section — assert counts and labels after every round trip.

How to Run It on .bloom*#

The platform runs Scanpy as a job, not as a notebook you maintain. .bloom* does not run R, so there is no Seurat path here — if Seurat is the right tool for your team, run it in your own R session and bring the .h5ad across.

What the Scanpy job actually computes. You supply one AnnData .h5ad (S3 URI or HTTPS URL, and the pipeline validates the .h5ad extension) and choose one workflow:

Normalisation, HVG selection, PCA, the neighbour graph and UMAP are steps

inside those workflows, not separately selectable stages — there is no workflow that returns only a PCA. Trajectory inference and gene-regulatory-network inference are not part of this job; fate and pseudotime are a separate CellRank tool and RNA velocity is the scvelo_velocity workflow.

The defaults this job ships with, which are the ones that make your results reproducible and which you should override deliberately: n_top_genes = 2000, hvg_flavor = seurat_v3, n_pcs = 50, n_neighbors = 15, resolution = 0.5, clustering_method = leiden (louvain is available), de_method = wilcoxon, min_genes = 200, min_cells = 3. Note that the resolution default is lower than either package's upstream default, and that seurat_v3 is the flavor that matches a default Seurat run — both choices are deliberate, and both are exactly the kind of thing this article argues you should state rather than inherit. The full parameter list, including the QC cutoffs, is in the Scanpy inputs reference.

What comes back: a processed .h5ad carrying the embeddings and cluster assignments, a UMAP coordinate CSV, a differential-expression CSV, a marker-gene table per cluster, and the figures — detailed in the Scanpy outputs reference. The processed .h5ad is the handoff: feed it to scvi-tools for a batch-corrected latent space, or convert it for an R session.

Access. scanpy_analyze is a free-tier tool and runs without an account. The scverse extension workflows on the same job — Harmony, Scanorama and ComBat batch correction, CellTypist annotation, AUCell and COSG scoring, Hotspot modules, CytoTrace2, scVelo — and scvi_analyze require a signed-in account. Start with the Scanpy quickstart, and the troubleshooting page if QC empties your object.

For neighbouring work: genomics tools for gene, pathway and variant context on the markers you find, spatial omics when the cells have coordinates, and which tool? to match a goal to a job.

Key Takeaways#

• There is no accuracy winner. Scanpy and Seurat implement the same pipeline; the choice is language, object model and downstream ecosystem.
• The same data gives different clusters because the defaults differ — normalisation target, HVG flavor, PC count, k, graph kernel, clustering algorithm, resolution, DE test, and p-value correction all diverge.
• Leiden is the better-behaved clustering algorithm, which is why Scanpy defaults to it and why Seurat exposes it as algorithm = 4.
• Scanpy's scale advantage has narrowed: Seurat v5 added on-disk layers and sketching. Check the version before repeating the old claim.
• Interop works, with assertions. SeuratDisk, zellkonverter and sceasy move data across; graphs and nested slots are where it frays.
• The platform runs the Python side. Scanpy and scvi-tools run as jobs here; Seurat is an external package, and that is a statement about runtimes, not about quality.

Pin your parameters, record your versions, and compare partitions rather than labels. Then run the Scanpy pipeline — QC through Leiden and marker genes, on your own .h5ad, with the defaults written down — on the .bloom* platform.

References#

1. Wolf, F.A., Angerer, P. & Theis, F.J. "SCANPY: large-scale single-cell gene expression data analysis." Genome Biology 19, 15 (2018). doi:10.1186/s13059-017-1382-0 · code: scverse/scanpy
2. Hao, Y., Hao, S., Andersen-Nissen, E. et al. "Integrated analysis of multimodal single-cell data." Cell 184, 3573–3587 (2021). doi:10.1016/j.cell.2021.04.048 · code: satijalab/seurat
3. Hao, Y., Stuart, T., Kowalski, M.H. et al. "Dictionary learning for integrative, multimodal and scalable single-cell analysis." Nature Biotechnology 42, 293–304 (2023). doi:10.1038/s41587-023-01767-y
4. Stuart, T., Butler, A., Hoffman, P. et al. "Comprehensive Integration of Single-Cell Data." Cell 177, 1888–1902 (2019). doi:10.1016/j.cell.2019.05.031
5. Satija, R., Farrell, J.A., Gennert, D. et al. "Spatial reconstruction of single-cell gene expression data." Nature Biotechnology 33, 495–502 (2015). doi:10.1038/nbt.3192
6. Traag, V.A., Waltman, L. & van Eck, N.J. "From Louvain to Leiden: guaranteeing well-connected communities." Scientific Reports 9, 5233 (2019). doi:10.1038/s41598-019-41695-z
7. Blondel, V.D., Guillaume, J.-L., Lambiotte, R. & Lefebvre, E. "Fast unfolding of communities in large networks." Journal of Statistical Mechanics: Theory and Experiment 2008, P10008 (2008). doi:10.1088/1742-5468/2008/10/P10008
8. McInnes, L., Healy, J. & Melville, J. "UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction." arXiv (2018). arXiv:1802.03426 · code: lmcinnes/umap
9. Gayoso, A., Lopez, R., Xing, G. et al. "A Python library for probabilistic analysis of single-cell omics data." Nature Biotechnology 40, 163–166 (2022). doi:10.1038/s41587-021-01206-w · code: scverse/scvi-tools
10. Lopez, R., Regier, J., Cole, M.B. et al. "Deep generative modeling for single-cell transcriptomics." Nature Methods 15, 1053–1058 (2018). doi:10.1038/s41592-018-0229-2
11. Luecken, M.D., Büttner, M., Chaichoompu, K. et al. "Benchmarking atlas-level data integration in single-cell genomics." Nature Methods 19, 41–50 (2021). doi:10.1038/s41592-021-01336-8
12. Soneson, C. & Robinson, M.D. "Bias, robustness and scalability in single-cell differential expression analysis." Nature Methods 15, 255–261 (2018). doi:10.1038/nmeth.4612
13. Bergen, V., Lange, M., Peidli, S. et al. "Generalizing RNA velocity to transient cell states through dynamical modeling." Nature Biotechnology 38, 1408–1414 (2020). doi:10.1038/s41587-020-0591-3 · code: theislab/scvelo