Features in `splatter` simulation will be inconsistent

[lazappi]

Currently the splatter method estimates parameters from each cell type and simulates them separately. The different simulated datasets are them combined.

task_spatial_simulators/src/methods/splatter/script.R

Lines 26 to 43 in cd39dd4

	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)))
	# combine the cell types
	if (is.null(simulated_result)) {
	simulated_result <- sim_spatial_cluster
	} else {
	simulated_result <- SingleCellExperiment::cbind(simulated_result, sim_spatial_cluster)
	}
	})
	}

This seems like a sensible approach but it will result in an inconsistent final dataset because the simulations for each cell type are completely independent. This means that the simulated "Gene1" for cell type A has no relationship to "Gene1" for cell type B and the two cell type can not be combined.

[rcannood]

If I let GitHub Copilot iterate on the component based on your feedback, I get this:

suppressMessages(library(SingleCellExperiment, quietly = TRUE))
suppressMessages(library(splatter, quietly = TRUE))
suppressMessages(library(Matrix, quietly = TRUE))

## VIASH START
par <- list(
  input = "resources_test/spatialsimbench_mobnew/dataset_sp.h5ad",
  base = "domain"
)
meta <- list(
  name = "splatter"
)
## VIASH END

cat("Reading input files\n")
input <- anndata::read_h5ad(par$input)

cat("Splatter simulation start\n")

if (par$base != "domain") {
  stop("ONLY domain base")
}

ordered_indices <- order(input$obs$spatial_cluster)
input_ordered <- input[ordered_indices]

# Estimate parameters from the entire dataset to maintain gene relationships across cell types
cat("Estimating splatter parameters from entire dataset\n")
global_params <- splatter::splatEstimate(as.matrix(t(input_ordered$layers[["counts"]])))

# Get the number of genes from input to ensure consistency
n_genes <- nrow(input_ordered$var)
global_params <- splatter::setParam(global_params, "nGenes", n_genes)

simulated_result <- NULL
cell_cluster_mapping <- c()  # Track which cluster each simulated cell belongs to

for (spatial_cluster in unique(input_ordered$obs[["spatial_cluster"]])) {
  res <- try({
    input_spatial_cluster <- input_ordered[input_ordered$obs[["spatial_cluster"]] == spatial_cluster]
    n_cells <- ncol(input_spatial_cluster)
    
    # Use global parameters but simulate only the number of cells for this cluster
    sim_spatial_cluster <- splatter::splatSimulate(global_params, batchCells = n_cells)
    sim_spatial_cluster$spatial_cluster <- spatial_cluster
    
    # Track cluster assignment for each cell
    cell_cluster_mapping <- c(cell_cluster_mapping, rep(spatial_cluster, ncol(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)
    }
  })
}

cat("Input cells:", nrow(input_ordered$obs), ", Simulated cells:", ncol(simulated_result), "\n")

# If the number of simulated cells differs from input, we need to handle this
if (ncol(simulated_result) != nrow(input_ordered$obs)) {
  # Sample or duplicate cells to match the input size
  n_input_cells <- nrow(input_ordered$obs)
  n_sim_cells <- ncol(simulated_result)
  
  if (n_sim_cells > n_input_cells) {
    # Too many simulated cells - sample down
    selected_cells <- sample(seq_len(n_sim_cells), n_input_cells)
    simulated_result <- simulated_result[, selected_cells]
    cell_cluster_mapping <- cell_cluster_mapping[selected_cells]
    cat("Sampled down from", n_sim_cells, "to", n_input_cells, "cells\n")
  } else {
    # Too few simulated cells - duplicate some
    additional_needed <- n_input_cells - n_sim_cells
    duplicated_cells <- sample(seq_len(n_sim_cells), additional_needed, replace = TRUE)
    simulated_result <- SingleCellExperiment::cbind(simulated_result, simulated_result[, duplicated_cells])
    cell_cluster_mapping <- c(cell_cluster_mapping, cell_cluster_mapping[duplicated_cells])
    cat("Duplicated cells from", n_sim_cells, "to", n_input_cells, "cells\n")
  }
}

# Now assign proper column and row names - making sure dimensions match
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))
]

