Keep simulated spots in the order they came in

CHANGELOG.md

@@ -12,6 +12,13 @@ Bug fixes:
     the counts they were attached to. Also make `cell_type_sel` and
     `cell_type_prop` agree in order, which they did not for 10 or more clusters.
     The simulated spots are returned in the order they came in.
+  - `splatter`, `sparsim`, `symsim` and `zinbwave`: return the simulated spots
+    in the order they were given, rather than sorted by `spatial_cluster`.
+    `clustering_ari` and `clustering_nmi` compare the real and simulated
+    clusterings element-wise, so those four were being scored against
+    mismatched spots.
+  - `generate_sim_spatialcluster`: check that the simulated dataset still lines
+    up with the real one, through a new `check_alignment()` helper.
   - `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.

src/helpers/utils.R

@@ -23,6 +23,43 @@ compute_logcounts <- function(adata) {
   Matrix::t(logcounts)
 }
 
+# A simulated dataset describes the same spots and genes as the dataset it was
+# simulated from, in the same order. The metrics compare the two element-wise --
+# ARI and NMI survive a permutation of the cluster labels, but not one of the
+# spots -- so a method that sorts its input has to sort the result back before
+# writing it out.
+check_alignment <- function(sim_adata, real_adata) {
+  checks <- list(
+    spots = list(sim = rownames(sim_adata$obs), real = rownames(real_adata$obs)),
+    genes = list(sim = rownames(sim_adata$var), real = rownames(real_adata$var))
+  )
+
+  for (what in names(checks)) {
+    sim_names <- checks[[what]]$sim
+    real_names <- checks[[what]]$real
+
+    if (identical(sim_names, real_names)) {
+      next
+    }
+
+    reason <- if (setequal(sim_names, real_names)) {
+      "they are in a different order"
+    } else {
+      paste0(
+        "they differ: ", length(setdiff(real_names, sim_names)), " missing, ",
+        length(setdiff(sim_names, real_names)), " unexpected"
+      )
+    }
+
+    stop(
+      "The simulated dataset's ", what, " do not match the real dataset's: ",
+      reason, "."
+    )
+  }
+
+  invisible(TRUE)
+}
+
 # spatial autocorrelation
 generate_moransI <- function(adata) {
   requireNamespace("spots", quietly = TRUE)

src/methods/sparsim/script.R

@@ -62,16 +62,18 @@ SPARSim_sim_param <- SPARSim::SPARSim_estimate_parameter_from_data(
 sim_result <- SPARSim::SPARSim_simulation(dataset_parameter = SPARSim_sim_param)
 
 # Reorder simulated results
-simulated_result_ordered <- sim_result$count_matrix[
-  match(rownames(sim_result$count_matrix), rownames(input_ordered$var)),
-  match(colnames(sim_result$count_matrix), rownames(input_ordered$obs))
-]
+simulated_result_ordered <- sim_result$count_matrix
+
+# put the spots back in the order they came in. The metrics compare the
+# simulated dataset against the real one spot for spot -- ARI and NMI are
+# invariant to a permutation of the labels, but not of the samples.
+restore_order <- order(ordered_indices)
 
 cat("Generating output\n")
 output <- anndataR::AnnData(
-  layers = list(counts = t(simulated_result_ordered)),
-  obs = input_ordered$obs[c("row", "col")],
-  var = input_ordered$var,
+  layers = list(counts = t(simulated_result_ordered)[restore_order, , drop = FALSE]),
+  obs = input$obs[c("row", "col")],
+  var = input$var,
   uns = c(
     input$uns,
     list(

src/methods/splatter/script.R

@@ -43,16 +43,18 @@ for (spatial_cluster in unique(input_ordered$obs[["spatial_cluster"]])) {
 colnames(simulated_result) <- rownames(input_ordered$obs)
 rownames(simulated_result) <- rownames(input_ordered$var)
 
-simulated_result_ordered <- counts(simulated_result)[
-  match(rownames(counts(simulated_result)), rownames(input_ordered$var)),
-  match(colnames(counts(simulated_result)), rownames(input_ordered$obs))
-]
+simulated_result_ordered <- counts(simulated_result)
+
+# put the spots back in the order they came in. The metrics compare the
+# simulated dataset against the real one spot for spot -- ARI and NMI are
+# invariant to a permutation of the labels, but not of the samples.
+restore_order <- order(ordered_indices)
 
 cat("Generating output\n")
 output <- anndataR::AnnData(
-  layers = list(counts = t(simulated_result_ordered)),
-  obs = input_ordered$obs[c("row", "col")],
-  var = input_ordered$var,
+  layers = list(counts = t(simulated_result_ordered)[restore_order, , drop = FALSE]),
+  obs = input$obs[c("row", "col")],
+  var = input$var,
   uns = c(
     input$uns,
     list(

src/methods/symsim/script.R

@@ -81,18 +81,20 @@ for (thisSpatialCluster in unique(input_ordered$obs[["spatial_cluster"]])) {
 colnames(simulated_result) <- rownames(input_ordered$obs)
 rownames(simulated_result) <- rownames(input_ordered$var)
 
-simulated_result_ordered <- counts(simulated_result)[
-  match(rownames(counts(simulated_result)), rownames(input_ordered$var)),
-  match(colnames(counts(simulated_result)), rownames(input_ordered$obs))
-]
+simulated_result_ordered <- counts(simulated_result)
+
+# put the spots back in the order they came in. The metrics compare the
+# simulated dataset against the real one spot for spot -- ARI and NMI are
+# invariant to a permutation of the labels, but not of the samples.
+restore_order <- order(ordered_indices)
 
 cat("Generating output\n")
 output <- anndataR::AnnData(
   layers = list(
-    counts = Matrix::t(simulated_result_ordered)
+    counts = Matrix::t(simulated_result_ordered)[restore_order, , drop = FALSE]
   ),
-  obs = input_ordered$obs[c("row", "col")],
-  var = input_ordered$var,
+  obs = input$obs[c("row", "col")],
+  var = input$var,
   uns = c(
     input$uns,
     list(

src/methods/zinbwave/script.R

@@ -39,18 +39,20 @@ simulated_result <- splatter::zinbSimulate(params)
 colnames(simulated_result) <- rownames(input_ordered$obs)
 rownames(simulated_result) <- rownames(input_ordered$var)
 
-simulated_result_ordered <- counts(simulated_result)[
-  match(rownames(counts(simulated_result)), rownames(input_ordered$var)),
-  match(colnames(counts(simulated_result)), rownames(input_ordered$obs))
-]
+simulated_result_ordered <- counts(simulated_result)
+
+# put the spots back in the order they came in. The metrics compare the
+# simulated dataset against the real one spot for spot -- ARI and NMI are
+# invariant to a permutation of the labels, but not of the samples.
+restore_order <- order(ordered_indices)
 
 cat("Generating output\n")
 output <- anndataR::AnnData(
   layers = list(
-    counts = Matrix::t(simulated_result_ordered)
+    counts = Matrix::t(simulated_result_ordered)[restore_order, , drop = FALSE]
   ),
-  obs = input_ordered$obs[c("row", "col")],
-  var = input_ordered$var,
+  obs = input$obs[c("row", "col")],
+  var = input$var,
   uns = c(
     input$uns,
     list(

src/process_datasets/generate_sim_spatialcluster/script.R

@@ -19,6 +19,10 @@ cat("Read input files\n")
 input_real_sp <- anndataR::read_h5ad(par$input_sp)
 input_simulated_sp <- anndataR::read_h5ad(par$input_sp_sim)
 
+# every method's output passes through here on its way to the metrics, so this
+# is the one place worth checking that it still lines up with the real dataset
+check_alignment(input_simulated_sp, input_real_sp)
+
 cat("add spatial cluster in simulated dataset:\n")
 sim_cluster <- generate_sim_spatialCluster(input_real_sp, input_simulated_sp)