Truecell Tutorials

Eighteen end-to-end tutorials covering increasingly complex single-cell analysis workflows, each pairing R Seurat code side-by-side with the equivalent Python Truecell code.

---

Tutorial Overview

#	Tutorial	Dataset	Key Concepts	Complexity
1	PBMC 3k — Guided Clustering	3,000 PBMCs · 10x Genomics (2016)	QC · Normalization · HVG/VST · PCA · Louvain · UMAP · Markers. Compared end to end, both sides running their own pipeline: the same 2,638 barcodes survive QC, 1,998/2,000 variable features shared, PCA matched |r| 0.9988, clusters at ARI 0.899 (8 vs 9 — one 32-cell DC population), and on the two clusters whose cells match exactly the marker tables are identical gene sets agreeing to 4.6e-14	Beginner
2	PBMC 8k — Advanced Subclustering	8,400 PBMCs · GRCh38 · 10x Genomics	All of Tutorial 1 + subclustering, hierarchical cell-type gating, T/NK annotation. Both stages compared by barcode: same 7,475 cells after QC, global clusters at ARI 0.977, and the T/NK compartment handed to stage 2 matches at Jaccard 0.9991 (4,631 of 4,635 cells) — subclusters then at ARI 0.916 and subset labels 98.2% concordant	Intermediate
3	CBMC CITE-seq — Multimodal	8,600 CBMCs · RNA + 13 surface proteins	Multi-assay objects · CLR normalization · Protein feature plots · RNA-protein comparison · WNN joint clustering. Compared per protein and per cell: CLR to 4.2e-15, WNN modality weights at Pearson 0.9847 over 8,617 shared barcodes, cell-type labels 99.29% concordant, all cluster counts identical. Settled the long-open progenitor question — it was a labelling difference, not a WNN one	Advanced
4	PBMC 3k — SCTransform	3,000 PBMCs · 10x Genomics (2016)	Regularized NB normalization · Pearson residuals · vars.to.regress · 30-PC workflow · SCT-vs-LogNormalize. The fitted model is compared per gene against Seurat's SCTModel feature attributes — detection_rate/gmean to machine precision, intercept and theta at Spearman 1.0000, the 3,848 non-overdispersed genes exactly the same set, residual variance at 0.9986; both arms now agree on cluster count (12 and 11)	Advanced
5	Xenium — Spatial (R vs Python)	36,602 cells · 10x Xenium mouse brain (CTX+HP)	load_xenium · ImageDimPlot/ImageFeaturePlot · nearest-neighbour distance · local density · BuildNicheAssay · composition_test — verified to 8 s.f. vs R Seurat	Spatial
6	Cell Hashing — Demultiplexing	39,842 cells · 8 HTOs · GSE108313 (human+mouse)	Hashtag assay · CLR (margin 1) · HTODemux ↔ hto_demux · MULTIseqDemux ↔ multiseq_demux · cross-species doublet ground truth — 99.81 % call-concordant with R	Advanced
7	Mixscape — Pooled CRISPR Screen	20,729 cells · 25 guides + NT · GSE153056 (THP-1 ECCITE-seq)	Perturbation signature (CalcPerturbSig) · KO-vs-escaper mixture (RunMixscape) · guide-separating LDA (MixscapeLDA) · PlotPerturbScore / MixscapeHeatmap — 97.45 % per-cell call-concordant with R	Advanced
8	Batch Integration — Harmony/CCA/RPCA	13,999 cells · CTRL/STIM · ifnb (Kang 2018)	Batch correction (RunHarmony ↔ run_harmony) · CCA/RPCA anchors (IntegrateLayers ↔ integrate_layers) · silhouette + cluster-ARI scoring — Harmony/CCA/RPCA all reach batch-mix 0.991; caught eight bugs, all fixed (a crash, an under-integration, and — found chasing what looked like the remaining implementation gap — IntegrateLayers silently running the v4 algorithm, a randomized-SVD PC drift in run_pca, and four in the neighbour graphs). The last gap standing turned out not to be a defect: Seurat's deeper modularity search buys 0.17 % by splitting CD14 Mono on batch, so truecell's coarser partition scores ARI 0.92 to the annotations against Seurat's 0.74	Advanced
9	Reference Mapping — Label Transfer	4,679 cells · celseq2→smartseq2 · panc8 (Baron 2016)	Cross-technology annotation transfer (FindTransferAnchors ↔ find_transfer_anchors) · TransferData ↔ transfer_data · MapQuery/ProjectUMAP ↔ map_query/project_umap — 98.71 % per-cell label-concordant with R, both ~98.5 % accurate vs ground truth	Advanced
10	Cell-cycle & Module Scoring	20,729 cells · THP-1 · GSE153056 (Papalexi 2021)	Gene-program scoring (AddModuleScore ↔ add_module_score) · cell-cycle phase (CellCycleScoring ↔ cell_cycle_scoring) · S/G2M scores + discrete phase — 96.6 % per-cell Phase-concordant with R, scores correlate at Pearson ≥ 0.998	Advanced
11	Dimensional-Reduction Extras	2,700 PBMCs · 10x Genomics (2016)	PC significance (JackStraw/ScoreJackStraw ↔ jack_straw/score_jackstraw) · RunICA ↔ run_ica · RunTSNE ↔ run_tsne — both tools keep 13 PCs; ICA matched |r| 0.982; caught two JackStraw bugs, both fixed (a too-tight null + the wrong aggregation test)	Advanced
12	Leverage-Score Sketching	13,999 cells · CTRL/STIM · ifnb (Kang 2018)	Scaling to atlas size (LeverageScore ↔ leverage_score · SketchData ↔ sketch_data · ProjectData ↔ project_data) · both of Seurat's regimes · uniform-sampling control · on-disk LazyMatrix — exact-regime Spearman 1.000000; leverage tracks rarity at −0.929 in both tools; caught two bugs, both fixed (full-rank leverage + anchor-based label transfer)	Advanced
13	The Object Model Itself	2,700 PBMCs · 10x Genomics (2016)	The container, not an algorithm: Cells/Features · the v5 layered assay (Layers/LayerData/split/JoinLayers) · Key · Embeddings/Loadings/Stdev · Graphs · FetchData · Idents/WhichCells/RenameIdents/subset · Command — 91 of 91 anchors match exactly, no tolerance (the two neighbour-graph anchors closed by PR #55); caught eleven bugs, all fixed (a split/join round trip that silently misordered columns, FetchData returning sparse objects instead of numbers, an inert command log)	Advanced
14	Spatial Statistics & the Spatial Container	36,602 cells · 248 genes · 10x Xenium mouse brain	The spatial container and the one spatial statistic never checked against R: LoadXenium ↔ load_xenium · CreateFOV/CreateCentroids/CreateSegmentation ↔ create_fov/create_centroids/create_segmentation · GetTissueCoordinates · Radius · FindSpatiallyVariableFeatures ↔ find_spatially_variable_features — 38 of 39 anchors match exactly; Moran's I to 1.6e-14 and 10/10 of Seurat's top genes, on a slide R cannot hold in memory; caught three bugs, all fixed (Moran's I on a kNN graph instead of R's inverse-square weights, centroids with no radius, unclosed polygons)	Advanced
15	The Differential-Expression Test Suite	2,700 PBMCs · 10x Genomics (2016)	All eight find_markers tests against FindMarkers — wilcox · t · bimod · LR · negbinom · roc · MAST · DESeq2 — on a shared cell assignment so no clustering difference can pose as a DE difference. All seven per-cell tests reproduce Seurat's top 50 exactly; avg_log2FC to 7.1e-15; caught two bugs, both fixed (Seurat's pseudocount on the group mean instead of the sum, which also changed which genes logfc_threshold returned; and a moment-dispersion LRT where Seurat runs an ML-dispersion Wald test)	Advanced
16	Out of Core — LazyMatrix vs BPCells	2,700 PBMCs · 10x Genomics (2016)	truecell's on-disk layer against BPCells, Seurat's: write_lazy_matrix/open_lazy_matrix ↔ write_matrix_dir/open_matrix_dir · streaming LogNormalize/VST/ScaleData ↔ Seurat's IterableMatrix methods · .CalcN · storage formats. 14 of 14 anchors match, 1998/2000 variable features shared with Seurat. The finding is each tool against itself: truecell's on-disk and in-memory paths are bit-identical, Seurat's differ by 1.0e-06 and pick a different variable feature. Caught seven bugs, all fixed (five functions that densified the whole store — so going on disk raised peak memory 4.6× — a constructor that ended laziness before analysis began, and a LOESS whose fit moved 28.8 % under a 1e-15 nudge)	Advanced
17	Visium — the Spatial Container	2,695 spots · 10x V1_Mouse_Brain_Sagittal_Anterior	The Visium loader and container: Load10X_Spatial/Read10X_Image ↔ load_visium · VisiumV2 · ScaleFactors · Radius · GetTissueCoordinates · SpatialDimPlot ↔ spatial_dim_plot — 24 of 24 anchors match, 17 exactly, coordinates to max\|dx\| = 0. The first tutorial where Seurat is the one that's wrong: it stores spot_diameter_fullres in the FOV's radius, which the slide's fixed 100 µm pitch shows is a diameter, and Radius() on its own VisiumV2 returns NULL. Caught one truecell bug (the tissue image came back 255× apart depending on whether matplotlib or Pillow was installed) and aligned three defaults to Seurat	Spatial
18	Anchor Internals — CCA & RPCA	2,400 cells · CTRL/STIM · ifnb (Kang 2018)	Both anchor paths at the level of the anchors and embeddings themselves, not the clustering they produce: v4 FindIntegrationAnchors ↔ find_integration_anchors · IntegrateData ↔ integrate_data · RunCCA · ReciprocalProject · ScoreAnchors · FilterAnchors · FindWeights; v5 IntegrateLayers ↔ integrate_layers · IntegrateEmbeddings — RPCA agrees on 100 % of Seurat's v4 anchors (649/649) and 30/30 PCs on the v5 embedding, up from 1/30; CCA v4 anchors 99.9 %, up from 70.0 %. Caught eighteen bugs, all fixed — CCA standardizing vs L2-normalizing, sklearn's randomized SVD drifting reciprocal-PCA's trailing PCs, integrate_layers silently running v4's IntegrateData behind the v5 IntegrateEmbeddings name, and a symmetrized KNN graph plus a missing SNN diagonal	Advanced

---

Quick Start

Clone the repo and install dependencies once, then pick any tutorial script:

git clone https://github.com/GenomicAI/truecell.git
cd truecell
uv venv && source .venv/bin/activate   # Windows: .venv\Scripts\activate
uv pip install -e ".[analysis]"

Each tutorial has a Python script that runs the analysis and prints validation output, and a figure-generation script that writes plots to a figures_*/ subfolder.

The datasets download automatically on first run, to ~/.truecell_data/ (~200 MB ~770 MB with every dataset cached).

Checking the tutorials still run

Two layers, because the fast one cannot reach everything:

pytest tests/test_tutorial_marker_tables.py      # fast, no data, runs in CI
TRUECELL_TUTORIAL_SMOKE=1 pytest tests/test_tutorial_smoke.py -v   # real runs, needs the cached data

The first covers the presentation helpers with synthetic frames. The second actually executes each tutorial and takes a few minutes; it is opt-in rather than "skip when the data is missing", so a skip always means nobody asked for it and never that it passed. Run the second before cutting a release — the unit suite passing tells you nothing about whether these scripts still work end to end.

Dataset loaders and the R export bridge

Every tutorial pairs R and Python on the same counts. Two ways that happens:

• Original public files — truecell.datasets downloads the dataset's own GEO/10x files, which R and Python both read directly. This covers pbmc3k, pbmc8k, cbmc_citeseq, xenium_mouse_brain, and the two hashing/perturbation datasets, pbmc_hashing (GSE108313) and thp1_eccite (GSE153056). Nothing extra to run — they cache to ~/.truecell_data/ on first use.