cat("Generating output\n")
output <- anndata::AnnData(
  layers = list(counts = t(simulated_result_ordered)),
  obs = input_ordered$obs[c("row", "col")],
  var = input_ordered$var,
  uns = c(
    input$uns,
    list(
      method_id = meta$name
    )
  )
)

cat("Write output files\n")
output$write_h5ad(par$output, compression = "gzip")

Which returns:

>> Running test 'run'
>> Checking whether input files exist
>> Running script as test
Reading input files
Splatter simulation start
Estimating splatter parameters from entire dataset
Getting parameters...
Creating simulation object...
Simulating library sizes...
Simulating gene means...
Simulating BCV...
Simulating counts...
Simulating dropout (if needed)...
Sparsifying assays...
Automatically converting to sparse matrices, threshold = 0.95
Skipping 'BatchCellMeans': estimated sparse size 1.5 * dense matrix
Skipping 'BaseCellMeans': estimated sparse size 1.5 * dense matrix
Skipping 'BCV': estimated sparse size 1.5 * dense matrix
Skipping 'CellMeans': estimated sparse size 1.5 * dense matrix
Skipping 'TrueCounts': estimated sparse size 1.76 * dense matrix
Skipping 'counts': estimated sparse size 1.76 * dense matrix
Done!
Getting parameters...
Creating simulation object...
Simulating library sizes...
Simulating gene means...
Simulating BCV...
Simulating counts...
Simulating dropout (if needed)...
Sparsifying assays...
Automatically converting to sparse matrices, threshold = 0.95
Skipping 'BatchCellMeans': estimated sparse size 1.5 * dense matrix
Skipping 'BaseCellMeans': estimated sparse size 1.5 * dense matrix
Skipping 'BCV': estimated sparse size 1.5 * dense matrix
Skipping 'CellMeans': estimated sparse size 1.5 * dense matrix
Skipping 'TrueCounts': estimated sparse size 1.76 * dense matrix
Skipping 'counts': estimated sparse size 1.76 * dense matrix
Done!
Getting parameters...
Creating simulation object...
Simulating library sizes...
Simulating gene means...
Simulating BCV...
Simulating counts...
Simulating dropout (if needed)...
Sparsifying assays...
Automatically converting to sparse matrices, threshold = 0.95
Skipping 'BatchCellMeans': estimated sparse size 1.5 * dense matrix
Skipping 'BaseCellMeans': estimated sparse size 1.5 * dense matrix
Skipping 'BCV': estimated sparse size 1.5 * dense matrix
Skipping 'CellMeans': estimated sparse size 1.5 * dense matrix
Skipping 'TrueCounts': estimated sparse size 1.76 * dense matrix
Skipping 'counts': estimated sparse size 1.76 * dense matrix
Done!
Getting parameters...
Creating simulation object...
Simulating library sizes...
Simulating gene means...
Simulating BCV...
Simulating counts...
Simulating dropout (if needed)...
Sparsifying assays...
Automatically converting to sparse matrices, threshold = 0.95
Skipping 'BatchCellMeans': estimated sparse size 1.5 * dense matrix
Skipping 'BaseCellMeans': estimated sparse size 1.5 * dense matrix
Skipping 'BCV': estimated sparse size 1.5 * dense matrix
Skipping 'CellMeans': estimated sparse size 1.5 * dense matrix
Skipping 'TrueCounts': estimated sparse size 1.76 * dense matrix
Skipping 'counts': estimated sparse size 1.76 * dense matrix
Done!
Input cells: 278 , Simulated cells: 728 
Sampled down from 728 to 278 cells
Generating output
Write output files
>> Checking whether output file exists
>> Reading h5ad files and checking formats
All checks succeeded!

I'm sure that there are quite some things that need to be improved -- rather than me investigating which approach would be best, maybe @lazappi you could already give some general feedback?

[lazappi]

This will still be a problem as long as you have multiple calls to splatSimulate(). The only way to make sure the genes match up across cell types/batches is to provide all parameters to a single call. Some of these are not estimated so will need to be calculated/set in the script.