Deconvolving Visium spots with RCTD¶
Runs the real RCTD (spacexr, Cable et al. 2022) on a
simulated Visium-style section and scores its per-spot cell-type proportions against the composition
iscc actually placed in each spot.
R notebook, kernel R (iscc-rctd), no simulation. Generated once by
python validation/make_analysis_data.py --only rctd:
| file | what it is | who sees it |
|---|---|---|
sp_counts.csv, sp_coords.csv |
the spatial section: genes × spots, and spot coordinates | the tool |
ref_counts.csv, ref_labels.csv |
a paired scRNA reference and its cell-type labels | the tool |
truth_type_composition.csv |
true per-spot cell-type proportions | scoring only |
truth_clone_composition.csv |
true per-spot CNA-clone proportions | scoring only |
Why this benchmark is worth having. Deconvolution methods are normally validated against a
reference that came from the same tissue but a different assay, so the "truth" is itself an estimate.
Here the truth is exact: iscc knows which cells it put in each spot. The reference used below is the
oracle one (built from the section's own cells) — an upper bound, not a realistic reference.
What the spots resolve. The section is cut from the same realistic ductal field as every other
dataset here, and a deme on that field is ~50 µm — so a Visium v1 spot is about one deme, and the
mixture inside a spot is a within-deme mixture rather than sub-spot structure. Deconvolving below
the spot would need demes smaller than cells, which iscc deliberately does not do. The field is
grown well past the default target for this dataset: at the default density the section averages
under two cells per spot, and RCTD's task degenerates from deconvolution into classification.
source("r_preamble.R")
suppressMessages({ library(spacexr); library(Matrix) })
data_dir <- analysis_dir("rctd")
read_counts <- function(f) as(as.matrix(read.csv(file.path(data_dir, f), row.names = 1,
check.names = FALSE)), "dgCMatrix")
ref_counts <- read_counts("ref_counts.csv") # genes x reference cells
sp_counts <- read_counts("sp_counts.csv") # genes x spots
labels <- read.csv(file.path(data_dir, "ref_labels.csv"), row.names = 1, check.names = FALSE)
coords <- read.csv(file.path(data_dir, "sp_coords.csv"), row.names = 1, check.names = FALSE)
cat(sprintf("reference: %d genes x %d cells\nsection: %d genes x %d spots\n",
nrow(ref_counts), ncol(ref_counts), nrow(sp_counts), ncol(sp_counts)))
print(table(labels[, 1]))
reference: 6000 genes x 450 cells section: 6000 genes x 2542 spots
cancer epithelial immune stromal
163 8 94 185
Fit¶
RCTD's own pipeline: a Reference from the annotated scRNA, a SpatialRNA puck from the section,
then run.RCTD in full mode — which allows any mixture of types per spot, the right setting when a
spot genuinely contains several cells.
cell_types <- factor(labels[colnames(ref_counts), 1])
names(cell_types) <- colnames(ref_counts)
# RCTD needs enough cells per type to estimate its profile; drop types below the floor rather than
# letting it fail deep inside the fit.
min_per_type <- 25
keep <- names(which(table(cell_types) >= min_per_type))
cells <- names(cell_types)[cell_types %in% keep]
ref_counts <- ref_counts[, cells]; cell_types <- droplevels(cell_types[cells])
cat(sprintf("keeping %d types with >= %d cells: %s\n",
length(keep), min_per_type, paste(keep, collapse = ", ")))
reference <- Reference(ref_counts, cell_types, nUMI = colSums(ref_counts), min_UMI = 1)
coords_df <- data.frame(x = coords$x, y = coords$y); rownames(coords_df) <- rownames(coords)
puck <- SpatialRNA(coords_df, sp_counts, nUMI = colSums(sp_counts), require_int = TRUE)
myRCTD <- create.RCTD(puck, reference, max_cores = 1, CELL_MIN_INSTANCE = min_per_type,
UMI_min = 1, counts_MIN = 1, test_mode = FALSE)
myRCTD <- run.RCTD(myRCTD, doublet_mode = "full")
W <- as.matrix(myRCTD@results$weights)
W <- sweep(W, 1, pmax(rowSums(W), 1e-9), "/") # weights -> proportions
cat(sprintf("\nfitted %d spots x %d cell types\n", nrow(W), ncol(W)))
round(head(W, 4), 3)
keeping 3 types with >= 25 cells: cancer, immune, stromal
Begin: process_cell_type_info
process_cell_type_info: number of cells in reference: 442
process_cell_type_info: number of genes in reference: 6000
cancer immune stromal
163 94 185
End: process_cell_type_info
create.RCTD: getting regression differentially expressed genes:
get_de_genes: cancer found DE genes: 713
get_de_genes: immune found DE genes: 715
get_de_genes: stromal found DE genes: 782
get_de_genes: total DE genes: 2210
create.RCTD: getting platform effect normalization differentially expressed genes:
get_de_genes: cancer found DE genes: 848
get_de_genes: immune found DE genes: 835
get_de_genes: stromal found DE genes: 893
get_de_genes: total DE genes: 2518
fitBulk: decomposing bulk
chooseSigma: using initial Q_mat with sigma = 1
Likelihood value: 1329903.25282265
Sigma value: 1.16
Likelihood value: 1299420.70373146
Sigma value: 1.32
Likelihood value: 1272478.84306214
Sigma value: 1.48
Likelihood value: 1249108.1850487
Sigma value: 1.64
Likelihood value: 1229241.59476611
Sigma value: 1.8
Likelihood value: 1212702.50173826
Sigma value: 1.96
Likelihood value: 1199248.07116178
Sigma value: 2
Likelihood value: 1196330.30270048
Sigma value: 2
fitted 2542 spots x 3 cell types
| cancer | immune | stromal | |
|---|---|---|---|
| S0 | 0.254 | 0.045 | 0.701 |
| S1 | 0.872 | 0.042 | 0.086 |
| S2 | 0.779 | 0.133 | 0.088 |
| S3 | 0.749 | 0.073 | 0.178 |
Scoring against the true composition¶
Now open the truth. Two measures, because they answer different questions: per-spot correlation asks whether the shape of each spot's mixture is right, and mean absolute error asks how far off the proportions are in absolute terms. A method can score well on one and poorly on the other.
truth <- read.csv(file.path(data_dir, "truth_type_composition.csv"),
row.names = 1, check.names = FALSE)
common_spots <- intersect(rownames(W), rownames(truth))
common_types <- intersect(colnames(W), colnames(truth))
P <- W[common_spots, common_types, drop = FALSE]
T <- as.matrix(truth[common_spots, common_types, drop = FALSE])
cat(sprintf("scoring %d spots over %d shared types: %s\n",
length(common_spots), length(common_types), paste(common_types, collapse = ", ")))
per_type_r <- sapply(common_types, function(k)
if (sd(T[, k]) < 1e-9 || sd(P[, k]) < 1e-9) NA_real_ else cor(P[, k], T[, k]))
mae <- mean(abs(P - T))
# A flat "every type in equal proportion" guess is the baseline worth beating on MAE.
flat <- matrix(1 / ncol(T), nrow(T), ncol(T))
cat(sprintf("\nRCTD vs iscc ground truth (ORACLE reference — an upper bound)\n"))
cat(sprintf(" mean |error| %.3f (flat-guess baseline %.3f)\n", mae, mean(abs(flat - T))))
cat(" per-type correlation:\n")
print(round(per_type_r, 3))
cat(sprintf("\n mean per-type r: %.3f\n", mean(per_type_r, na.rm = TRUE)))
scoring 2542 spots over 3 shared types: cancer, immune, stromal
RCTD vs iscc ground truth (ORACLE reference — an upper bound)
mean |error| 0.058 (flat-guess baseline 0.277)
per-type correlation:
cancer immune stromal 0.979 0.967 0.977
mean per-type r: 0.974
The map, drawn by RCTD's own plotting¶
plot_puck_continuous is spacexr's spatial plotter — the function its vignette uses to show a value
across a puck. spacexr also ships plot_weights(), which draws exactly these panels for every
fitted type, but writes them to a PDF on disk; inline, that is the same native plotter in a loop.
The point is not the correlation coefficient, it is that the map is legible as tissue. The true per-spot composition to hold these against is in The analysis dataset and its ground truth.
suppressMessages(library(ggplot2))
options(repr.plot.width = 5.5, repr.plot.height = 4.6)
spots <- rownames(W)
for (k in colnames(W)) {
suppressWarnings(print(
plot_puck_continuous(puck, spots, W[, k], ylimit = c(0, 1),
title = sprintf("RCTD weight: %s", k), size = 1.4)))
}