• The R bridge — a few datasets exist only as curated SeuratData .rda objects with no clean cross-language raw source (ifnb, panc8). For those, run the one-time export first, which writes a gzipped 10x folder both languages read so the counts are byte-identical:

Rscript tutorials/export_seuratdata.R ifnb     # ~394 MB, first run only
Rscript tutorials/export_seuratdata.R panc8    # ~117 MB

Then truecell.datasets.ifnb() / panc8() load that export in Python. (This is the mirror of the *_verify.R scripts, where R instead depends on the Python download.) Needs R with Seurat + SeuratData.

---

Tutorial 1 — PBMC 3k Guided Clustering

Walkthrough: pbmc3k_tutorial.md

A step-by-step port of the official Seurat PBMC 3k tutorial. Covers every step from raw counts to annotated cell types — the entry point for new users.

python tutorials/pbmc3k_tutorial.py    # runs analysis + validation
python tutorials/generate_plots.py     # writes tutorials/figures/

What you'll learn: - Load 10x Genomics data and create a Truecell object - Compute QC metrics (percent mitochondrial reads) and filter low-quality cells - Normalize counts (LogNormalize), select highly variable genes (VST), and scale data - Run PCA, build the KNN/SNN neighbor graph, and cluster with Louvain - Embed with UMAP and visualize clusters - Find cluster markers with Wilcoxon rank-sum and annotate 8 cell types (Seurat's nine minus DC, which merges into CD14+ Mono at this resolution)

Key output figures (in tutorials/figures/):

Figure	Description
01_qc_violin.png	QC metrics violin plot
03_variable_features.png	HVG mean–variance plot
07_umap_clusters.png	UMAP coloured by cluster
08_feature_plots.png	Canonical marker feature plots
10_marker_heatmap.png	Top-10 markers per cluster heatmap
11_umap_labeled.png	UMAP with cell-type labels

Measured against R Seurat 5.5.1 — run Rscript tutorials/pbmc3k_verify.R then python tutorials/pbmc3k_tutorial.py --report. Both sides run their own pipeline from the same 10x bytes; nothing is pinned across them.

Step	truecell vs Seurat
Cells surviving QC	the same 2,638 barcodes, nCount/nFeature exact, percent.mt to 5.3e-15
VST per gene, all 13,714	mean 4.8e-14, variance 1.6e-11, variance.expected 2.5e-2 and variance.standardized 2.6e-2 relative — the whole gap is the LOESS fit
The 2,000 variable features	1,998 shared; the two swaps are genes 0.03 apart in a LOESS fit, at ranks 1,982–2,000
PCA, the 10 dims clustering uses	matched |r| mean 0.9988, min 0.9946, no reordering
kNN graph	52,760 = 2,638 × 20 on both
Clusters at resolution 0.5	truecell 8, Seurat 9 — ARI 0.899, 2,519/2,638 cells agree
Clusters across the sweep	0.4 → 9 vs 9 (ARI 0.896) · 0.8 → 11 vs 11 (0.826) · 1.2 → 12 vs 12 (0.800). The cluster count matches at every resolution except 0.5
The one cluster Seurat has and truecell does not	its 32 DC cells land, all 32, in truecell's CD14+ Mono cluster
Markers on the two clusters whose cells match exactly	identical gene sets (151/151 and 242/242), avg_log2FC to 4.9e-15 and 4.6e-14 respectively

The DC split is the honest caveat: at resolution 0.5 a 32-cell dendritic-cell population sits on the boundary, and the two runs land on opposite sides of it. Tutorial 2 shows the same boundary going the other way — there truecell resolves DC from Platelet and Seurat merges them.

---

Tutorial 2 — PBMC 8k Advanced Subclustering

Walkthrough: advanced_pbmc8k_subclustering.md

A larger dataset (~8,400 cells) that demonstrates the standard Seurat subclustering workflow. After global clustering and broad annotation, the T/NK lymphoid compartment is isolated and re-analysed to resolve four functionally distinct subsets that global clustering merges into one cluster.

python tutorials/pbmc8k_subclustering_tutorial.py   # analysis + validation
python tutorials/generate_advanced_plots.py         # writes tutorials/figures_advanced/

What you'll learn: - Run the full pipeline on a larger dataset and use z-scored relative enrichment for automated cluster annotation (avoids bias from high-magnitude housekeeping genes) - Subset a cell lineage and re-run HVG → PCA → neighbors → clusters → UMAP within it - Apply hierarchical gating to annotate T/NK subsets: NK cells (CD3⁻ NKG7/GNLY⁺) → CD8 T (CD8A/B⁺ or cytotoxic CD3⁺, incl. γδ T and MAIT) → CD4 Naive (CCR7/SELL⁺) → CD4 Memory

Key output figures (in tutorials/figures_advanced/):

Figure	Description
03_umap_global_clusters.png	UMAP — global clusters
04_umap_global_celltypes.png	UMAP — broad cell types
07_umap_tnk_subclusters.png	T/NK subcluster UMAP
08_umap_tnk_subsets.png	T/NK subsets annotated
09_tnk_subset_featureplots.png	CD3, CD4, CD8A, CCR7, NKG7 feature plots
11_tnk_markers_heatmap.png	T/NK subset marker heatmap

Measured against R Seurat 5.5.1 — run Rscript tutorials/pbmc8k_subclustering_verify.R then python tutorials/pbmc8k_subclustering_tutorial.py --report.

Step	truecell vs Seurat
Cells surviving QC	the same 7,475 barcodes, metrics exact to 5.3e-15
Stage 1 — global clusters	truecell 13, Seurat 12; ARI 0.977, 7,341/7,475 cells agree
Broad lineage label, per cell	0.9858
The compartment handed to stage 2	Jaccard 0.9991 — 4,631 of 4,635 T/NK cells are the same barcodes
Stage 2 — T/NK subclusters	truecell 12, Seurat 11; ARI 0.916 on the shared cells
T/NK subset label, per cell	0.9821; subset sizes agree within 25 cells of 4,631
kNN graph	149,500 = 7,475 × 20 on both

The compartment Jaccard is the load-bearing number. Everything in stage 2 is conditioned on which cells stage 1 handed it, so a subclustering fed from the wrong global clusters would still produce a compartment, still produce subclusters, and still pass a count check — which is why the handoff compares barcodes rather than sizes.

The extra truecell cluster is the mirror image of Tutorial 1's caveat: Seurat's 100-cell cluster 11 holds both Platelet and DC, and truecell splits it in two (54 DC + 53 Platelet). Here truecell resolves what Seurat merges; on PBMC 3k it is the other way round. Neither run is uniformly finer than the other.

---

Tutorial 3 — CBMC CITE-seq Multimodal

Walkthrough: multimodal_citeseq.md

A Python port of Seurat's multimodal vignette on the CBMC CITE-seq dataset (GSE100866, Stoeckius et al. 2017): ~8,600 cord-blood mononuclear cells profiled simultaneously for the transcriptome and 13 surface proteins.

python tutorials/cbmc_citeseq_tutorial.py     # analysis + validation
python tutorials/generate_multimodal_plots.py # writes tutorials/figures_multimodal/

What you'll learn: - Load and align RNA + ADT matrices (handles the mouse spike-in gene prefix) - Attach a second ADT assay to the Truecell object - CLR-normalize the surface protein counts (margin=2 centers each cell across the protein panel, matching Seurat's flag exactly — axis included) - Overlay protein expression on the RNA UMAP — protein gives smooth, high-SNR signal where the encoding mRNA is sparse - Annotate cells with a protein-priority gating strategy + RNA fallback for populations the 13-protein panel cannot resolve (platelets, erythroid, pDC, cycling cells) - Cluster both modalities jointly with WNN, and read the per-cell modality weights to see which lineages are decided by protein and which fall back to RNA

Key output figures (in tutorials/figures_multimodal/). Each has an r_-prefixed R Seurat counterpart for the side-by-side comparisons in the walkthrough:

Figure	Description
01_rna_umap_clusters.png	UMAP — RNA clusters
02_rna_umap_celltypes.png	UMAP — protein + RNA cell types
03_adt_featureplots.png	8 surface protein feature plots on RNA UMAP
04_protein_vs_rna.png	Protein vs RNA side-by-side (CD19, CD3, CD8, CD14)
05_adt_ridgeplots.png	ADT ridge plots by cell type
06_adt_scatter_CD4_CD8.png	CD4 vs CD8 protein scatter
07_adt_scatter_CD19_CD3.png	CD19 vs CD3 protein scatter
08_wnn_umap_clusters.png	WNN joint clusters on the joint embedding
09_wnn_vs_rna_umap.png	Same cells and labels — RNA-only vs WNN joint UMAP
10_adt_weight_by_celltype.png	Learned per-cell ADT weight, by cell type

---

Tutorial 4 — PBMC 3k SCTransform

Walkthrough: sctransform_vignette.md

A Python port of Seurat's sctransform vignette. SCTransform replaces the NormalizeData → FindVariableFeatures → ScaleData trio with a single regularized negative-binomial model whose Pearson residuals are the normalized values, then clusters over 30 PCs. The tutorial runs SCTransform and the standard log-normalization workflow on the same cells to compare their resolution.

python tutorials/pbmc3k_sctransform_tutorial.py   # analysis + SCT-vs-std comparison
python tutorials/generate_sctransform_plots.py    # writes tutorials/figures_sctransform/

What you'll learn: - Run sctransform(vars_to_regress=["percent.mt"]) — one call replacing three, producing a new SCT assay with corrected counts and residual scale.data - Cluster over 30 PCs (the vignette's deeper embedding) and read the resulting cytotoxicity gradient: Naive CD4 → Memory CD4 → CD8 Effector → NK - Compare SCTransform against log-normalization on identical cells

Key output figures (in tutorials/figures_sctransform/):

Figure	Description
01_sct_umap_clusters.png	UMAP — SCT clusters
02_sct_umap_celltypes.png	UMAP — annotated cell types
03_sct_featureplots_1.png	CD8A, GZMK, CCL5, S100A4, ANXA1, CCR7
04_sct_featureplots_2.png	CD3D, ISG15, TCL1A, FCER2, XCL1, FCGR3A
05_sct_violins.png	Vignette marker violins per cluster
06_sct_vs_std_umap.png	SCTransform vs LogNormalize, side by side

Accuracy vs R: verified against a live Seurat 5.5.1 / sctransform 0.4.3 run on the same cells — the regularized intercept matches at Spearman 1.0000, theta at 0.96, residual variance at 0.9986, and 2,913 of the 3,000 variable features agree. Like Seurat 5, sctransform defaults to vst_flavor="v2"; under vst_flavor="v1" Python and R land on 13 clusters each. At the v2 default Truecell resolves 13 where R resolves 12 — a genuine one-cluster difference (R is stable at 12 across seeds) from vst's random step-1 gene sample and the different clustering libraries, the same ±1 caveat as the other tutorials.

---

Tutorial 5 — Xenium Spatial (R vs Python)

Walkthrough: xenium_spatial_tutorial.md

