We would like to assess the statistical significance of module scores obtained using Seurat’s AddModuleScore, particularly to define a threshold to classify cells as positive or negative. We are working with several gene signatures that we want to visualize spatially. After running AddModuleScore, we observe different score ranges across signatures: some reach values around 2, while others remain closer to 0.5. To address this, we are considering a permutation-based approach:

• For each gene signature, we generate random gene sets matched by expression level (binning genes into 24 bins, as in Seurat).

• Module scores are computed both for the real signature and for the permuted (random) gene sets.

• Permutation multiple times and calculate a summary statistic is calculated (e.g., mean score) for each cell.Based on the null distribution from vthe permutations, we compute a p-value for each cell’s observed module score.

• We then apply multiple testing correction and count how many cells have adjusted p-values < 0.05.

all adjusted p-values are 0.

gene_means <- rowMeans(obj@assays$SCT@data) bins <- cut(gene_means, breaks = n_bins, labels = FALSE) names(bins) <- names(gene_means)

  all_results <- list()

  for (set_name in names(gene_sets)) {

    cat("Processing:", set_name, "\n")
    gene_set <- gene_sets[[set_name]]
    gene_set <- gene_set[gene_set %in% names(bins)]


    random_gene_sets <- replicate(n_perm, {
      sapply(gene_set, function(gene) {
        bin_id <- bins[gene]
        pool <- names(bins[bins == bin_id])
        sample(pool, 1)
      })
    }, simplify = FALSE)


    obj <- AddModuleScore(obj, features = list(gene_set), name = paste0(set_name, "_real"))


    obj <- AddModuleScore(obj, features = random_gene_sets, name = paste0(set_name, "_perm"))


    real_values <- obj[[paste0(set_name, "_real1")]]

    perm_matrix <- as.matrix(obj@meta.data[, grep(paste0(set_name, "_perm"), colnames(obj@meta.data))])

    # p-values
    pvals <- sapply(1:length(real_values), function(i) {
      mean(perm_matrix[i, ] >= real_values[i])
    })

    pvals_adj <- p.adjust(pvals, method = "BH")


    obj[[paste0(set_name, "_pval")]] <- pvals
    obj[[paste0(set_name, "_fdr")]] <- pvals_adj
    obj[[paste0(set_name, "_signif")]] <- pvals_adj < 0.05


    all_results[[set_name]] <- data.frame(
      signature = set_name,
      n_signif = sum(pvals_adj < 0.05),
      prop_signif = mean(pvals_adj < 0.05)
    )
  }

  summary_df <- do.call(rbind, all_results)

  return(list(obj = obj, summary = summary_df))
}