Skip to content

Truecell Roadmap

This document tracks features planned for future releases, organized by milestone. Each item includes the R Seurat equivalent, implementation notes, and dependencies so any item can be picked up and scoped independently.

What v0.1.0 already covers (not listed below):
LogNormalize · CLR · VST · ScaleData (+ covariate regression) · SCTransform · AddModuleScore · CellCycleScoring · PCA · UMAP · KNN/SNN · Louvain/Leiden · FindMarkers/FindAllMarkers (wilcox/t/LR/negbinom/roc) · JackStraw · DimPlot · FeaturePlot · VlnPlot · DotPlot · ElbowPlot · DoHeatmap · DimHeatmap · FeatureScatter · VariableFeaturePlot · RidgePlot · CITE-seq storage (Assay5 multi-layer) · AnnData interop · Spatial data structures (FOV/Centroids/Segmentation/Molecules)

v0.2.0 — Batch Correction & Integration — ✅ complete

Why first: batch effects are unavoidable in real datasets; Harmony is the most widely-used, has a pip-installable Python package, and has a well-defined scope. CCA/RPCA follow naturally.

Status: all three integration paths delivered — Harmony (run_harmony), CCA/RPCA anchors (find_integration_anchors / integrate_data), and the integrate_layers dispatcher (method="harmony"|"cca"|"rpca").

Harmony integration — ✅ delivered

  • Implemented in truecell/integration.py as run_harmony(...); stores a DimReduc("harmony") and is verified to lower per-batch silhouette while preserving cell-type separation (tests/test_integration.py). Enable with pip install truecell[integration] (adds the harmonypy dep).
  • R: RunHarmony(obj, group.by.vars = "batch") (via harmony package)
  • Python dep: harmonypy (pip)
  • Plan:
  • run_harmony(seurat, group_by, theta, lambda_, sigma, nclust, max_iter, random_seed) → stores "harmony" in obj.reductions
  • Input: PCA embeddings from obj.reductions["pca"].cell_embeddings
  • Output: DimReduc("harmony", embeddings) — same shape as PCA
  • Downstream: pass reduction="harmony" to find_neighbors / run_umap
  • Tests: corrected embeddings have lower silhouette separation by batch than raw PCA