The spatial counterpart to the others: a 10x Xenium section (mouse brain coronal CTX+HP subset — the dataset in Seurat's spatial vignette, 36,602 cells × 248 genes) analysed side by side in R Seurat and Truecell. It reproduces the style of an internal Xenium mast-cell/neighbourhood workflow on fully public data, exercising the spatial Seurat-parity layer end to end.

python tutorials/generate_spatial_plots.py     # auto-downloads ~20 MB → figures_spatial/
Rscript tutorials/xenium_spatial_verify.R      # R reference figures + r_reference.json
python tutorials/compare_xenium_anchors.py     # prints the R-vs-Python parity table

What you'll learn: - load_xenium — build an object with expression and centroids, keeping only Gene Expression features (like LoadXenium) - Deterministic marker-panel cell typing (the KIT+ TPSAB1+-style rule) - image_dim_plot / image_feature_plot — plot cells in tissue space (immune to the ggplot2 4.x ImageDimPlot blank-render bug) - nearest_neighbor_distance / local_neighborhood — FNN::get.knn idioms - build_niche_assay — neighbourhood-composition niches (BuildNicheAssay) - composition_test — Fisher + BH enrichment across a spatial split

Key output figures (in tutorials/figures_spatial/, r_* = R Seurat):

Figure	Description
03_image_celltype.png	Marker cell types in tissue space
05_image_feature_Slc17a7.png	Excitatory-neuron marker in space
06_image_niches.png	Neighbourhood niches
07_image_focal.png	Focal (Vascular) cells highlighted

Accuracy vs R: every deterministic anchor matches to 8 significant figures — cell counts, all cell-type counts, nearest-neighbour distances, local density, and the composition test (log2 ratios, Fisher p, BH padj, χ² p). Clustering / UMAP / niche layout are stochastic and agree in structure only, the same caveat as the other tutorials.

---

Tutorial 6 — Cell Hashing Demultiplexing (R vs Python)

Walkthrough: hashing_vignette.md

A port of Seurat's hashing vignette: eight samples, each tagged with an antibody hashtag (HTO), pooled on one lane. Both of Truecell's demultiplexers run against their R references on byte-identical GEO input (GSE108313, Stoeckius et al. 2018 — the original plain-text matrices, not the vignette's binary .rds). It is the first real-data check of the CLR-margin fix (#32) and the clara default (#34).

python tutorials/pbmc_hashing_tutorial.py     # downloads ~34 MB, prints the report
Rscript tutorials/pbmc_hashing_verify.R       # Seurat reference → r_calls.csv + r_*.png
python tutorials/pbmc_hashing_tutorial.py     # re-run → prints the R-vs-Python concordance
python tutorials/generate_hashing_plots.py    # Truecell figures → figures_hashing/py_*.png

What you'll learn: - Attach a hashtag panel as an "HTO" assay and CLR-normalise it (margin 1 — per hashtag across cells, Seurat's default for hashing) - hto_demux — Seurat's HTODemux: a negative-binomial cutoff per hashtag, learned from its background cluster - multiseq_demux — Seurat's MULTIseqDemux: a kernel-density cutoff instead, and where the two methods diverge - Reading a combined human+mouse alignment as an independent cross-species doublet signal the hashtags never saw

Key output figures (in tutorials/figures_hashing/, r_* = R Seurat):

Figure	Description
py_01_ridge.png	Hashtag CLR enrichment per assigned sample (the canonical QC read)
py_02_scatter.png	Two hashtags, coloured by HTODemux global class
py_03_ncount_violin.png	Total hashtag counts per global class
py_04_species_scatter.png	Human vs mouse UMIs — the barnyard ground truth

Accuracy vs R (call-for-call on identical input, 39,842 cells):

Comparison	Agreement
HTODemux global class & sample assignment	99.81 % (77 cells differ)
MULTIseqDemux call (MULTI_ID)	94.67 %

HTODemux reproduces Seurat almost exactly — the 0.2 % residual is clara's un-seedable sampling. MULTIseqDemux's larger gap is a genuine method-level difference in the KDE step (scipy gaussian_kde vs R density() — bandwidth and grid, not a single knob), documented rather than smoothed over. Because these are raw unfiltered barcodes, the Negative rate is high (~42 %) and the totals do not match the vignette's pre-filtered PBMC headline — by design.

---

Tutorial 7 — Mixscape Pooled CRISPR Screen (R vs Python)

Walkthrough: mixscape_vignette.md

A port of Seurat's Mixscape vignette: the THP-1 ECCITE-seq screen (GSE153056, Papalexi et al. 2021), 20,729 cells each carrying one guide against one of 25 immune-regulatory genes. Mixscape solves the problem that carrying a guide is not the same as being perturbed — it separates true knockouts (KO) from non-perturbed escapers (NP) that look like controls. Run against the R reference on the same GEO bytes and a shared variable-feature basis.

python tutorials/thp1_mixscape_tutorial.py    # downloads ~66 MB, writes HVGs, prints the report
Rscript tutorials/thp1_mixscape_verify.R      # Seurat reference → r_calls.csv + r_*.png
python tutorials/thp1_mixscape_tutorial.py    # re-run → prints the R-vs-Python concordance
python tutorials/generate_mixscape_plots.py   # Truecell figures → figures_mixscape/py_*.png

What you'll learn: - calc_perturb_sig — Seurat's CalcPerturbSig: subtract each cell's nearest-NT mean (over 40 PCs, within replicate) to leave the local perturbation signature - run_mixscape — Seurat's RunMixscape: an iterative two-component Gaussian mixture per guide, splitting KO from NP by perturbation score - mixscape_lda — Seurat's MixscapeLDA: a supervised map of the guide populations - Reading a per-guide knockout rate — which edits actually took, and why the checkpoint genes read as zero at the RNA level (their effect is on protein)

Key output figures (in tutorials/figures_mixscape/, r_* = R Seurat):

Figure	Description
py_01_perturb_score.png	IFNGR2 perturbation-score density — the KO/NP/NT split
py_02_lda.png	MixscapeLDA map, coloured by global class
py_03_heatmap.png	DE heatmap, NT vs IFNGR2 KO, cells ordered by KO probability
py_04_ko_rate.png	Per-guide knockout rate — the headline result

Accuracy vs R (call-for-call on shared input + shared HVG basis, 20,729 cells):

Comparison	Agreement
Mixscape global class (KO/NP/NT)	97.45 % (528 cells differ)
Mixscape full class (<gene> KO/NP)	97.45 %

Both tools agree on all 2,386 NT cells, sort the same 14 guides to a flat-zero knockout rate, and rank the strong interferon-γ hits identically (STAT1 98.6 %, JAK2 97.8 %, IFNGR2 97.7 % per-cell concordant). The disagreement concentrates in the weak, boundary guides (MYC, SPI1, BRD4, CUL3) whose score sits close to the NT mode — a genuine method-level difference in the EM mixture and the per-gene DE, not a bug. No defect found.

---

Tutorial 8 — Batch Integration (R vs Python)

Walkthrough: integration_vignette.md

A port of Seurat's integration vignette: the ifnb dataset (Kang et al. 2018), ~14,000 PBMCs split CTRL vs interferon-β-stimulated. Interferon drives a global response, so uncorrected the cells cluster by condition rather than cell type — the standard batch-effect benchmark. Runs all three of Truecell's integration paths against Seurat on identical counts and a shared variable-feature basis.

Rscript tutorials/export_seuratdata.R ifnb        # one-time counts export (~394 MB pkg)
python  tutorials/ifnb_integration_tutorial.py    # writes HVGs, prints the scoreboard
Rscript tutorials/ifnb_integration_verify.R       # Seurat reference → r_calls.csv + r_*.png
python  tutorials/ifnb_integration_tutorial.py    # re-run → prints the R-vs-Python concordance
python  tutorials/generate_integration_plots.py   # Truecell figures → figures_integration/py_*.png

What you'll learn: - run_harmony — Seurat's RunHarmony: correct the PCA embedding so batches mix - integrate_layers(method="cca"/"rpca") — Seurat's IntegrateLayers: anchor-based correction in a shared CCA space, or reciprocal-PCA - Scoring an integration without coordinate-comparable embeddings: silhouette by batch (mixing) and by cell type (preservation), and cluster ARI against the annotations — all partition-based, so R and Python compare on identical terms

Key output figures (in tutorials/figures_integration/, r_* = R Seurat):

Figure	Description
py_01_uncorrected_stim.png	UMAP of raw PCA, by condition — the batch effect
py_02_harmony_stim.png	UMAP after Harmony, by condition — now mixed
py_03_harmony_celltype.png	Same map by cell type — the biology survived
py_04_scoreboard.png	Batch mixing vs cell-type recovery, per method

Accuracy vs R (partition-based, 13,999 cells; mix = batch-mixing entropy):

method	py mix	R mix	py ARI→celltype	R ARI→celltype
Harmony	0.991	0.991	0.917	0.930
CCA	0.990	0.991	0.884	0.873
RPCA	0.867	0.914	0.677	0.735

Harmony and CCA reproduce Seurat's integration to three decimals — the first real-data confirmation of run_harmony / integrate_layers. RPCA is where the tutorial earned its keep: it surfaced two defects — a crash on unequal batch sizes (fixed in #41) and a deeper under-integration bug (fixed here: per-object scaling + Seurat's reciprocal-embedding normalization lift batch mixing from 0.222 to 0.867). Both came with regression tests; neither reproduces on any synthetic fixture, only on the real Seurat comparison.

---

Tutorial 9 — Reference Mapping (R vs Python)

Walkthrough: refmap_vignette.md

A port of Seurat's reference mapping vignette: the panc8 pancreatic-islet dataset (Baron et al. 2016), ~14,900 cells across five sequencing technologies. Reference mapping is integration's asymmetric cousin — keep an annotated atlas fixed and borrow its labels for a new, unlabelled dataset. The tutorial annotates a SMART-seq2 query (2,394 cells) from a CEL-seq2 reference (2,285 cells) — a genuine cross-chemistry transfer — and, because the query's true celltype ships with the data, scores the result directly for accuracy. A single-technology reference is used on purpose, to isolate the transfer machinery from the integration path Tutorial 8 covers.

Rscript tutorials/export_seuratdata.R panc8            # one-time counts export (~117 MB pkg)
python  tutorials/panc8_reference_mapping_tutorial.py  # writes HVGs, prints accuracy + per-class recall
Rscript tutorials/panc8_reference_mapping_verify.R     # Seurat reference → r_calls.csv + r_*.png
python  tutorials/panc8_reference_mapping_tutorial.py  # re-run → prints the R-vs-Python concordance
python  tutorials/generate_refmap_plots.py             # Truecell figures → figures_refmap/py_*.png

What you'll learn: - find_transfer_anchors — Seurat's FindTransferAnchors(reduction="pcaproject"): project the query through the reference's PCA loadings and score the anchors - transfer_data — Seurat's TransferData: a distance-weighted, score-scaled vote over the reference labels → predicted.id + per-class prediction.score.* - map_query / project_umap — Seurat's MapQuery / ProjectUMAP: run the reference's fitted UMAP model in transform-only mode so the query lands on the atlas you already know how to read - Reading a per-class recall table: where a small single-tech reference annotates confidently, and where its rare types are too thin to anchor

Key output figures (in tutorials/figures_refmap/, r_* = R Seurat):

Figure	Description
py_01_reference_umap_celltype.png	The reference atlas UMAP, by cell type
py_02_query_projected_predicted.png	The query projected into the reference UMAP, by transferred label
py_03_query_projected_truth.png	The same projection, by the held-out true label
py_04_perclass_recall.png	Per-cell-type transfer recall vs support

Accuracy vs R (per-cell label transfer on shared input + shared HVG basis, 2,394 query cells):

Comparison	Agreement
Same predicted.id per cell (truecell vs Seurat)	98.71 % (31 cells differ)
truecell accuracy vs ground-truth celltype	98.45 %
Seurat accuracy vs ground-truth celltype	98.79 %

Both tools annotate the query almost identically — 2,363 of 2,394 cells get the same label — and each is ~98.5 % accurate against the held-out truth. Every abundant cell type is recovered at ≥98 %; the entire error budget is the rare types (epsilon, quiescent_stellate, schwann, <10 reference cells each), where both tools stumble the same way — the honest limit of a small single-technology reference, not a divergence. No defect found — the transfer stack ports faithfully on its first real-data benchmark.

