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 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 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) } 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)