diff --git a/testdata/downsampling/seurat_objects/filtered_seurat_tenx_test_medium.rds b/testdata/downsampling/seurat_objects/filtered_seurat_tenx_test_medium.rds index 84b89e4..a47bfc6 100644 Binary files a/testdata/downsampling/seurat_objects/filtered_seurat_tenx_test_medium.rds and b/testdata/downsampling/seurat_objects/filtered_seurat_tenx_test_medium.rds differ diff --git a/tests/downsampling/make_medium_fixture.R b/tests/downsampling/make_medium_fixture.R index 103f7de..cf4c575 100644 --- a/tests/downsampling/make_medium_fixture.R +++ b/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) diff --git a/tests/downsampling/test_downsample_clusters_rule.py b/tests/downsampling/test_downsample_clusters_rule.py index d9e87d4..2fd259f 100644 --- a/tests/downsampling/test_downsample_clusters_rule.py +++ b/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) diff --git a/tests/downsampling/test_downsample_scaling.py b/tests/downsampling/test_downsample_scaling.py index aeac43e..7fe5fb3 100644 --- a/tests/downsampling/test_downsample_scaling.py +++ b/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, ( diff --git a/tests/downsampling/utils.py b/tests/downsampling/utils.py index ce3223b..4a5bb26 100644 --- a/tests/downsampling/utils.py +++ b/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"]) diff --git a/workflow/scripts/downsample_clusters.R b/workflow/scripts/downsample_clusters.R index 8167624..0d37790 100644 --- a/workflow/scripts/downsample_clusters.R +++ b/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))