---

Tutorial 10 — Cell-cycle & Module Scoring (R vs Python)

Walkthrough: cellcycle_vignette.md

A port of Seurat's cell-cycle vignette and AddModuleScore, run on the THP-1 ECCITE-seq dataset (GSE153056, Papalexi et al. 2021) — a proliferating leukemia line, so unlike the resting PBMCs of earlier tutorials it has real S and G2/M populations (the right substrate for cell-cycle scoring). Reuses the Mixscape counts; only the raw RNA matters.

python  tutorials/thp1_cellcycle_tutorial.py    # downloads ~66 MB (shared with Mixscape), writes gene lists
Rscript tutorials/thp1_cellcycle_verify.R       # Seurat reference → r_calls.csv + r_*.png
python  tutorials/thp1_cellcycle_tutorial.py    # re-run → prints the R-vs-Python concordance
python  tutorials/generate_cellcycle_plots.py   # Truecell figures → figures_cellcycle/py_*.png

What you'll learn: - add_module_score — Seurat's AddModuleScore: score a gene program as its mean expression minus a control set drawn from the same expression bins - cell_cycle_scoring — Seurat's CellCycleScoring: module scoring on the Tirosh S / G2M sets (CC_GENES = cc.genes.updated.2019), then a discrete Phase - Why the scores are not expected to be bit-identical across tools (both sample control genes at random, and NumPy's RNG is not R's) — and why the discrete Phase is robust to it anyway

Key output figures (in tutorials/figures_cellcycle/, r_* = R Seurat):

Figure	Description
py_01_score_scatter.png	S.Score vs G2M.Score, coloured by assigned phase
py_02_phase_bar.png	Cell count per phase (G1/S/G2M)
py_03_ifn_hist.png	Interferon-response module-score distribution (AddModuleScore)

Accuracy vs R (identical counts + shared resolved gene lists, 20,729 cells):

Comparison	Agreement
Per-cell Phase concordance (truecell vs Seurat)	96.62 % (701 cells differ)
S.Score / G2M.Score correlation (Pearson)	0.998 / 0.999
IFN.Response module score correlation (Pearson)	0.9995

The continuous scores correlate at Pearson ≥ 0.998 — the algorithm is faithfully ported; the only reason they are not bit-identical is the random control-gene set (NumPy vs R). 96.62 % of cells get the same phase, with the disagreements sitting on the S=0 / G2M=0 boundary where the small score wobble tips the call — the same boundary-sensitivity as Mixscape's weak guides. No defect found — the same "don't chase the RNG" residual as clara (hashing) and the MULTI-seq KDE.

---

Tutorial 11 — Dimensional-Reduction Extras (R vs Python)

Walkthrough: dimreduc_vignette.md

The three reductions that sit beside the standard PCA → UMAP path, and the question every guided-clustering run has to answer first: how many PCs are real? Same PBMC 3k object and HVG basis as Tutorial 1.

python  tutorials/pbmc3k_dimreduc_tutorial.py   # downloads ~8 MB, writes the shared lists
Rscript tutorials/pbmc3k_dimreduc_verify.R      # slow: 100 PCA refits, ~2 min
python  tutorials/generate_dimreduc_plots.py    # figures + the side-by-side numbers

What you'll learn: - jack_straw / score_jackstraw — Seurat's JackStraw/ScoreJackStraw: a permutation test for how many PCs carry real signal, not just the elbow guess - run_ica — Seurat's RunICA, matched to R one-to-one by |Pearson r| via the Hungarian algorithm, since a component's sign and order are not defined - run_tsne — Seurat's RunTSNE, compared on neighbourhood structure preserved from the PCA space rather than on raw coordinates, which never agree across Barnes-Hut and scikit-learn - Why a stochastic PC-count answer needs a declared band, not a point comparison — and how a 60-seed sweep turns "close enough" into something --report can fail on

Key output figures (in tutorials/figures_dimreduc/, r_* = R Seurat):

Figure	Description
py_01_jackstraw_scores.png	JackStraw significance score per PC
py_02_nsig_compare.png	Features below threshold per PC, truecell vs R
py_03_elbow.png	The cheap first-look elbow plot
py_04_tsne.png	t-SNE embedding, coloured by LYZ
py_05_ica.png	ICA embedding, coloured by LYZ

Accuracy vs R (PCA basis matched through PC 20, |r| ≥ 1.0000):