CCA / RPCA integration (IntegrateData v4 API) — ✅ delivered

  • Implemented in truecell/anchors.py as find_integration_anchors(...) + integrate_data(...), with the IntegrationAnchors result container. Anchors are mutual nearest neighbours in a shared CCA (SVD of the cross-covariance AᵀB) or reciprocal-PCA space, scored by neighbourhood consistency and filtered against the raw expression space; integrate_data corrects each query onto the reference with a Gaussian-weighted sum of anchor correction vectors and returns a merged object carrying an active "integrated" assay. Verified (tests/test_integration.py) to cluster by cell type, not batch.
  • RPCA (found in T6, integration_vignette.md): the reduction="rpca" path had two real defects the synthetic tests missed, both now fixed — a crash on unequal batch sizes (fixed in #41) and a 4× under-integration (batch-mix 0.222 → 0.867 vs Seurat's 0.914; fixed by per-object scaling + Seurat's reciprocal-embedding SD/L2 normalization + disabling the RPCA anchor filter, all with regression tests). reduction="cca" and run_harmony match Seurat to three decimals; RPCA now integrates well, ~5% off Seurat on the residual (exact-vs-annoy NN, sklearn-vs-irlba PCA).
  • R: FindIntegrationAnchors(list, reduction="cca") → IntegrateData(anchors)
  • Implemented:
  • find_integration_anchors(objects, reduction, dims, k_anchor, k_filter, k_score, reference) → IntegrationAnchors
  • integrate_data(anchors, new_assay, k_weight, sd_weight) → corrected "integrated" assay on merged object
  • CCA: numpy.linalg.svd on the cross-covariance; RPCA: reciprocal PCA projections, MNN via sklearn.neighbors
  • Note: reference-based — anchors link each dataset to reference=0, and every other dataset is corrected onto it (one of Seurat's supported modes). A full guide-tree over all pairwise anchors is a later refinement.
  • Reuse: the same IntegrationAnchors object is what v0.3.0's reference mapping (FindTransferAnchors / TransferData) is built to consume.

IntegrateLayers (Seurat v5 API) — ✅ delivered (harmony + cca + rpca)

  • Implemented as integrate_layers(obj, method=..., group_by=...) in truecell/integration.py. method="harmony" corrects an existing reduction; method="cca" / "rpca" split the object by group_by, run the anchor pipeline above, and store the batch-corrected embedding as a new reduction.
  • R: IntegrateLayers(obj, method = HarmonyIntegration, orig.reduction = "pca")
  • Dep: Harmony and CCA/RPCA functions above

v0.3.0 — Reference Mapping & Label Transfer — ✅ complete

Why: label transfer from a curated atlas to a query dataset is a standard first-annotation step; all machinery (KNN, PCA, anchors) already exists.

Status: all three pieces delivered — find_transfer_anchors (pcaproject / cca) and transfer_data (classification + imputation) in truecell/transfer.py, plus project_umap / map_query in truecell/mapping.py (placing the query in the reference UMAP). Reuses the anchors.py machinery end-to-end and adds no dependency.

FindTransferAnchors — ✅ delivered

  • Implemented as find_transfer_anchors(reference, query, reduction="pcaproject") in truecell/transfer.py, returning a TransferAnchors container (anchor pairs cell1/cell2/score + the query embedding used for weighting). Default pcaproject projects the query through the reference's PCA loadings (computed on the shared anchor features, so the reference need not already carry a pca reduction); cca builds a jointly-learned space. MNN → score → filter reuse the same helpers as find_integration_anchors (tests/test_transfer.py).
  • R: FindTransferAnchors(reference, query, dims = 1:30)
  • Dep: CCA/RPCA anchors from v0.2.0 (truecell/anchors.py)

TransferData — ✅ delivered

  • Implemented as transfer_data(anchors, refdata) in truecell/transfer.py. Builds a per-query-cell weight over the anchors (the same distance-weighted, score-scaled Gaussian kernel integrate_data uses), then either classifies (a metadata column or 1-D label array → predicted.id + prediction.score.<class> per class, rows summing to 1, + prediction.score.max) or imputes (a 2-D features × reference-cells matrix → predicted query expression). Verified to recover a query's true cell types across an injected batch block (tests/test_transfer.py).
  • R: TransferData(anchors, refdata = reference$celltype)
  • Tests: query cells get the correct transferred label at >85% accuracy; per-class scores form a distribution; imputation recovers the marker split.

MapQuery + ProjectUMAP — ✅ delivered

  • Implemented in truecell/mapping.py as project_umap(query, reference) and the map_query(anchors, refdata="celltype") convenience. project_umap is a two-step map: project the query through the reference's PCA loadings into the reference's PC space, then run the reference's fitted UMAP model in transform-only mode (umap-learn's UMAP.transform), so the query lands in the reference's existing embedding. It aligns the loadings to only the reference-PCA features the query actually carries scaled, so a query missing a few variable features still projects. Result stored as query.reductions["ref.umap"].
  • map_query composes the whole workflow from a TransferAnchors: transfer_data the labels onto query.meta_data (predicted.id / prediction.score.*), then project_umap the query into the reference UMAP — one call from anchors to an annotated, atlas-placed query. Verified (tests/test_mapping.py) that projected query cells land nearest the matching reference type's UMAP centroid (>85%) despite an injected batch block, and that map_query annotates + places in one call.
  • R: MapQuery(anchorset, query, reference, refdata = ...) then ProjectUMAP(...)
  • Note: the fitted model lives in reference.reductions["umap"].misc["umap_model"], stored by run_umap when embedding from a reduction (not a graph).

v0.4.0 — Weighted Nearest Neighbor (WNN)

Why: truecell already stores RNA + ADT assays; WNN is the natural joint analysis step for CITE-seq data and is well-scoped.

FindMultiModalNeighbors — ✅ delivered (full port)

  • Implemented as find_multi_modal_neighbors(...) in truecell/multimodal.py. Both WNN stages are ported from the R/C++ source; stores wknn/wsnn graphs and <assay>.weight columns. Verified on synthetic complementary-modality data to recover structure RNA alone cannot (tests/test_multimodal_wnn.py).
  • R: FindMultiModalNeighbors(obj, reduction.list = list("pca","apca"), dims.list = list(1:30, 1:18))
  • Stage 1 — FindModalityWeights (_modality_weights): L2-normalise each embedding; impute each cell from its own modality's neighbours and from the other's; d = ||x - x_hat|| - d(nearest neighbour), ReLU'd; per-cell kernel bandwidth from the SNN graph; exp(-d / sigma); score within / (cross + 1e-4) clipped to [0, 200]; softmax across modalities.
  • Stage 2 — MultiModalNN (_multi_modal_nn): each modality nominates knn_range = 200 candidates, they are unioned per cell, and each is scored by sum_r exp(-d_r / sigma_r) * weight_r. The top k_nn become the cell's joint neighbours; wknn/wsnn are built from that ranking.
  • Tests: WNN clusters CBMC data closer to protein-defined ground truth than RNA alone

Superseded: the earlier weighting approximation

The first implementation approximated stage 1 with a linear distance ratio, theta_m = d_cross / (d_same + d_cross), normalised across modalities — and skipped stage 2 entirely, blending per-modality SNN graphs as SNN_wnn = w * SNN_RNA + (1-w) * SNN_ADT instead. It was monotone in the right quantity, so the direction was right, but the linear ratio has no dynamic range: on CBMC every cell landed in 0.46–0.53, versus R's 0.21–0.65. A weight pinned near 0.5 cannot say "this cell is decided by protein", which is the one thing WNN exists to say. The exponential kernel plus the clipped softmax is what supplies the range, and the joint neighbour search is what turns it into a graph. Caught by putting the truecell and Seurat weight violins side by side in Tutorial 3 — the tutorial figure was the test.

Four details in the R source that the natural reading gets wrong, all load-bearing: - FindMultiModalNeighbors never passes prune.SNN down to FindModalityWeights, so the bandwidth's SNN graph uses prune = 0, not 1/15. - SNN_SmallestNonzero_Dist (src/snn.cpp) averages distances to the least-similar SNN partners, and on ties at the k-th weight keeps the k largest distances. - PredictAssay silently drops the self column, so imputation averages k_nn - 1 neighbours. - R requests k.nn neighbours including self, not k.nn + 1.

  • Known scaling limit: _multi_modal_nn loops per cell over the unioned candidate pool. Fine at CBMC's ~8.6k cells (Tutorial 3 runs end to end in ~58s), but it is the bottleneck if WNN is pointed at a much larger dataset. Vectorising it is the obvious follow-on if that comes up.

Fixed: CLR's margin flag was inverted against Seurat

Not a WNN bug, but found through WNN and worth recording next to it. Seurat's CustomNormalize applies the CLR kernel via apply(data, MARGIN = margin, clr_function), and R's apply treats MARGIN = 1 as rows and MARGIN = 2 as columns. With counts stored features × cells that makes Seurat's margin=1 per-feature across cells (its default) and margin=2 per-cell across features (what ADT panels want). _clr_normalize had the two swapped, so normalize_data(..., margin=2) computed what Seurat's margin=1 computes. The per-vector kernel was always exact, including Seurat's quirk of summing log1p over non-zero entries while dividing by the full length — only the axis was wrong.

Consequences, all now resolved: the ADT matrix feeding apca was the wrong transform, so every ADT.weight was off (erythroid 0.62 vs Seurat's 0.40, platelet 0.53 vs 0.30, and the whole table biased high); and the annotation thresholds in cbmc_citeseq_verify.R had been retuned around the discrepancy and documented as a CLR "scale" difference between the languages, which it was not. Both scripts now share one set of thresholds.

hto_demux and multiseq_demux were accidentally correct — they passed margin=2 and, under the inverted code, got per-hashtag-across-cells, which is numerically what Seurat's hashing vignette does at its default margin=1. Their defaults moved 2 → 1 in the same change so their behaviour is unchanged; that coupling is the trap in this fix.

It survived because test_clr_matches_seurat_formula verified the kernel but derived the axis mapping in Python, restating the same inverted assumption it was meant to check — right formula, assumed axis, exactly the failure mode of the superseded weighting above. test_clr_margin_matches_r_ground_truth replaces that with fixed output captured from a real Seurat 5.5.1 run; the old implementation fails it on both margins.

Breaking: any caller passing margin explicitly to normalize_data, hto_demux or multiseq_demux gets different output than before. Callers who relied on the defaults are unaffected.

WNN UMAP + clustering — ✅ delivered

  • find_clusters(graph_name="wsnn") already routed correctly; run_umap now accepts a graph= kwarg that embeds a precomputed graph via UMAP's simplicial_set_embedding (tests/test_reductions_extra.py).
  • CBMC CITE-seq tutorial (Tutorial 3) Step 8 covers WNN end to end, with R-side figures for the joint embedding and the modality weights.
  • Agreement with R: both sides find 16 RNA clusters and 21 WNN clusters at resolution = 0.6 and resolve the same nine lineages with the same multiplicities. Eight of nine per-cell-type ADT.weight means match Seurat to 0.02 or better; progenitor is the exception at 0.06, on a 146-cell population the two sides do not cut identically. truecell's neighbour search is exact where R's is approximate (annoy) and the Louvain implementations differ, so cluster boundaries still move slightly — small populations feel it most.

v0.5.0 — Additional Dimensionality Reductions — ✅ complete

t-SNE and ICA are single functions wrapping a scikit-learn call. sPCA and GLM-PCA are not — both are implemented directly against NumPy/SciPy, and neither added a dependency.

run_tsne — ✅ delivered

  • Implemented in truecell/reduction.py (run_tsne), mirrors run_umap; stores DimReduc("tsne"). Tested in tests/test_reductions_extra.py.
  • R: RunTSNE(obj, dims = 1:10)
  • Python dep: scikit-learn (TSNE) — already a dep

run_ica — ✅ delivered

  • Implemented in truecell/reduction.py (run_ica); stores DimReduc("ica", embeddings + loadings); find_neighbors/run_umap accept reduction="ica". Tested in tests/test_reductions_extra.py.
  • R: RunICA(obj, nics = 30)
  • Python dep: sklearn.decomposition.FastICA

run_spca (supervised PCA) — ✅ delivered

  • Implemented in truecell/reduction.py as run_spca(obj, graph="wsnn") (tests/test_spca_glmpca.py).
  • The plan previously written here described the wrong algorithm — a gene-graph Laplacian smoothing a gene × gene matrix. Seurat's RunSPCA takes a cell × cell graph (the documented call is RunSPCA(reference, assay = "SCT", graph = "wsnn")) and eigendecomposes XᵀGX, which is features × features. Corrected here.
  • What it does: ordinary PCA maximises vᵀXᵀXv and knows nothing about which cells you consider neighbours. sPCA swaps the identity for a graph you already trust and maximises vᵀXᵀGXv — the gene axes that best reproduce that graph. Pass G = I and PCA falls back out exactly, which is how the implementation is tested (loadings match run_pca to a cosine > 0.999).
  • Why it matters: the output is a linear map from genes to components, so a query dataset can be pushed into a reference's graph-defined space with one matrix multiply. That is why Azimuth maps onto sPCA rather than PCA, and it is the reduction v0.3.0's reference mapping will want.
  • One departure from R: Seurat runs irlba (an SVD) on XᵀGX, ranking components by |λ|; we take the largest eigenvalues themselves, since vᵀXᵀGXv is the quantity being maximised and a graph can push eigenvalues negative. With non-negative edge weights the leading eigenvalues are positive, so the two orderings differ only in the tail.
  • R: RunSPCA(obj, assay, graph)

glm_pca (GLM-PCA) — ✅ delivered (Poisson + negative binomial)

  • Implemented in truecell/glmpca.py as glm_pca(obj, n_components=10), following Townes et al. (2019). Pure NumPy/SciPy — no glmpca-py dependency, in keeping with how MAST and bimod were done (tests/test_spca_glmpca.py).
  • What it does: log-normalise-then-PCA assumes the transformed counts are Gaussian with constant variance. They are not, and the pseudocount needed to survive log(0) distorts exactly the low-expression genes where the zeros live. GLM-PCA drops the transform and fits a low-rank model on the count scale: Y[g,c] ~ Poisson(μ), log μ = a[g] + o[c] + Σ_l U[g,l]·V[c,l], with the log library size as a fixed offset o so sequencing depth is a known quantity rather than a factor to be rediscovered. Factors land in cell_embeddings, loadings in feature_loadings, so find_neighbors(reduction="glmpca") and run_umap work downstream unchanged.
  • Fitting: Fisher scoring, alternating over intercept → loadings → factors. Each block gets a diagonal Newton step (score ÷ Fisher information); under a log link and Poisson noise both are one matrix product. Any step that fails to lower the deviance is rejected and retried at half the step size, so the deviance falls monotonically by construction; the trace is kept in misc["deviance"].
  • Initialisation is load-bearing, not a detail. U = V = 0 is an exact saddle of the log-likelihood — each block's score is a product with the other, so both vanish there. Starting near zero (the obvious choice, and glmpca's own default) leaves the fit inching away from the saddle, and any relative-improvement stopping rule then declares convergence on a model that has fitted nothing. It fails convincingly: one step is enough to orient the factors, so clusters separate cleanly in a plot while the deviance sits at its null value. So the factors are seeded from the SVD of the intercept-only model's residuals instead, as Townes recommends. test_glmpca_actually_fits_rather_than_stalling guards it.
  • Negative binomial: family="nb" fits Y ~ NB(μ, θ), Var = μ + μ²/θ. Poisson understates the overdispersion in most scRNA-seq, and a few noisy genes then dominate a Poisson fit; NB down-weights them. The Fisher scoring loop is the Poisson one with a single divisor 1 + μ/θ on both the residual and the working weight — as θ → ∞ it collapses back onto Poisson exactly. The shared dispersion θ is estimated by maximum likelihood between factor updates (optimize_theta, MASS theta.ml-style Newton seeded from a method-of-moments estimate) or pinned at a value you pass. Stored in misc["theta"] (inf for a Poisson fit). A moving θ re-scales the deviance, so the monotone-deviance guarantee holds only when θ is held fixed (optimize_theta=False).
  • Scale: the fit is dense in genes × cells. Pass a few thousand variable features, as you would to run_pca.
  • R: RunGLMPCA(obj, L = 10) (via SeuratWrappers + glmpca)

v0.6.0 — Pseudobulk DE & Advanced Marker Methods — ✅ complete

AggregateExpression (pseudobulk) — ✅ delivered

  • Implemented as aggregate_expression(...) in truecell/aggregate.py. Sums raw counts per group via a single sparse counts @ indicator matmul; group_by accepts one or more metadata columns (joined with "_", as Seurat does) or "ident". Returns a features×groups pd.DataFrame (a dict for multiple assays), or a Truecell object with one "cell" per group when return_object=True (tests/test_pseudobulk_conserved.py).
  • R: AggregateExpression(obj, group.by = c("celltype","donor"))
  • Intended input for DESeq2-style testing (see below).

DESeq2-style pseudobulk DE — ✅ delivered

  • Implemented as the test_use="deseq2" branch of find_markers (_deseq2_pseudobulk in truecell/markers.py). Sums counts to one pseudobulk profile per (group × sample_col) — the AggregateExpression operation — then fits pydeseq2.DeseqDataSet(design="~condition") and contrasts group 1 vs group 2. Returns Seurat-shaped columns (p_val/avg_log2FC/pct.1/pct.2/p_val_adj); warns below 2 replicates per group. Enable with pip install truecell[deseq2] (tests/test_deseq2_pseudobulk.py).
  • R: FindMarkers(obj, test.use = "DESeq2") (via DESeq2 R package)
  • Python dep: pydeseq2 (pip, optional [deseq2] extra)

MAST — ✅ delivered

  • Implemented as the test_use="mast" branch of find_markers (_mast_pvalue in truecell/markers.py): a pure-Python two-part hurdle LRT — a logistic model of detection (expr > 0) plus a Gaussian model of magnitude among detected cells, each ~ group (+ latent). The combined statistic is the sum of the two components' LR statistics on the sum of their df (components with no signal drop out). No R dep — statsmodels is already present. Pass the cellular detection rate via latent_vars to match Seurat's CDR covariate (tests/test_mast_de.py).
  • R: FindMarkers(obj, test.use = "MAST")

FindConservedMarkers — ✅ delivered

  • Implemented as find_conserved_markers(...) in truecell/markers.py. Runs find_markers independently within each level of grouping_var, keeps genes that are markers in every level, and combines their per-level p-values with Fisher's method (scipy.stats.combine_pvalues). Output has per-level prefixed stats plus max_pval and combined_p_val (sorted by the latter); levels lacking a comparison group are skipped with a warning (tests/test_pseudobulk_conserved.py).
  • R: FindConservedMarkers(obj, ident.1, grouping.var)

bimod test (likelihood-ratio on bimodal model) — ✅ delivered

  • Implemented as the test_use="bimod" branch of find_markers (_bimod_pvalue / _bimod_likelihood in truecell/markers.py), a faithful port of Seurat's DifferentialLRT/bimodLikData: each group's expression is modelled as a point mass at zero (Bernoulli detection rate) plus a Gaussian on the detected values, and 2·(logLik₁ + logLik₂ − logLik_pooled) is tested as χ²(df=3). Pure Python, no new dep (tests/test_bimod_de.py).
  • R: FindMarkers(obj, test.use = "bimod")

v0.7.0 — Spatial Transcriptomics

Delivered. The data structures (FOV, Centroids, Segmentation, Molecules) plus these loaders and analysis functions are done and validated end-to-end against R Seurat in Tutorial 5 (deterministic anchors match to 8 significant figures):

  • Loaders: load_xenium, load_visium, load_cosmx (each returns a Truecell object with coordinates in an FOV slot); spatial-aware from_anndata (rebuilds images from obsm['spatial'])
  • Analysis: get_tissue_coordinates, spatial_knn, nearest_neighbor_distance, local_neighborhood, build_niche_assay (BuildNicheAssay), composition_test, add_module_score(search=), find_spatially_variable_features (FindSpatiallyVariableFeatures) with both the moransi and markvariogram methods
  • Plots: image_dim_plot (ImageDimPlot), image_feature_plot (ImageFeaturePlot) — sub-cellular Xenium/CosMx centroids; spatial_dim_plot (SpatialDimPlot), spatial_feature_plot (SpatialFeaturePlot) — Visium spots over the H&E tissue image

This milestone is complete.

load_merscope — ✅ delivered

  • Implemented in truecell/spatial/loaders.py as load_merscope(...), mirroring load_cosmx: reads cell_by_gene.csv + cell_metadata.csv (center_x / center_y) into a Truecell object with populated images (one per fov). Drops Blank-* control barcodes by default (as LoadVizgen does; override with keep_controls=True) and tolerates both the named and unnamed cell-id column layouts Vizgen emits. Verified to flow through the whole spatial stack (spatial_knn → nearest_neighbor_distance → build_niche_assay) (tests/test_merscope_loader.py).
  • R: LoadVizgen(data.dir) (Vizgen MERSCOPE)
  • File format: cell_by_gene.csv, cell_metadata.csv

FindSpatiallyVariableFeatures — ✅ delivered (both methods)

  • Implemented as find_spatially_variable_features(...) in truecell/spatial/variable_features.py, dispatching on method=. Pure NumPy/SciPy — no libpysal and no spatstat equivalent needed.

method="moransi" (default) builds a row-standardised sparse KNN weight matrix from spatial_knn, then computes I = (N/S0)·(zᵀWz)/(zᵀz) vectorised across all genes in one sparse matmul. Significance uses the closed-form E[I] = −1/(N−1) and normality-assumption variance (computed once, since it depends only on W) → z-score → two-sided p, plus BH adjustment. Writes moransi / moransi_pval / moransi_padj / moransi_rank into the assay's feature metadata (as find_variable_features does). Validated against a brute-force double sum.

method="markvariogram" computes the normalised mark variogram γ(r) = E[½·(m_i − m_j)² | d_ij ≈ r] / Var(m) — the expression difference between cells about r apart, relative to the gene's own variance. γ ≈ 1 means two cells r apart differ as much as two picked at random (no structure); γ < 1 means they still resemble each other. Writes markvariogram / markvariogram_rank; rank 1 = lowest γ. No p-value — the variogram has no closed-form null, and R does not offer one either. Also validated against a brute-force loop over every cell pair (tests/test_markvariogram.py).

  • Two deliberate departures from R, both documented in the docstring:
  • r_metric is in nearest-neighbour spacings, not raw coordinate units. R passes r.metric straight through to spatstat, so the same script answers differently on a slide in pixels and in microns, and the default of 5 is only meaningful if you know your coordinate scale. Here r_metric=5 means "five cells apart" on any slide.
  • γ is a kernel-weighted (Nadaraya-Watson) ratio estimator, not spatstat's translation-corrected one, so absolute γ values are close to but not identical with R's. The gene ranking — what the function is for — carries over.
  • Performance: the pairwise differences are never materialised (that array would be genes × pairs). Because the kernel matrix K is symmetric with a zero diagonal, Σ_{i<j} K_ij·(m_i − m_j)² = Σ_i s_i·m_i² − mᵀKm with s = K·1, which is two sparse products regardless of how many pairs land in the band. A cKDTree range query means only pairs near r are ever built.
  • Caveat documented in the docstring: when a few strongly spatial genes dominate library size, log-normalisation leaks their structure into flat genes and inflates their score — a property of compositional normalisation, not of either statistic.
  • R: FindSpatiallyVariableFeatures(obj, method = "moransi" | "markvariogram")

Visium tissue image (VisiumV2) — ✅ delivered (data layer)

  • truecell/spatial/visium.py adds VisiumV2 (an FOV subclass mirroring Seurat v5's) carrying the H&E image, the ScaleFactors from scalefactors_json.json, and which resolution is stored. load_visium(..., image=True) (the default) now reads spatial/tissue_{hires,lowres}_image.png + spatial/scalefactors_json.json and returns VisiumV2 images; bundles with neither still load as a plain FOV, so this is backwards compatible. Also adds filter_by_tissue= (drops in_tissue == 0 spots from the matrix as well as the coordinates), and registers the previously unimplemented get_image generic on SpatialImage (tests/test_visium_image.py).
  • Coordinate convention: spot coordinates stay in full-resolution pixels (the space tissue_positions.csv uses), so spatial_knn / nearest_neighbor_distance / Moran's I keep seeing real, image-independent distances. VisiumV2.scale_coordinates() and .spot_radius() convert to the stored image's pixel space on demand — one multiply, at draw time.
  • Deps: none added. The PNG is read with matplotlib (or Pillow) lazily; with neither installed the loader warns and skips the image rather than failing.
  • R: Load10X_Spatial(data.dir), GetImage(obj[["slice1"]]), ScaleFactors()

Tissue-image spatial plots (Visium H&E) — ✅ delivered

Function R equivalent Notes
spatial_dim_plot SpatialDimPlot clusters/ident over the H&E tissue image
spatial_feature_plot SpatialFeaturePlot gene expression over the tissue image
  • Both in truecell/plotting.py. Each panel imshows the photo held by the VisiumV2 image, then overlays the spots on top of it (tests/test_spatial_plots.py).
  • Spots are drawn at their true diameter, not as scatter points: they are an EllipseCollection in units="xy", so a spot's size is expressed in data units and stays registered against the tissue when the axes are zoomed or the resolution changes. A scatter(s=…) size is in points² and would drift. pt_size_factor= (default 1.6, as in Seurat) scales relative to the real spot.
  • The image anchors the coordinate space. Coordinates come from VisiumV2.scale_coordinates() (fullres → image pixels) and the diameter from .spot_radius(); imshow then supplies the frame. crop= (default True) zooms to the spots rather than the whole slide, image_alpha= fades the tissue, and resolution= overrides which PNG is drawn.
  • Two independent fallbacks, because the bundle has two independent holes. No PNG (a plain FOV, or load_visium(image=False)) → a bare scatter of the same spots, y-axis still pointing down so it looks the same. No scalefactors_json.json → the photo still draws, but with no spot_diameter_fullres there is nothing to size the spots to, so they degrade to fixed-size scatter points.
  • R: SpatialDimPlot(obj), SpatialFeaturePlot(obj, features = "Gad1")

v0.8.0 — Scale & Performance

Status: ✅ complete. Leverage-score sketching (sketch_data / project_data + leverage_score, truecell/sketch.py) draws an information-dense subset of a huge dataset and extends the sketch's analysis back to every cell; BPCells-style lazy on-disk matrices (LazyMatrix, truecell/lazy.py) keep a matrix out-of-core and stream over it. Both add no new dependency — sketching reuses the mapping.py / transfer.py machinery, and LazyMatrix is built on NumPy's memory-mapping alone.

BPCells-style lazy/on-disk matrices — ✅ delivered

  • Implemented as LazyMatrix in truecell/lazy.py, with write_lazy_matrix(matrix, path) to persist a matrix out-of-core and open_lazy_matrix(path) to map it back (is_lazy(x) to test). A matrix is stored as a directory of three memory-mapped .npy arrays — the data / indices / indptr triple of a compressed-sparse-column matrix, exactly scipy's csc_matrix layout — plus a JSON header. Opening maps the arrays without reading them; peak memory is the slices you touch, not the whole matrix (tests/test_lazy.py).
  • A faithful drop-in for a sparse layer. Indexing (m[:, cells], m[np.ix_(genes, cells)], slices, boolean masks) reads only the touched columns off disk and returns an ordinary scipy.sparse.csc_matrix, so every hot path that already accepts a sparse layer keeps working. Assay5 layers hold any array-like, so assay.set_layer_data("counts", lazy) just works, and the full-densify helpers (as_dense / np.asarray) route through __array__ — the deliberate "materialise everything" escape hatch you avoid on the huge path.
  • Streaming reductions. sum / mean over axis 0 (per-cell) or 1 (per-feature) run in a single pass over the on-disk arrays; nnz_per_col() (the nFeature count) comes straight from indptr; col_blocks(block_size) walks the matrix in cell-blocks — the primitive for processing a million cells at bounded RAM.
  • R: BPCells package enables out-of-core analysis on millions of cells.
  • Design note (CSC, not CSR): truecell matrices are features × cells and the scale-out ops (sketching, cell subsetting, per-cell normalisation) select columns, which CSC makes cheap. Row-only (feature) subsetting still scans the touched columns; keeping both orientations on disk (as BPCells does) to make that cheap too, and a full .toarray() audit of the hot paths to gate densification behind shape checks, are the natural follow-ons.

SketchData (leverage-score sub-sampling) — ✅ delivered

  • Implemented as sketch_data(obj, ncells=5000, method="LeverageScore") in truecell/sketch.py, with the leverage computation exposed on its own as leverage_score(obj). Each cell is sampled without replacement with probability proportional to its statistical leverage, so the rare states a uniform sample would drop are kept (indeed over-represented, which the tests check directly). Returns a standalone subset Truecell object whose active assay is renamed to "sketch"; the scores are written back onto the source object's leverage.score metadata (tests/test_sketch.py).
  • The leverage computation is the point, not a detail. The exact leverage of cell i is ℓ_i = ‖U_i‖² (the row norm of the left singular vectors), an O(n·d²) SVD — the very cost sketching exists to avoid. So (as Seurat's LeverageScore does) a sparse CountSketch S embeds the n rows into a small nsketch × d matrix B with BᵀB ≈ AᵀA in a single pass over the non-zeros; whitening A with B's right singular vectors gives an approximately orthonormal Z and ℓ_i ≈ ‖Z_i‖². When nsketch ≥ n no sketch is taken and the scores are exact — validated against the textbook ‖U_i‖² to 1e-6.
  • One caveat, documented in the docstring and a test: a single-hash CountSketch is a faithful subspace embedding only once it has ~d² rows, so nsketch should comfortably exceed the feature count (the default 5000 does, for the few-thousand variable features a real sketch runs on).
  • R: SketchData(obj, ncells = 5000, method = "LeverageScore")
  • Divergence from R: Seurat stores the sketch as an extra assay on the same object; here it is a separate object (as the plan called for), which fits truecell's subset model and keeps the full/sketch data cleanly apart.

ProjectData — ✅ delivered

  • Implemented as project_data(full, sketch, ...) in truecell/sketch.py, the inverse of sketch_data. Projects every full-dataset cell through the sketch's PCA loadings (stored as full.reductions["pca.full"]) — the same transform-only linear map project_umap uses — and, when the sketch carries a fitted UMAP, through that model too (full.reductions["ref.umap"], via project_umap). Optional refdata ({new_col: sketch_col} or a column name) carries the sketch's labels to the full data via find_transfer_anchors + transfer_data. Verified that projected sketch cells reproduce their sketch-PCA coordinates and that transferred labels recover the full data's cell types at >85% (tests/test_sketch.py).
  • R: ProjectData(obj, reference = sketch, reduction = "pca")
  • Reuse: no new machinery — the PCA/UMAP projection is mapping.py's and the label transfer is transfer.py's; project_data is the composition.

v0.9.0 — Specialized Assay Methods

HTODemux / MULTIseqDemux (cell hashing)

  • R: HTODemux(obj), MULTIseqDemux(obj)
  • hto_demux — ✅ delivered (truecell/hto.py): clustering into k = n_hashtags + 1 groups (kfunc="clara" by default, as in Seurat, or "kmeans" — see below) on CLR-normalised HTO counts, then a per-hashtag negative-binomial fit (maximum likelihood) to the tag's lowest-expressing cluster — its background — thresholded at positive_quantile (0.99) to call positive cells. Cells positive for zero / one / many hashtags are classified singlet / doublet / negative. Writes the Seurat columns HTO_maxID, HTO_secondID, HTO_margin, HTO_classification, HTO_classification.global plus a convenient hash.ID (also set as the active identity); learned cutoffs are stashed in obj.misc["hto_demux"]. normalize=False reuses a prior normalize_data(method="CLR") layer. Built on the existing CLR helper, with clara in-tree and sklearn k-means optional — no new dependency.
  • multiseq_demux — ✅ delivered (truecell/multiseq.py): the MULTI-seq chemistry (McGinnis et al.). Reuses hto_demux's CLR extraction, then thresholds each barcode straight off its distribution shape rather than a background fit — a Gaussian KDE over a 100-point grid exposes the background and positive modes as the two tallest local maxima and the cutoff sits a fraction quantile (0.7) between them. Cells clearing zero / one / many barcodes are singlet / doublet / negative. autothresh=True runs McGinnis's iterative deMULTIplex sweep — pick the q maximising the singlet rate, peel off negatives, re-threshold the remainder up to maxiter rounds. Writes MULTI_ID (also the active identity) and MULTI_classification; thresholds stashed in obj.misc["multiseq_demux"]. KDE via scipy.stats.gaussian_kde — no new dependency.
  • kfunc="clara" — ✅ delivered (truecell/_clara.py): Seurat's default clustering for HTODemux, a port of the clara C routine in R's cluster package (2.1.8.2). CLARA is k-medoids for data too big for PAM: it draws nsamples (100, Seurat's default) sub-samples of min(n, 40 + 2k) cells, runs PAM on each, assigns every cell to the resulting medoids, and keeps the sub-sample with the lowest total dissimilarity. Ported rather than taken from a library because the details that decide the answer are all non-standard: clara draws from its own 16-bit LCG (rngR = FALSE), so set.seed cannot reach it and clara correctly takes no seed argument; at pamLike = FALSE the swap rule is the pre-2011 one the C source itself calls "a bit illogical", not pam()'s; ties break last-wins in BUILD but first-wins in SWAP; and cluster numbering follows first appearance. kfunc="clara" is the default, matching Seurat; "kmeans" remains available. hto_demux first shipped defaulting to "kmeans", so this is a behaviour change for callers who never passed kfunc — one for the CHANGELOG when v0.10.0 adds it. The switch moves ~1% of calls on synthetic panels (rising with tag count — ~3.5% at 12 tags, where clara is also the more accurate of the two) because only each tag's least-expressing cluster feeds the background fit. Both scale linearly in cells; clara costs a roughly constant 4× (~1.3 s vs ~0.3 s at 100k cells). No new dependency; clara needs no sklearn.

Note: clara's answer depends on the CPU it was compiled for

clara accepts a swap on any improvement below zero — R genuinely takes swaps worth -2.2e-16 — so a one-ulp difference in a single distance flips a swap, then the winning sub-sample, then the whole clustering. That ulp is a compiler's choice: clara.c built for arm64 contracts clk += d*d into a fused multiply-add (one rounding), while the same source built for baseline x86_64, whose ISA has no FMA, rounds twice. Compiling the real clara.c for both targets and running them against each other confirms it — same source, same input, materially different clusterings on ~2% of random inputs and ~7% of realistic hashtag panels. R's clara is therefore not reproducible across architectures, and "bit-exact to R" is not a well-defined target.

truecell follows plain IEEE double arithmetic, which is what numpy gives on every platform and what clara.c gives on x86_64. Against that reference the port is exact — 200/200 pathological and 40/40 realistic HTO cases — and the only inputs where it differs from an arm64 R are exactly the ones where clara.c disagrees with itself across architectures. Emulating FMA to chase the arm64 answer would simply break fidelity the other way; there is no choice that satisfies both. Three details in _clara.py exist solely to hold that reference — np.cumsum rather than np.sum in _selec, the per-j accumulation of dz in _bswap2, and dz's h-major layout — and provably have no observable effect on any architecture-stable input (no fixture among 1490 can distinguish them). They are not redundant: without them the port drifts from the IEEE reference in the chaotic regime. Don't "simplify" them.

Mixscape (pooled CRISPR screen analysis) — ✅ delivered

  • R: CalcPerturbSig(obj, ...) → RunMixscape(obj, labels, nt.class.name) → MixscapeLDA(obj, labels, nt.label); plots PlotPerturbScore(obj, target.gene.ident) and MixscapeHeatmap(obj, ident.1, ident.2)
  • calc_perturb_sig — ✅ delivered (truecell/mixscape.py): Seurat's CalcPerturbSig. For every cell, the mean expression of its num_neighbors (20) nearest non-targeting (NT) control cells — in the first ndims of a reduction (default pca), optionally within each split_by batch — is subtracted from its own expression. The residual local perturbation signature (shared technical variation cancelled) is stored as a new assay (default "PRTB"). Neighbours via sklearn.neighbors.NearestNeighbors — no new dep.
  • run_mixscape — ✅ delivered (truecell/mixscape.py): Seurat's RunMixscape. Per target gene: (1) DE vs NT on de_assay (reuses find_markers) picks the perturbation-response genes; a gene with fewer than min_de_genes (5) is all NP. (2) On the signature over those genes, a perturbation vector (mean KO − mean NT) projects every cell to a single perturbation score, and an iterative 2-component sklearn.mixture.GaussianMixture splits the knockout (KO) mode from the non-perturbed (NP) mode — the vector is recomputed from each new KO set and re-fit until it stabilises (up to iter_num rounds). Writes mixscape_class ("<gene> KO" / "<gene> NP" / "NT", also the active identity), mixscape_class.global, and mixscape_class_p_<type>; per-gene bookkeeping stashed in obj.misc["mixscape"].
  • Two documented departures from R: the mixture is fit on the pooled NT + gene-cell scores (so the NP mode is anchored by the controls), and the signature is read from the assay's data layer directly — Seurat's pre-ScaleData centring is a provable no-op for the KO/NP calls (it leaves the perturbation vector unchanged and only globally shifts the scores the re-fit mixture absorbs).
  • mixscape_lda — ✅ delivered (truecell/mixscape.py): Seurat's MixscapeLDA (its PrepLDA + RunLDA). The complementary question to run_mixscape — not which cells are perturbed but how do the guide populations differ from each other and from control — answered with one supervised map. Per guide: DE vs NT picks its response genes (a guide needs npcs + 1, default 10, to contribute); on the signature assay over those genes a PCA is fit on that guide's cells plus the NT cells, and every cell is projected onto that guide's npcs-dim subspace. The per-guide blocks are concatenated and a single sklearn.discriminant_analysis.LinearDiscriminantAnalysis is fit with the guide label as the class. Stores an lda reduction (key LDA_, n_classes − 1 dims, loadings + assignments/posterior/genes_used in misc) and the lda_assignments / LDAP_<class> metadata columns.
  • Prerequisite is only calc_perturb_sig: faithful to R, the LDA groups cells by their raw guide label, not mixscape_class — the KO/NP calls are not consumed and no cells are dropped, so NP escapers stay in their guide's group (and, being control-like, generally land on NT).
  • Two more documented departures from R: Seurat's ScaleData → RunPCA → ProjectCellEmbeddings chain is composed directly (centre/scale each response gene against the guide-plus-NT reference, project through the reference PCA loadings — the same map, scale_max=10 clip included); and MASS's leave-one-out CV posterior (lda(..., CV = TRUE), stashed in misc by R and read by nothing) is not computed, only the resubstitution assignment and posterior.
  • plot_perturb_score — ✅ delivered (truecell/plotting.py): Seurat's PlotPerturbScore. The diagnostic for why mixscape split a guide the way it did — the perturbation score is the one axis the mixture is fit on, and the plot overlays the NT control density against the guide's own. A guide with a real effect is bimodal: one lobe on the NT curve (the escapers), one shifted away (the knockouts). before_mixscape=False (default) colours by mixscape_class (NT / "<gene> NP" / "<gene> KO", i.e. where mixscape drew the line); before_mixscape=True colours by the raw guide label, the pre-mixscape view. Cells are also drawn as a jittered strip (controls above the axis, target gene below), and split_by facets a multi-cell-type screen. Densities via scipy.stats.gaussian_kde (R uses geom_density) — no new dep.
  • The perturbation score is now persisted (truecell/mixscape.py): run_mixscape previously computed the score and discarded it. R keeps it in the Tool(object, "RunMixscape") slot precisely so PlotPerturbScore can read it back, so it is now stored per gene under obj.misc["mixscape"][assay]["genes"] as a scores frame (a pvec column plus the guide-label column, indexed by cell — R's gv data.frame). Two details are faithful to R's source rather than to the obvious reading: the score is kept from the first iteration (R's if (n.iter == 0)), before any KO/NP split feeds back into the perturbation vector — a later round would shift the axis; and it is normalised by vec · vec (R's ProjectVec divides by v2 %*% v2), which the classifier never needed (a positive constant cannot move the split) but a plotted axis does. Genes that short-circuit to NP without a mixture fit store scores = None.
  • mixscape_heatmap — ✅ delivered (truecell/plotting.py): Seurat's MixscapeHeatmap. The genes underneath the score: DE genes between two mixscape_class levels (e.g. "NT" vs "IFNGR2 KO"), with every cell ordered left-to-right by its knockout posterior, so the expression block is seen turning on in step with the probability. balanced=True takes up to max_genes from each fold-change direction; max_cells_group downsamples per class; order_by_prob=False shuffles instead (as R does). Delegates to do_heatmap exactly as R delegates to DoHeatmap.
  • do_heatmap gained R's cells argument to support that ordering (restrict to those cells and plot them in that exact order). Two latent bugs surfaced and were fixed with it: a scale.data layer scaled over a feature subset (scale_data(obj, features = [...])) has fewer rows than the assay has features, so indexing rows by the assay's full feature list read the wrong genes — layer-aware resolution now goes through _resolve_layer; and the group-label pass assumed each group was one contiguous block, which an explicit cells order interleaves, so it now walks runs.
  • Departure from R: the R colour names in the signatures (col = "orange2", and the fixed grey49 / grey79 for controls and escapers) are kept so the call reads like the R one, and translated to hex for matplotlib.
  • Still open: DEenrichRPlot (mixscape's third plot) is not ported — it calls the enrichR web service for pathway enrichment, which would put a live network dependency in the test path.

v0.10.0 — Package Infrastructure

PyPI publication — ✅ delivered

  • pip install truecell works: published as truecell, currently 0.9.0 (build + twine added to [dev] extras; published to TestPyPI then PyPI; verified with a clean-venv install + import + mini-pipeline smoke test)
  • Still open: replace the hard-coded __version__ string — ✅ delivered: truecell/__init__.py reads importlib.metadata.version("truecell"), with a PEP 440-valid 0.0.0+unknown fallback for a source tree that has no installed distribution. pyproject.toml is now the only place the version is written. Tradeoff worth knowing: that lookup resolves through the install-time dist-info snapshot, so an editable install keeps reporting the old version after a bump until it is reinstalled — a way to be wrong that the hard-coded string did not have. tests/test_packaging.py::test_version_matches_pyproject exists to make it loud rather than silent.
  • Still open — decide what the >= dependency floors mean. pyproject declares only lower bounds (pandas>=2.0, numpy>=1.24, anndata>=0.10, umap-learn>=0.5, scikit-learn>=1.3, …), so a fresh install resolves to whatever shipped this week and no two installs need agree. This has now caused four distinct failures, none of them theoretical:
  • tutorial figures drifting against the committed R panels;
  • mypy aborting on anndata 0.13.2's PEP 695 syntax and silently checking nothing (fixed by dropping python_version);
  • the mypy baseline spanning 81–85 across the CI matrix, because each leg resolves different versions;

  • pandas 3 crashing pbmc3k_tutorial.py on every fresh install while a developer venv pinned at pandas 2 stayed green (#36). Note the asymmetry that makes this expensive: the maintainer never sees it. Existing venvs hold old versions; only new users get the break. Pinning upper bounds is not obviously right either — it locks users out of new releases and ages badly. The decision to make is which of {upper bounds, a tested lockfile for CI, a scheduled "resolve latest" CI leg} this project wants; the status quo is "find out from a user".

  • Cut a release — ✅ delivered. The tags had stopped at 0.2.0 (2026-07-05) while milestones v0.3.0–v0.9.0 had all landed on main, so pip install truecell shipped almost none of the README's feature list: of 22 advertised entry points sampled, 19 were absent from the published wheel (reference mapping, sketching, LazyMatrix, hashing, Mixscape, run_spca/glm_pca, pseudobulk DE, MERSCOPE/Visium). 0.9.0 (2026-07-27) closes that gap in one jump — every milestone through v0.9.0 is now on PyPI, which is also why the version number lines up with the milestone it completes (a coincidence of this one release; see the note atop CHANGELOG.md).

GitHub Actions CI — ✅ delivered

  • File: .github/workflows/ci.yml
  • Matrix: Python 3.12, 3.13 on ubuntu-latest, via astral-sh/setup-uv (was 3.10–3.12; moved to track SPEC 0, which had already retired 3.10 in Oct 2024 and 3.11 in Oct 2025).
  • Open — add the 3.14 leg when harmonypy ships a cp314 wheel. 3.14 was in this change and came out again: it is the only package in the set without one (manylinux cp39–cp313 only), and both ways around it are worse than waiting. Building from source needs BLAS plus a CMake-fetched armadillo that ubuntu-latest lacks — that is the failure we hit. Forcing a wheels-only resolve instead backtracks to harmonypy 0.2.0, which depends on torch and pulls triton and 24 nvidia-* packages. Everything else resolves clean: 95 packages, wheels only, on x86_64-manylinux_2_28/3.14. Recheck with uv pip compile pyproject.toml --extra all --python-platform x86_64-manylinux_2_28 --python-version 3.14 --only-binary :all:.
  • Verify wheels on the target platform, not the dev machine. The 3.14 leg was signed off locally on macOS arm64, where the harmonypy source build succeeds because Accelerate supplies BLAS — so a clean local install proved nothing about Linux, and CI found it in 23 seconds. --python-platform + --only-binary :all: is the check that would have caught it; the hand-picked PyPI wheel survey that preceded it missed harmonypy entirely.
  • Open — the matrix no longer varies dependencies. Dropping 3.10 removed something nobody had written down: because the floors are >=, that leg was resolving numpy 2.2.6 / scipy 1.15.3 / pandas 2.3.3 where the others got 2.4.6 / 1.18.0 / 3.0.3, so it was the only coverage of an older scientific stack. All three legs now resolve one identical set — three interpreters, one dependency version. This surfaced through test_leverage_scores_are_not_flat, whose docstring had attributed its 3.10 failure to CPython when the real variable was the resolved numpy. The fix is a dedicated oldest-supported-deps leg (install against the declared floors, not the latest), not keeping an EOL Python around as an accidental proxy for it.
  • Triggers: push to main, all PRs
  • test job: ruff check truecell and mypy, both advisory (pre-existing debt not yet cleared, so neither gates the build), then pytest tests/ -q with --cov=truecell (80%); the coverage XML is uploaded as an artifact from the 3.12 leg rather than sent to a third-party service
  • Still open: a dedicated build job and coverage reporting — ✅ delivered: the build job runs uv build (sdist + wheel), twine check on both, then installs the wheel into a clean venv and imports it from outside the source tree, asserting __version__ matches pyproject.toml. That last step covers two things the test job structurally cannot: the test job installs -e .[all], so it never proves the wheel is complete or that a plain pip install truecell works without the optional scientific stack; and __version__ now comes from installed metadata, which only a real non-editable install exercises.

Type annotations

  • Add mypy to the CI lint job — ✅ delivered: advisory, on the same footing as ruff. [tool.mypy] in pyproject.toml pins files = ["truecell"] so a bare mypy checks exactly what CI checks (ruff takes its scope from the command line, where ruff check truecell = 75 and ruff check . = 163 are both correct and easy to confuse), and ignore_missing_imports silences the stub-less scientific stack — without it the run is 222 errors, 139 of them purely "this third party has no stubs".
  • Baseline to work down: ~80 errors in 15 modules on default settings, ~540 in 47 under --strict. Neither is a fixed target: both drift with the interpreter and with the dependency versions each resolves against the >= floors in pyproject.toml (the CI matrix spans 81–85 on default settings), so read them as a scale rather than a score. The package already ships py.typed, so these annotations are what a downstream mypy trusts.
  • Start here — 21 of the 81 are [name-defined], and they are real rather than pedantic: string annotations naming symbols that are never in module scope (-> "plt.Figure" across 17 plotting functions where plt is only imported inside function bodies; graph.py/neighbor.py annotating each other's classes to dodge a circular import). Runtime is unaffected — string annotations are never evaluated, which is why 444 tests pass — but typing.get_type_hints() raises NameError on them. Fixed with if TYPE_CHECKING: imports, keeping the lazy runtime imports exactly as they are (they keep matplotlib optional and the import graph acyclic). This was recorded as blocking the documentation site, and it was not. get_type_hints() still raises on all 19 of them, because a TYPE_CHECKING import does not execute — that is its purpose. mkdocstrings never calls it: griffe reads annotations statically, without importing the module, so the signatures render and cross-link precisely because those imports are there. The site was built on top of exactly this state; see below.
  • Add from __future__ import annotations to all modules (already done on some)
  • Annotate all public function signatures (mypy --strict clean)

Tutorial coverage — the R-fidelity net for everything after PR #10

  • Why this matters most. The two real defects ever found in this port — the CLR margin inversion (#32) and the SCTransform model (#37) — were both caught by a tutorial with an R side-by-side, not by the test suite, which was green through both. Meanwhile 24 feature PRs landed after #10 (integration, reference mapping, sketching, lazy matrices, HTO/MULTI-seq, Mixscape, spatial) and none of them has ever been compared to real Seurat — their tests assert self-consistency on synthetic default_rng fixtures. At the start of the initiative 36 of 103 public exports appeared in a runnable tutorial; 67 did not. As of T-de it is 82 of 104, counted as exports named anywhere in tutorials/*.py — the runnable scripts, not the prose. That method is stated because the earlier figures in this bullet were not reproducible from it (the T-obj entry said 81 where this count gives 74), so treat the series as indicative and the method as the thing to keep. Closing the rest is the highest-leverage correctness work left.
  • Plan: ~12 new side-by-side tutorials in four waves, one tutorial per PR, in the existing shape (<name>_tutorial.py + <name>_verify.R + <name>.md + figures_<name>/). Wave 1 = integration (ifnb), cell hashing (GSE108313), reference mapping (panc8), Mixscape (GSE153056). Wave 2 = cell-cycle / module scores (THP-1), dim-reduction extras (pbmc3k), leverage sketching (ifnb), object internals (pbmc3k) and spatial statistics / the spatial container (Xenium). Wave 3 = the DE-test suite (pbmc3k) and out-of-core LazyMatrix (pbmc3k). All four waves are now delivered — 16 tutorials, 29 defects, every one fixed. BPCells was installed for the final one (GitHub, not CRAN; needs libhdf5) after fingerprinting ten Seurat references to prove it was inert for ordinary input. T-vis (Visium) was added afterwards as a seventeenth, on the same method but with the reference question reopened, and T-int (anchor internals) as an eighteenth — bringing the totals to 18 tutorials and 48 defects, every one fixed — see Beyond the four waves below. The visualization gallery never became its own tutorial: dot_plot was the last plotting export uncovered and was folded into the pbmc3k gallery in T-sp.
  • Wave 0 — ✅ delivered (#38). The data plumbing every side-by-side needs: truecell.datasets loaders for the raw-source datasets, tutorials/export_seuratdata.R for the two SeuratData-only ones (ifnb/panc8, verified to round-trip R's counts exactly), and the R deps (SeuratData + harmony).
  • Wave 1 — ✅ complete.
  • T7 cell hashing — ✅ delivered (#39). hto_demux / multiseq_demux on GSE108313 (hashing_vignette.md). Result: HTODemux is 99.81 % call-concordant with Seurat on identical input — the first real-data confirmation of the CLR fix (#32) and clara default (#34). MULTIseqDemux lands at 94.67 %; the gap is a real KDE-implementation difference (scipy gaussian_kde vs R density() — bandwidth and grid), logged not papered over. No defect found — the demuxers hold up.
  • T9 Mixscape — ✅ delivered (#40). calc_perturb_sig / run_mixscape / mixscape_lda on GSE153056 (mixscape_vignette.md). On a shared variable-feature basis, per-cell class concordance is 97.45 % (KO/NP/NT and the full <gene> KO/NP class) — all NT cells agree, the same 14 guides read zero-effect on both sides, strong IFN-γ hits ≥97 %. Divergence is isolated to the weak boundary guides (MYC/SPI1/BRD4/CUL3) where the EM mixture is init-sensitive — a real method-level residual on a far more stochastic pipeline than the demuxers, not a bug. No defect found.
  • T6 integration — ✅ delivered (#41). First defects of the initiative. run_harmony / integrate_layers on ifnb (integration_vignette.md). Harmony and CCA reproduce Seurat's batch mixing and cell-type recovery to three decimals (batch-mixing entropy py/R 0.991 & 0.990/0.991). But investigating RPCA surfaced two real bugs, both now fixed: (1) a crash on unequal batch sizes — the reciprocal-PCA MNN args were mis-ordered, IndexError whenever n_query > n_ref, masked by the balanced synthetic test fixtures — fixed in #41; (2) a 4× under-integration (batch-mix 0.222 vs 0.914) — fixed in its own follow-up PR by matching Seurat's reciprocal-PCA construction: per-object scaling, the ReciprocalProject SD/L2 embedding normalization, and disabling the RPCA anchor filter, lifting batch-mix to 0.867 (py ARI→celltype 0.444 → 0.677). The old "quality, not count" read was wrong — it was mostly count (global scaling under-found anchors) plus the missing embedding normalization.
  • T8 reference mapping — ✅ delivered (#43). Wave 1's last tutorial. find_transfer_anchors / transfer_data / map_query / project_umap on panc8 (refmap_vignette.md). Transferring celltype from a CEL-seq2 reference to a SMART-seq2 query on a shared HVG basis: per-cell label concordance with Seurat is 98.71 % (2,363 of 2,394 cells), and each tool is ~98.5 % accurate against the held-out truth (truecell 0.9845, Seurat 0.9879). Every abundant type is ≥98 %; the whole error budget is the rare types (<10 reference cells) where both tools stumble the same way — a small single-tech reference's honest limit, not a divergence. No defect found — the transfer stack ports faithfully on its first real-data benchmark, the initiative's other valid outcome.
  • Wave 2 — ✅ complete. Five tutorials, sixteen defects, against Wave 1's four tutorials and two. The wave's lesson is that the quietest bugs are not in the algorithms: T-obj audited the container and found eleven, T-sp audited a statistic nobody had checked and found three more, and in both cases the broken code was returning plausible output that no test and no plot flagged.
  • T-cc cell-cycle & module scoring — ✅ delivered (#44). Wave 2's first. add_module_score / cell_cycle_scoring on THP-1 (cellcycle_vignette.md), a proliferating line with real S/G2/M populations (unlike resting PBMCs). On identical counts + shared resolved gene lists, per-cell Phase concordance is 96.62 % (20,028/20,729 cells) and the S/G2M/module scores correlate at Pearson ≥ 0.998. Both functions sample control genes at random and NumPy's RNG is not R's, so the scores are not bit-identical by construction — the residual is purely that RNG (the discrete Phase is robust to it), the same "don't chase the RNG" story as clara (hashing) and the MULTI-seq KDE. No defect found.
  • T-dr dim-reduction extras — ✅ delivered (#45). Two more defects. jack_straw / score_jackstraw / run_ica / run_tsne on pbmc3k (dimreduc_vignette.md). ICA reproduces R's subspace (Hungarian-matched |r| 0.982) and t-SNE preserves its input's neighbourhoods as well as R's does (30-NN retention 0.470 vs 0.477) — both clean. JackStraw was not. Two independent bugs, both fixed here: (1) the permutation null was built by projecting the scrambled rows onto the fixed embedding instead of re-running the PCA as R's JackRandom does, making the null far too tight — on the pure-noise PCs 14-20 that put 109-203 of 2000 features below p ≤ 1e-5 where R finds 0-5; (2) score_jackstraw aggregated with a one-sided KS test rather than R's prop.test, so its largest score across all 20 PCs was 8.1e-112 and no PC ever failed. Together they made truecell keep all 20 PCs where Seurat keeps 13. After the fix both keep 13, and the residual is permutation scatter (13-15 across seeds; R fixes its per-replicate seeds and is deterministic). JackStrawData.fake_reduction_scores, declared but never populated, is filled too.
  • T-sk leverage-score sketching — ✅ delivered. Two more defects. leverage_score / sketch_data / project_data on ifnb (sketch_vignette.md), exercising both of Seurat's regimes by moving nsketch rather than the dataset. After the fixes the exact regime is a per-cell match (Spearman 1.000000, max abs diff 3.4e-6 — below R's own unseeded-irlba noise), leverage tracks cell-type rarity at Spearman −0.929 in both tools, and project_data reaches Seurat's accuracy exactly (0.9050 each). Two bugs, both fixed here: (1) leverage was whitened against the full rank rather than Seurat's rank-50 truncation, so on 2000 HVGs the scores were crushed to a max/median of 1.34 against R's 6.48 — uniform sampling scores 1.00, so leverage sampling had become an expensive way to sample uniformly; (2) project_data transferred labels through the integration anchors, where ProjectData uses a weighted k-NN vote in the projected reduction — the anchor route scored better (0.936 vs 0.905), which is why it survived, but it costs exactly what sketching exists to remove and is unusable at the scale the API targets. Fixing it moved the headline number down and per-cell agreement up to 98.1 %.
  • T-obj object internals — ✅ delivered. Eleven more defects, the largest haul of the initiative. Cells/Features, the layered assay, Key, Embeddings/Loadings/Stdev, Graphs, FetchData, Idents/WhichCells/ RenameIdents/subset, and the command log, on pbmc3k (objects_vignette.md). The container rather than an algorithm, which makes it the sharpest net in the series: nothing here is stochastic, so 89 of the 91 anchors are compared with no tolerance at all — orders, names, dimensions, keys and non-zero counts either match or they do not. Coverage went from 36/103 exports at the start of the initiative to 81/104 as counted at the time — see the note on counting method in the plan bullet above. Eleven bugs, all fixed here. Five in the layered assay: split/JoinLayers was not a round trip, returning a layer named joined whose columns were in split order while the assay's own cell vector never moved — the right numbers in the wrong columns, silently misaligned against the metadata that indexes them, with every shape and checksum intact. The no-argument join_layers() additionally raised ValueError on any prepared assay, and generics.split_layers was declared but never registered for any type. Three in FetchData: np.asarray on a sparse matrix wraps it rather than densifying, so every one of 2,700 rows held a copy of the whole matrix — on the most-called accessor in Seurat, on the default assay class, guarded by a test that asserted only the column name and row count. Plus PC_1 was unaddressable and an unqualified fetch read counts where R reads data. Three in the bookkeeping: log_truecell_command had zero call sites so obj.commands was always empty against Seurat's five entries; no orig.ident; and add_meta_data rejected the plain vector R documents. join_layers/split_layers had zero call sites and zero tests before this — the defining feature of the v5 object model, never once run. Deferred here, closed in PR #55: find_neighbors symmetrized the kNN graph where Seurat's is directed (nnz 75,740 against exactly 2,700 × 20 = 54,000; degree 20-83, mean 28.1), and _build_snn dropped the self-edge Seurat keeps (~2,700 of the 4,044-edge gap). Both were fixed with the rest of the graph work, and the tutorial now matches 91 of 91 anchors — once its R script pins nn.method = "rann", because Seurat's default annoy is approximate and was costing 182 SNN edges on its own.
  • T-sp spatial statistics & the spatial container — ✅ delivered. Wave 2's last, three more defects. load_xenium / create_fov / create_centroids / create_segmentation / find_spatially_variable_features on the Xenium mouse brain (svf_vignette.md). 38 of 39 anchors match Seurat exactly. The existing Xenium tutorial built its R side with Read10X plus a coordinate data.frame, so R never constructed an FOV and the entire boundary layer had gone uncompared; FindSpatiallyVariableFeatures had never been run against R at all. Three bugs, all fixed here. (1) Moran's I was computed on a k-nearest-neighbour graph where Seurat row-standardises 1/d² over every pair — a good approximation, which is why it survived (Pearson 0.986, 46 of R's top 50), but a median 1.23× high and agreeing on only 7 of R's top 10, which is the part of the output anyone reads. Now 1.6e-14 and 10/10, evaluated in row blocks so R's O(n²) weighting runs on the full 36,602-cell slide in 5.3 s at 0.95 GB where RunMoransI needs a 10.7 GB dense allocation; weights="knn" keeps the old path as a documented approximation. (2) Centroids never received a radius — .AutoRadius gives 42.83 here — and _spatial_panel read it off the FOV, where R keeps NULL, so every true-to-scale spot renderer silently fell back to a fixed-size scatter on every non-Visium FOV. (3) Segmentation stored polygons open where R closes each ring. Left standing on purpose: the Moran's I p-value (R's 999-permutation test returns 14 distinct values and ties 233 of 248 genes at its floor — matching it would cost information), and the FOV Key, which appears only in __repr__.
  • Also fixed here, found in passing: the pbmc3k figure generator's hardcoded cluster→cell-type map had 1↔2 and 3↔4 transposed, so the monocyte compartment carried T-cell names in every labelled figure — including 11_umap_labeled.png, the tutorial's headline image, published since 93f33a2. The R block printed beside it in pbmc3k_tutorial.md had the correct order all along. Verified against markers (cluster 1 is LYZ 5.07/CD14 1.51, cluster 3 is MS4A1 2.14), corrected in both files, figures regenerated, and guarded by test_cell_type_map_matches_the_markers, which asserts each label leads on its own discriminative marker rather than trusting the map. That guard did not hold, and the reason is worth recording: the map broke again a release later — carrying nine entries for eight clusters, so the platelet cluster was captioned "DC" throughout 1.0.0 — and the test caught none of it, because test_tutorial_smoke.py is opt-in and excluded from CI. A correct assertion nobody runs is not a guard. It is now paired with test_cell_type_map_covers_exactly_the_clusters_produced, and the gate itself is closed for this dataset: the tutorials CI job caches PBMC 3k and runs nine of its eleven smoke tests on each pull request, treating a skip as a failure so it cannot go green having checked nothing. lazy_bpcells and pbmc3k_de are held out for runtime (65% of it) and are covered only by the full opt-in run. The other datasets (~200 MB) stay developer-run, so TRUECELL_TUTORIAL_SMOKE=1 pytest tests/test_tutorial_smoke.py before a release is still the rule. dot_plot — the last plotting export with no tutorial coverage — was folded into the same gallery, and drawing it is what exposed the mislabelling.
  • Wave 3 — ✅ complete.
  • T-de the differential-expression test suite — ✅ delivered. Two defects, one of them the most consequential of the initiative. All eight find_markers tests against FindMarkers on pbmc3k clusters 0 vs 1 (695 vs 477 cells, 13,714 genes, de_vignette.md), on a shared cell assignment so no clustering difference can masquerade as a DE difference. Eight tests and none had ever been compared to R. Result: all seven per-cell tests reproduce Seurat's top 50 genes exactly, with wilcox/t/bimod/LR at p-value Spearman ≥ 0.99997. (1) avg_log2FC put Seurat's pseudocount on the group mean instead of the group sum. Seurat 5's log1pdata.mean.fxn is log2((sum(expm1(x)) + 1) / n) — worth 1/n on the mean scale, not 1; truecell used Seurat 4's formula. This floored every fold change: a gene detected in 0 % of one cluster and 24 % of the other read −1.26 against Seurat's −9.92. It is not a display bug, because logfc_threshold filters on this value — at Seurat's own default of 0.1 truecell returned 4,903 genes where Seurat returned 13,009, and at 0.25 it was 2,298 against 11,931, Jaccard 0.193. The error concentrated in sparse marker-like genes (where both groups express a gene the two formulas agree at Spearman 0.990), which is exactly what DE is looking for. Now 7.11e-15 across all 13,712 genes. test_avg_log2fc_matches_seurat_formula already existed, encoded the same wrong formula, and was green throughout — corrected here. (2) negbinom ran a different test. Seurat's GLMDETest uses MASS::glm.nb (ML-estimated dispersion) and the Wald p-value; truecell used a fixed moment dispersion and a likelihood-ratio test, reading HLA-DRA at 5.5e-128 against R's 1.1e-321. After the fix the p-values agree exactly (median |log10 ratio| 0.000) for every gene detected above 5 %; the residual sits below that, where the GLM fits near-empty rows and Seurat's own min.cells.feature drops the genes (all 2,354 R dropped here were below 1 % detection in both groups). Left standing: deseq2 is pseudobulk where Seurat's DESeq2DETest tests cells as replicates — the practice Squair et al. 2021 showed inflates false positives — and since it requires sample_col it raises rather than silently substituting; mast is a hand-rolled hurdle model (Spearman 0.9993 on detected genes) because MAST has no Python equivalent; Seurat rounds myAUC to 3 dp inside DifferentialAUC, so ROC agrees only to 5e-4 by construction. Also corrected: the docstring advised passing CDR "to match Seurat's default CDR covariate", but MASTDETest fits ~ condition alone and adds no CDR term.
  • R-side deps: MAST 1.38.0 + DESeq2 1.52.0, and deliberately NOT glmGamPoi. It is a Suggests of DESeq2 rather than an Imports, FindMarkers never calls it, and its presence flips sctransform's vst onto a different backend — which would move the SCTransform R reference the SCT tutorial is pinned against. Installing with default dependencies leaves Suggests alone. Verified rather than assumed: an SCTransform fingerprint taken before and after the install is byte-identical.
  • T-lazy out-of-core — ✅ delivered. Seven defects, and the initiative's last tutorial. LazyMatrix against BPCells 0.3.1, the same pbmc3k matrix down both tools' on-disk and in-memory paths (lazy_vignette.md). 14 of 14 anchors match, 1998/2000 variable features shared with Seurat. The finding is each tool measured against itself: truecell's on-disk and in-memory paths are bit-identical, and Seurat's are not — BPCells computes in single precision, so its out-of-core run differs by 1.0e-06 on normalisation, 2.1e-02 on variance.standardized, and picks a different variable feature (1999/2000). truecell also runs all eight DE tests on a lazy layer where FindMarkers takes wilcox alone. The defects: (1) five functions densified the whole store — _log_normalize, _vst_hvg, scale_data, find_markers, add_module_score — so going on disk raised peak memory 4.6× and left a dense layer 16× larger; col_blocks had no callers at all. (2) percentage_feature_set crashed on a lazy layer. (3) create_assay5_object and (4) calc_n ended laziness in the constructor, making the obvious usage the one path that could not work — found only because a memory measurement came out backwards. (5) _loess2 was chaotically sort-dependent: 85.5 % of pbmc3k genes share a log10(mean), argsort is unstable, and windows were chosen by position, so a 1e-15 nudge moved fitted values 28.8 %; the fix moved toward R (HVG overlap 99.65 % → 99.90 %). (6) both HVG selectors broke ties opposite to R's order. (7) the sparse and lazy paths were two implementations agreeing to 1e-14, which a tie-break turned into 147 reordered features and 9 clusters vs 8; now one implementation. Left standing: truecell's store is uncompressed — 26.88 MB against BPCells' 4.50 MB bitpacked uint32, 6.0×. Storing counts as uint32 would take it to 17.92 MB; the remaining 4× is BP128 delta encoding, a format rather than a tolerance, and is not attempted here. LazyMatrix is also CSC-only, and Seurat warns that column-major is the wrong orientation for DE.
  • R-side dep: BPCells 0.3.1, and it is safe to install. Not on CRAN — remotes::install_github("bnprks/BPCells/r") — and it hard-requires libhdf5 (brew install hdf5 pkg-config); its configure script aborts without it. Unlike glmGamPoi, every reference to it inside Seurat is gated on inherits(x, "IterableMatrix") rather than requireNamespace, so a dgCMatrix run cannot reach a BPCells branch. Verified rather than assumed: ten Seurat references fingerprinted before and after — none moved.
  • Beyond the four waves — T-vis, and the reference question reopened.
  • T-vis Visium — ✅ delivered. One truecell defect, two Seurat ones. load_visium / VisiumV2 against Read10X_Image / Load10X_Spatial on the 10x V1_Mouse_Brain_Sagittal_Anterior Space Ranger 1.1.0 bundle, 2,695 in-tissue spots × 32,285 genes (visium_vignette.md). 24 of 24 anchors match, 17 exactly; coordinates to max|dx| = max|dy| = 0 over every spot; all four scale factors exact; the tissue image to 1.68e-08. This is the first tutorial where the difference was Seurat's. Every earlier one treated R as the reference and a mismatch as truecell's defect, which was right 29 times. Here it would have introduced one: Seurat builds the FOV with radius = scale.factors[["spot"]], and that field holds spot_diameter_fullres — a diameter in a slot named radius. Decided from the slide rather than from either tool: Visium spots sit on a fixed 100 µm centre-to-centre grid, the measured spacing is 137.000 px so 1.3700 px/µm, and read as a radius the field describes 130.62 µm capture spots on a 100 µm pitch — overlapping wells, which cannot be. Read as a diameter it is 65.31 µm, 10x's 65 µm reference spot to 1.0047×. truecell keeps diameter / 2 and a test pins the 2× divergence so a later "parity fix" cannot reintroduce it. Second Seurat finding: Radius() on a VisiumV2 returns NULL — methods("Radius") has Centroids, STARmap, SlideSeq, SpatialImage and VisiumV1, but no VisiumV2. The truecell defect: _imread fell back from matplotlib to Pillow, which return float32 in [0,1] and uint8 in [0,255] — 255× apart, from the same file, and neither library is a declared dependency, so get_image() was a function of the environment. Plotting hid it because imshow takes both. Three defaults aligned to Seurat (breaking): filter_by_tissue False → True, resolution hires → lowres, image key "spatial" → "slice1". And GetTissueCoordinates now returns Seurat's frame — x, y, cell with the cells on the index too, not x, y alone. Left standing: the LOESS residual reaches here as well — 1995/2000 shared variable features and a 1.9e-3 PCA gap. Isolated rather than tolerated: on Seurat's own feature list the decompositions agree to 2.49e-05, so the gap is feature selection, not the PCA.
  • T-int anchor internals — ✅ delivered. Eighteen defects, all fixed. find_integration_anchors / integrate_data against FindIntegrationAnchors / IntegrateData on a 2,400-cell ifnb subsample, 2,000 shared anchor features, both reductions (anchors_vignette.md). RPCA agrees on 100 % of Seurat's anchors (649/649, 99.5 % of scores identical to the last bit); CCA on 99.9 %, up from 70.0 % at 60.3 % precision. Corrected expression over the query half: mean|Δ| 0.0993 against Seurat's 0.0992. Why T6 missed all of it: T6 compares clusterings. A partition survives a great deal of damage before an adjusted Rand index notices, and the anchor set underneath it was barely two-thirds right the whole time. Anchors are also where a compiled reference has to be read rather than guessed — FindWeightsC and IntegrateDataC were pinned by calling them directly with controlled inputs (w = 1 − exp(−d̃·score/(2/sd)²), corrected = query − Wᵀ(query − ref), both reproducing to max|diff| = 0). The dominant defect: RunCCA standardizes each cell where truecell L2-normalized it — a correlation matrix between cells versus a cosine-similarity one, hence different singular vectors and a different anchor for every cell. The subtlest: _pca_loadings used sklearn's randomized SVD, fine in the leading components and drifting in the trailing ones (12–14 of 30 PCs matched irlba above 0.99). Reciprocal PCA standardizes each projected dimension by its own SD, which is not rotation-invariant, so a drifted trailing axis silently became a different reciprocal space — RPCA recall 44.9 % → 100 %. A metric got worse before it got better. Fixing the weight kernel dropped RPCA's batch mixing from 0.867 to 0.689 — the correct kernel had stopped masking the bad anchors, because a broad undiscriminating Gaussian smears batches together whether or not the anchors are right, and batch-mix entropy rewards exactly that. Restoring it to keep 0.867 would have meant reinstating a bug to flatter a metric; chasing the anchors instead is what found the SVD defect. Second time in this port — see project_data. The "residual" on the v5 IntegrateLayers path was two more bugs, not an implementation gap. RPCAIntegration/CCAIntegration don't call IntegrateData — they call IntegrateEmbeddings, which corrects the PCA embedding directly by transposing it into a fake assay; truecell's integrate_layers was running the v4 workflow (correct expression, re-scale, re-PCA) behind that name, agreeing with Seurat's actual output on 1 of 30 dimensions. Alongside it, run_pca itself carried the same randomized-SVD drift as defect 8 above, invisible while only leading PCs were read downstream but not once IntegrateEmbeddings corrects the embedding itself. Fixing both took the v5 embedding to 30/30 PCs above |r| = 0.99 (full 13,999-cell, unequal-batch ifnb) and RPCA batch mixing from 0.883 to 0.991 — above Seurat's own 0.917. And the last remnant was not a bug at all — the third framing of this same gap, and the one that finally holds. Clustering Seurat's own RPCA embedding through truecell reproduces truecell's numbers, not Seurat's, so the divergence is downstream of the embedding. Split into three stages it resolves cleanly: the KNN indices are identical, the SNN graphs agree off-diagonal to 2.8e-08, and only the community detection differs. Seurat's n.start = 10 restarts find a partition 0.17 % higher in modularity than truecell's single igraph pass (0.899903 vs 0.898336, and 20 truecell seeds never reach it) — but that extra modularity is spent splitting CD14 Mono into a 73.8 %-CTRL and an 83.3 %-STIM cluster, re-discovering the batch effect integration had just removed. truecell scores ARI 0.9195 to seurat_annotations against Seurat's 0.7368, and 85 % of 20 seeds beat Seurat. Searching less thoroughly is the advantage here, so the Louvain search was left alone and no n.start equivalent added. Four real defects did fall out of proving that: the nn graph was symmetrised (Seurat's is directed, nnz = n·k, column sums carrying the in-degree), the SNN diagonal was dropped (Seurat stores SNN[i,i] = 1 on all 13,999), run_umap didn't strip it the way RunUMAP.Graph does, and Jaccard ran in float32. All four were invisible to find_clusters, whose igraph conversion takes the strict upper triangle and discarded exactly the entries that were missing. GroupSingletons is ported alongside them. The guide tree (BuildSampleTree, three or more datasets) remains out of scope, as does a faithful port of RunModularityClusteringCpp.
  • Expect bugs, and read a mismatch as a bug report. Wave 1 went T7, T9 and T8 clean, while T6 found the first two defects, T-dr the next two, T-sk two more, T-obj eleven, T-sp three, T-de two, T-lazy seven, T-vis one (plus two in Seurat itself) and T-int eighteen — exactly the point: a green synthetic suite (balanced batches, self-consistent fixtures) hid a crash, a 4× under-integration, a mis-specified permutation null, the wrong significance test, a flattened sampling weight and a label transfer through the wrong machinery — all of which one real Seurat comparison exposed immediately. JackStraw and leverage sketching are the sharpest cases: JackStraw's tests asserted only that signal genes score lower than noise genes, which stayed true while the function was recommending every PC; leverage_score's asserted its own full-rank definition to six decimals, and its sampling fixture was too small to be in the algorithm's regime at all. The known-good tolerance is narrow — deterministic values match exactly, Louvain cluster counts drift ±1 — so anything outside that band gets investigated, not written up as an expected difference. This repo has twice let a real defect hide behind a documented "language difference" caveat.
  • A fix can make the headline number worse, and still be the fix. T-sk's project_data bug scored above Seurat before it was corrected. Fidelity to the reference implementation is the goal; when a divergence flatters us, that is a reason to look harder, not to keep it. Where a residual really is RNG, prove it distribution-against-distribution over matched seeds — single-run pairs were actively misleading on T-sk's sketch composition.

Documentation site — ✅ delivered

  • Tool: MkDocs + Material + mkdocstrings, config in mkdocs.yml, pages in docs/.
  • Structure: as planned, except that docs/tutorials is a symlink to the whole tutorials/ directory rather than one symlink per vignette — the vignettes reference their figures by relative path, so the figures have to come with them. docs/CHANGELOG.md and docs/ROADMAP.md are symlinks under their repo names, which keeps the relative links inside those two files working unchanged.
  • Deploy: .github/workflows/docs.yml — mkdocs build --strict on every pull request, deploy to GitHub Pages on push to main. Needs Pages set to "GitHub Actions" as its source in the repo settings; until then the build runs and the deploy step is the only part that fails.
  • The blocker this item led with was wrong, and worth recording as wrong. It said mkdocstrings resolves annotations and that typing.get_type_hints() raising NameError on the plotting and graph/neighbor signatures would stop the site on day one. get_type_hints() still raises, on 19 callables. mkdocstrings never calls it — griffe reads annotations statically, without importing the module — so those signatures render and cross-link because the if TYPE_CHECKING: imports are there. Clearing mypy (#69, #70) was worth doing on its own merits; it was not a prerequisite for this.
  • What the site found. Three Returns sections parsed as parameters, all 370 documented parameters putting their description in the type slot, 129 Sphinx roles rendering as literal text, fifteen Slots blocks collapsing into one line, and two dead heading anchors. Details in CHANGELOG.md; the parameter and role fixes live in tools/griffe_sphinx_roles.py as a build-time translation rather than as source churn.
  • Guarded by tests/test_docs.py: every public export reaches an API page, every ::: directive resolves, no symbol is documented twice, every nav entry and figure exists, and the site builds under --strict.

Changelog — ✅ delivered

  • File: CHANGELOG.md at repo root, in Keep a Changelog format
  • Populated retroactively for 0.1.0, 0.1.1 and 0.2.0 — but deliberately not straight from the git log as this item originally said: there is a chore: bump version to 0.1.2 commit, yet 0.1.2 was never tagged and never reached PyPI, so it is not a release and gets no entry. Taking the log at face value would have documented a version that never existed. (0.1.0 was tagged but never published; 0.1.1 was the first release on PyPI.)
  • Carried a standing note that the milestones on this page are not releases, because the two sequences had diverged far enough to mislead: the tags stopped at 0.2.0 while the milestones ran to v0.9.0, and a milestone can straddle releases (v0.7.0's spatial loaders shipped in 0.1.1; the rest of v0.7.0 stayed unreleased until 0.9.0). 0.9.0 (2026-07-27) closed that gap — everything through v0.9.0, including the three entries queued by earlier PRs (#32's BREAKING CLR margin fix, #33's clara cross-architecture caveat, and #34's hto_demux default change), moved from [Unreleased] into a dated [0.9.0] section. The standing note stays on the page: the two sequences will drift again the moment a milestone lands on main without a release to match it.

Dependency budget

Each milestone's new pip deps:

Milestone New deps
v0.2.0 harmonypy
v0.3.0 (none — uses sklearn already present)
v0.4.0 (none)
v0.5.0 (none — sPCA and GLM-PCA are pure NumPy/SciPy; glmpca-py proved unnecessary)
v0.6.0 pydeseq2 (optional)
v0.7.0 (none — Moran's I is pure NumPy/SciPy; the Visium PNG uses matplotlib, already in [analysis])
v0.8.0 (none — sketching reuses existing machinery; LazyMatrix is built on NumPy memory-mapping alone)
v0.9.0 (none — sklearn already present)
v0.10.0 build, twine, mkdocs, mkdocstrings, ruff, mypy (all dev-only)

Optional deps go in a new [spatial], [integration], or [all] extra in pyproject.toml so the base install stays lightweight.

Priority order

If milestones are too large, these are the highest-value individual items:

  1. Harmony (v0.2.0) — ✅ delivered (run_harmony / integrate_layers)
  2. WNN (v0.4.0) — ✅ delivered (find_multi_modal_neighbors + run_umap(graph=)); full two-stage port, CBMC tutorial section complete
  3. GitHub Actions CI (v0.10.0) — ✅ delivered
  4. FindTransferAnchors / TransferData / MapQuery / ProjectUMAP ✅ (v0.3.0) — truecell/transfer.py (find_transfer_anchors pcaproject/cca + transfer_data classification/imputation) and truecell/mapping.py (project_umap + map_query) deliver atlas-based annotation and place the query in the reference UMAP. Built on the v0.2.0 anchor machinery — v0.3.0 is complete.
  5. FindSpatiallyVariableFeatures + SpatialFeaturePlot ✅ (v0.7.0) — all four loaders, niche/neighbourhood analysis, both spatially-variable-feature methods (Moran's I and markvariogram), the VisiumV2 tissue-image data layer and the spatial_* H&E plots delivered; v0.7.0 is complete
  6. AggregateExpression + DESeq2 ✅ (v0.6.0) — aggregate_expression, find_conserved_markers, and pseudobulk DESeq2 (test_use="deseq2") delivered; MAST (test_use="mast") and bimod (test_use="bimod") too — v0.6.0 complete
  7. SketchData + BPCells-style lazy matrices (v0.8.0) — ✅ delivered. sketch_data / project_data + leverage_score (truecell/sketch.py): leverage-weighted subsampling for million-cell datasets, and projection of the sketch's PCA/UMAP/labels back to the full data. LazyMatrix (truecell/lazy.py): out-of-core, memory-mapped compressed-sparse-column storage with lazy column reads, streaming block reductions, and a clean Assay5-layer drop-in — no new dependency. v0.8.0 is complete.
  8. run_spca + glm_pca (Poisson + negative binomial) ✅ (v0.5.0) — v0.5.0 is complete; GLM-PCA now fits both family="poisson" and family="nb" (dispersion estimated by ML), closing the last gap in it