diff --git a/CHANGELOG.md b/CHANGELOG.md
index 83fa29ba..b50eabbf 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,6 +1,12 @@
 # task_spatial_simulators dev
 
 Bug fixes:
+  - `scdesign2`: order the input by `spatial_cluster` before simulating.
+    `simulate_count_scDesign2()` returns cells grouped per cell type, so the
+    coordinates taken from the unsorted input belonged to different spots than
+    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.
   - `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.
diff --git a/src/methods/scdesign2/script.R b/src/methods/scdesign2/script.R
index aac648d2..ee061fd5 100644
--- a/src/methods/scdesign2/script.R
+++ b/src/methods/scdesign2/script.R
@@ -18,13 +18,24 @@ if (par$base != "domain") {
 }
 
 cat("scDesign2 simulation start\n")
-counts <- Matrix::t(input$layers[["counts"]])
-colnames(counts) <- as.character(input$obs$spatial_cluster)
-spatial_cluster_prop <- table(input$obs$spatial_cluster)
+# simulate_count_scDesign2() returns the cells grouped per cell type, so order
+# the spots by spatial cluster first -- otherwise the coordinates we attach
+# below belong to different spots than the counts we attach them to
+ordered_indices <- order(input$obs$spatial_cluster)
+input_ordered <- input[ordered_indices]
+
+counts <- Matrix::t(input_ordered$layers[["counts"]])
+colnames(counts) <- as.character(input_ordered$obs$spatial_cluster)
+
+# cell_type_prop is matched to the fitted models by position, so both have to be
+# in the same order -- and that order has to be the one the spots are sorted in,
+# which unique() gives us since the columns are already grouped per cluster
+cell_types <- unique(colnames(counts))
+spatial_cluster_prop <- prop.table(table(colnames(counts))[cell_types])
 
 copula_result <- scDesign2::fit_model_scDesign2(
   data_mat = counts,
-  cell_type_sel = unique(colnames(counts)),
+  cell_type_sel = cell_types,
   sim_method = "copula",
   ncores = ifelse(!is.null(meta$cpus), meta$cpus, 8L)
 )
@@ -33,17 +44,22 @@ sim_out <- scDesign2::simulate_count_scDesign2(
   copula_result,
   ncol(counts),
   sim_method = "copula",
-  cell_type_prop = prop.table(spatial_cluster_prop)
+  cell_type_prop = spatial_cluster_prop
 )
 
-rownames(sim_out) <- input$var_names
-colnames(sim_out) <- input$obs_names
+rownames(sim_out) <- input_ordered$var_names
+colnames(sim_out) <- input_ordered$obs_names
+
+# 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(sim_out)
+    counts = Matrix::t(sim_out)[restore_order, , drop = FALSE]
   ),
   obs = input$obs[c("row", "col")],
   var = input$var,