Stop losing downstream scores to NA, and let splatter/symsim fail properly

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
 

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

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

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)

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)