Simplify dev

tests/downsampling/make_medium_fixture.R

@@ -35,12 +35,25 @@ sampled_source_idx <- unlist(lapply(names(per_cluster_target), function(cl) {
 sampled_counts <- counts[, sampled_source_idx, drop = FALSE]
 sampled_clusters <- unname(clusters[sampled_source_idx])
 
-jitter <- matrix(
-  rpois(length(sampled_counts), lambda = 0.02),
-  nrow = nrow(sampled_counts),
-  ncol = ncol(sampled_counts)
-)
-sampled_counts <- sampled_counts + jitter
+sampled_counts@x <- sampled_counts@x + rpois(length(sampled_counts@x), lambda = 0.02)
+
+# Sparse jitter alone doesn't reliably separate resampled cells: a cell with few
+# nonzero genes has a non-negligible chance every jitter draw is 0 (P(X=0) ~= 0.98
+# per entry with lambda=0.02), leaving repeat draws of the same source cell exactly
+# identical. Deterministically perturb each repeat draw so it diverges from both the
+# first (unbumped) draw of that source and every other draw of it: bump the Nth draw
+# of a source at row N-1. Keying the row on the per-source occurrence rank (not a
+# global counter) keeps rows distinct within each source group without wrapping, so
+# repeat draws of one source never land on the same row.
+occ_rank <- ave(sampled_source_idx, sampled_source_idx, FUN = seq_along)
+repeat_pos <- which(occ_rank > 1L)
+if (length(repeat_pos) > 0L) {
+  bump_row <- occ_rank[repeat_pos] - 1L
+  stopifnot(max(bump_row) <= nrow(sampled_counts))
+  bump_idx <- cbind(bump_row, repeat_pos)
+  sampled_counts[bump_idx] <- sampled_counts[bump_idx] + 1
+}
+
 colnames(sampled_counts) <- paste0("cell_", seq_len(ncol(sampled_counts)))
 
 new_obj <- CreateSeuratObject(counts = sampled_counts)

tests/downsampling/test_downsample_clusters_rule.py

@@ -1,12 +1,6 @@
 import pytest
 
-from utils import (
-    assert_stability_table,
-    combined_output,
-    repo_root,
-    run_command,
-    snakemake_executable,
-)
+from utils import assert_stability_table, repo_root, run_downsample_clusters_rule
 
 
 SAMPLE = "filtered_seurat_tenx_test"
@@ -17,44 +11,5 @@ def test_downsample_clusters_rule_produces_stability_table(tmp_path, pytestconfi
         pytest.skip("use --run-downsample-rule to execute the downsample_clusters rule")
 
     root = repo_root()
-    seurat_dir = root / "testdata" / "downsampling" / "seurat_objects"
-    results_dir = tmp_path / "results"
-    target = results_dir / f"{SAMPLE}_clusterdownsampling.tsv"
-
-    cmd = [
-        snakemake_executable(),
-        str(target),
-        "--snakefile",
-        "workflow/Snakefile",
-        "--configfile",
-        "config/config.yaml",
-        "--config",
-        "workflow_mode=downsample_only",
-        f"downsampleSeuratObjectDir={seurat_dir.as_posix()}",
-        f"downsampleResultsDir={results_dir.as_posix()}",
-        "nDownsampleReplicates=2",
-        "downsampleRate=0.5",
-        "workflowSeed=12345",
-        "--profile",
-        "none",
-        "--workflow-profile",
-        "none",
-        "--executor",
-        "local",
-        "--cores",
-        "1",
-        "--jobs",
-        "1",
-        "--latency-wait",
-        "30",
-        "--rerun-incomplete",
-        "--use-conda",
-    ]
-    conda_prefix = pytestconfig.getoption("--snakemake-conda-prefix")
-    if conda_prefix:
-        cmd.extend(["--conda-prefix", conda_prefix])
-
-    result = run_command(cmd, root, timeout=1800)
-    assert result.returncode == 0, combined_output(result)
-    assert target.exists(), f"missing rule output: {target}"
+    target = run_downsample_clusters_rule(SAMPLE, tmp_path, pytestconfig, root)
     assert_stability_table(target)

tests/downsampling/test_downsample_scaling.py

@@ -2,14 +2,7 @@
 
 import pytest
 
-from utils import (
-    EXPECTED_COLUMNS,
-    combined_output,
-    read_tsv,
-    repo_root,
-    run_command,
-    snakemake_executable,
-)
+from utils import assert_stability_table, repo_root, run_downsample_clusters_rule
 
 
 SMALL_SAMPLE = "filtered_seurat_tenx_test"
@@ -28,64 +21,12 @@
 CELL_RATIO = MEDIUM_NCELLS / SMALL_NCELLS
 MAX_ALLOWED_RATIO = CELL_RATIO * SAFETY_MULTIPLIER
 
-# This test uses nDownsampleReplicates=1 (unlike the other downsampling tests,
-# which use 2) to keep the scaling comparison itself fast, so it can't reuse
-# utils.assert_stability_table (which hardcodes EXPECTED_BOOTSTRAPS={1, 2}).
-def _assert_basic_stability_table(path):
-    columns, rows = read_tsv(path)
-    assert rows, f"{path}: downsample output is empty"
-    assert columns == EXPECTED_COLUMNS
-    for row in rows:
-        assert row["clusterid"] != ""
-        assert int(row["bootstrap_number"]) == 1
-        assert 0 <= float(row["max_jaccard"]) <= 1
-
-
-def _run_downsample_clusters(sample, tmp_path, pytestconfig, root):
-    seurat_dir = root / "testdata" / "downsampling" / "seurat_objects"
-    results_dir = tmp_path / sample / "results"
-    target = results_dir / f"{sample}_clusterdownsampling.tsv"
-
-    cmd = [
-        snakemake_executable(),
-        str(target),
-        "--snakefile",
-        "workflow/Snakefile",
-        "--configfile",
-        "config/config.yaml",
-        "--config",
-        "workflow_mode=downsample_only",
-        f"downsampleSeuratObjectDir={seurat_dir.as_posix()}",
-        f"downsampleResultsDir={results_dir.as_posix()}",
-        "nDownsampleReplicates=1",
-        "downsampleRate=0.5",
-        "workflowSeed=12345",
-        "--profile",
-        "none",
-        "--workflow-profile",
-        "none",
-        "--executor",
-        "local",
-        "--cores",
-        "1",
-        "--jobs",
-        "1",
-        "--latency-wait",
-        "30",
-        "--rerun-incomplete",
-        "--use-conda",
-    ]
-    conda_prefix = pytestconfig.getoption("--snakemake-conda-prefix")
-    if conda_prefix:
-        cmd.extend(["--conda-prefix", conda_prefix])
 
+def _timed_run(sample, tmp_path, pytestconfig, root):
     t0 = time.monotonic()
-    result = run_command(cmd, root, timeout=1800)
+    target = run_downsample_clusters_rule(sample, tmp_path, pytestconfig, root, n_replicates=1)
     elapsed = time.monotonic() - t0
-
-    assert result.returncode == 0, combined_output(result)
-    assert target.exists(), f"missing rule output: {target}"
-    _assert_basic_stability_table(target)
+    assert_stability_table(target, expected_bootstraps={1})
     return elapsed
 
 
@@ -95,8 +36,8 @@ def test_downsample_clusters_runtime_scales_reasonably_with_ncells(tmp_path, pyt
 
     root = repo_root()
 
-    small_elapsed = _run_downsample_clusters(SMALL_SAMPLE, tmp_path, pytestconfig, root)
-    medium_elapsed = _run_downsample_clusters(MEDIUM_SAMPLE, tmp_path, pytestconfig, root)
+    small_elapsed = _timed_run(SMALL_SAMPLE, tmp_path, pytestconfig, root)
+    medium_elapsed = _timed_run(MEDIUM_SAMPLE, tmp_path, pytestconfig, root)
 
     ratio = medium_elapsed / small_elapsed
     assert ratio <= MAX_ALLOWED_RATIO, (

tests/downsampling/utils.py

@@ -86,6 +86,59 @@ def expected_output_for_sample(sample, results_dir=TEST_RESULTS_DIR):
     return Path(results_dir) / f"{sample}_clusterdownsampling.tsv"
 
 
+def run_downsample_clusters_rule(
+    sample,
+    tmp_path,
+    pytestconfig,
+    root,
+    n_replicates=2,
+    downsample_rate=0.5,
+    workflow_seed=12345,
+    timeout=1800,
+):
+    seurat_dir = root / "testdata" / "downsampling" / "seurat_objects"
+    results_dir = tmp_path / sample / "results"
+    target = results_dir / f"{sample}_clusterdownsampling.tsv"
+
+    cmd = [
+        snakemake_executable(),
+        str(target),
+        "--snakefile",
+        "workflow/Snakefile",
+        "--configfile",
+        "config/config.yaml",
+        "--config",
+        "workflow_mode=downsample_only",
+        f"downsampleSeuratObjectDir={seurat_dir.as_posix()}",
+        f"downsampleResultsDir={results_dir.as_posix()}",
+        f"nDownsampleReplicates={n_replicates}",
+        f"downsampleRate={downsample_rate}",
+        f"workflowSeed={workflow_seed}",
+        "--profile",
+        "none",
+        "--workflow-profile",
+        "none",
+        "--executor",
+        "local",
+        "--cores",
+        "1",
+        "--jobs",
+        "1",
+        "--latency-wait",
+        "30",
+        "--rerun-incomplete",
+        "--use-conda",
+    ]
+    conda_prefix = pytestconfig.getoption("--snakemake-conda-prefix")
+    if conda_prefix:
+        cmd.extend(["--conda-prefix", conda_prefix])
+
+    result = run_command(cmd, root, timeout=timeout)
+    assert result.returncode == 0, combined_output(result)
+    assert target.exists(), f"missing rule output: {target}"
+    return target
+
+
 def reference_output_for_sample(sample):
     return REFERENCE_ROOT / REFERENCE_RESULTS_DIR / f"{sample}_clusterdownsampling.tsv"
 
@@ -97,13 +150,13 @@ def read_tsv(path):
         return reader.fieldnames or [], rows
 
 
-def assert_stability_table(path):
+def assert_stability_table(path, expected_bootstraps=EXPECTED_BOOTSTRAPS):
     columns, rows = read_tsv(path)
     assert rows, f"{path}: downsample output is empty"
     assert columns == EXPECTED_COLUMNS
 
     bootstraps = {int(row["bootstrap_number"]) for row in rows}
-    assert bootstraps == EXPECTED_BOOTSTRAPS
+    assert bootstraps == expected_bootstraps
     for row in rows:
         assert row["clusterid"] != ""
         max_jaccard = float(row["max_jaccard"])

workflow/scripts/downsample_clusters.R

@@ -57,11 +57,8 @@ SubSampleReSCTSeuratObject <- function(seurat_obj, subrate, replicate_seed) {
     drop = TRUE
   ]
 
-  if ("percent.mt" %in% colnames(subsampled_seurat@meta.data)) {
-    subsampled_seurat <- SCTransform(subsampled_seurat, vars.to.regress = "percent.mt", verbose = FALSE)
-  } else {
-    subsampled_seurat <- SCTransform(subsampled_seurat, verbose = FALSE)
-  }
+  vars_to_regress <- if ("percent.mt" %in% colnames(subsampled_seurat@meta.data)) "percent.mt" else NULL
+  subsampled_seurat <- SCTransform(subsampled_seurat, vars.to.regress = vars_to_regress, verbose = FALSE)
   subsampled_seurat <- RunPCA(subsampled_seurat, verbose = FALSE)
   pca_dims <- seq_len(min(30, ncol(Embeddings(subsampled_seurat, "pca"))))
   subsampled_seurat <- FindNeighbors(subsampled_seurat, dims = pca_dims, verbose = FALSE)
@@ -80,18 +77,18 @@ GetJaccardMaxByCluster <- function(seurat_obj, bootstrap) {
     bootstrap_number = integer()
   )
 
+  dat <- tibble::tibble(
+    cell_id = names(seurat_obj@active.ident),
+    cluster = seurat_obj$seurat_clusters
+  ) %>%
+    tidyr::nest(data = -cluster) %>%
+    dplyr::arrange(cluster)
+
   for (original_cluster in unique(seurat_obj$presub_clusters)) {
     barcodes <- rownames(
       subset(seurat_obj@meta.data, presub_clusters == original_cluster)
     )
 
-    dat <- tibble::tibble(
-      cell_id = names(seurat_obj@active.ident),
-      cluster = seurat_obj$seurat_clusters
-    ) %>%
-      tidyr::nest(data = -cluster) %>%
-      dplyr::arrange(cluster)
-
     maxstat <- dat %>%
       dplyr::mutate(
         jaccard = purrr::map(data, ~ JaccardSimilarity(barcodes, .x$cell_id))