From 05d0887c986b3a7bc229a940e4afa40079cbe853 Mon Sep 17 00:00:00 2001 From: Robrecht Cannoodt Date: Tue, 28 Jul 2026 13:57:30 +0200 Subject: [PATCH 1/4] stop losing downstream scores to NA * generate_cosine(): filter the Moran's I values rather than the objects holding them. `real[!is.na(real) & !is.na(sim)]` operates on a list, so it never dropped anything and one NaN took the whole metric with it -- 32 of 99 runs came out NA, while generate_mantel() passes na.rm and lost none. * calculate_precision(): report 0 when the simulation has no spatially variable genes at all. That is a bad simulation rather than a missing result, and reporting NA left both negative controls without a score on any dataset -- 8 of 30 expected control scores -- so nothing anchored the bottom of the scale. * calculate_recall(): left as NA on an empty denominator, and said why. There it is the real dataset that has no spatially variable genes, so there is nothing for any simulator to recover. --- src/helpers/utils.R | 25 +++++++++++++++++++++---- 1 file changed, 21 insertions(+), 4 deletions(-) diff --git a/src/helpers/utils.R b/src/helpers/utils.R index c4b46ea7..1063f51a 100644 --- a/src/helpers/utils.R +++ b/src/helpers/utils.R @@ -46,9 +46,20 @@ generate_moransI <- function(adata) { generate_cosine <- function(real, sim) { requireNamespace("lsa", quietly = TRUE) - real_new <- real[!is.na(real) & !is.na(sim)] - sim_new <- sim[!is.na(real) & !is.na(sim)] - similarity <- lsa::cosine(lsa::as.textmatrix(cbind(as.vector(real_new$Morans.I), as.vector(sim_new$Morans.I)))) + + real_i <- as.vector(real$Morans.I) + sim_i <- as.vector(sim$Morans.I) + + # drop the pairs where either side is missing. This used to filter `real` and + # `sim` themselves, which are lists, so it never dropped anything and a single + # NaN took the whole metric with it -- 32 of 99 runs came out NA, while + # generate_mantel() passes na.rm and lost none. + keep <- is.finite(real_i) & is.finite(sim_i) + if (sum(keep) < 2) { + return(NA_real_) + } + + similarity <- lsa::cosine(lsa::as.textmatrix(cbind(real_i[keep], sim_i[keep]))) return(mean(similarity)) } @@ -82,7 +93,10 @@ calculate_precision <- function(real_svg, sim_svg) { filtered_compared_data <- dplyr::filter(sim_svg$res_mtest, adjustedPval < 0.05) tp <- length(intersect(row.names(filtered_real_data), row.names(filtered_compared_data))) fp <- length(setdiff(row.names(filtered_compared_data), row.names(filtered_real_data))) - precision <- ifelse((tp + fp) > 0, tp / (tp + fp), NA) + # a simulation that produces no spatially variable genes at all is a bad + # simulation, not a missing result. Reporting NA meant both negative controls + # had no score on any dataset, so nothing anchored the bottom of the scale. + precision <- ifelse((tp + fp) > 0, tp / (tp + fp), 0) return(precision) } @@ -91,6 +105,9 @@ calculate_recall <- function(real_svg, sim_svg) { filtered_compared_data <- dplyr::filter(sim_svg$res_mtest, adjustedPval < 0.05) tp <- length(intersect(row.names(filtered_real_data), row.names(filtered_compared_data))) fn <- length(setdiff(row.names(filtered_real_data), row.names(filtered_compared_data))) + # unlike precision, an empty denominator here says something about the real + # dataset rather than about the simulation: there are no spatially variable + # genes to recover, so no simulator can be scored on it. NA is right. recall <- ifelse((tp + fn) > 0, tp / (tp + fn), NA) return(recall) } From 1595780db1d19032f4e8fbef9f75045dd35ad96b Mon Sep 17 00:00:00 2001 From: Robrecht Cannoodt Date: Tue, 28 Jul 2026 13:59:45 +0200 Subject: [PATCH 2/4] let splatter and symsim fail properly * Drop the `try()` wrapping the per-cluster loop body in both. A cluster that failed left `simulated_result` a cluster short and the loop carried on; the script only fell over later, on a length mismatch, and `res` was never looked at anyway. * symsim: sample gene lengths with replacement. `gene_len_pool` is finite, so `replace = FALSE` errors outright on any dataset with more genes than the pool holds. --- src/methods/splatter/script.R | 26 +++++----- src/methods/symsim/script.R | 90 +++++++++++++++++------------------ 2 files changed, 56 insertions(+), 60 deletions(-) diff --git a/src/methods/splatter/script.R b/src/methods/splatter/script.R index b63735f9..e1662936 100644 --- a/src/methods/splatter/script.R +++ b/src/methods/splatter/script.R @@ -25,21 +25,19 @@ input_ordered <- input[ordered_indices] simulated_result <- NULL for (spatial_cluster in unique(input_ordered$obs[["spatial_cluster"]])) { - res <- try({ - input_spatial_cluster <- input_ordered[input_ordered$obs[["spatial_cluster"]] == spatial_cluster] - params <- splatter::splatEstimate(as.matrix(t(input_spatial_cluster$layers[["counts"]]))) - sim_spatial_cluster <- splatter::splatSimulate(params) - sim_spatial_cluster$spatial_cluster <- spatial_cluster - colnames(sim_spatial_cluster) <- paste0(spatial_cluster, colnames(sim_spatial_cluster)) - names(rowData(sim_spatial_cluster)) <- paste(spatial_cluster, names(rowData(sim_spatial_cluster))) + input_spatial_cluster <- input_ordered[input_ordered$obs[["spatial_cluster"]] == spatial_cluster] + params <- splatter::splatEstimate(as.matrix(t(input_spatial_cluster$layers[["counts"]]))) + sim_spatial_cluster <- splatter::splatSimulate(params) + sim_spatial_cluster$spatial_cluster <- spatial_cluster + colnames(sim_spatial_cluster) <- paste0(spatial_cluster, colnames(sim_spatial_cluster)) + names(rowData(sim_spatial_cluster)) <- paste(spatial_cluster, names(rowData(sim_spatial_cluster))) - # combine the cell types - if (is.null(simulated_result)) { - simulated_result <- sim_spatial_cluster - } else { - simulated_result <- SingleCellExperiment::cbind(simulated_result, sim_spatial_cluster) - } - }) + # combine the cell types + if (is.null(simulated_result)) { + simulated_result <- sim_spatial_cluster + } else { + simulated_result <- SingleCellExperiment::cbind(simulated_result, sim_spatial_cluster) + } } colnames(simulated_result) <- rownames(input_ordered$obs) diff --git a/src/methods/symsim/script.R b/src/methods/symsim/script.R index e921b5a2..b6e67681 100644 --- a/src/methods/symsim/script.R +++ b/src/methods/symsim/script.R @@ -27,57 +27,55 @@ ordered_indices <- order(input$obs$spatial_cluster) input_ordered <- input[ordered_indices] for (thisSpatialCluster in unique(input_ordered$obs[["spatial_cluster"]])) { - res <- try({ - input_thiscelltype <- input_ordered[input_ordered$obs[["spatial_cluster"]] == thisSpatialCluster] - - # this is because if some genes are 0 , this will cause error in simulation - keep_feature <- colSums(input_thiscelltype$layers[["counts"]] > 0) > 0 - input_thiscelltype_f <- input_thiscelltype[, keep_feature] - - best_matches_UMI <- BestMatchParams( - tech = "UMI", - counts = as.matrix(t(input_thiscelltype_f$layers[["counts"]])), - plotfilename = "best_params.umi.qqplot", - n_optimal = 1 - ) + input_thiscelltype <- input_ordered[input_ordered$obs[["spatial_cluster"]] == thisSpatialCluster] - sim_thiscelltype <- SimulateTrueCounts( - ncells_total = dim(input_thiscelltype)[1], - ngenes = dim(input_thiscelltype)[2], - evf_type = "one.population", - randseed = 1, - Sigma = best_matches_UMI$Sigma[1], - gene_effects_sd = best_matches_UMI$gene_effects_sd[1], - scale_s = best_matches_UMI$scale_s[1], - gene_effect_prob = best_matches_UMI$gene_effect_prob[1], - prop_hge = best_matches_UMI$prop_hge[1], - mean_hge = best_matches_UMI$mean_hge[1] - ) + # this is because if some genes are 0 , this will cause error in simulation + keep_feature <- colSums(input_thiscelltype$layers[["counts"]] > 0) > 0 + input_thiscelltype_f <- input_thiscelltype[, keep_feature] - gene_len <- sample(gene_len_pool, dim(input_thiscelltype)[2], replace = FALSE) - sim_thiscelltype <- True2ObservedCounts( - true_counts = sim_thiscelltype[[1]], - meta_cell = sim_thiscelltype[[3]], - protocol = tech, - alpha_mean = best_matches_UMI$alpha_mean[1], - alpha_sd = best_matches_UMI$alpha_sd[1], - gene_len = gene_len, - depth_mean = best_matches_UMI$depth_mean[1], - depth_sd = best_matches_UMI$depth_sd[1] - ) + best_matches_UMI <- BestMatchParams( + tech = "UMI", + counts = as.matrix(t(input_thiscelltype_f$layers[["counts"]])), + plotfilename = "best_params.umi.qqplot", + n_optimal = 1 + ) + + sim_thiscelltype <- SimulateTrueCounts( + ncells_total = dim(input_thiscelltype)[1], + ngenes = dim(input_thiscelltype)[2], + evf_type = "one.population", + randseed = 1, + Sigma = best_matches_UMI$Sigma[1], + gene_effects_sd = best_matches_UMI$gene_effects_sd[1], + scale_s = best_matches_UMI$scale_s[1], + gene_effect_prob = best_matches_UMI$gene_effect_prob[1], + prop_hge = best_matches_UMI$prop_hge[1], + mean_hge = best_matches_UMI$mean_hge[1] + ) + + gene_len <- sample(gene_len_pool, dim(input_thiscelltype)[2], replace = TRUE) + sim_thiscelltype <- True2ObservedCounts( + true_counts = sim_thiscelltype[[1]], + meta_cell = sim_thiscelltype[[3]], + protocol = tech, + alpha_mean = best_matches_UMI$alpha_mean[1], + alpha_sd = best_matches_UMI$alpha_sd[1], + gene_len = gene_len, + depth_mean = best_matches_UMI$depth_mean[1], + depth_sd = best_matches_UMI$depth_sd[1] + ) - # tidy up the names - sim_thiscelltype <- SingleCellExperiment(list(counts = as.matrix(sim_thiscelltype$counts))) - sim_thiscelltype$spatial_cluster <- thisSpatialCluster + # tidy up the names + sim_thiscelltype <- SingleCellExperiment(list(counts = as.matrix(sim_thiscelltype$counts))) + sim_thiscelltype$spatial_cluster <- thisSpatialCluster - # combine the cell types - if (is.null(simulated_result)) { - simulated_result <- sim_thiscelltype - } else { - simulated_result <- SingleCellExperiment::cbind(simulated_result, sim_thiscelltype) - } + # combine the cell types + if (is.null(simulated_result)) { + simulated_result <- sim_thiscelltype + } else { + simulated_result <- SingleCellExperiment::cbind(simulated_result, sim_thiscelltype) + } - }) } colnames(simulated_result) <- rownames(input_ordered$obs) From 2983410c66a35313eb49c79a0ba5ee91e55d123d Mon Sep 17 00:00:00 2001 From: Robrecht Cannoodt Date: Tue, 28 Jul 2026 14:03:16 +0200 Subject: [PATCH 3/4] correct the logcounts type file_dataset_sp.yaml declared `logcounts` as integer. Log-transformed counts are doubles, and the shipped datasets store them as such. --- src/api/file_dataset_sp.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/api/file_dataset_sp.yaml b/src/api/file_dataset_sp.yaml index ffb6a077..c97c9811 100644 --- a/src/api/file_dataset_sp.yaml +++ b/src/api/file_dataset_sp.yaml @@ -13,7 +13,7 @@ info: name: counts description: Raw counts required: true - - type: integer + - type: double name: logcounts description: Log-transformed counts required: true From 1d0a6fd2b34230d3ddba4b2adbeffffcdf478f0c Mon Sep 17 00:00:00 2001 From: Robrecht Cannoodt Date: Tue, 28 Jul 2026 14:14:14 +0200 Subject: [PATCH 4/4] update changelog --- CHANGELOG.md | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index cff62be8..756c30a8 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -25,6 +25,18 @@ Bug fixes: `process_datasets/generate_sim_spatialcluster` instead of running BayesSpace a second time. The two runs disagreed, which put the positive control at `clustering_ari` 0.55 while the best simulator scored 0.22. + - `generate_cosine()`: filter the Moran's I values rather than the objects + holding them, so that a single NaN no longer takes the whole metric with + it. `crosscor_cosine` was NA for 32 of 99 runs. + - `calculate_precision()`: report 0 rather than NA when a simulation has no + spatially variable genes at all. Both negative controls were left without a + score on any dataset, so nothing anchored the bottom of the scale. + - `splatter` and `symsim`: drop the `try()` around the per-cluster loop, which + let a failed cluster pass silently and only surfaced later as a length + mismatch. + - `symsim`: sample gene lengths with replacement, which errored outright on + datasets holding more genes than `gene_len_pool`. + - `file_dataset_sp.yaml`: `logcounts` is a double, not an integer. # task_spatial_simulators 0.1.0