Calling copy number from scRNA with Numbat¶
Runs the real Numbat (Gao et al. 2023) — allele-aware
CNA inference from single-cell RNA — on iscc data, and scores it against the copy number iscc
actually gave each cell.
R notebook, kernel R (iscc-numbat), no simulation. Generated once by
python validation/make_analysis_data.py --only numbat:
| file | what it is | who sees it |
|---|---|---|
count_mat.csv |
genes × cells expression counts | the tool |
ref_counts.csv |
genes × normal cells, the expression reference | the tool |
df_allele.csv |
phased allele counts per cell and SNP | the tool |
gtf.csv |
gene → segment annotation | the tool |
truth.csv |
per cell: malignant or not, and its clone | scoring only |
Why Numbat needs more than the others. It is allele-aware: it reads per-homolog signal, not just
total expression. That is why this dataset is the only one grown with allele-specific expression
params — without them the tumour has no cell_rna_baf layer and Numbat cannot be built at all.
Scope note. validation/numbat_runner.R wraps the call below in fallback handling for the case
where Numbat's LLR filter leaves no CNV at all, so that the benchmark still produces scoreable
neutral output. This notebook shows the pipeline itself; if Numbat recovers nothing here, the cell
below will say so rather than pretending to a result.
source("r_preamble.R")
suppressMessages({ library(numbat); library(dplyr); library(data.table); library(Matrix) })
data_dir <- analysis_dir("numbat")
count_mat <- as(as.matrix(read.csv(file.path(data_dir, "count_mat.csv"), row.names = 1,
check.names = FALSE)), "dgCMatrix") # genes x cells
ref_mat <- as(as.matrix(read.csv(file.path(data_dir, "ref_counts.csv"), row.names = 1,
check.names = FALSE)), "dgCMatrix") # genes x normal cells
df_allele <- as.data.frame(fread(file.path(data_dir, "df_allele.csv")))
df_allele$CHROM <- as.integer(df_allele$CHROM)
df_allele$cell <- as.character(df_allele$cell)
gtf <- as.data.frame(fread(file.path(data_dir, "gtf.csv")))
cat(sprintf("expression: %d genes x %d cells\nreference: %d genes x %d normal cells\n",
nrow(count_mat), ncol(count_mat), nrow(ref_mat), ncol(ref_mat)))
cat(sprintf("alleles: %d rows over %d SNPs\nannotation: %d genes over %d segments\n",
nrow(df_allele), length(unique(df_allele$snp_id)),
nrow(gtf), length(unique(gtf$CHROM))))
expression: 6000 genes x 1677 cells reference: 6000 genes x 150 normal cells
alleles: 1129161 rows over 4458 SNPs annotation: 6000 genes over 12 segments
Fit¶
aggregate_counts turns the normal cells into Numbat's expression reference (lambdas_ref), then
run_numbat does the joint expression + allele inference. max_iter = 1 keeps this to a single pass
— enough to recover the CNA structure, and the phylogeny refinement is not what is being scored.
annot <- data.frame(cell = colnames(ref_mat), group = "normal")
lambdas_ref <- aggregate_counts(ref_mat, annot)
out_dir <- file.path(tempdir(), "numbat_out"); dir.create(out_dir, showWarnings = FALSE)
ok <- tryCatch({
run_numbat(count_mat, lambdas_ref, df_allele, gtf = gtf, genome = "hg38", out_dir = out_dir,
min_cells = 20, ncores = 1, ncores_nni = 1, max_iter = 1,
t = 1e-5, min_LLR = 5.0, plot = FALSE, verbose = FALSE)
TRUE
}, error = function(e) { cat(sprintf("run_numbat error: %s\n", conditionMessage(e))); FALSE })
iters <- as.integer(gsub(".*joint_post_(\\d+)\\.tsv$", "\\1",
list.files(out_dir, pattern = "joint_post_\\d+\\.tsv$")))
cat(sprintf("\nrun completed: %s | iterations on disk: %s\n",
ok, ifelse(length(iters) == 0, "none", paste(iters, collapse = ","))))
cell_dict normal 150
Mem used: 0.704Gb
Mem used: 0.704Gb
number of genes left: 3292
running hclust...
Mem used: 0.64Gb
Expression noise level (MSE): high (4.9). Consider using a custom expression reference profile.
less than 5% of genome is in neutral region - including LOH in baseline
Testing for multi-allelic CNVs ..
0 multi-allelic CNVs found:
Mem used: 0.691Gb
All cells succeeded
No multi-allelic CNVs, skipping ..
No multi-allelic CNVs, skipping ..
No multi-allelic CNVs, skipping ..
Mem used: 0.741Gb
Mem used: 0.754Gb
Iter 2 -93.3595213189637, 0.015s
Iter 3 -93.3595213189637, 0.016s
Warning message: “There was 1 warning in `mutate()`. ℹ In argument: `GT = unlist(...)`. Caused by warning: ! The `father` argument of `bfs()` is deprecated as of igraph 2.2.0. ℹ Please use the `parent` argument instead. ℹ The deprecated feature was likely used in the tidygraph package. Please report the issue at <https://github.com/thomasp85/tidygraph/issues>.”
less than 5% of genome is in neutral region - including LOH in baseline
All done!
run completed: TRUE | iterations on disk: 1
Scoring¶
The headline question for a CNA-from-expression method is whether it separates malignant from
normal cells. Numbat reports a per-cell aneuploidy probability, which we score as an AUC against
iscc's ground truth — a threshold-free measure, so it does not depend on where a cut-off is placed.
truth <- read.csv(file.path(data_dir, "truth.csv"), stringsAsFactors = FALSE)
if (length(iters) == 0) {
cat("Numbat's LLR filter left no CNV, so there is no per-cell call to score.\n")
cat("That is a real outcome, not an error — it is what numbat_runner.R's neutral fallback covers.\n")
} else {
nb <- Numbat$new(out_dir, i = max(iters))
cp <- as.data.frame(nb$clone_post)
p_aneu <- if ("p_cnv" %in% names(cp)) cp$p_cnv else 1 - cp$p_1
is_mal <- truth$is_malignant[match(cp$cell, truth$cell)]
keep <- !is.na(is_mal) & !is.na(p_aneu)
auc1 <- function(score, pos) { # Mann-Whitney AUC, threshold-free
r <- rank(score); np <- sum(pos); nn <- sum(!pos)
if (np == 0 || nn == 0) return(NA_real_)
(sum(r[pos]) - np * (np + 1) / 2) / (np * nn)
}
auc <- auc1(p_aneu[keep], is_mal[keep] == 1)
cat(sprintf("Numbat vs iscc ground truth\n"))
cat(sprintf(" scored %d cells (%d malignant, %d normal)\n",
sum(keep), sum(is_mal[keep] == 1), sum(is_mal[keep] == 0)))
cat(sprintf(" malignant-vs-normal AUC %.3f (0.5 = no separation)\n", auc))
}
Numbat vs iscc ground truth scored 1677 cells (1527 malignant, 150 normal) malignant-vs-normal AUC 0.977 (0.5 = no separation)
Numbat's own view of the tumour¶
plot_phylo_heatmap is the figure from Numbat's vignette: the single-cell phylogeny beside the
per-cell copy-number posterior, so clones appear as blocks of shared CNA. This is what a user would
actually look at after running Numbat — and it is being drawn from iscc-simulated expression.
suppressMessages(library(ggplot2))
if (length(iters) > 0) {
options(repr.plot.width = 12, repr.plot.height = 7)
suppressWarnings(print(
nb$plot_phylo_heatmap(clone_bar = TRUE, p_min = 0.9)
))
} else cat("no CNV recovered, so there is nothing to draw\n")
And the consensus copy-number segments Numbat called across the tumour.
if (length(iters) > 0) {
options(repr.plot.width = 11, repr.plot.height = 3.2)
suppressWarnings(print(plot_consensus(nb$segs_consensus)))
}
Numbat's phylogeny-and-heatmap and its consensus segments are its own two headline figures, and both are above. Everything they should be compared against — the true malignant/normal labels and the copy number iscc actually gave each cell — is in
The analysis dataset and its ground truth.