Integrating scDNA and scRNA with clonealign¶
Runs the real clonealign (Campbell et al. 2019)
on data simulated by iscc, and scores its clone assignments against ground truth the tool never
sees.
R notebook, kernel R (iscc-clonealign) — the code below is the R you would write yourself, not
a Python wrapper. It contains no simulation; the dataset was generated once by
python validation/make_analysis_data.py --only clonealign:
| file | what it is | who sees it |
|---|---|---|
Y.csv.gz |
scRNA counts, cells × genes | the tool |
L.csv.gz |
copy number, genes × clones, called by HMMcopy | the tool |
truth.csv |
each cell's true clone | scoring only |
What this shows. clonealign assigns single cells to copy-number clones by assuming a gene's
expression scales with the number of copies its clone carries. The question here is whether iscc's
simulated expression carries that dosage relationship faithfully enough for the real tool to exploit
it — i.e. whether simulated data is amenable to the analysis a real study would run.
Where the copy number comes from. Not from iscc. L holds the per-clone profiles
HMMcopy called from single-cell read depth on the same cells — the order a
real study runs things in. Handing clonealign the true copy number would answer half of its own
question for it.
That costs accuracy, and the notebook is worth reading with that in mind. On the true profiles the same fit reaches 0.81; on called profiles it does not, because this tumour is whole-genome doubled and ploidy is not identifiable from depth, so the two large clones — which truly differ in one segment out of twelve — come back with nearly the same called profile. Watch the per-clone AUCs: the two small clones, which differ substantially, are recovered almost perfectly.
The tumour. A realistic breach-gated ductal field (grid 96, 5 glands, 6,000 genes), WGD+, with four subclones.
# reticulate must be pointed at THIS env's Python, or it searches a cached uv interpreter and
# reports "Valid installation of TensorFlow not found" — clonealign's backend is TensorFlow via
# reticulate, so this preamble is load-bearing, not boilerplate.
local({
env_py <- file.path(dirname(dirname(R.home())), "bin", "python")
if (file.exists(env_py)) Sys.setenv(RETICULATE_PYTHON = env_py)
})
suppressWarnings(suppressMessages({
library(clonealign)
library(tensorflow)
library(jsonlite)
}))
data_dir <- file.path("..", "analysis_data", "clonealign")
stopifnot(dir.exists(data_dir))
Y <- as.matrix(read.csv(file.path(data_dir, "Y.csv.gz"), row.names = 1, check.names = FALSE))
L <- as.matrix(read.csv(file.path(data_dir, "L.csv.gz"), row.names = 1, check.names = FALSE))
meta <- fromJSON(file.path(data_dir, "meta.json"))
cat(sprintf("Y: %d cells x %d genes\nL: %d genes x %d clones\n",
nrow(Y), ncol(Y), nrow(L), ncol(L)))
cat(sprintf("clone sizes: %s\n", paste(meta$clone_sizes, collapse = ", ")))
cat(sprintf("baselines — chance %.2f, majority %.2f\n",
meta$chance_baseline, meta$majority_baseline))
cat(sprintf("copy number source: %s\n", meta$cn_source))
cat("\nclone copy-number profiles clonealign is given (rows = clones, cols = segments):\n")
print(meta$cn_called)
Y: 1501 cells x 6000 genes L: 6000 genes x 4 clones
clone sizes: 575, 730, 87, 109
baselines — chance 0.25, majority 0.49
clone copy-number consensus (rows = clones, cols = segments):
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [1,] 4 4 2 4 2 4 3 4 2 2 4 4 [2,] 4 4 2 4 2 4 4 4 2 2 4 4 [3,] 4 4 2 5 2 4 3 5 2 2 2 4 [4,] 2 2 1 2 1 3 2 2 2 1 1 2
What to expect¶
Most of the 6,000 genes sit on a segment whose called copy number varies between clones, so there is real dosage signal for the model to use — but less of it than the truth carries.
Two baselines matter when reading the score. Chance is 1/4 = 0.25. The majority baseline — always guessing the biggest clone — is 0.49 here, and that is the bar a useful method must clear. Accuracy alone can flatter or libel a method under class imbalance, so ARI and one-vs-rest AUC are reported beside it.
# Align genes across the two modalities, then fit.
common <- intersect(colnames(Y), rownames(L))
Yc <- Y[, common, drop = FALSE]; storage.mode(Yc) <- "double"
Lc <- L[common, , drop = FALSE]; storage.mode(Lc) <- "double"
set.seed(1); tf$compat$v1$set_random_seed(1L)
# clonealign's own multi-restart wrapper. The variational objective is sensitive to initialisation,
# so a single fit can settle in a poor local optimum; this is the package's intended entry point.
fit <- run_clonealign(Yc, Lc,
initial_shrinks = c(0, 5, 10),
n_repeats = 3,
max_iter = 200,
print_elbos = FALSE,
verbose = FALSE)
probs <- fit$ml_params$clone_probs
colnames(probs) <- colnames(Lc); rownames(probs) <- rownames(Yc)
cat(sprintf("fitted: %d cells assigned over %d clones\n", nrow(probs), ncol(probs)))
table(fit$clone)
Warning message in (function (gene_expression_data, copy_number_data, max_iter = 200, : “Less than 75% of genes positively correlated with expression - assignment may have failed ”
Warning message in (function (gene_expression_data, copy_number_data, max_iter = 200, : “Less than 75% of genes positively correlated with expression - assignment may have failed ”
Warning message in (function (gene_expression_data, copy_number_data, max_iter = 200, : “Less than 75% of genes positively correlated with expression - assignment may have failed ”
Warning message in (function (gene_expression_data, copy_number_data, max_iter = 200, : “Less than 75% of genes positively correlated with expression - assignment may have failed ”
Warning message in (function (gene_expression_data, copy_number_data, max_iter = 200, : “Less than 75% of genes positively correlated with expression - assignment may have failed ”
Warning message in (function (gene_expression_data, copy_number_data, max_iter = 200, : “Less than 75% of genes positively correlated with expression - assignment may have failed ”
Warning message in (function (gene_expression_data, copy_number_data, max_iter = 200, : “Less than 75% of genes positively correlated with expression - assignment may have failed ”
Warning message in (function (gene_expression_data, copy_number_data, max_iter = 200, : “Less than 75% of genes positively correlated with expression - assignment may have failed ”
fitted: 1501 cells assigned over 4 clones
clone0 clone1 clone2 clone3 unassigned
1134 111 102 151 3
clonealign's own figure¶
plot_clonealign is the package's own diagnostic — the same figure as in the paper. The lower panel is
copy number per clone along the genome (scDNA); the upper panel is the mean expression of the cells
assigned to each clone (scRNA). The model's entire assumption is that the two track each other: where
a clone carries extra copies, the cells assigned to it should express more.
iscc's genome is a run of contiguous segments rather than named chromosomes, and the figure draws one
contig at a time, so the whole genome is handed over as a single contig — every segment then appears in
one panel, in genome order.
One label to disambiguate: the lower panel's legend reads Ground truth clone, which is plot_clonealign's name for the scDNA copy-number profiles it was given as input — L, not iscc's held-back answer key. Nothing on this figure has seen truth.csv; that is opened in the next section.
suppressWarnings(suppressMessages({
library(SingleCellExperiment); library(cowplot); library(tidyr)
}))
options(repr.plot.width = 10, repr.plot.height = 6.5)
# 400 cells is plenty for the per-gene z-scores the figure averages, and keeps the gene x cell
# reshape inside it small.
set.seed(1); cells <- sort(sample(nrow(Yc), 400))
cts <- t(Yc[cells, , drop = FALSE]) # genes x cells, as SingleCellExperiment wants
sce <- SingleCellExperiment(assays = list(counts = cts))
sf <- colSums(cts); sf <- sf / mean(sf) # library-size normalise, then log
logcounts(sce) <- log2(sweep(cts, 2, sf, "/") + 1)
rowData(sce)$chr <- "1"
rowData(sce)$start_position <- seq_len(nrow(sce))
rowData(sce)$end_position <- seq_len(nrow(sce))
rowData(sce)$ensembl_gene_id <- rownames(cts)
suppressWarnings(suppressMessages(
print(plot_clonealign(sce, fit$clone[cells], Lc, chromosome = "1"))
))
Scoring against ground truth¶
Only now do we open truth.csv. Predicted clone labels are arbitrary, so we take the best
permutation of predicted-to-true labels (with three clones there are only six, so we enumerate them
rather than reaching for the Hungarian algorithm).
truth <- read.csv(file.path(data_dir, "truth.csv"), stringsAsFactors = FALSE)
true_clone <- truth$true_clone[match(rownames(probs), truth$cell)]
pred <- max.col(probs) - 1L # 0-based, to match the truth file
# best label permutation (3 clones -> 6 permutations, so enumerate rather than use Hungarian)
K <- ncol(probs)
perms <- as.matrix(expand.grid(rep(list(0:(K - 1)), K)))
perms <- perms[apply(perms, 1, function(r) length(unique(r)) == K), , drop = FALSE]
best <- which.max(apply(perms, 1, function(p) mean(p[pred + 1L] == true_clone)))
acc <- mean(perms[best, ][pred + 1L] == true_clone)
maj <- max(table(true_clone)) / length(true_clone)
# ACCURACY ALONE IS MISLEADING HERE: the clones are 1062/211/254, so "always guess the biggest"
# already scores ~0.70. Rank-based measures say whether the tool found signal at all, independent
# of that imbalance.
adj_rand <- function(a, b) { # adjusted Rand index, from the contingency table
tab <- table(a, b); n <- length(a)
ch2 <- function(x) sum(x * (x - 1) / 2)
idx <- ch2(tab); ea <- ch2(rowSums(tab)); eb <- ch2(colSums(tab)); tot <- n * (n - 1) / 2
(idx - ea * eb / tot) / (0.5 * (ea + eb) - ea * eb / tot)
}
auc1 <- function(score, pos) { # one-vs-rest AUC (Mann-Whitney)
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)
}
mapped_true <- match(true_clone, perms[best, ]) - 1L # true labels in predicted-column space
aucs <- sapply(seq_len(K), function(k) auc1(probs[, k], mapped_true == (k - 1L)))
cat("clonealign vs iscc ground truth\n")
cat(sprintf(" accuracy %.2f (chance %.2f, majority %.2f) -> %s the majority baseline\n",
acc, 1 / K, maj, ifelse(acc > maj + 0.02, "ABOVE", "AT OR BELOW")))
cat(sprintf(" ARI %.2f (0 = no better than random grouping)\n", adj_rand(pred, true_clone)))
cat(sprintf(" mean AUC %.2f per clone: %s\n", mean(aucs, na.rm = TRUE),
paste(sprintf("%.2f", aucs), collapse = ", ")))
# Make the verdict follow the numbers rather than hard-coding a conclusion: a hard-coded closing
# sentence silently goes stale the moment the dataset changes.
if (acc > maj + 0.02) {
cat(sprintf("\n Above the majority baseline on hard assignment, with mean AUC %.2f: iscc's simulated\n", mean(aucs, na.rm = TRUE)))
cat(" expression carries the copy-number dosage signal clonealign relies on.\n")
} else {
cat(sprintf("\n AUC %.2f above 0.5 with accuracy at the majority baseline: signal was found, but not\n", mean(aucs, na.rm = TRUE)))
cat(" enough to win hard assignments against the largest clone.\n")
}
clonealign vs iscc ground truth
accuracy 0.54 (chance 0.25, majority 0.49) -> ABOVE the majority baseline
ARI 0.26 (0 = no better than random grouping)
mean AUC 0.88 per clone: 0.74, 0.81, 0.98, 0.99
Above the majority baseline on hard assignment, with mean AUC 0.88: iscc's simulated expression carries the copy-number dosage signal clonealign relies on.
Reading the result¶
If the accuracy sits at the majority baseline, clonealign has not found the clone structure — and
on this tumour that is the expected outcome rather than a failure of the tool. Copy-number dosage is
almost uninformative when clones share copy number in ten of twelve segments.
The companion notebook Integrating scDNA and scRNA with clonealign shows what does separate these clones: the allele layer. Two clones can both sit at total copy number 4 while differing in their allelic composition (4+0 versus 2+2), and B-allele frequency sees a difference that total copy number cannot.
That contrast is the point. A benchmark on a substrate where every tool succeeds tells you nothing about which modality carries the signal.
The truth this is scored against — and the panels to compare these figures with — are in The analysis dataset and its ground truth.