Comparison	Agreement
PCs kept (truecell vs R's 13)	14, band |Δ| ≤ 2 (modal answer over 60 seeds: 13)
ICA, mean matched |Pearson r| over 20 components	0.9991 (worst pair 0.9960)
t-SNE, 30-NN retained from PCA (truecell / R)	0.474 / 0.477
t-SNE, 30-NN shared between the two tools	0.862

Found and fixed two defects — jack_straw built its permutation null against a fixed PCA basis instead of refitting per replicate, and score_jackstraw ran a KS test where R runs prop.test, which together kept all 20 PCs where Seurat keeps 13. Both are ported exactly and now match R to nine significant figures.

---

Tutorial 12 — Leverage-Score Sketching (R vs Python)

Walkthrough: sketch_vignette.md

The Seurat v5 answer to "this atlas has two million cells and I have a laptop": subsample by statistical leverage rather than uniformly, analyse the small subset, then push the answers back onto every cell. Reuses the ifnb object from Tutorial 8.

Rscript tutorials/export_seuratdata.R ifnb      # one-time counts export
python  tutorials/ifnb_sketch_tutorial.py       # writes the shared lists
Rscript tutorials/ifnb_sketch_verify.R          # ~3 min
python  tutorials/generate_sketch_plots.py

What you'll learn: - leverage_score — Seurat's LeverageScore: a rank-50 truncated-SVD hat-matrix score that decides which cells are worth keeping in a sketch - sketch_data — Seurat's SketchData, both method="LeverageScore" and the uniform-sampling control that proves the leverage sketch is doing something a random sample would not - project_data — Seurat's ProjectData: extend a sketch's clustering, UMAP and labels back onto every cell the sketch never analysed - Why leverage has to be checked against something the algorithm was never told — rarity by cell type — not just against R's numbers

Key output figures (in tutorials/figures_sketch/, r_* = R Seurat):

Figure	Description
py_01_leverage_by_type.png	Mean leverage per cell type, rarest to most common
py_02_sketch_enrichment.png	Leverage sketch vs uniform draw, enrichment by type
py_03_leverage_vs_r.png	Per-cell leverage score, truecell vs R
py_04_rarity.png	Leverage vs cell-type abundance
py_05_projected_umap.png	The sketch's UMAP and labels, projected onto all 13,999 cells

Accuracy vs R (13,999 cells, exact-regime leverage):

Comparison	Agreement
Leverage, exact regime — per-cell Spearman vs R	1.000000 (max abs diff 3.4e-6)
Does leverage track rarity? Spearman(mean leverage, type size)	truecell −0.929 · R −0.929
project_data labels — per-cell agreement with R	94.9 % (98.1 % on a shared sketch)
project_data accuracy vs held-out annotations	truecell 0.9050 · R 0.9050

Found and fixed two defects — leverage_score whitened against the full matrix rank instead of Seurat's rank-50 truncation, which flattened the sampling weights until leverage sampling was indistinguishable from uniform sampling. And project_data transferred labels through integration anchors, which is not what Seurat does and costs exactly what sketching exists to remove — the wrong version even scored higher, which is why it survived review until this tutorial checked the mechanism, not just the number.

---

Tutorial 13 — The Object Model Itself (R vs Python)

Walkthrough: objects_vignette.md

Every other tutorial compares an algorithm. This one compares the container — the accessors, the v5 layer machinery, and the bookkeeping that Seurat's own command cheat sheet is made of. Same PBMC 3k object as Tutorial 1.

python tutorials/pbmc3k_objects_tutorial.py     # downloads ~24 MB on first run
Rscript tutorials/pbmc3k_objects_verify.R       # writes r_anchors.json
python tutorials/pbmc3k_objects_tutorial.py     # now prints the side-by-side table
python tutorials/generate_objects_plots.py      # figures

What you'll learn: - layers / layer_data / split_layers / join_layers — Seurat's v5 layered Assay: Layers, LayerData, split, JoinLayers - fetch_data — Seurat's FetchData, including addressing an embedding column by its reduction Key (PC_1) - idents / which_cells / rename_idents / subset — Seurat's identity machinery, compared cell-by-cell and in order, not just by count - obj.commands — Seurat's Command() log, and why it is a lookup table users query rather than decoration - Why almost nothing here is stochastic — 89 of 91 anchors are exact matches, no tolerance at all, which is what makes this the sharpest net in the series

Key output figures (in tutorials/figures_objects/):

Figure	Description
py_01_split_join.png	The split → JoinLayers round trip, before and after
py_02_concordance.png	Anchors matching R, by section
py_03_nn_degree.png	kNN graph degree distribution, truecell vs R

Accuracy vs R (2,700 cells × 13,714 genes, 91 anchors total):

Comparison	Agreement
Anchors matching Seurat exactly	91 / 91
Cells — order-sensitive md5 of the barcode vector	identical
nCount_RNA total	6,386,518 on both
FetchData — CD3E column total	2,831.901591 on both
split → JoinLayers round trip	identity — name, cell order and matrix restored
kNN / SNN graph edges	RNA_nn 54,000, RNA_snn 199,616 on both

Found and fixed eleven defects — the layered assay's split/JoinLayers pair returned the right numbers in the wrong columns, silently; FetchData returned sparse-matrix objects instead of expression values on the most-called accessor in the library; and the command log and orig.ident had never been wired up at all. All eleven are fixed in the same pull request as this tutorial; two further neighbour-graph differences it found but did not fix were closed later, in PR #55.

---

Tutorial 14 — Spatial Statistics & the Spatial Container (R vs Python)

Walkthrough: svf_vignette.md

The spatial container — FOV, Centroids, Segmentation — and the one spatial statistic in the library, FindSpatiallyVariableFeatures, neither of which Tutorial 5 checked against R. Same Xenium mouse-brain slide.

python  tutorials/xenium_svf_tutorial.py    # writes anchors + cells.txt
Rscript tutorials/xenium_svf_verify.R       # needs Rfast2
python  tutorials/xenium_svf_tutorial.py --report
python  tutorials/generate_svf_plots.py

What you'll learn: - create_centroids / create_fov / create_segmentation — Seurat's CreateCentroids/CreateFOV/CreateSegmentation, down to auto-radius and closed-ring polygon behaviour - find_spatially_variable_features — Seurat's FindSpatiallyVariableFeatures(selection.method="moransi"), evaluated in row blocks so the full n×n weight matrix — 10.7 GB on this slide — is never built - Why a statistic has to match its weights, not just look plausible: a k-nearest-neighbour approximation scored Pearson 0.986 against R and still missed 3 of R's top 10 genes

Key output figures (in tutorials/figures_svf/):

Figure	Description
py_01_moransi_vs_r.png	Moran's I per gene, both weightings against Seurat
py_02_top_n_overlap.png	Recovery of Seurat's top-N gene ranking
py_03_spot_radius.png	The centroid spot radius, before and after

Accuracy vs R (36,602 cells × 248 genes, 39 anchors total):

Comparison	Agreement
Anchors matching Seurat exactly	38 / 39
Moran's I per gene vs R (2,000-cell shared subset)	max abs diff 1.6e-14, Pearson 1.0000000000
Moran's I — Seurat's top 10 recovered	10 / 10, same order
Centroids.radius() — auto radius vs .AutoRadius	42.82543 vs 42.82543
Full-slide Moran's I, R's exact weights	5.3 s / 0.95 GB, vs a 10.7 GB dense matrix in R

Found and fixed three defects — Moran's I was computed on a k-nearest-neighbour graph instead of Seurat's row-standardised inverse-square weights, answering a different question under the same name; Centroids never received a radius, which silently disabled every true-to-scale spot renderer; and Segmentation stored polygons open where R closes them. Two further divergences — the permutation p-value and the FOV Key string — are left standing on purpose, and the reasoning is in the vignette.

---

Tutorial 15 — The Differential-Expression Test Suite (R vs Python)

Walkthrough: de_vignette.md

All eight find_markers tests against FindMarkers, on a shared cell assignment so no clustering difference can pose as a DE difference. Runs on two clusters (692 and 515 cells) from PBMC 3k.

python  tutorials/pbmc3k_de_tutorial.py     # writes groups.csv and py_<test>.csv
Rscript tutorials/pbmc3k_de_verify.R        # needs MAST + DESeq2
python  tutorials/pbmc3k_de_tutorial.py --report
python  tutorials/generate_de_plots.py

What you'll learn: - find_markers(test_use=...) — all eight of Seurat's tests: wilcox · t · bimod · LR · negbinom · roc · mast · deseq2 - Why avg_log2FC's pseudocount placement matters beyond the number it reports — it feeds logfc_threshold, so a formula error changes which genes come back, not just how they're described - Why mast and deseq2 are deliberate reimplementations rather than calls to the R packages, and where each is expected to diverge from them

Key output figures (in tutorials/figures_de/):

Figure	Description
py_01_log2fc_vs_r.png	avg_log2FC, before and after, vs R
py_02_threshold_impact.png	Genes returned at each logfc_threshold, truecell vs R
py_03_test_concordance.png	All eight tests against Seurat

Accuracy vs R (13,712 shared genes):

Comparison	Agreement
avg_log2FC vs Seurat, all genes	max abs diff 6.44e-15
Tests reproducing Seurat's top 50 genes	7 of 7 per-cell tests
wilcox / t / bimod / LR — p-value Spearman	1.000000 / 0.999980 / 0.999994 / 0.999975
mast — Spearman, detected >5 %	0.9979
negbinom — Spearman, detected >5 %	0.9165

Found and fixed two defects — avg_log2FC put Seurat's pseudocount on the group mean rather than the group sum (Seurat 4's formula, not Seurat 5's), which floored every fold change and changed which genes logfc_threshold returned; and negbinom ran a moment-dispersion likelihood-ratio test where Seurat runs an ML-dispersion Wald test.

---

Tutorial 16 — Out of Core: LazyMatrix against BPCells (R vs Python)

Walkthrough: lazy_vignette.md

Not a port comparison — there is no LazyMatrix inside Seurat to match value-for-value. Instead: does each tool's out-of-core path agree with its own in-memory path, and what does each format cost. Same PBMC 3k counts.

python tutorials/lazy_bpcells_tutorial.py       # writes anchors + cell list
Rscript tutorials/lazy_bpcells_verify.R         # needs BPCells (not on CRAN)
python tutorials/lazy_bpcells_tutorial.py --report
python tutorials/generate_lazy_plots.py

What you'll learn: - write_lazy_matrix / open_lazy_matrix — truecell's on-disk backend, against Seurat's BPCells (write_matrix_dir/open_matrix_dir) - Streaming LogNormalize / VST / ScaleData over an on-disk layer without ever densifying the whole matrix - Why "agrees with R to 1e-14" is not the same claim as "agrees with itself" — BPCells computes in single precision, so Seurat's own on-disk and in-memory runs disagree and select a different variable feature; truecell's two paths are bit-identical

Key output figures (in tutorials/figures_lazy/):

Figure	Description
py_01_self_consistency.png	Each tool's on-disk run against its own in-memory run
py_02_memory.png	Peak memory, sparse vs lazy, before and after the fix
py_03_storage.png	On-disk size: truecell LazyMatrix vs BPCells

Accuracy vs R (13,714 genes × 2,638 cells):

Comparison	Agreement
Anchors matching Seurat	14 / 14 (7 exact)
Variable features shared with Seurat's in-memory run	1998 / 2000
truecell on-disk vs in-memory, normalised values	0 — bit-identical (Seurat's own: 1.02e-06)
truecell on-disk vs in-memory, variable features agreeing	2000 / 2000 (Seurat's own: 1999 / 2000)
DE tests supported out of core	all 8 (Seurat's IterableMatrix: wilcox only)

Found and fixed seven defects — five functions densified an entire on-disk layer on read, so going out of core made peak memory worse, not better; object construction ended laziness before any analysis could run; and a NumPy LOESS fit was chaotically dependent on sort order, moving fitted values by up to 28.8 % under a 1e-15 input nudge. All are pre-existing bugs this tutorial was the first to notice, not consequences of going out of core.

---

Tutorial 17 — Visium: the Spatial Container (R vs Python)

Walkthrough: visium_vignette.md

The Visium loader and container — and the first tutorial in the series where Seurat is the one that's wrong. 2,695 in-tissue spots, 10x V1_Mouse_Brain_Sagittal_Anterior.

python  tutorials/visium_tutorial.py    # downloads the bundle (~64 MB)
Rscript tutorials/visium_verify.R
python  tutorials/visium_tutorial.py --report
python  tutorials/generate_visium_plots.py

What you'll learn: - load_visium — Seurat's Load10X_Spatial/Read10X_Image, including the headerless Space Ranger 1.1.0 tissue_positions_list.csv layout - VisiumV2.radius() — reachable where Seurat's own Radius() returns NULL on its current-generation Visium class - Why the Visium slide's fixed 100 µm spot pitch settles, from geometry alone and without consulting either tool, that Seurat stores a diameter in a slot named radius — and why truecell deliberately does not match that

Key output figures (in tutorials/figures_visium/):

Figure	Description
py_01_radius_geometry.png	Spot spacing and the diameter-vs-radius geometry
py_02_spots_on_tissue.png	Spots plotted on the tissue image
py_03_pca_isolation.png	PCA gap isolated to the upstream feature-selection difference

Accuracy vs R (2,695 spots × 32,285 genes):

Comparison	Agreement
Anchors matching Seurat exactly	24 / 24 (17 exactly)
Spot coordinates	max|dx| = max|dy| = 0 across all 2,695 spots
Tissue image, both backends vs R	1.68e-08 (float32 epsilon)
PCA stdev, on Seurat's own 2,000 features	max relative diff 2.49e-05
radius() ratio, truecell vs Seurat (deliberate)	2.000000

Found one truecell bug and aligned three defaults to Seurat — the tissue image came back 255× apart depending on whether matplotlib or Pillow was installed, from the same file; and the in-tissue filter, image resolution and image key defaults were changed to match Seurat's (filter_by_tissue=True, lowres, "slice1"). Seurat's own radius bug — a diameter stored in a slot named radius, with Radius() returning NULL on VisiumV2 — is reported, not matched.

---

Tutorial 18 — Anchor Internals: CCA & RPCA (R vs Python)

Walkthrough: anchors_vignette.md

Tutorial 8 compares the clusterings integration produces. This one compares the anchors underneath them — which mutual-nearest-neighbour pairs each tool calls an anchor, what score it gives them, and what the corrected expression comes out as — on both the v4 (FindIntegrationAnchors/IntegrateData) and v5 (IntegrateLayers) paths. ifnb, 2,400 cells (v4) and the full 13,999 (v5).

Rscript tutorials/export_seuratdata.R ifnb     # one-time counts export
python  tutorials/anchors_tutorial.py          # writes the cell list + HVGs
Rscript tutorials/anchors_verify.R             # writes the Seurat anchors
python  tutorials/anchors_tutorial.py --report
python  tutorials/generate_anchors_plots.py

What you'll learn: - find_integration_anchors / integrate_data — Seurat's v4 FindIntegrationAnchors/IntegrateData, checked anchor-by-anchor rather than only by the clustering that follows - integrate_layers — Seurat's v5 IntegrateLayers, and why it dispatches to a genuinely different algorithm (IntegrateEmbeddings) than the v4 path, not the same code under a new name - Reading compiled Seurat internals (FindWeightsC, IntegrateDataC) by calling them directly on controlled input, rather than inferring their behaviour from the R wrapper

Key output figures (in tutorials/figures_anchors/):

Figure	Description
py_01_anchor_agreement.png	Anchor recall and precision, truecell vs Seurat
py_02_correction.png	Corrected expression over the query half, truecell vs Seurat

Accuracy vs R (v4 path, 2,400-cell subsample, 2,000 shared anchor features):

Comparison	truecell	Seurat
RPCA anchors recovered	100.0 % (649/649)	—
RPCA anchor score, Pearson r	0.99997	—
CCA anchors recovered	99.9 % (2,815 vs 2,814)	—
RPCA embedding, v5 path, dims |r| > 0.99 (full 13,999 cells)	30 / 30	—

Found and fixed eighteen defects across two rounds — the dominant one was RunCCA standardizing each cell where truecell L2-normalized it, which rotates the whole anchor space; on the v5 path, integrate_layers was silently running v4's IntegrateData behind the v5 IntegrateEmbeddings name, and run_pca's randomized SVD drifted enough in its trailing components to corrupt reciprocal PCA, which standardizes per-dimension and is not rotation invariant. RPCA now agrees with Seurat on every v4 anchor and all 30 v5 embedding dimensions, up from 1 of 30.

---

R Seurat → Truecell API Quick Reference

Task	R (Seurat)	Python (Truecell)
Create object	CreateSeuratObject(counts, min.cells, min.features)	create_truecell_object(counts, min_cells, min_features)
% mito genes	PercentageFeatureSet(pbmc, pattern="^MT-")	percentage_feature_set(pbmc, pattern=r"^MT-")
Normalize	NormalizeData(pbmc, method, scale.factor)	normalize_data(pbmc, normalization_method, scale_factor)
CLR (ADT)	NormalizeData(pbmc, method="CLR", margin=2)	normalize_data(pbmc, normalization_method="CLR", margin=2)
SCTransform	SCTransform(pbmc, vars.to.regress="percent.mt")	sctransform(pbmc, vars_to_regress=["percent.mt"])
HVGs	FindVariableFeatures(pbmc, selection.method, nfeatures)	find_variable_features(pbmc, selection_method, nfeatures)
Scale	ScaleData(pbmc, features)	scale_data(pbmc, features)
PCA	RunPCA(pbmc, features, npcs)	run_pca(pbmc, features, n_pcs)
Supervised PCA	RunSPCA(pbmc, graph="wsnn")	run_spca(pbmc, graph="wsnn")
GLM-PCA	RunGLMPCA(pbmc, L=10)	glm_pca(pbmc, n_components=10) — Poisson; family="nb" for negative binomial
Neighbors	FindNeighbors(pbmc, dims)	find_neighbors(pbmc, dims, k_param)
Cluster	FindClusters(pbmc, resolution)	find_clusters(pbmc, resolution, algorithm)
UMAP	RunUMAP(pbmc, dims)	run_umap(pbmc, dims)
Harmony integration	RunHarmony(pbmc, "batch")	run_harmony(pbmc, "batch") / integrate_layers(pbmc, method="harmony", group_by="batch")
CCA/RPCA anchors	FindIntegrationAnchors(list, reduction="cca") → IntegrateData(anchors)	find_integration_anchors(objs, reduction="cca") → integrate_data(anchors)
CCA/RPCA layers	IntegrateLayers(obj, method=CCAIntegration)	integrate_layers(obj, method="cca", group_by="batch") (or "rpca")
Transfer anchors	FindTransferAnchors(reference, query, reduction="pcaproject")	find_transfer_anchors(reference, query, reduction="pcaproject") (or "cca")
Transfer labels	TransferData(anchors, refdata=reference$celltype)	transfer_data(anchors, refdata="celltype")
Project into reference UMAP	ProjectUMAP(query, reference, reduction.model="umap")	project_umap(query, reference)
Map query (annotate + place)	MapQuery(anchors, query, reference, refdata=list(id="celltype"))	map_query(anchors, refdata="celltype")
Leverage scores	LeverageScore(obj)	leverage_score(obj)
Sketch a large dataset	SketchData(obj, ncells=5000, method="LeverageScore")	sketch_data(obj, ncells=5000)
Project sketch → full data	ProjectData(obj, sketched.assay="sketch", reduction="pca")	project_data(full, sketch, refdata={"cluster_full": "seurat_clusters"})
Matrix to disk (out-of-core)	write_matrix_dir(mat, "counts.mat") (BPCells)	write_lazy_matrix(mat, "counts.mat")
Open on-disk matrix	open_matrix_dir("counts.mat") (BPCells)	open_lazy_matrix("counts.mat")
Demultiplex hashtags	HTODemux(obj, assay="HTO")	hto_demux(obj, assay="HTO") — both default to kfunc="clara"
Demultiplex (MULTI-seq)	MULTIseqDemux(obj, assay="HTO")	multiseq_demux(obj, assay="HTO")
Perturbation signature	CalcPerturbSig(obj, gd.class="gene", nt.cell.class="NT")	calc_perturb_sig(obj, labels="gene", nt_class="NT")
Mixscape (CRISPR KO calls)	RunMixscape(obj, labels="gene", nt.class.name="NT")	run_mixscape(obj, labels="gene", nt_class="NT")
Mixscape LDA (guide separation)	MixscapeLDA(obj, labels="gene", nt.label="NT")	mixscape_lda(obj, labels="gene", nt_class="NT")
Mixscape perturbation score	PlotPerturbScore(obj, target.gene.ident="IFNGR2")	plot_perturb_score(obj, target_gene_ident="IFNGR2")
Mixscape DE heatmap	MixscapeHeatmap(obj, ident.1="NT", ident.2="IFNGR2 KO")	mixscape_heatmap(obj, ident_1="NT", ident_2="IFNGR2 KO")
Markers	FindMarkers(pbmc, ident.1)	find_markers(pbmc, ident_1)
All markers	FindAllMarkers(pbmc, only.pos, logfc.threshold)	find_all_markers(pbmc, only_pos, logfc_threshold)
Conserved markers	FindConservedMarkers(pbmc, ident.1, grouping.var)	find_conserved_markers(pbmc, ident_1, grouping_var)
Pseudobulk	AggregateExpression(pbmc, group.by)	aggregate_expression(pbmc, group_by)
Pseudobulk DESeq2	FindMarkers(pbmc, test.use="DESeq2")	find_markers(pbmc, ident_1, test_use="deseq2", sample_col=...)
MAST hurdle DE	FindMarkers(pbmc, test.use="MAST")	find_markers(pbmc, ident_1, test_use="mast")
Bimodal LRT DE	FindMarkers(pbmc, test.use="bimod")	find_markers(pbmc, ident_1, test_use="bimod")
Rename idents	RenameIdents(pbmc, new.ids)	pbmc.rename_idents(mapping_dict)
Subset cells	subset(pbmc, subset = condition)	pbmc.subset(cells=keep_list)
Add assay	pbmc[["ADT"]] <- CreateAssayObject(counts)	obj.assays["ADT"] = create_assay5_object(counts, key="adt_")
Switch assay	DefaultAssay(cbmc) <- "ADT"	feature_plot(..., assay="ADT")
Access metadata	pbmc@meta.data	pbmc.meta_data
Access assay	pbmc[["RNA"]]	pbmc.assays["RNA"]
HVF statistics	HVFInfo(pbmc[["RNA"]])	pbmc.assays["RNA"].meta_data
Active idents	Idents(pbmc)	pbmc.idents

Beyond PCA — supervised PCA and GLM-PCA

Two reductions that answer questions PCA cannot. Both store a DimReduc exactly as run_pca does, so find_neighbors, find_clusters and run_umap take them with a reduction= argument and nothing else changes.

from truecell import run_spca, glm_pca, find_neighbors, run_umap

# Supervised PCA: the gene axes that best explain a graph you already trust.
find_multi_modal_neighbors(cbmc, reduction_list=["pca", "apca"])   # builds "wsnn"
run_spca(cbmc, graph="wsnn", npcs=50)
run_umap(cbmc, reduction="spca", dims=range(30))

# GLM-PCA: a Poisson fit straight on the raw counts — no log, no pseudocount.
glm_pca(pbmc, n_components=10)
find_neighbors(pbmc, reduction="glmpca", dims=range(10))
print(pbmc.reductions["glmpca"].misc["converged"])

# Negative binomial for over-dispersed counts; θ is estimated by ML.
glm_pca(pbmc, n_components=10, family="nb")
print(pbmc.reductions["glmpca"].misc["theta"])          # the fitted dispersion

run_spca maximises vᵀXᵀGXv where PCA maximises vᵀXᵀXv — the same problem with the identity swapped for a cell-cell graph G. So it finds the gene directions that reproduce a neighbourhood structure you have already decided is the right one (typically the WNN graph, which knows about protein as well as RNA). Hand it G = I and you get PCA back exactly. Its value is that the result is a plain linear map from genes to components, so a query dataset can be projected into a reference's space with one matrix multiply — which is why Azimuth maps onto sPCA and not PCA.

glm_pca models counts as counts. The standard pipeline log-normalises and then runs PCA, which quietly assumes constant-variance Gaussian data; a gene averaging 0.1 UMIs and one averaging 100 do not have remotely comparable noise, and the pseudocount you add to survive log(0) distorts exactly the low-expression genes where most of the zeros are. GLM-PCA drops the transform and fits log μ = a[g] + o[c] + U·Vᵀ with a Poisson likelihood, holding the log library size as a fixed offset so sequencing depth is a known quantity rather than something the factors have to spend themselves discovering.

Real UMI counts are usually noisier than Poisson allows — the same gene in two copies of one cell state still varies more than the mean — and a handful of such genes will dominate a Poisson fit. family="nb" swaps in a negative binomial, Var = μ + μ²/θ, with a single shared dispersion θ estimated by maximum likelihood alongside the factors (pass optimize_theta=False and a theta= to pin it). As θ → ∞ it collapses back onto Poisson, so NB never fits worse — only more forgivingly. The fitted θ lands in misc["theta"].

Two things to know. glm_pca fits densely in genes × cells, so pass a few thousand variable features, as you would to run_pca. And check misc["converged"]; if it is False, or misc["deviance"] is still falling steeply at the end, raise max_iter. (With family="nb" and θ being estimated, the deviance is re-scaled as θ moves, so it is only strictly monotone when you pin θ with optimize_theta=False.)

Integrating datasets — Harmony, CCA, RPCA

Two batches of the same tissue rarely line up: a shared cell type sits in a different place in each dataset's PCA, so cells cluster by batch before they cluster by biology. Truecell offers two remedies.

Harmony (run_harmony) corrects an embedding you already have — it takes the joint PCA and iteratively pulls the batches together while keeping cell types apart. It is fast and needs only a group_by column.

Anchor integration (Seurat's CCA/RPCA) works from scratch, without assuming the datasets even share a coordinate system. It builds a shared space for a pair of datasets — by canonical correlation (reduction="cca", the SVD of the cross-covariance) or reciprocal PCA (reduction="rpca") — then keeps only the mutual nearest neighbours across datasets as anchors: cell i here and cell j there, each among the other's closest matches, are almost certainly the same state seen twice. Anchors are scored for neighbourhood consistency, filtered against the raw expression, and then used to pull every query dataset onto the reference.

from truecell import (
    find_integration_anchors, integrate_data, integrate_layers,
    run_harmony, scale_data, run_pca, find_clusters, run_umap,
)

# --- Harmony: correct an existing joint PCA in place -------------------------
run_harmony(pbmc, group_by="batch")             # stores reductions["harmony"]
run_umap(pbmc, reduction="harmony", dims=range(30))

# --- CCA/RPCA anchors: a list of per-batch objects → a corrected assay -------
anchors = find_integration_anchors([ref, query], reduction="cca")   # or "rpca"
merged = integrate_data(anchors)                # active assay is now "integrated"
scale_data(merged)
run_pca(merged)                                 # clusters by cell type, not batch
find_clusters(merged, resolution=0.5)

# --- One-call Seurat v5 path: split one object by batch, integrate, embed ----
integrate_layers(pbmc, method="cca", group_by="batch", new_reduction="integrated")
run_umap(pbmc, reduction="integrated", dims=range(30))

find_integration_anchors is reference-based: objects[reference] (index 0 by default) is the anchor every other dataset is corrected onto. CCA shines when the datasets share structure but differ globally (cross-species, cross- technology); RPCA is stricter and faster, a better fit when the batches are already similar or very large. Both return the same IntegrationAnchors, which is also what v0.3.0's reference mapping is built to consume.

The correction is applied to the log-normalised data of the shared anchor features and stored as an "integrated" assay — so scale_data + run_pca on it is the natural next step. The reference dataset is left untouched.

Reference mapping — annotating a query from an atlas

Integration mixes several datasets into one shared space. Reference mapping is the asymmetric cousin: you keep an annotated atlas fixed and borrow its labels for a new, unlabelled dataset. The anchor machinery is the same — build a shared space, find mutual nearest neighbours, score and filter them — but the reference is never moved and the anchors carry information reference → query.

from truecell import find_transfer_anchors, transfer_data

# reference: annotated atlas (has a "celltype" column). query: new, unlabelled.
# Both normalized + find_variable_features + scale_data, as usual.
anchors = find_transfer_anchors(reference, query, reduction="pcaproject")

# Classification: predict the query's cell types from the reference labels.
pred = transfer_data(anchors, refdata="celltype")
query.add_meta_data(pred["predicted.id"], col_name="predicted.celltype")
query.add_meta_data(pred["prediction.score.max"], col_name="prediction.score")
# pred also has one prediction.score.<class> column per reference class.

# Imputation: carry reference expression (features × ref-cells) onto the query.
ref_expr = reference.get_assay().layer_data("data", features=["CD3D", "MS4A1"])
imputed = transfer_data(anchors, refdata=ref_expr, refdata_features=["CD3D", "MS4A1"])

reduction="pcaproject" (the default) projects the query through the reference's PCA loadings — the principal axes are learned once on the reference and the query is pushed through the same map. Because those axes never saw the query, batch-specific structure the reference lacks simply lands nowhere, which is what makes projection robust for annotation. reduction="cca" learns a joint space instead, for the harder cross-modality / cross-species cases.

transfer_data weights each query cell's anchors with the same distance-weighted, score-scaled Gaussian kernel integrate_data uses, so a query cell surrounded by confident, consistent anchors of one type gets a high prediction.score.max; an ambiguous one gets a low score you can filter on.

Placing the query in the reference UMAP

Annotating the query is half the job; the other half is seeing it on the atlas you already know how to read. project_umap runs the reference's fitted UMAP model in transform-only mode, so the query cells land in the reference's existing embedding rather than in a fresh, unrelated one. map_query composes the whole workflow — transfer the labels and project the UMAP — in a single call.

from truecell import run_pca, run_umap, find_transfer_anchors, map_query, project_umap

# The reference needs a fitted PCA + UMAP; run_umap stashes the umap-learn model
# in reference.reductions["umap"].misc["umap_model"] for transform-only projection.
run_pca(reference)
run_umap(reference)                       # embeds from "pca", keeps the model

anchors = find_transfer_anchors(reference, query, reduction="pcaproject")

# One call: transfer_data writes predicted.id / prediction.score.* onto the
# query's metadata, and project_umap places it in the reference UMAP.
pred = map_query(anchors, refdata="celltype")
query.reductions["ref.umap"]              # the query, on the reference's UMAP

# Or just the projection, without label transfer:
project_umap(query, reference)            # -> query.reductions["ref.umap"]

project_umap is itself a two-step map: it projects the query through the reference's PCA loadings into the reference's PC space (the same "project into a space the query never helped define" logic as pcaproject above), then runs the reference's UMAP model's .transform. A query cell that resembles reference T cells is pinned near the reference's T-cell island and optimised to sit there — so the query overlays the atlas, batch block and all.

Sketching — analysing a million cells through a small subset

A huge atlas is mostly redundant: the common states are thousands of near-identical cells stacked on top of each other, while the rare states — usually the interesting ones — are a handful of points each. Sketching picks a small, information-dense subset, does the expensive clustering / UMAP on that, and projects the answers back onto every cell. The subset is drawn by leverage, not uniformly: a cell in a dense redundant cloud has low leverage (drop it and nothing changes), a cell in a sparse distinctive corner has high leverage (it is the only evidence that corner exists) — so leverage sampling keeps the rare states a uniform sample would lose.

from truecell import sketch_data, project_data, run_pca, run_umap
from truecell.neighbors import find_neighbors
from truecell.clustering import find_clusters

# full: a large, normalized + scaled object. Draw a leverage-weighted sketch.
sketch = sketch_data(full, ncells=5000)   # a standalone object; assay -> "sketch"
# full.meta_data["leverage.score"] now holds the per-cell scores.

# Do the heavy analysis on the small sketch.
run_pca(sketch, n_pcs=30)
find_neighbors(sketch, dims=range(30))
find_clusters(sketch, resolution=0.8)     # -> sketch.meta_data["seurat_clusters"]
run_umap(sketch, dims=range(30))          # keeps the umap-learn model

# Extend the sketch's PCA, UMAP and cluster labels back to *every* cell.
project_data(full, sketch, refdata={"cluster_full": "seurat_clusters"})
full.reductions["pca.full"]               # every cell in the sketch's PC space
full.reductions["ref.umap"]               # every cell on the sketch's UMAP
full.meta_data["cluster_full"]            # every cell's transferred cluster

leverage_score computes the scores with a sparse CountSketch rather than a full SVD of the data — the whole reason sketching scales — so it stays cheap on millions of cells. project_data is the mirror image of reference mapping: the sketch plays the role of the reference, and each full cell is pushed through the sketch's PCA loadings (and its fitted UMAP model) exactly as project_umap places a query on an atlas, with find_transfer_anchors + transfer_data carrying the cluster labels across.

leverage_score and sketch_data default to nsketch=5000; the CountSketch is a faithful embedding once it has more rows than roughly the squared feature count, which the default comfortably clears for the few-thousand variable features a sketch runs on.

Lazy on-disk matrices — keeping a million cells out of RAM

Sketching shrinks how many cells you analyse; a lazy matrix shrinks how much of the matrix is in memory at once. A dense million-by-twenty-thousand float64 matrix is 160 GB; even the sparse counts leave no room for the copies each step makes. Seurat's answer is the BPCells package, which keeps the matrix on disk and streams over it. LazyMatrix is the truecell analogue, built on NumPy's memory-mapping — no new dependency.

A matrix is written to a directory as the three memory-mapped arrays of a compressed-sparse-column matrix (scipy's csc_matrix layout). Opening it maps those arrays without reading them; a slice pulls only the touched cells off disk and hands back an ordinary scipy.sparse block, so it drops straight into an assay layer:

from truecell import write_lazy_matrix, open_lazy_matrix

assay = obj.get_assay()

# Persist the counts layer out-of-core, then map it back in and swap it in place.
write_lazy_matrix(assay.layers["counts"], "counts.mat")
lazy = open_lazy_matrix("counts.mat")
assay.set_layer_data("counts", lazy)          # a LazyMatrix is a valid layer

# Slicing reads only the selected cells' non-zeros off disk...
block = assay.layer_data("counts", cells=obj.cell_names()[:1000])

# ...and reductions stream in a single pass without materialising the matrix.
per_cell = lazy.sum(axis=0)                    # nCount, one pass over the store
per_gene = lazy.mean(axis=1)                   # per-feature mean

# col_blocks is the streaming primitive: walk a million cells at bounded RAM.
for start, stop, chunk in lazy.col_blocks(block_size=50_000):
    ...                                        # chunk is a csc_matrix of those cells

LazyMatrix stores columns (cells) contiguously, because the operations that dominate at scale — sketching, cell subsetting, per-cell normalisation — select cells, and CSC makes reading an arbitrary set of columns cost only their own non-zeros. as_dense(lazy) / np.asarray(lazy) still materialise the whole thing when you genuinely need it — the escape hatch you keep for the small datasets and avoid on the million-cell path.

Cell hashing — demultiplexing pooled samples

Cell Hashing (Stoeckius et al.) tags each sample with a distinct antibody-oligo hashtag before pooling the samples on one lane. Every droplet then carries a little vector of hashtag counts saying which sample it came from — and droplets that caught two cells carry two tags. hto_demux (Seurat's HTODemux) turns that hashtag matrix back into per-cell calls: singlet (one tag), doublet (two or more), or negative (none).

The hashtags live in their own assay alongside the RNA. hto_demux learns each tag's positive cutoff from the data — it CLR-normalizes the counts, clusters the cells into k = n_hashtags + 1 groups, fits a negative binomial to each tag's background (the cluster where it is least expressed), and thresholds at the 0.99 quantile:

from truecell import hto_demux

# `obj` has an "HTO" assay of hashtag counts (features = hashtags, columns = cells).
hto_demux(obj, assay="HTO")                    # positive_quantile=0.99 by default

# Per-cell results land in meta_data under the Seurat column names:
obj.meta_data["HTO_classification.global"]     # "Singlet" / "Doublet" / "Negative"
obj.meta_data["HTO_maxID"]                      # the top hashtag per cell
obj.meta_data["HTO_classification"]            # tag name (singlet) or "HTOa_HTOb" (doublet)
obj.meta_data["hash.ID"]                        # tag name / "Doublet" / "Negative"

# hash.ID is also set as the active identity, so keeping only the clean singlets
# of one sample is a one-liner:
singlets = obj.meta_data["HTO_classification.global"] == "Singlet"
sample3 = obj.subset(cells=obj.meta_data.index[
    singlets & (obj.meta_data["hash.ID"] == "HTO-3")
].tolist())

# The learned per-hashtag cutoffs are kept for inspection:
obj.misc["hto_demux"]["HTO"]["cutoffs"]         # {"HTO-1": 6.0, "HTO-2": 5.0, ...}

By default hto_demux CLR-normalizes internally at margin=1 — per hashtag, across cells — which is what Seurat's hashing vignette does and what HTODemux assumes. If you already ran normalize_data(obj, normalization_method="CLR", margin=1, assay="HTO"), pass normalize=False to reuse that data layer. (Note this is the opposite margin to the ADT/CITE-seq recipe above, where margin=2 centers each cell across a small protein panel; hashing wants each tag centered across cells.) The negative binomial — not a fixed threshold — is what makes the call robust to each antibody's own staining background and each run's own depth.

Choosing the clustering. The clustering only decides which cluster is each tag's background, so it rarely changes the calls. truecell ships both and defaults to clara (k-medoids), matching Seurat:

hto_demux(obj, assay="HTO")                    # clara, nsamples=100 -- as in Seurat
hto_demux(obj, assay="HTO", kfunc="kmeans")    # the alternative

On synthetic panels the two put ~1% of cells in different classes, rising with tag count (~3.5% at 12 tags, where clara is also the more accurate). Both scale linearly in cells; clara costs a roughly constant 4× (~1.3 s vs ~0.3 s at 100k cells), so pick on fidelity, not speed.

clara takes no seed — its sampling runs off a generator R's set.seed cannot reach either, so it is deterministic in the data alone. One caveat worth knowing if you diff against R: clara accepts a swap on any improvement below zero (R really does take swaps worth -2.2e-16), so its answer is decided at the last bit, and R's own clara returns different clusterings on arm64 vs x86_64. truecell follows plain IEEE arithmetic — identical to R's on x86_64, and to R everywhere the two agree. See ROADMAP.md for the measurements.

MULTI-seq — the other demultiplexer. MULTI-seq (McGinnis et al.) uses lipid-anchored barcodes instead of antibodies, and Seurat's MULTIseqDemux calls it a different way: rather than a background fit, it reads each barcode's cutoff straight off the shape of its distribution. multiseq_demux is a drop-in alternative that often disagrees with hto_demux at the margins, so it is handy as a second opinion:

from truecell import multiseq_demux

# Same "HTO" assay of barcode counts. For each barcode a Gaussian KDE exposes its
# background and positive modes; the cutoff sits a fraction `quantile` between them.
multiseq_demux(obj, assay="HTO", quantile=0.7)     # 0.7 is the Seurat default

obj.meta_data["MULTI_ID"]                          # barcode name / "Doublet" / "Negative"
obj.meta_data["MULTI_classification"]              # a character copy of MULTI_ID
obj.misc["multiseq_demux"]["HTO"]["thresholds"]    # learned per-barcode cutoffs

# Don't want to guess `quantile`? Let it sweep for the value that maximises the
# singlet rate, peeling off negatives and re-thresholding until it settles:
multiseq_demux(obj, assay="HTO", autothresh=True)  # ignores `quantile`

MULTI_ID is set as the active identity, so subsetting one sample's singlets works exactly as with hash.ID above. As with hto_demux, pass normalize=False to reuse an existing CLR data layer instead of recomputing it.

Pooled CRISPR screens — Mixscape

In a pooled CRISPR screen every cell carries a guide RNA, but carrying a guide is not the same as being perturbed: some cells escape the knockout and look just like controls. Mixscape (Papalexi, Mimitou et al.) separates the true knockouts (KO) from those non-perturbed escapers (NP) so downstream analysis runs on genuinely perturbed cells. It is a two-step workflow — build a per-cell perturbation signature, then classify against it — mirroring Seurat's CalcPerturbSig + RunMixscape. It expects a guide-assignment column (each cell's target gene, with the controls labelled "NT") and a computed reduction (pca):

from truecell import calc_perturb_sig, run_mixscape

# 1. Local perturbation signature: subtract each cell's 20 nearest NT controls
#    (in PCA space) to cancel cell-cycle / depth / batch variation. Stored as a
#    new "PRTB" assay. `split_by` keeps neighbours within a replicate.
calc_perturb_sig(obj, assay="RNA", labels="gene", nt_class="NT",
                 reduction="pca", ndims=15, num_neighbors=20)   # → obj.assays["PRTB"]

# 2. Per guide: DE vs NT picks the response genes, then an iterative 2-component
#    Gaussian mixture over the perturbation score splits KO from NP.
run_mixscape(obj, assay="PRTB", labels="gene", nt_class="NT", de_assay="RNA")

obj.meta_data["mixscape_class"]           # "IFNGR2 KO" / "IFNGR2 NP" / "NT"
obj.meta_data["mixscape_class.global"]    # "KO" / "NP" / "NT"
obj.meta_data["mixscape_class_p_ko"]      # KO posterior per guide cell (NaN for NT)
obj.misc["mixscape"]["PRTB"]["genes"]     # per-gene DE-gene / iteration / KO counts

mixscape_class is set as the active identity, so obj.subset(idents="IFNGR2 KO") pulls just the confirmed knockouts. A guide with too few cells, or too few DE genes to show a phenotype (min_de_genes, default 5), has all its cells called NP. For a knock-down rather than a knockout screen, pass prtb_type="KD" — the class labels and the mixscape_class_p_kd posterior column follow the name.

Separating the guide populations — mixscape_lda

run_mixscape asks which cells are perturbed. The complementary question is how the guide populations differ from each other and from control, and mixscape_lda (Seurat's MixscapeLDA) answers it with a single supervised map on which each guide class forms its own cloud. It needs only the PRTB signature — it groups cells by their raw guide label, so the KO/NP calls are not used and it can follow calc_perturb_sig directly:

from truecell import mixscape_lda

# Per guide: DE vs NT picks its response genes, a PCA is fit on that guide's cells
# plus the NT cells, and every cell is projected onto the guide's npcs-dim subspace.
# The blocks are concatenated and one LDA is fit with the guide as the class.
mixscape_lda(obj, assay="PRTB", labels="gene", nt_class="NT", npcs=10)

obj.reductions["lda"]                      # n_classes - 1 discriminant dimensions
obj.reductions["lda"].misc["genes_used"]   # guides that contributed a block
obj.meta_data["lda_assignments"]           # predicted guide class per cell
obj.meta_data["LDAP_IFNGR2"]               # posterior for that class, per cell

Plot it like any other reduction — dim_plot(obj, reduction="lda", group_by="gene"). A guide needs at least npcs + 1 DE genes to contribute a block (it cannot support npcs components otherwise) and is skipped if it falls short; if no guide clears the bar, mixscape_lda raises and you should lower npcs. Because the grouping is by guide rather than by mixscape class, NP escapers stay in their guide's group — and, being control-like, they generally land on top of the NT cloud, which is itself a useful read on how much of a guide escaped.

Checking the calls — plot_perturb_score and mixscape_heatmap

Mixscape's KO/NP split is a threshold on one number per cell: the perturbation score, each cell's projection onto its guide's response axis. plot_perturb_score (Seurat's PlotPerturbScore) draws that axis, overlaying the NT control density against the guide's own — the single most useful check that a guide worked.

from truecell import plot_perturb_score, mixscape_heatmap

fig = plot_perturb_score(obj, target_gene_ident="IFNGR2")
fig.savefig("perturb_score.png", dpi=150, bbox_inches="tight")

A guide with a real effect is bimodal: one lobe sitting on the NT curve (the escapers) and one shifted away from it (the knockouts) — exactly the structure the mixture model is asked to find. A guide that simply did not work is one curve on top of the controls, and no threshold will rescue it. By default the curves are coloured by mixscape_class, so you see where mixscape actually drew the line; pass before_mixscape=True for the raw guide label instead, which is the view you would have had without mixscape at all. For a screen spanning several cell types, split_by="celltype" facets it.

mixscape_heatmap (Seurat's MixscapeHeatmap) then shows the genes underneath that score — the DE genes between two classes, with every cell ordered by its knockout probability:

fig = mixscape_heatmap(obj, ident_1="NT", ident_2="IFNGR2 KO",
                       max_genes=20, balanced=True)

Read with the class colour bar along the top, a clean screen shows the expression block turning on in step with the probability, the low-probability escapers at one end still looking like control. ident_1 / ident_2 are mixscape_class levels ("NT", "IFNGR2 KO", "IFNGR2 NP"), which run_mixscape also leaves as the active identity. balanced=True takes up to max_genes from each direction of the fold change rather than only the up-regulated ones, and max_cells_group downsamples each class for a large screen.

Both plots read the perturbation score that run_mixscape stores in obj.misc["mixscape"]["PRTB"]["genes"], so they need run_mixscape to have run — unlike mixscape_lda, which needs only calc_perturb_sig. A guide whose cells were called NP without a mixture fit (too few cells, or too few DE genes) has no score axis, and plot_perturb_score says so rather than drawing an empty panel.

Pseudobulk & conserved markers

Two multi-sample DE helpers (mirroring Seurat's AggregateExpression and FindConservedMarkers):

from truecell import aggregate_expression, find_conserved_markers

# Pseudobulk: sum raw counts per (cell type × donor) → features × groups DataFrame.
# Pass return_object=True to get a Truecell object with one "cell" per group instead
# (the standard input for pyDESeq2-style sample-level testing).
pb = aggregate_expression(obj, group_by=["cell_type", "donor"])

# Conserved markers: genes up in cluster "B" in *every* condition. Runs FindMarkers
# per level of grouping_var, keeps genes significant in all, and combines their
# p-values with Fisher's method (the `combined_p_val` column; `max_pval` is the
# worst single-condition p-value).
cons = find_conserved_markers(obj, ident_1="B", grouping_var="condition",
                              only_pos=True)
cons.head()   # per-condition stats + max_pval + combined_p_val, sorted by combined_p_val

Pseudobulk DESeq2 (find_markers(test_use="deseq2")) tests between conditions rather than between clusters: set obj.idents to the two conditions, aggregate to one profile per replicate (sample_col), and run DESeq2 on those samples. Needs pip install truecell[deseq2]:

obj.idents = obj.meta_data["condition"]              # e.g. "stim" vs "ctrl"
de = find_markers(obj, ident_1="stim", ident_2="ctrl",
                  test_use="deseq2", sample_col="donor")
de.head()   # p_val / avg_log2FC (DESeq2 log2FoldChange, +ve = up in stim) / pct.1 / pct.2 / p_val_adj

Plotting

R (Seurat)	Python (Truecell)
VlnPlot(pbmc, features, slot)	vln_plot(pbmc, features, layer)
FeaturePlot(pbmc, features)	feature_plot(pbmc, features, assay)
DimPlot(pbmc, reduction, label, pt.size)	dim_plot(pbmc, reduction, label, pt_size)
ElbowPlot(pbmc)	elbow_plot(pbmc)
FeatureScatter(pbmc, feature1, feature2)	feature_scatter(pbmc, feature1, feature2)
VariableFeaturePlot(pbmc)	variable_feature_plot(pbmc)
VizDimLoadings(pbmc, dims, reduction)	viz_dim_loadings(pbmc, dims, reduction)
DimHeatmap(pbmc, dims, cells)	dim_heatmap(pbmc, dims, cells)
DoHeatmap(pbmc, features)	do_heatmap(pbmc, features)
RidgePlot(pbmc, features, ncol)	ridge_plot(pbmc, features, ncol)
ImageDimPlot(obj, group.by)	image_dim_plot(obj, group_by)
ImageFeaturePlot(obj, features)	image_feature_plot(obj, feature)
SpatialDimPlot(obj, group.by)	spatial_dim_plot(obj, group_by) — spots over the H&E image
SpatialFeaturePlot(obj, features)	spatial_feature_plot(obj, feature)

Spatial

R (Seurat)	Python (Truecell)
LoadXenium(dir) / Load10X_Spatial / LoadNanostring	load_xenium(dir) / load_visium(dir) / load_cosmx(dir)
LoadVizgen(dir) (MERSCOPE)	load_merscope(dir) — drops Blank-* controls by default
GetTissueCoordinates(obj)	get_tissue_coordinates(obj)
FNN::get.knn(coords, k) / get.knnx	spatial_knn(coords, k, query)
FNN::get.knn (nearest same-type)	nearest_neighbor_distance(obj, group_by, reference)
(hand-rolled neighbourhood counts)	local_neighborhood(obj, group_by, reference, k)
BuildNicheAssay(obj, fov, group.by, niches.k)	build_niche_assay(obj, group_by, k, niches)
FindSpatiallyVariableFeatures(obj, method="moransi")	find_spatially_variable_features(obj, k=10) — Moran's I
FindSpatiallyVariableFeatures(obj, method="markvariogram", r.metric=5)	find_spatially_variable_features(obj, method="markvariogram", r_metric=5) — r_metric is in cell spacings, not pixels
(hand-rolled Fisher + p.adjust)	composition_test(obj, group_by, split_by)
GetImage(obj[["slice1"]])	obj.images["spatial"].get_image() — the Visium H&E image
ScaleFactors(obj[["slice1"]])	obj.images["spatial"].scale_factors
SpatialDimPlot(obj)	spatial_dim_plot(obj)
SpatialFeaturePlot(obj, features)	spatial_feature_plot(obj, feature)

Plot output: R renders to the graphics device automatically. Truecell functions return a matplotlib.Figure — call fig.savefig("out.png") to save or display inline in Jupyter.

Visium tissue images

load_visium reads the H&E tissue image and scalefactors_json.json by default, giving each image slot a VisiumV2 (Seurat v5's class) instead of a bare FOV:

from truecell import load_visium

obj = load_visium("visium_out/")            # image=True by default
fov = obj.images["spatial"]

fov.get_image().shape                       # (H, W, 3) — tissue_hires_image.png
fov.scale_factors.spot                      # spot diameter, full-resolution pixels
fov.radius()                                # spot radius, full-resolution pixels

# Spot coordinates stay in FULL-RESOLUTION pixels, so spatial_knn /
# nearest_neighbor_distance / Moran's I keep measuring real distances.
# Convert to the stored image's pixel space only when you draw:
xy = fov.scale_coordinates()                # x, y scaled onto the PNG
r  = fov.spot_radius()                      # matching spot radius, in image pixels

Pass image_resolution="lowres" for the smaller PNG, image=False to skip it entirely, or filter_by_tissue=True to keep only spots with in_tissue == 1. A bundle with no PNG still loads — you get a plain FOV, exactly as before.

Plotting spots on the tissue

spatial_dim_plot / spatial_feature_plot (SpatialDimPlot / SpatialFeaturePlot) draw the H&E photo and overlay the spots on top of it — the scaling above happens for you:

from truecell import spatial_dim_plot, spatial_feature_plot

fig = spatial_dim_plot(obj, group_by="seurat_clusters")
fig = spatial_feature_plot(obj, "Gad1", image_alpha=0.4)   # fade the tissue
fig.savefig("visium.png", dpi=150, bbox_inches="tight")

Spots are drawn at their true diameter (from spot_diameter_fullres), not as fixed-size points, so they stay registered against the tissue at any zoom. pt_size_factor= (default 1.6, as in Seurat) scales them; crop=False shows the whole slide instead of zooming to the spots; resolution="lowres" draws the smaller PNG.

Neither function needs an image to work. Plot an object loaded with image=False and you get a plain scatter of the same spots — useful for Xenium/CosMx, or for a Visium bundle whose PNG is missing.

Finding spatially variable genes

find_spatially_variable_features (FindSpatiallyVariableFeatures) ranks genes by how strongly their expression is organised in space. Both of R's methods are here:

from truecell import find_spatially_variable_features

# Moran's I — is this gene autocorrelated at all? Comes with a p-value.
svf = find_spatially_variable_features(obj, k=10)
svf.head()          # moransi, moransi_pval, moransi_padj, moransi_rank

# Mark variogram — how much of the gene's variance has decayed by distance r?
mv = find_spatially_variable_features(obj, method="markvariogram", r_metric=5)
mv.head()           # markvariogram, markvariogram_rank

fig = spatial_feature_plot(obj, mv.index[0])   # the most spatially variable gene

Both return a table sorted so that rank 1 is the most spatially variable, and both also write their columns into the assay's feature metadata, the way find_variable_features does. On a large panel, pass features= to score only the variable genes — the mark variogram is the heavier of the two.

The two statistics ask different questions. Moran's I gives one number for the whole slide and a p-value with it. The mark variogram is read at a distance: markvariogram is the average squared expression difference between cells about r_metric apart, divided by the gene's own variance. So ≈ 1 means no spatial structure (two cells that far apart differ as much as any two cells picked at random) and below 1 means they still resemble each other. There is no p-value — the variogram has no closed-form null, and R does not offer one either.

r_metric is not in R's units. R passes r.metric straight through to spatstat in raw coordinate units, so the same script answers differently on a slide measured in pixels and one in microns. Here r_metric is measured in nearest-neighbour spacings, so the default of 5 means "five cells apart" on any slide. Widen bandwidth= (also in spacings) if the tissue is sparse and too few cell pairs land near r_metric to average over.