| name | rctd-py |
| description | GPU-accelerated Python RCTD cell type deconvolution for spatial transcriptomics (Visium, VisiumHD, Xenium, MERFISH, Slide-seq). Use when running rctd-py, when a Python alternative to R spacexr is needed, when deconvolving spatial data on GPU, when preparing h5ad references from Seurat/scanpy/CELLxGENE, or when integrating rctd-py results back into a Seurat/AnnData workflow. Covers the CLI (`rctd run`, `rctd validate`, `rctd info`), the Python API (`run_rctd`, `Reference`, `RCTDConfig`), all three modes (doublet / multi / full), GPU dtype selection, per-arch eigh tuning, and result interpretation. Triggers on "rctd-py", "python rctd", "GPU RCTD", "spacexr python", "h5ad deconvolution", "spatial cell type deconvolution python", "convert spacexr reference to h5ad". |
rctd-py: GPU-Accelerated Spatial Deconvolution
Python reimplementation of RCTD (Cable et al., Nature Biotechnology 2022) with PyTorch GPU acceleration. ~4-41× faster than R spacexr with >99% cell-type-call concordance on identical pixel sets. Pure GPU implementation -- no R runtime required.
When to use
- Cell-type deconvolution on spatial transcriptomics (Visium, VisiumHD, Xenium, MERFISH, Slide-seq) where R
spacexr is too slow or impractical to integrate into a Python pipeline.
- Large datasets (>50k spots) on GPU.
- Workflows that read/write h5ad and want to keep the deconvolution in-process with scanpy/anndata.
- Validating R
spacexr results against a faster Python reference.
For the R spacexr workflow itself, use that package directly; this skill assumes you have already decided on the Python path.
Install
uv pip install torch
uv pip install torch --index-url https://download.pytorch.org/whl/cu126
uv pip install torch --index-url https://download.pytorch.org/whl/cu128
uv pip install rctd-py
rctd info
rctd info prints torch version + CUDA availability + GPU name. If cuda: False is reported on a machine with a GPU, the wrong torch wheel was installed -- reinstall from the CUDA-matched index URL.
Quick start: CLI
rctd validate spatial.h5ad reference.h5ad --cell-type-col cell_type
rctd run spatial.h5ad reference.h5ad \
--mode doublet \
--device cuda \
--umi-min 20 \
--output results.h5ad \
--json
--json prints a one-line summary on completion (mode, pixel count, runtime, output path). Default output is <input_stem>_rctd.h5ad.
Quick start: Python API
import anndata
from rctd import Reference, RCTDConfig, run_rctd
spatial = anndata.read_h5ad("spatial.h5ad")
ref = Reference(anndata.read_h5ad("reference.h5ad"), cell_type_col="cell_type")
config = RCTDConfig(device="cuda", UMI_min=20, dtype="float32")
result = run_rctd(spatial, ref, mode="doublet", config=config, batch_size=10000)
For results stored back into the spatial AnnData, use the CLI -- it does the merge for you. From Python, see "Result interpretation" below.
Mode selection
| Mode | Best for | What it returns |
|---|
doublet (default) | Xenium, MERFISH, Slide-seq, any single-cell-resolution data | 1-2 cell types per spot with confidence class |
multi | Dense spatial with several contributing types (Xenium with heavy contamination) | Up to MAX_MULTI_TYPES (default 4) types per spot via greedy forward selection |
full | Visium, VisiumHD, anything with multi-cell spots | Continuous weight vector over all K types (unconstrained LS) |
For single-cell-resolution platforms, prefer doublet. Reach for multi only when a domain-specific reason makes >2 types per spot biologically expected.
Reference preparation
Reference inputs are AnnData objects with:
.X = raw integer counts (cells × genes). Not normalized, not log-transformed..obs[cell_type_col] = string cell-type label.
Per-type subsampling to ~500 cells is recommended (more cells rarely improves the fit and slows reference loading).
From CZ CELLxGENE via gget
import anndata
import gget
gget.setup("cellxgene")
gget.cellxgene(
species="mus_musculus", tissue="cortex", disease="normal",
out="reference.h5ad",
)
ref = anndata.read_h5ad("reference.h5ad")
ref.var_names = ref.var["feature_name"].values
ref.var_names_make_unique()
ref.write_h5ad("reference.h5ad")
From Seurat (R → h5ad)
library(anndataR)
seurat_ref <- readRDS("reference.rds")
Idents(seurat_ref) <- "cell_type"
seurat_ref <- subset(seurat_ref, downsample = 500)
adata <- as_anndata(seurat_ref, x_layer = "counts")
adata$obs$cell_type <- as.character(seurat_ref$cell_type)
write_h5ad(adata, "reference.h5ad")
From a pre-built spacexr::Reference .rds
The .rds may be either a bare spacexr::Reference S4 object or a list containing one as $reference. Handle both:
library(spacexr)
library(anndata)
obj <- readRDS("spacexr_reference.rds")
ref <- if (is.list(obj) && !is.null(obj$reference)) obj$reference else obj
counts <- t(ref@counts)
obs <- data.frame(
cell_type = as.character(ref@cell_types),
nUMI = ref@nUMI,
row.names = names(ref@cell_types)
)
write_h5ad(AnnData(X = counts, obs = obs), "reference.h5ad")
From existing scanpy / h5ad
import anndata
adata = anndata.read_h5ad("reference.h5ad")
assert "cell_type" in adata.obs.columns
CLI flag reference (rctd run)
| Flag | Default | Notes |
|---|
--mode | doublet | full, doublet, multi |
--output / -o | <stem>_rctd.h5ad | Output h5ad path |
--device | auto | auto, cpu, cuda |
--batch-size | auto | Pixels per IRWLS launch; auto sizes from free VRAM |
--dtype | float64 | float32 is ~2× faster on consumer GPUs (see GPU note below) with identical cell-type calls |
--cell-type-col | cell_type | Reference obs column with labels |
--umi-min | 100 | Use 20 for Xenium; 100 is the R spacexr default |
--umi-max | 20000000 | Per-pixel UMI ceiling |
--sigma-override | (auto) | Skip sigma estimation, force a fixed value (use for exact R concordance testing) |
--gene-cutoff | 0.000125 | Bulk gene expression threshold |
--fc-cutoff | 0.5 | Bulk fold-change threshold |
--gene-cutoff-reg | 0.0002 | DE gene expression threshold |
--fc-cutoff-reg | 0.75 | DE fold-change threshold |
--max-multi-types | 4 | Max types per spot (multi mode) |
--confidence-threshold | 5.0 | Singlet confidence margin |
--doublet-threshold | 20.0 | Doublet certainty margin |
--cell-min | 25 | Min cells per reference type |
--n-max-cells | 10000 | Max cells per reference type (downsampling) |
--eigh-threshold | (arch default) | Force the K cutoff for the GPU-vs-CPU eigh dispatch (issue #22 diagnostic) |
--json | off | Machine-readable summary |
--quiet / -q | off | Suppress progress output |
RCTDConfig exposes the same fields under the Python API; consult _types.py for the canonical list.
GPU notes
dtype="float32" is usually the better default on GPU
The default dtype="float64" matches R spacexr arithmetic and is the safest choice when validating against R. On a data-center GPU with strong fp64 (A100, H100) the cost is small. On consumer / workstation GPUs (RTX, L40S, anything Blackwell sm_120, anything pre-Hopper) fp64 throughput is 1/32-1/64 of fp32, and dtype="float32" is roughly 2× faster on doublet mode with the same cell-type calls. Empirically (L40S, sm_89, N=20000 spots, K=30 types): 76 s → 39 s; spot_class agreement with the fp64 run is 100%, weights diff <1e-6. See tests/test_fp32_concordance.py for the locked-in assertion.
If you need bit-for-bit R spacexr concordance for validation, stay on fp64.
Per-architecture K cutoff for the eigh dispatch
_psd_batch decides between batched GPU eigh (cuSOLVER) and CPU offload (LAPACK syevd) based on K and the GPU's compute capability:
sm < 9 (Volta, Turing, Ampere, Ada / L40S): K ≤ 16 stays on GPU; K > 16 offloads to CPU.sm ≥ 9 (Hopper, Blackwell): K ≤ 128 stays on GPU.
This was tuned empirically and is correct for the configurations we've tested. If you hit a perf cliff or a CUSOLVER_STATUS_INVALID_VALUE crash, override it:
--eigh-threshold 0 (or RCTDConfig(eigh_threshold=0)): force CPU eigh on every architecture. Cheap and robust if a particular cuSOLVER batch shape is misbehaving.--eigh-threshold 64 etc.: force GPU eigh up to a higher K. Useful on V100/L20 if the CPU offload turns out to be slower than GPU eigh for your specific workload.
CPU eigh thread oversubscription (fixed in v0.3.5)
On hosts with many cores running on Volta / Turing / Ampere / Ada (i.e. forced into the CPU eigh branch at moderate K), versions ≤ 0.3.4 could stall for thousands of seconds because OpenBLAS spawned one thread per core under batched syevd. v0.3.5 wraps the CPU eigh call with torch.set_num_threads(1) and restores the caller's count on exit. No user action required from v0.3.5 onwards. If you're on an older release and see Step 1 of doublet mode stuck, upgrade or set OMP_NUM_THREADS=1.
VRAM guidance
| Available VRAM | Recommended --batch-size | Peak VRAM (K=45) |
|---|
| ≥ 24 GB | auto (often 10000+) | ~4 GB |
| 8-16 GB | 5000 | ~2 GB |
| < 8 GB | 2000 | ~1 GB |
--batch-size auto reads free VRAM at startup. On a shared GPU where another job's allocation can grow during the run, pin an explicit --batch-size to avoid mid-run OOM, and consider PYTORCH_ALLOC_CONF=expandable_segments:True.
Generic SLURM template
#!/bin/bash
set -euo pipefail
echo "Started: $(date) on $(hostname)"
nvidia-smi --query-gpu=name,memory.total --format=csv,noheader || true
rctd run spatial.h5ad reference.h5ad \
--cell-type-col cell_type \
--mode doublet \
--device cuda \
--dtype float32 \
--umi-min 20 \
--output results_rctd.h5ad \
--json
Result interpretation
The CLI writes results into a copy of the spatial AnnData:
| Slot | Content | Modes |
|---|
obs["rctd_dominant_type"] | Top cell type per spot | all |
obs["rctd_spot_class"] | singlet / doublet_certain / doublet_uncertain / reject / filtered | doublet |
obs["rctd_first_type"] | Primary type (string) | doublet |
obs["rctd_second_type"] | Secondary type (string) | doublet |
obs["rctd_converged"] | IRWLS convergence flag | full |
obs["rctd_n_types"] | Types assigned per spot | multi |
obsm["rctd_weights"] | (N × K) weight matrix | all |
obsm["rctd_weights_doublet"] | (N × 2) top-2 weights | doublet |
obsm["rctd_sub_weights"] | Selected-type weights | multi |
uns["rctd_cell_type_names"] | Ordered K labels for the weight matrix | all |
uns["rctd_mode"] | Mode used | all |
uns["rctd_config"] | Configuration dict | all |
uns["rctd_version"] | rctd-py version that produced the file | all |
Spot class is 0-indexed in the Python API (unlike R spacexr's 1-indexed spot_class):
| Code | Label | Meaning |
|---|
| 0 | reject | Confidence too low |
| 1 | singlet | Single dominant type |
| 2 | doublet_certain | Two types, high confidence |
| 3 | doublet_uncertain | Two types, lower confidence |
Filtered spots (below UMI_min) carry NaN weights and "filtered" labels.
import anndata
adata = anndata.read_h5ad("results_rctd.h5ad")
print(adata.obs["rctd_spot_class"].value_counts())
singlets = adata[adata.obs["rctd_spot_class"] == "singlet"]
print(singlets.obs["rctd_first_type"].value_counts().head(10))
print(adata.obsm["rctd_weights"].shape)
print(adata.uns["rctd_cell_type_names"])
Integration back to Seurat
library(anndataR)
library(Seurat)
result <- read_h5ad("results_rctd.h5ad")
obs <- as.data.frame(result$obs)
spatial_obj <- readRDS("spatial.rds")
common <- intersect(colnames(spatial_obj), rownames(obs))
spatial_obj$rctd_first_type <- obs[common, "rctd_first_type"]
spatial_obj$rctd_spot_class <- obs[common, "rctd_spot_class"]
spatial_obj$rctd_dominant_type <- obs[common, "rctd_dominant_type"]
W <- as.matrix(result$obsm[["rctd_weights"]])
colnames(W) <- result$uns[["rctd_cell_type_names"]]
for (ct in colnames(W)) {
safe <- paste0("rctd.weight.", make.names(ct))
spatial_obj[[safe]] <- W[match(colnames(spatial_obj), rownames(obs)), ct]
}
saveRDS(spatial_obj, "spatial_rctd.rds")
Best practices
- Validate first.
rctd validate checks gene overlap, count integrality, and label coverage in seconds without a GPU.
--umi-min 20 for Xenium. The 100 default is inherited from R spacexr and is set for Visium/Slide-seq UMI distributions.- Match reference to sample condition. If your cohort mixes diagnoses (e.g. malignancy A vs B), do not run every sample against one combined reference containing disease-specific malignant prototypes -- RCTD will misassign healthy/reactive cells of the wrong-diagnosis samples to those prototypes. Branch on the sample's metadata column and pick the reference per sample.
- Avoid
/ in cell-type names. zarr/h5 may convert it to _. Pick consistent ASCII labels.
- Downsample large references to ~500 cells/type; more rarely improves the fit.
- Prefer
dtype="float32" on consumer GPUs. Concordance is verified; speedup is ~2×. Stay fp64 only for exact R spacexr validation runs.
Gotchas
spot_class is 0-indexed in Python; R spacexr is 1-indexed. When cross-comparing, add 1 to the Python values (or subtract 1 from the R values).- MT-gene case sensitivity vs R
spacexr. The DE filter is case-insensitive (.lower().startswith("mt-")), which drops human MT-CO1-style symbols from the candidate pool. R spacexr's case-sensitive grep("mt-", ...) keeps them. Net effect on real data: rctd-py keeps fewer pixels because the differing gene_list_bulk shifts the per-pixel counts_MIN filter. Cell-type calls on the pixels both runs keep agree at ~100%.
- Concordance metrics that count filtered pixels as mismatches understate agreement. When the pixel-filtering decision differs between rctd-py and R
spacexr, a concordance = same_first_type / total_pixels will appear low even when the kept-pixel calls agree exactly. Compare assignments only on the pixel intersection.
- Doublet mode under-assigns fine subtypes on small panels. When a reference fine subtype has markers that the spatial panel can't reliably detect (e.g. CD4 markers on a 300-500-gene Xenium panel), doublet mode falls back to a coarser label or
doublet_uncertain rather than over-calling the subtype. Always value_counts() your subtype × group counts before slicing by subtype for downstream analysis -- sparsity at the call layer can look like a biological finding. To re-split a coarse pool into subtypes, use a reference classifier (SingleR, ProjecTILs) rather than a per-cell marker argmax (low capture biases argmax toward high-capture markers).
batch_size="auto" on a shared GPU can OOM mid-run. It sizes from free VRAM at startup; if another job grows afterwards, the IRWLS launch can fail. Pin --batch-size 5000-8000 on contested nodes.- Determinism. Python runs are bit-reproducible (
np.random.default_rng(42) throughout). R spacexr with max_cores > 1 is not (mclapply advances the parent's global RNG as a side effect, perturbing the later sample() in chooseSigma). Use --sigma-override for run-to-run determinism matching a fixed R baseline.
Multi-modal (RNA + protein)
Joint RNA + protein deconvolution (Gaussian/WLS protein block sharing the per-spot cell-type weights) is in development on the feat/multimodal-protein branch. The objective is L(w) = L_RNA(w) + λ · L_protein(w); protein enters the IRWLS gradient/Hessian and the doublet pair/singlet scores. Per-type protein profiles are bootstrapped from RNA-confident singlets, optionally combined with curated signed gates (positive + negative markers per type) to break RNA-degenerate types apart. Check the PR list on https://github.com/p-gueguen/rctd-py/pulls for current status before assuming the API is on main.
Gotchas observed on the branch (API: run_rctd_multimodal(...), RCTDConfig(protein_weight=..., protein_norm=...)):
- Never route continuous protein through count-fusion.
run_rctd_multimodal([RNA, protein], fusion="count") treats every modality as Poisson counts; fed continuous intensity (~1e6-1e7) the protein likelihood is flat and the joint collapses to a single constant weight vector across all pixels (corr ~0, dominant-type accuracy ~chance) -- confirmed on a ~62-plex CosMx dataset and a synthetic RNA-degenerate-pair test. Use fusion="protein_wls" (protein enters as a Gaussian/WLS block), never fusion="count", for any continuous protein modality.
- Raise the UMI ceiling for protein-intensity values. The RNA-count ceiling (
--umi-max default 2e7) silently filters out protein-scale per-cell sums (e.g. CosMx protein ~3e7, ~6e7 for a 2-cell pseudospot) -- nearly all pixels get dropped and scoring crashes on a shape mismatch. When feeding protein-scale values, set the ceiling to 1e9-1e12; RNA totals are hundreds-to-thousands, so a high ceiling leaves RNA membership unchanged.
protein_weight="auto" underweights a rich, informative panel. The per-feature heuristic balances gradient magnitude at init, not marker informativeness -- on a deep, trustworthy panel (60+ markers) it lands low (~0.1) and only ties RNA-only, while a fixed weight ~4 can be best-in-class. Sweep it up to 4+ even for sharpening (not just to override a wrong RNA call) when the panel is deep. protein_norm="arcsinh_robust" handles both IF intensity and protein counts; prefer it over per-cell CLR for IF (CLR couples the few channels via the per-cell geometric mean).class_df (doublet mode, also on this branch) relaxes the singlet/doublet/reject confidence threshold for within-class ambiguity -- it does NOT change which type wins. Grouping same-lineage subtypes (e.g. several fibroblast subtypes into one Fibroblast class) recovers more usable singlets with no accuracy downside, but it will NOT fix a wrong-type-assignment problem (that lives upstream in candidate scoring). Keep ungrouped any ambiguity you actually need to resolve.
Resources
