diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 2cfa65d..2502cf4 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -31,17 +31,17 @@ jobs: SNAKEMAKE_CONDA_FRONTEND: mamba steps: - name: Check out repository - uses: actions/checkout@v4 + uses: actions/checkout@v7 - name: Create test environment - uses: mamba-org/setup-micromamba@v2 + uses: mamba-org/setup-micromamba@v3 with: environment-file: environment.yml cache-environment: true cache-downloads: true - name: Cache Snakemake rule environments - uses: actions/cache@v4 + uses: actions/cache@v6 with: path: .snakemake/conda key: snakemake-conda-${{ runner.os }}-${{ hashFiles('workflow/envs/*.yml') }} @@ -49,7 +49,65 @@ jobs: snakemake-conda-${{ runner.os }}- - name: Run tests - run: pytest tests -q + run: pytest tests -q --run-downsample-scaling + + downsample-rule: + runs-on: ubuntu-latest + defaults: + run: + shell: bash -el {0} + env: + SNAKEMAKE_CONDA_FRONTEND: mamba + steps: + - name: Check out repository + uses: actions/checkout@v7 + + - name: Create test environment + uses: mamba-org/setup-micromamba@v3 + with: + environment-file: environment.yml + cache-environment: true + cache-downloads: true + + - name: Cache Snakemake rule environments + uses: actions/cache@v6 + with: + path: .snakemake/conda + key: snakemake-conda-${{ runner.os }}-${{ hashFiles('workflow/envs/*.yml') }} + restore-keys: | + snakemake-conda-${{ runner.os }}- + + - name: Run downsampling rule smoke test + run: pytest tests/downsampling --run-downsample-rule -q + + downsample-workflow: + runs-on: ubuntu-latest + defaults: + run: + shell: bash -el {0} + env: + SNAKEMAKE_CONDA_FRONTEND: mamba + steps: + - name: Check out repository + uses: actions/checkout@v7 + + - name: Create test environment + uses: mamba-org/setup-micromamba@v3 + with: + environment-file: environment.yml + cache-environment: true + cache-downloads: true + + - name: Cache Snakemake rule environments + uses: actions/cache@v6 + with: + path: .snakemake/conda + key: snakemake-conda-${{ runner.os }}-${{ hashFiles('workflow/envs/*.yml') }} + restore-keys: | + snakemake-conda-${{ runner.os }}- + + - name: Run full downsampling workflow test + run: pytest tests/downsampling --run-downsample-workflow -q conda-validation: if: github.event_name == 'workflow_dispatch' && inputs.run_conda_validation @@ -60,6 +118,7 @@ jobs: env_name: - cellbender.yml - doubletfinder.yml + - downsample_clusters.yml - emptydrops.yml - posthocfilter.yml - scdblfinder.yml @@ -72,10 +131,10 @@ jobs: SNAKEMAKE_CONDA_FRONTEND: mamba steps: - name: Check out repository - uses: actions/checkout@v4 + uses: actions/checkout@v7 - name: Create test environment - uses: mamba-org/setup-micromamba@v2 + uses: mamba-org/setup-micromamba@v3 with: environment-file: environment.yml cache-environment: true @@ -92,10 +151,10 @@ jobs: shell: bash -el {0} steps: - name: Check out repository - uses: actions/checkout@v4 + uses: actions/checkout@v7 - name: Create test environment - uses: mamba-org/setup-micromamba@v2 + uses: mamba-org/setup-micromamba@v3 with: environment-file: environment.yml cache-environment: true @@ -117,10 +176,10 @@ jobs: GITHUB_PAT: ${{ github.token }} steps: - name: Check out repository - uses: actions/checkout@v4 + uses: actions/checkout@v7 - name: Create test environment - uses: mamba-org/setup-micromamba@v2 + uses: mamba-org/setup-micromamba@v3 with: environment-file: environment.yml cache-environment: true diff --git a/.gitignore b/.gitignore index cb76481..e2afe8f 100644 --- a/.gitignore +++ b/.gitignore @@ -13,9 +13,16 @@ scratch/ *out results/ testdata/results/ +testdata/downsampling/results/ # Keep golden test outputs trackable even though results/ is ignored. !tests/reference_outputs/ !tests/reference_outputs/** !tests/reference_outputs/testdata/ !tests/reference_outputs/testdata/results/ +testdata/downsampling/results/ !tests/reference_outputs/testdata/results/** +!tests/downsampling/reference_outputs/ +!tests/downsampling/reference_outputs/** +!tests/downsampling/reference_outputs/testdata/ +!tests/downsampling/reference_outputs/testdata/results/ +!tests/downsampling/reference_outputs/testdata/results/** diff --git a/README.md b/README.md index 137145f..a5bf4bb 100644 --- a/README.md +++ b/README.md @@ -20,6 +20,17 @@ In order to assess different tools and options for pre-processing scRNA-seq data | Mouse | C57BL/6 | L8TX_190509_01_E09 | striatum, striatal amygdala | cells | Chromium 10x 3' Gene Expression v3 | NovaSeq6000 | 13,475 | 82,801 | 3,400 | [nemo](https://assets.nemoarchive.org/dat-qg7n1b0) | [run1](https://data.nemoarchive.org/biccn/grant/u19_zeng/zeng/transcriptome/scell/10x_v2/mouse/raw/STR/NW_TX0007-8_S01_L003.fastq.tar);[run2](https://data.nemoarchive.org/biccn/grant/u19_zeng/zeng/transcriptome/scell/10x_v2/mouse/raw/STR/NW_TX0010-7_S01_L003.fastq.tar) | 2 runs on same library | | Mouse | C57BL/6 | L8TX_210204_01_H05 | olfactory region: main and accessory olfactory bulbs | cells | Chromium 10x 3' Gene Expression v3 | NovaSeq6000 | 10,895 | 136,593 | 3,971 | [nemo](https://assets.nemoarchive.org/dat-qg7n1b0) | [fastq]( https://data.nemoarchive.org/biccn/grant/u19_zeng/zeng/transcriptome/scell/10x_v3/mouse/raw/NW_TX0109-4_S01_L003-001.fastq.tar) | | + +## Workflow modes + +The main Snakemake entrypoint supports three `workflow_mode` values in `config/config.yaml` or via `--config`: + +- `preprocess`: run the preprocessing workflow only. This is the default and preserves the original behavior. +- `preprocess_and_downsample`: run preprocessing and then downsample the generated Seurat `.rds` outputs with `workflow/rules/downsample_clusters.smk`. +- `downsample_only`: skip preprocessing and run downsampling on existing Seurat `.rds` files from `downsampleSeuratObjectDir`. + +Downsampling outputs are written to `downsampleResultsDir`, defaulting to `results/downsampling`. Use `downsampleTargets` to restrict downsampling to selected Seurat object basenames, or leave it as `all` to use every available input for the selected mode. + ## Tests For information on how to run the test suite, or run the workflow in test mode, see tests/README.md. diff --git a/config/config.yaml b/config/config.yaml index 6ac8b47..741101a 100644 --- a/config/config.yaml +++ b/config/config.yaml @@ -1,6 +1,7 @@ conda-channel-priority: "strict" sampleTable: "samplesheet.tsv" +workflow_mode: "preprocess_and_downsample" # {"preprocess", "preprocess_and_downsample", "downsample_only"} workflow_seed: 12345 emptydrop_removal_methods: ["tenx","emptydrops"] @@ -10,3 +11,9 @@ posthoc_methods: ["threshold", "mad"] min_nfeature: 200 min_ncount: 500 max_mtdna: 5 + +downsampleSeuratObjectDir: "seurat_objects" +downsampleResultsDir: "results/downsampling" +downsampleTargets: "all" +nDownsampleReplicates: 100 +downsampleRate: 0.8 diff --git a/profiles/slurm/config.yaml b/profiles/slurm/config.yaml index c81df73..14957f9 100644 --- a/profiles/slurm/config.yaml +++ b/profiles/slurm/config.yaml @@ -7,6 +7,8 @@ set-threads: soupx: 1 scdblfinder: 1 emptydrops: 1 + downsample_cluster_replicate: 1 + downsample_clusters: 1 set-resources: tenx2seuratrds: @@ -43,3 +45,13 @@ set-resources: nodes: 1 tasks: 1 cpus_per_task: 1 + + downsample_cluster_replicate: + nodes: 1 + tasks: 1 + cpus_per_task: 1 + + downsample_clusters: + nodes: 1 + tasks: 1 + cpus_per_task: 1 diff --git a/scrnaseq_preprocess_slurmrunner.sh b/scrnaseq_preprocess_slurmrunner.sh index c8eedce..e11056a 100755 --- a/scrnaseq_preprocess_slurmrunner.sh +++ b/scrnaseq_preprocess_slurmrunner.sh @@ -35,4 +35,4 @@ PATH_TO_MY_CONDA_ENVS=$1 snakemake --unlock --snakefile workflow/Snakefile --configfile config/config.yaml --use-conda --workflow-profile profiles/slurm --profile cannon -snakemake --conda-prefix $PATH_TO_MY_CONDA_ENVS --snakefile workflow/Snakefile --rerun-incomplete --retries 2 --jobs 200 --latency-wait 120 --configfile config/config.yaml --use-conda --workflow-profile profiles/slurm --profile cannon +snakemake --conda-prefix $PATH_TO_MY_CONDA_ENVS --snakefile workflow/Snakefile --rerun-incomplete --retries 2 --keep-going --jobs 1500 --max-jobs-per-timespan "10/1s" --max-status-checks-per-second 5 --latency-wait 120 --configfile config/config.yaml --use-conda --workflow-profile profiles/slurm --profile cannon diff --git a/testdata/downsampling/seurat_objects/filtered_seurat_emptydrops_test.rds b/testdata/downsampling/seurat_objects/filtered_seurat_emptydrops_test.rds new file mode 100644 index 0000000..99c3a33 Binary files /dev/null and b/testdata/downsampling/seurat_objects/filtered_seurat_emptydrops_test.rds differ diff --git a/testdata/downsampling/seurat_objects/filtered_seurat_tenx_test.rds b/testdata/downsampling/seurat_objects/filtered_seurat_tenx_test.rds new file mode 100644 index 0000000..99d5516 Binary files /dev/null and b/testdata/downsampling/seurat_objects/filtered_seurat_tenx_test.rds differ diff --git a/testdata/downsampling/seurat_objects/filtered_seurat_tenx_test_medium.rds b/testdata/downsampling/seurat_objects/filtered_seurat_tenx_test_medium.rds new file mode 100644 index 0000000..84b89e4 Binary files /dev/null and b/testdata/downsampling/seurat_objects/filtered_seurat_tenx_test_medium.rds differ diff --git a/tests/README.md b/tests/README.md index bffdefe..c7c07af 100644 --- a/tests/README.md +++ b/tests/README.md @@ -22,6 +22,8 @@ The dry-run DAG test uses the small CellRanger-style input data in `testdata/` a GitHub Actions runs the default test suite via `.github/workflows/tests.yml`, using the same top-level `environment.yml` test runner environment and caching Snakemake-created rule environments under `.snakemake/conda`. +The default suite also includes the downsampling fixture/reference checks and a dry-run of `workflow_mode=downsample_only`. Optional downsampling execution tests can be run with `pytest tests/downsampling --run-downsample-rule -q` or `pytest tests/downsampling --run-downsample-workflow -q`, and both are available as manual GitHub Actions dispatch jobs. + The default test suite also includes a focused local rule-execution smoke test. It runs the real `tenx2seuratrds`, `find_markers`, and `combine_markers` rule chain against `testdata/`, using Snakemake's `--use-conda` support and writing outputs under pytest's temporary directory. This catches broken R package imports, script argument drift, invalid Seurat object creation, and marker CSV schema changes without submitting to SLURM. The R output validator also checks that the Seurat object has at least 100 features and 100 cells; metadata rows match the cell count; barcode row names are present, unique, and nonempty; `orig.ident`, `nCount_RNA`, `nFeature_RNA`, `percent.mt`, and `seurat_clusters` metadata columns exist; RNA count and feature-count metadata values are finite and positive; mitochondrial percentages are finite and within `[0, 100]`; at least two clusters are present; PCA and UMAP reductions exist; the marker table is nonempty and has the expected columns; marker gene symbols are present and nonempty; marker numeric columns are finite; marker p-value and percent columns are within `[0, 1]`; marker clusters are present in the Seurat metadata; markers are reported for at least two clusters; and the marker `workflow` column matches the expected test workflow label. The test runner environment is defined in the repository-level `environment.yml`; the rule-specific R environment is still created by Snakemake from `workflow/envs/tenx2seuratrds.yml`. A separate lightweight checkpoint-expansion test uses a fake `Rscript` to materialize the `marker_manifest` checkpoint, verify that dynamic `find_markers` jobs are generated for each cluster id, and confirm that `combine_markers` receives the expected marker chunks. @@ -48,7 +50,7 @@ The default tests build and inspect the DAG and run a focused local R-rule smoke pytest tests --run-workflow ``` -The full-run test calls `tests/run_test_workflow.sh`, which uses `testdata/samplesheet_test.tsv` and overrides the workflow output directory with `resultsDir=testdata/results`. The manifest in `tests/test_sample_rule_output_files.txt` is therefore written with paths under `testdata/results/`. +The full-run test calls `tests/run_test_workflow.sh`, which uses `testdata/samplesheet_test.tsv` and overrides the workflow output directories with `resultsDir=testdata/results` and `downsampleResultsDir=testdata/results/downsampling`. The manifest in `tests/test_sample_rule_output_files.txt` is therefore written with paths under `testdata/results/`. For testing, omit `--snakemake-conda-prefix` so Snakemake uses its default `.snakemake/conda` location under the repository root. The runner assumes that the current environment already provides `snakemake` on `PATH`. diff --git a/tests/conftest.py b/tests/conftest.py index 4e69f78..89e6d63 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -10,6 +10,24 @@ def pytest_addoption(parser): default=None, help="optional value to pass to Snakemake --conda-prefix for --run-workflow", ) + parser.addoption( + "--run-downsample-rule", + action="store_true", + default=False, + help="execute the downsample_clusters Snakemake rule on one downsampling test fixture", + ) + parser.addoption( + "--run-downsample-workflow", + action="store_true", + default=False, + help="execute the full downsampling Snakemake workflow on downsampling testdata", + ) + parser.addoption( + "--run-downsample-scaling", + action="store_true", + default=False, + help="run downsample_clusters on small and medium fixtures and check runtime scales reasonably", + ) parser.addoption( "--run-conda-validation", action="store_true", diff --git a/tests/downsampling/README.md b/tests/downsampling/README.md new file mode 100644 index 0000000..e60af6f --- /dev/null +++ b/tests/downsampling/README.md @@ -0,0 +1,53 @@ +# Downsampling Tests + +The default pytest suite uses the Seurat RDS fixtures in `testdata/downsampling/seurat_objects/`, validates the saved reference TSVs, and runs a Snakemake dry run in `workflow_mode=downsample_only` against `testdata/downsampling/results/`: + +```bash +pytest tests/downsampling -q +``` + +The dry-run test checks that the DAG includes the `downsample_clusters` rule (which loops over all replicates for a target in a single job) and expected testdata outputs without executing R code. + +To execute the single workflow rule on one test fixture, run the optional rule test: + +```bash +pytest tests/downsampling --run-downsample-rule -q +``` + +To execute the full downsampling test workflow and write outputs under `testdata/downsampling/results/`, run: + +```bash +pytest tests/downsampling --run-downsample-workflow -q +``` + +The full-workflow test removes `testdata/downsampling/results/` before running so that both expected test fixture outputs are regenerated from scratch. It validates each output table and compares it to the reference snapshot under `tests/downsampling/reference_outputs/testdata/results/`. + +All three of `--run-downsample-rule`, `--run-downsample-workflow`, and `--run-downsample-scaling` run automatically in CI on every push/PR (as separate jobs in `.github/workflows/tests.yml`) rather than being opt-in-only — these flags mainly exist so they can also be run manually/locally without the rest of the suite. + +```bash +pytest tests/downsampling --run-downsample-scaling -q +``` + +It runs `downsample_clusters` once on the small fixture (`filtered_seurat_tenx_test`, 300 cells) and once on a synthetic mid-scale fixture (`filtered_seurat_tenx_test_medium`, 2500 cells, generated by stratified-resampling the small fixture via `tests/downsampling/make_medium_fixture.R`), then asserts the wall-clock ratio between the two doesn't wildly exceed the ~8.3x cell-count ratio. This is the kind of check that would have caught the `do.call(SCTransform, ...)` performance regression (see `workflow/scripts/downsample_clusters.R`) — that bug produced no incorrect output and no error, only a data-size-dependent slowdown invisible to the tiny, single-size fixtures the other tests use. + +By default Snakemake creates rule conda environments under `.snakemake/conda` in the repository root. To use an explicit Snakemake conda prefix: + +```bash +pytest tests/downsampling --run-downsample-workflow --snakemake-conda-prefix /path/to/snakemake-conda-envs -q +``` + +Snakemake uses the explicit Linux lock file `workflow/envs/downsample_clusters.linux-64.pin.txt` when creating the rule environment on Linux. The workflow test uses two downsampling replicates, a 50% downsampling rate, and `workflowSeed=12345`. + +Reference outputs are compared by `(bootstrap_number, clusterid)`. The `max_jaccard` column is numeric-compared with `DOWNSAMPLE_REFERENCE_REL_TOLERANCE` defaulting to `1e-6` and `DOWNSAMPLE_REFERENCE_ABS_TOLERANCE` defaulting to `1e-8`. + +To refresh the reference outputs after an intentional workflow change, first run: + +```bash +pytest tests/downsampling --run-downsample-workflow -q +``` + +Then copy the regenerated outputs into `tests/downsampling/reference_outputs/`: + +```bash +python tests/downsampling/update_reference_outputs.py +``` diff --git a/tests/downsampling/conftest.py b/tests/downsampling/conftest.py new file mode 100644 index 0000000..f5754e3 --- /dev/null +++ b/tests/downsampling/conftest.py @@ -0,0 +1 @@ +# Downsampling-specific options are registered in tests/conftest.py. diff --git a/tests/downsampling/make_medium_fixture.R b/tests/downsampling/make_medium_fixture.R new file mode 100644 index 0000000..103f7de --- /dev/null +++ b/tests/downsampling/make_medium_fixture.R @@ -0,0 +1,58 @@ +# Generates a mid-scale downsampling test fixture by stratified resampling +# (with light Poisson jitter) of the existing tiny fixture, preserving +# per-cluster proportions. This exists purely to give the scaling test two +# data points of meaningfully different size, so quadratic-or-worse +# performance regressions (e.g. the do.call()/SCTransform hang) are +# detectable in CI without needing production-scale data. +# +# Usage: Rscript make_medium_fixture.R [seed] + +args <- commandArgs(trailingOnly = TRUE) +if (length(args) < 3) { + stop("Usage: make_medium_fixture.R [seed]", call. = FALSE) +} + +input_rds <- args[1] +output_rds <- args[2] +target_ncells <- as.integer(args[3]) +seed <- if (length(args) >= 4) as.integer(args[4]) else 12345L + +suppressPackageStartupMessages(library(Seurat)) +set.seed(seed) + +obj <- readRDS(input_rds) +counts <- obj[["RNA"]]$counts +clusters <- obj$seurat_clusters + +cluster_props <- table(clusters) / length(clusters) +per_cluster_target <- round(cluster_props * target_ncells) + +sampled_source_idx <- unlist(lapply(names(per_cluster_target), function(cl) { + pool <- which(clusters == cl) + sample(pool, size = per_cluster_target[[cl]], replace = TRUE) +})) + +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 +colnames(sampled_counts) <- paste0("cell_", seq_len(ncol(sampled_counts))) + +new_obj <- CreateSeuratObject(counts = sampled_counts) +pct_mt <- PercentageFeatureSet(new_obj, pattern = "^MT-") +new_obj$percent.mt <- if (is.data.frame(pct_mt)) pct_mt[[1]] else unname(pct_mt) +new_obj$seurat_clusters <- factor(sampled_clusters, levels = levels(clusters)) +Idents(new_obj) <- new_obj$seurat_clusters + +dir.create(dirname(output_rds), showWarnings = FALSE, recursive = TRUE) +saveRDS(new_obj, output_rds) + +message(sprintf( + "Wrote %s: %d cells x %d genes (%d clusters)", + output_rds, ncol(new_obj), nrow(new_obj), length(levels(new_obj$seurat_clusters)) +)) diff --git a/tests/downsampling/reference_outputs/testdata/results/filtered_seurat_emptydrops_test_clusterdownsampling.tsv b/tests/downsampling/reference_outputs/testdata/results/filtered_seurat_emptydrops_test_clusterdownsampling.tsv new file mode 100644 index 0000000..1c374d0 --- /dev/null +++ b/tests/downsampling/reference_outputs/testdata/results/filtered_seurat_emptydrops_test_clusterdownsampling.tsv @@ -0,0 +1,9 @@ +clusterid max_jaccard bootstrap_number +3 0.8421052631578947 1 +2 0.29896907216494845 1 +0 0.3431372549019608 1 +1 0.32989690721649484 1 +3 0.9393939393939394 2 +2 0.64 2 +1 0.5263157894736842 2 +0 0.3384615384615385 2 diff --git a/tests/downsampling/reference_outputs/testdata/results/filtered_seurat_tenx_test_clusterdownsampling.tsv b/tests/downsampling/reference_outputs/testdata/results/filtered_seurat_tenx_test_clusterdownsampling.tsv new file mode 100644 index 0000000..227270d --- /dev/null +++ b/tests/downsampling/reference_outputs/testdata/results/filtered_seurat_tenx_test_clusterdownsampling.tsv @@ -0,0 +1,7 @@ +clusterid max_jaccard bootstrap_number +0 0.5966386554621849 1 +1 0.4492753623188406 1 +2 0.0847457627118644 1 +0 0.5957446808510638 2 +1 0.6666666666666666 2 +2 0.13157894736842105 2 diff --git a/tests/downsampling/test_downsample_clusters_rule.py b/tests/downsampling/test_downsample_clusters_rule.py new file mode 100644 index 0000000..d9e87d4 --- /dev/null +++ b/tests/downsampling/test_downsample_clusters_rule.py @@ -0,0 +1,60 @@ +import pytest + +from utils import ( + assert_stability_table, + combined_output, + repo_root, + run_command, + snakemake_executable, +) + + +SAMPLE = "filtered_seurat_tenx_test" + + +def test_downsample_clusters_rule_produces_stability_table(tmp_path, pytestconfig): + if not pytestconfig.getoption("--run-downsample-rule"): + 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}" + assert_stability_table(target) diff --git a/tests/downsampling/test_downsample_scaling.py b/tests/downsampling/test_downsample_scaling.py new file mode 100644 index 0000000..aeac43e --- /dev/null +++ b/tests/downsampling/test_downsample_scaling.py @@ -0,0 +1,107 @@ +import time + +import pytest + +from utils import ( + EXPECTED_COLUMNS, + combined_output, + read_tsv, + repo_root, + run_command, + snakemake_executable, +) + + +SMALL_SAMPLE = "filtered_seurat_tenx_test" +SMALL_NCELLS = 300 +MEDIUM_SAMPLE = "filtered_seurat_tenx_test_medium" +MEDIUM_NCELLS = 2500 + +# How much slower the medium fixture is allowed to run relative to the small +# fixture, expressed as a multiple of the plain cell-count ratio. Real Seurat +# steps (SCTransform, PCA, neighbor graphs, clustering) scale roughly linearly +# to mildly superlinearly with cell count, so a few-fold safety margin above +# the linear expectation comfortably covers normal variance while still +# catching a quadratic-or-worse regression (e.g. the do.call()/SCTransform +# hang this test suite failed to catch: see workflow/scripts/downsample_clusters.R). +SAFETY_MULTIPLIER = 4 +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]) + + t0 = time.monotonic() + result = run_command(cmd, root, timeout=1800) + elapsed = time.monotonic() - t0 + + assert result.returncode == 0, combined_output(result) + assert target.exists(), f"missing rule output: {target}" + _assert_basic_stability_table(target) + return elapsed + + +def test_downsample_clusters_runtime_scales_reasonably_with_ncells(tmp_path, pytestconfig): + if not pytestconfig.getoption("--run-downsample-scaling"): + pytest.skip("use --run-downsample-scaling to execute the downsample_clusters scaling check") + + 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) + + ratio = medium_elapsed / small_elapsed + assert ratio <= MAX_ALLOWED_RATIO, ( + f"downsample_clusters took {ratio:.1f}x longer on {MEDIUM_NCELLS} cells than on " + f"{SMALL_NCELLS} cells ({small_elapsed:.1f}s -> {medium_elapsed:.1f}s), exceeding the " + f"{MAX_ALLOWED_RATIO:.1f}x budget ({SAFETY_MULTIPLIER}x the {CELL_RATIO:.1f}x cell-count " + "ratio). This suggests a non-linear performance regression in the R script." + ) diff --git a/tests/downsampling/test_downsampling_reference_outputs.py b/tests/downsampling/test_downsampling_reference_outputs.py new file mode 100644 index 0000000..5551594 --- /dev/null +++ b/tests/downsampling/test_downsampling_reference_outputs.py @@ -0,0 +1,14 @@ +from utils import ( + TEST_SAMPLES, + assert_stability_table, + reference_output_for_sample, + repo_root, +) + + +def test_reference_outputs_exist_and_have_expected_schema(): + root = repo_root() + for sample in TEST_SAMPLES: + reference = root / reference_output_for_sample(sample) + assert reference.exists(), f"missing reference output: {reference}" + assert_stability_table(reference) diff --git a/tests/downsampling/test_downsampling_snakemake_dry_run.py b/tests/downsampling/test_downsampling_snakemake_dry_run.py new file mode 100644 index 0000000..f0acf7a --- /dev/null +++ b/tests/downsampling/test_downsampling_snakemake_dry_run.py @@ -0,0 +1,29 @@ +from utils import ( + TEST_RESULTS_DIR, + TEST_SAMPLES, + TEST_SEURAT_DIR, + base_snakemake_cmd, + combined_output, + repo_root, + run_command, +) + + +def test_testdata_fixtures_exist(): + root = repo_root() + for sample in TEST_SAMPLES: + assert (root / TEST_SEURAT_DIR / f"{sample}.rds").exists() + + +def test_testdata_dry_run_builds_expected_dag(): + root = repo_root() + cmd = base_snakemake_cmd() + ["-np"] + + result = run_command(cmd, root) + output = combined_output(result) + + assert result.returncode == 0, output + assert "downsample_clusters" in output + for sample in TEST_SAMPLES: + expected = TEST_RESULTS_DIR / f"{sample}_clusterdownsampling.tsv" + assert expected.as_posix() in output diff --git a/tests/downsampling/test_downsampling_workflow_run.py b/tests/downsampling/test_downsampling_workflow_run.py new file mode 100644 index 0000000..be5a07b --- /dev/null +++ b/tests/downsampling/test_downsampling_workflow_run.py @@ -0,0 +1,40 @@ +import shutil + +import pytest + +from utils import ( + TEST_RESULTS_DIR, + TEST_SAMPLES, + assert_stability_table, + assert_workflow_outputs_match_references, + combined_output, + expected_output_for_sample, + repo_root, + run_command, + workflow_run_cmd, +) + + +def test_testdata_workflow_run(pytestconfig): + if not pytestconfig.getoption("--run-downsample-workflow"): + pytest.skip("use --run-downsample-workflow to execute the full downsampling workflow") + + root = repo_root() + results_dir = root / TEST_RESULTS_DIR + if results_dir.exists(): + shutil.rmtree(results_dir) + + cmd = workflow_run_cmd() + conda_prefix = pytestconfig.getoption("--snakemake-conda-prefix") + if conda_prefix: + cmd.extend(["--conda-prefix", conda_prefix]) + + result = run_command(cmd, root, timeout=3600) + assert result.returncode == 0, combined_output(result) + + for sample in TEST_SAMPLES: + output = root / expected_output_for_sample(sample) + assert output.exists(), f"missing workflow output: {output}" + assert_stability_table(output) + + assert_workflow_outputs_match_references(root) diff --git a/tests/downsampling/update_reference_outputs.py b/tests/downsampling/update_reference_outputs.py new file mode 100644 index 0000000..5576fa0 --- /dev/null +++ b/tests/downsampling/update_reference_outputs.py @@ -0,0 +1,22 @@ +import shutil +from pathlib import Path + +from utils import TEST_SAMPLES, expected_output_for_sample, reference_output_for_sample, repo_root + + +def main(): + root = repo_root() + copied = 0 + for sample in TEST_SAMPLES: + source = root / expected_output_for_sample(sample) + if not source.exists(): + raise SystemExit(f"missing workflow output: {source}") + dest = root / reference_output_for_sample(sample) + dest.parent.mkdir(parents=True, exist_ok=True) + shutil.copy2(source, dest) + copied += 1 + print(f"Copied {copied} reference output files") + + +if __name__ == "__main__": + main() diff --git a/tests/downsampling/utils.py b/tests/downsampling/utils.py new file mode 100644 index 0000000..ce3223b --- /dev/null +++ b/tests/downsampling/utils.py @@ -0,0 +1,157 @@ +import csv +import math +import os +import shutil +import subprocess +from pathlib import Path + + +TEST_SEURAT_DIR = Path("testdata/downsampling/seurat_objects") +TEST_RESULTS_DIR = Path("testdata/downsampling/results") +REFERENCE_ROOT = Path("tests/downsampling/reference_outputs") +REFERENCE_RESULTS_DIR = Path("testdata/results") +TEST_SAMPLES = [ + "filtered_seurat_emptydrops_test", + "filtered_seurat_tenx_test", +] +TEST_CONFIG = [ + "workflow_mode=downsample_only", + f"downsampleSeuratObjectDir={TEST_SEURAT_DIR.as_posix()}", + f"downsampleResultsDir={TEST_RESULTS_DIR.as_posix()}", + "nDownsampleReplicates=2", + "downsampleRate=0.5", + "workflowSeed=12345", +] +EXPECTED_COLUMNS = ["clusterid", "max_jaccard", "bootstrap_number"] +EXPECTED_BOOTSTRAPS = {1, 2} + + +def repo_root(): + return Path(__file__).resolve().parents[2] + + +def snakemake_executable(): + snakemake = shutil.which("snakemake") + assert snakemake is not None, "snakemake is not available on PATH" + return snakemake + + +def run_command(cmd, cwd, timeout=300): + return subprocess.run( + cmd, + cwd=cwd, + text=True, + capture_output=True, + check=False, + timeout=timeout, + ) + + +def combined_output(result): + return result.stdout + result.stderr + + +def base_snakemake_cmd(): + return [ + snakemake_executable(), + "--snakefile", + "workflow/Snakefile", + "--configfile", + "config/config.yaml", + "--config", + *TEST_CONFIG, + "--profile", + "none", + "--workflow-profile", + "none", + "--executor", + "local", + ] + + +def workflow_run_cmd(): + return base_snakemake_cmd() + [ + "--cores", + "1", + "--jobs", + "1", + "--latency-wait", + "30", + "--rerun-incomplete", + "--use-conda", + ] + + +def expected_output_for_sample(sample, results_dir=TEST_RESULTS_DIR): + return Path(results_dir) / f"{sample}_clusterdownsampling.tsv" + + +def reference_output_for_sample(sample): + return REFERENCE_ROOT / REFERENCE_RESULTS_DIR / f"{sample}_clusterdownsampling.tsv" + + +def read_tsv(path): + with path.open(newline="") as handle: + reader = csv.DictReader(handle, delimiter="\t") + rows = list(reader) + return reader.fieldnames or [], rows + + +def assert_stability_table(path): + 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 + for row in rows: + assert row["clusterid"] != "" + max_jaccard = float(row["max_jaccard"]) + assert 0 <= max_jaccard <= 1 + + +def rows_by_bootstrap_cluster(path, rows): + indexed = {} + for row in rows: + key = (int(row["bootstrap_number"]), row["clusterid"]) + assert key not in indexed, f"{path}: duplicate bootstrap/cluster row: {key}" + indexed[key] = row + return indexed + + +def assert_stability_table_matches_reference(current_path, reference_path): + current_columns, current_rows = read_tsv(current_path) + reference_columns, reference_rows = read_tsv(reference_path) + assert current_columns == reference_columns == EXPECTED_COLUMNS + + current_index = rows_by_bootstrap_cluster(current_path, current_rows) + reference_index = rows_by_bootstrap_cluster(reference_path, reference_rows) + assert set(current_index) == set(reference_index), ( + f"{current_path}: bootstrap/cluster rows differ from reference" + ) + + rel_tol = float(os.getenv("DOWNSAMPLE_REFERENCE_REL_TOLERANCE", "1e-6")) + abs_tol = float(os.getenv("DOWNSAMPLE_REFERENCE_ABS_TOLERANCE", "1e-8")) + for key, reference_row in reference_index.items(): + current_row = current_index[key] + assert current_row["bootstrap_number"] == reference_row["bootstrap_number"] + assert current_row["clusterid"] == reference_row["clusterid"] + current_jaccard = float(current_row["max_jaccard"]) + reference_jaccard = float(reference_row["max_jaccard"]) + assert math.isclose( + current_jaccard, + reference_jaccard, + rel_tol=rel_tol, + abs_tol=abs_tol, + ), ( + f"{current_path}: max_jaccard differs for bootstrap/cluster {key}: " + f"{current_jaccard} != {reference_jaccard}" + ) + + +def assert_workflow_outputs_match_references(root): + for sample in TEST_SAMPLES: + current = root / expected_output_for_sample(sample) + reference = root / reference_output_for_sample(sample) + assert reference.exists(), f"missing reference output: {reference}" + assert_stability_table_matches_reference(current, reference) diff --git a/tests/run_test_workflow.sh b/tests/run_test_workflow.sh index f3136a4..ed11e1f 100755 --- a/tests/run_test_workflow.sh +++ b/tests/run_test_workflow.sh @@ -50,7 +50,7 @@ fi common_args=( --snakefile workflow/Snakefile --configfile config/config.yaml - --config sampleTable=testdata/samplesheet_test.tsv resultsDir=testdata/results + --config sampleTable=testdata/samplesheet_test.tsv resultsDir=testdata/results downsampleResultsDir=testdata/results/downsampling --use-conda --workflow-profile profiles/slurm --profile cannon diff --git a/tests/test_conda_container_validation.py b/tests/test_conda_container_validation.py index f3893bb..9084699 100644 --- a/tests/test_conda_container_validation.py +++ b/tests/test_conda_container_validation.py @@ -14,6 +14,7 @@ EXPECTED_ENV_FILES = { "cellbender.yml", "doubletfinder.yml", + "downsample_clusters.yml", "emptydrops.yml", "posthocfilter.yml", "scdblfinder.yml", @@ -26,6 +27,7 @@ R_IMPORTS_BY_ENV = { "cellbender.yml": ["Seurat", "tidyverse", "bluster"], "doubletfinder.yml": ["Seurat", "tidyverse", "remotes", "fields", "Matrix", "KernSmooth", "ROCR", "igraph", "glmGamPoi", "bluster"], + "downsample_clusters.yml": ["Seurat", "tidyverse", "glmGamPoi"], "emptydrops.yml": ["Seurat", "tidyverse", "DropletUtils", "scater", "glmGamPoi", "bluster", "Matrix", "R.utils"], "posthocfilter.yml": ["Seurat", "tidyverse", "glmGamPoi", "scater", "bluster"], "scdblfinder.yml": ["Seurat", "tidyverse", "scDblFinder", "glmGamPoi", "bluster"], diff --git a/tests/test_config_validation.py b/tests/test_config_validation.py index bef5076..dcf25e6 100644 --- a/tests/test_config_validation.py +++ b/tests/test_config_validation.py @@ -25,6 +25,7 @@ ALLOWED_DECON_METHODS = {"soupx", "cellbender_fromraw"} ALLOWED_DOUBLET_METHODS = {"doubletfinder", "scdblfinder"} ALLOWED_POSTHOC_METHODS = {"threshold", "mad"} +ALLOWED_WORKFLOW_MODES = {"preprocess", "preprocess_and_downsample", "downsample_only"} def load_default_config(repo_root): @@ -92,7 +93,12 @@ def test_default_config_has_required_keys_and_valid_values(): assert REQUIRED_CONFIG_KEYS <= set(config) assert config["conda-channel-priority"] == "strict" assert isinstance(config["sampleTable"], str) and config["sampleTable"] + assert config.get("workflow_mode") in ALLOWED_WORKFLOW_MODES assert isinstance(config["workflow_seed"], int) and not isinstance(config["workflow_seed"], bool) + assert isinstance(config["nDownsampleReplicates"], int) and config["nDownsampleReplicates"] > 0 + assert isinstance(config["downsampleRate"], (int, float)) and 0 < config["downsampleRate"] <= 1 + assert isinstance(config["downsampleSeuratObjectDir"], str) and config["downsampleSeuratObjectDir"] + assert isinstance(config["downsampleResultsDir"], str) and config["downsampleResultsDir"] for key, allowed_values in [ ("emptydrop_removal_methods", ALLOWED_EMPTYDROP_METHODS), @@ -120,6 +126,8 @@ def test_default_config_has_required_keys_and_valid_values(): (lambda cfg: cfg.update({"min_nfeature": 0}), "min_nfeature must be a positive integer", True), (lambda cfg: cfg.update({"max_mtdna": 101}), "max_mtdna must be a number between 0 and 100", True), (lambda cfg: cfg.update({"resultsDir": ""}), "resultsDir must be a non-empty string", False), + (lambda cfg: cfg.update({"workflow_mode": "bad_mode"}), "workflow_mode must be one of", True), + (lambda cfg: cfg.update({"workflow_mode": "downsample_only", "downsampleRate": 1.5}), "downsampleRate must be > 0 and <= 1", True), ], ) def test_invalid_config_fails_early_with_clear_message(tmp_path, mutate, expected_message, override_results_dir): @@ -180,3 +188,38 @@ def test_emptydrop_method_variant_excludes_emptydrops_specific_targets(tmp_path) assert str(results_dir / "seurat_filtered" / "filtered_seurat_tenx_test_markergenes.csv") in output assert str(results_dir / "soupx" / "seurat_soupx_tenx_test_markergenes.csv") in output assert "emptydrops" not in output + + +def test_downsample_only_mode_builds_dag_from_external_seurat_objects(tmp_path): + repo_root = Path(__file__).resolve().parents[1] + config = { + "conda-channel-priority": "strict", + "workflow_mode": "downsample_only", + "downsampleSeuratObjectDir": "testdata/downsampling/seurat_objects", + "downsampleResultsDir": (tmp_path / "downsampling").as_posix(), + "nDownsampleReplicates": 2, + "downsampleRate": 0.5, + "workflowSeed": 12345, + } + config_file = write_config(tmp_path / "downsample_only.yaml", config) + + result = run_snakemake(repo_root, config_file, None, "-np") + output = combined_output(result) + + assert result.returncode == 0, output + assert "downsample_clusters" in output + assert str(tmp_path / "downsampling" / "filtered_seurat_tenx_test_clusterdownsampling.tsv") in output + assert "tenx2seuratrds" not in output + + +def test_preprocess_and_downsample_mode_adds_downsample_targets(tmp_path): + output, results_dir = run_variant_dry_run( + tmp_path, + workflow_mode="preprocess_and_downsample", + downsampleTargets=["filtered_seurat_tenx_test"], + downsampleResultsDir=(tmp_path / "downsampling").as_posix(), + ) + + assert "tenx2seuratrds" in output + assert "downsample_clusters" in output + assert str(tmp_path / "downsampling" / "filtered_seurat_tenx_test_clusterdownsampling.tsv") in output diff --git a/tests/test_sample_rule_output_files.txt b/tests/test_sample_rule_output_files.txt index 460fad6..3a85591 100644 --- a/tests/test_sample_rule_output_files.txt +++ b/tests/test_sample_rule_output_files.txt @@ -8,6 +8,29 @@ testdata/results/doubletfinder/seurat_doubletfinder_soupx_emptydrops_test.rds testdata/results/doubletfinder/seurat_doubletfinder_soupx_emptydrops_test_markergenes.csv testdata/results/doubletfinder/seurat_doubletfinder_soupx_tenx_test.rds testdata/results/doubletfinder/seurat_doubletfinder_soupx_tenx_test_markergenes.csv +testdata/results/downsampling/filtered_seurat_emptydrops_test_clusterdownsampling.tsv +testdata/results/downsampling/filtered_seurat_tenx_test_clusterdownsampling.tsv +testdata/results/downsampling/seurat_cellbender_fromraw_test_clusterdownsampling.tsv +testdata/results/downsampling/seurat_doubletfinder_cellbender_fromraw_test_clusterdownsampling.tsv +testdata/results/downsampling/seurat_doubletfinder_soupx_emptydrops_test_clusterdownsampling.tsv +testdata/results/downsampling/seurat_doubletfinder_soupx_tenx_test_clusterdownsampling.tsv +testdata/results/downsampling/seurat_posthocfilt_mad_doubletfinder_cellbender_fromraw_test_clusterdownsampling.tsv +testdata/results/downsampling/seurat_posthocfilt_mad_doubletfinder_soupx_emptydrops_test_clusterdownsampling.tsv +testdata/results/downsampling/seurat_posthocfilt_mad_doubletfinder_soupx_tenx_test_clusterdownsampling.tsv +testdata/results/downsampling/seurat_posthocfilt_mad_scdblfinder_cellbender_fromraw_test_clusterdownsampling.tsv +testdata/results/downsampling/seurat_posthocfilt_mad_scdblfinder_soupx_emptydrops_test_clusterdownsampling.tsv +testdata/results/downsampling/seurat_posthocfilt_mad_scdblfinder_soupx_tenx_test_clusterdownsampling.tsv +testdata/results/downsampling/seurat_posthocfilt_threshold_doubletfinder_cellbender_fromraw_test_clusterdownsampling.tsv +testdata/results/downsampling/seurat_posthocfilt_threshold_doubletfinder_soupx_emptydrops_test_clusterdownsampling.tsv +testdata/results/downsampling/seurat_posthocfilt_threshold_doubletfinder_soupx_tenx_test_clusterdownsampling.tsv +testdata/results/downsampling/seurat_posthocfilt_threshold_scdblfinder_cellbender_fromraw_test_clusterdownsampling.tsv +testdata/results/downsampling/seurat_posthocfilt_threshold_scdblfinder_soupx_emptydrops_test_clusterdownsampling.tsv +testdata/results/downsampling/seurat_posthocfilt_threshold_scdblfinder_soupx_tenx_test_clusterdownsampling.tsv +testdata/results/downsampling/seurat_scdblfinder_cellbender_fromraw_test_clusterdownsampling.tsv +testdata/results/downsampling/seurat_scdblfinder_soupx_emptydrops_test_clusterdownsampling.tsv +testdata/results/downsampling/seurat_scdblfinder_soupx_tenx_test_clusterdownsampling.tsv +testdata/results/downsampling/seurat_soupx_emptydrops_test_clusterdownsampling.tsv +testdata/results/downsampling/seurat_soupx_tenx_test_clusterdownsampling.tsv testdata/results/emptydrops/filtered_seurat_emptydrops_test.rds testdata/results/emptydrops/filtered_seurat_emptydrops_test_markergenes.csv testdata/results/emptydrops/test_emptydrops_filtered_matrix/barcodes.tsv.gz diff --git a/tests/test_snakemake_dry_run.py b/tests/test_snakemake_dry_run.py index c03186c..3eb8ebe 100644 --- a/tests/test_snakemake_dry_run.py +++ b/tests/test_snakemake_dry_run.py @@ -8,6 +8,7 @@ SAMPLE_ID = "test" TEST_SAMPLE_SHEET = Path("testdata/samplesheet_test.tsv") RESULTS_DIR = "testdata/results" +DOWNSAMPLE_RESULTS_DIR = f"{RESULTS_DIR}/downsampling" EXPECTED_OUTPUTS = Path(__file__).with_name("test_sample_rule_output_files.txt") @@ -28,6 +29,7 @@ def run_snakemake(repo_root, *extra_args): "--config", f"sampleTable={TEST_SAMPLE_SHEET.as_posix()}", f"resultsDir={RESULTS_DIR}", + f"downsampleResultsDir={DOWNSAMPLE_RESULTS_DIR}", ] return subprocess.run( cmd, diff --git a/tests/test_snakemake_lint.py b/tests/test_snakemake_lint.py index ba252d8..d8fbb99 100644 --- a/tests/test_snakemake_lint.py +++ b/tests/test_snakemake_lint.py @@ -21,6 +21,7 @@ def test_snakemake_lint_has_no_findings(): "--config", f"sampleTable={TEST_SAMPLE_SHEET.as_posix()}", "resultsDir=testdata/results", + "downsampleResultsDir=testdata/results/downsampling", ] result = subprocess.run( cmd, diff --git a/workflow/Snakefile b/workflow/Snakefile index 44b8d5e..bf22523 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -1,13 +1,28 @@ +import re +from pathlib import Path + import pandas as pd import yaml include: "rules/common.smk" validated_config = validate_workflow_config(config) -sampleinfo = pd.read_table(config["sampleTable"], dtype={"sampleid": str}) -sampleinfo = validate_sample_sheet(sampleinfo, config["sampleTable"]) +WORKFLOW_MODE = validated_config["workflow_mode"] +RUN_PREPROCESS = WORKFLOW_MODE in {"preprocess", "preprocess_and_downsample"} +RUN_DOWNSAMPLE = WORKFLOW_MODE in {"preprocess_and_downsample", "downsample_only"} +RESULTS_DIR = config.get("resultsDir", "results").rstrip("/") +WORKFLOW_SEED = int(config.get("workflow_seed", config.get("workflowSeed", 12345))) +config.setdefault("min_nfeature", 1) +config.setdefault("min_ncount", 1) +config.setdefault("max_mtdna", 100) -SAMPLES = list(sampleinfo.loc[:, "sampleid"]) +if RUN_PREPROCESS: + sampleinfo = pd.read_table(config["sampleTable"], dtype={"sampleid": str}) + sampleinfo = validate_sample_sheet(sampleinfo, config["sampleTable"]) + SAMPLES = list(sampleinfo.loc[:, "sampleid"]) +else: + sampleinfo = pd.DataFrame({"sampleid": [], "tenx_datadir": []}) + SAMPLES = [] with open("profiles/slurm/config.yaml") as f: slurm_profile_config = yaml.safe_load(f) @@ -16,89 +31,121 @@ DECON_METHODS = validated_config["decon_methods"] DOUBLET_METHODS = validated_config["doublet_methods"] EMPTYDROP_METHODS = validated_config["emptydrop_methods"] POSTHOC_METHODS = validated_config["posthoc_methods"] -RESULTS_DIR = config.get("resultsDir", "results").rstrip("/") -WORKFLOW_SEED = int(config["workflow_seed"]) +DOWNSAMPLE_RESULTS_DIR = config.get("downsampleResultsDir", f"{RESULTS_DIR}/downsampling").rstrip("/") +DOWNSAMPLE_SEURAT_OBJECT_DIR = config.get( + "downsampleSeuratObjectDir", + config.get("seuratObjectDir", "seurat_objects"), +).rstrip("/") -TENX_TARGETS = expand( - f"{RESULTS_DIR}/seurat_filtered/filtered_seurat_tenx_{{sample}}_markergenes.csv", - sample=SAMPLES, +TENX_PREFIXES = ( + expand("seurat_filtered/filtered_seurat_tenx_{sample}", sample=SAMPLES) + if RUN_PREPROCESS + else [] ) -EMPTYDROPS_TARGETS = ( - expand( - f"{RESULTS_DIR}/emptydrops/filtered_seurat_emptydrops_{{sample}}_markergenes.csv", - sample=SAMPLES, - ) - if "emptydrops" in EMPTYDROP_METHODS +EMPTYDROPS_PREFIXES = ( + expand("emptydrops/filtered_seurat_emptydrops_{sample}", sample=SAMPLES) + if RUN_PREPROCESS and "emptydrops" in EMPTYDROP_METHODS else [] ) -SOUPX_TARGETS = ( +SOUPX_PREFIXES = ( expand( - f"{RESULTS_DIR}/soupx/seurat_soupx_{{empty_method}}_{{sample}}_markergenes.csv", + "soupx/seurat_soupx_{empty_method}_{sample}", sample=SAMPLES, empty_method=EMPTYDROP_METHODS, ) - if "soupx" in DECON_METHODS + if RUN_PREPROCESS and "soupx" in DECON_METHODS else [] ) -CELLBENDER_TARGETS = ( - expand( - f"{RESULTS_DIR}/cellbender_fromraw/seurat_cellbender_fromraw_{{sample}}_markergenes.csv", - sample=SAMPLES, - ) - if "cellbender_fromraw" in DECON_METHODS +CELLBENDER_PREFIXES = ( + expand("cellbender_fromraw/seurat_cellbender_fromraw_{sample}", sample=SAMPLES) + if RUN_PREPROCESS and "cellbender_fromraw" in DECON_METHODS else [] ) -SOUPX_DOUBLET_TARGETS = ( +SOUPX_DOUBLET_PREFIXES = ( expand( - f"{RESULTS_DIR}/{{doublet_method}}/seurat_{{doublet_method}}_soupx_{{empty_method}}_{{sample}}_markergenes.csv", + "{doublet_method}/seurat_{doublet_method}_soupx_{empty_method}_{sample}", sample=SAMPLES, doublet_method=DOUBLET_METHODS, empty_method=EMPTYDROP_METHODS, ) - if "soupx" in DECON_METHODS + if RUN_PREPROCESS and "soupx" in DECON_METHODS else [] ) -CELLBENDER_DOUBLET_TARGETS = ( +CELLBENDER_DOUBLET_PREFIXES = ( expand( - f"{RESULTS_DIR}/{{doublet_method}}/seurat_{{doublet_method}}_cellbender_fromraw_{{sample}}_markergenes.csv", + "{doublet_method}/seurat_{doublet_method}_cellbender_fromraw_{sample}", sample=SAMPLES, doublet_method=DOUBLET_METHODS, ) - if "cellbender_fromraw" in DECON_METHODS + if RUN_PREPROCESS and "cellbender_fromraw" in DECON_METHODS else [] ) -SOUPX_POSTHOC_TARGETS = ( +SOUPX_POSTHOC_PREFIXES = ( expand( - f"{RESULTS_DIR}/posthocfilter/seurat_posthocfilt_{{posthoc_method}}_{{doublet_method}}_soupx_{{empty_method}}_{{sample}}_markergenes.csv", + "posthocfilter/seurat_posthocfilt_{posthoc_method}_{doublet_method}_soupx_{empty_method}_{sample}", sample=SAMPLES, posthoc_method=POSTHOC_METHODS, doublet_method=DOUBLET_METHODS, empty_method=EMPTYDROP_METHODS, ) - if "soupx" in DECON_METHODS + if RUN_PREPROCESS and "soupx" in DECON_METHODS else [] ) -CELLBENDER_POSTHOC_TARGETS = ( +CELLBENDER_POSTHOC_PREFIXES = ( expand( - f"{RESULTS_DIR}/posthocfilter/seurat_posthocfilt_{{posthoc_method}}_{{doublet_method}}_cellbender_fromraw_{{sample}}_markergenes.csv", + "posthocfilter/seurat_posthocfilt_{posthoc_method}_{doublet_method}_cellbender_fromraw_{sample}", sample=SAMPLES, posthoc_method=POSTHOC_METHODS, doublet_method=DOUBLET_METHODS, ) - if "cellbender_fromraw" in DECON_METHODS + if RUN_PREPROCESS and "cellbender_fromraw" in DECON_METHODS else [] ) -ALL_TARGETS = ( - TENX_TARGETS - + EMPTYDROPS_TARGETS - + SOUPX_TARGETS - + CELLBENDER_TARGETS - + SOUPX_DOUBLET_TARGETS - + CELLBENDER_DOUBLET_TARGETS - + SOUPX_POSTHOC_TARGETS - + CELLBENDER_POSTHOC_TARGETS +PREPROCESS_PREFIXES = ( + TENX_PREFIXES + + EMPTYDROPS_PREFIXES + + SOUPX_PREFIXES + + CELLBENDER_PREFIXES + + SOUPX_DOUBLET_PREFIXES + + CELLBENDER_DOUBLET_PREFIXES + + SOUPX_POSTHOC_PREFIXES + + CELLBENDER_POSTHOC_PREFIXES +) +PREPROCESS_MARKER_TARGETS = marker_targets(PREPROCESS_PREFIXES) +PREPROCESS_SEURAT_TARGETS = rds_targets(PREPROCESS_PREFIXES) + + +if RUN_DOWNSAMPLE: + if WORKFLOW_MODE == "downsample_only": + DOWNSAMPLE_INPUTS_BY_TARGET = downsample_inputs_from_external_dir() + else: + DOWNSAMPLE_INPUTS_BY_TARGET = downsample_inputs_from_preprocess_outputs() + DOWNSAMPLE_INPUTS_BY_TARGET = select_downsample_inputs(DOWNSAMPLE_INPUTS_BY_TARGET) + if not DOWNSAMPLE_INPUTS_BY_TARGET: + raise ValueError( + "No Seurat RDS inputs were found for downsampling. " + "For downsample_only mode, set downsampleSeuratObjectDir to a directory containing *.rds files." + ) +else: + DOWNSAMPLE_INPUTS_BY_TARGET = {} + +DOWNSAMPLE_TARGETS = [ + f"{DOWNSAMPLE_RESULTS_DIR}/{target}_clusterdownsampling.tsv" + for target in sorted(DOWNSAMPLE_INPUTS_BY_TARGET) +] +DOWNSAMPLE_TARGET_REGEX = ( + "|".join(re.escape(target) for target in sorted(DOWNSAMPLE_INPUTS_BY_TARGET)) + if DOWNSAMPLE_INPUTS_BY_TARGET + else r"[^/]+" ) +ALL_TARGETS = [] +if RUN_PREPROCESS: + ALL_TARGETS.extend(PREPROCESS_MARKER_TARGETS) +if RUN_DOWNSAMPLE: + ALL_TARGETS.extend(DOWNSAMPLE_TARGETS) + + default_target: True @@ -119,3 +166,4 @@ include: "rules/soupx_emptydrops.smk" include: "rules/posthocfilter_mad.smk" include: "rules/posthocfilter_threshold.smk" include: "rules/markers.smk" +include: "rules/downsample_clusters.smk" diff --git a/workflow/Snakefile_downsample b/workflow/Snakefile_downsample new file mode 100644 index 0000000..c403095 --- /dev/null +++ b/workflow/Snakefile_downsample @@ -0,0 +1,57 @@ +import re + +RESULTS_DIR = config.get("resultsDir", "results").rstrip("/") +DOWNSAMPLE_RESULTS_DIR = config.get("downsampleResultsDir", RESULTS_DIR).rstrip("/") +DOWNSAMPLE_SEURAT_OBJECT_DIR = config.get( + "downsampleSeuratObjectDir", + config.get("seuratObjectDir", "seurat_objects"), +).rstrip("/") + + +def select_downsample_inputs(inputs_by_target): + requested = config.get("downsampleTargets", "all") + if requested == "all" or requested == ["all"]: + return dict(inputs_by_target) + + missing = sorted(set(requested) - set(inputs_by_target)) + if missing: + raise ValueError( + "downsampleTargets contains unknown target(s): " + + ", ".join(missing) + + "; available targets are: " + + ", ".join(sorted(inputs_by_target)) + ) + return {target: inputs_by_target[target] for target in requested} + + +DOWNSAMPLE_INPUTS_BY_TARGET = { + target: f"{DOWNSAMPLE_SEURAT_OBJECT_DIR}/{target}.rds" + for target in sorted( + glob_wildcards(f"{DOWNSAMPLE_SEURAT_OBJECT_DIR}/{{target}}.rds").target + ) +} +DOWNSAMPLE_INPUTS_BY_TARGET = select_downsample_inputs(DOWNSAMPLE_INPUTS_BY_TARGET) +if not DOWNSAMPLE_INPUTS_BY_TARGET: + raise ValueError( + "No Seurat RDS inputs were found for downsampling. " + "Set downsampleSeuratObjectDir to a directory containing *.rds files." + ) + +DOWNSAMPLE_TARGETS = [ + f"{DOWNSAMPLE_RESULTS_DIR}/{target}_clusterdownsampling.tsv" + for target in sorted(DOWNSAMPLE_INPUTS_BY_TARGET) +] +DOWNSAMPLE_TARGET_REGEX = "|".join( + re.escape(target) for target in sorted(DOWNSAMPLE_INPUTS_BY_TARGET) +) + + +default_target: True + + +rule all: + input: + DOWNSAMPLE_TARGETS + + +include: "rules/downsample_clusters.smk" diff --git a/workflow/envs/downsample_clusters.linux-64.pin.txt b/workflow/envs/downsample_clusters.linux-64.pin.txt new file mode 100644 index 0000000..9f3a361 --- /dev/null +++ b/workflow/envs/downsample_clusters.linux-64.pin.txt @@ -0,0 +1,315 @@ +# This file may be used to create an environment using: +# $ conda create --name --file +# platform: linux-64 +# created-by: conda 25.11.0 +@EXPLICIT +https://conda.anaconda.org/conda-forge/noarch/_r-mutex-1.0.1-anacondar_1.tar.bz2 +https://conda.anaconda.org/conda-forge/noarch/ca-certificates-2026.5.20-hbd8a1cb_0.conda +https://conda.anaconda.org/conda-forge/noarch/font-ttf-dejavu-sans-mono-2.37-hab24e00_0.tar.bz2 +https://conda.anaconda.org/conda-forge/noarch/font-ttf-inconsolata-3.000-h77eed37_0.tar.bz2 +https://conda.anaconda.org/conda-forge/noarch/font-ttf-source-code-pro-2.038-h77eed37_0.tar.bz2 +https://conda.anaconda.org/conda-forge/noarch/font-ttf-ubuntu-0.83-h77eed37_3.conda +https://conda.anaconda.org/conda-forge/noarch/kernel-headers_linux-64-4.18.0-he073ed8_9.conda +https://conda.anaconda.org/conda-forge/noarch/libgcc-devel_linux-64-15.2.0-hcc6f6b0_119.conda +https://conda.anaconda.org/conda-forge/linux-64/libgomp-15.2.0-he0feb66_19.conda +https://conda.anaconda.org/conda-forge/noarch/libstdcxx-devel_linux-64-15.2.0-hd446a21_119.conda +https://conda.anaconda.org/conda-forge/linux-64/libzlib-1.3.2-h25fd6f3_2.conda +https://conda.anaconda.org/conda-forge/linux-64/pandoc-3.10-ha770c72_0.conda +https://conda.anaconda.org/conda-forge/noarch/tzdata-2025c-hc9c84f9_1.conda +https://conda.anaconda.org/conda-forge/linux-64/_openmp_mutex-4.5-20_gnu.conda +https://conda.anaconda.org/conda-forge/noarch/fonts-conda-forge-1-hc364b38_1.conda +https://conda.anaconda.org/conda-forge/noarch/sysroot_linux-64-2.28-h4ee821c_9.conda +https://conda.anaconda.org/conda-forge/linux-64/zstd-1.5.7-hb78ec9c_6.conda +https://conda.anaconda.org/conda-forge/noarch/fonts-conda-ecosystem-1-0.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/ld_impl_linux-64-2.45.1-default_hbd61a6d_102.conda +https://conda.anaconda.org/conda-forge/linux-64/libgcc-15.2.0-he0feb66_19.conda +https://conda.anaconda.org/conda-forge/linux-64/binutils_impl_linux-64-2.45.1-default_hfdba357_102.conda +https://conda.anaconda.org/conda-forge/linux-64/bzip2-1.0.8-hda65f42_9.conda +https://conda.anaconda.org/conda-forge/linux-64/c-ares-1.34.6-hb03c661_0.conda +https://conda.anaconda.org/conda-forge/linux-64/fribidi-1.0.16-hb03c661_0.conda +https://conda.anaconda.org/conda-forge/linux-64/keyutils-1.6.3-hb9d3cd8_0.conda +https://conda.anaconda.org/conda-forge/linux-64/libdeflate-1.25-h17f619e_0.conda +https://conda.anaconda.org/conda-forge/linux-64/libexpat-2.8.1-hecca717_1.conda +https://conda.anaconda.org/conda-forge/linux-64/libffi-3.5.2-h3435931_0.conda +https://conda.anaconda.org/conda-forge/linux-64/libgcc-ng-15.2.0-h69a702a_19.conda +https://conda.anaconda.org/conda-forge/linux-64/libgfortran5-15.2.0-h68bc16d_19.conda +https://conda.anaconda.org/conda-forge/linux-64/libiconv-1.18-h3b78370_2.conda +https://conda.anaconda.org/conda-forge/linux-64/libjpeg-turbo-3.1.4.1-hb03c661_0.conda +https://conda.anaconda.org/conda-forge/linux-64/liblzma-5.8.3-hb03c661_0.conda +https://conda.anaconda.org/conda-forge/linux-64/libpng-1.6.58-h421ea60_0.conda +https://conda.anaconda.org/conda-forge/linux-64/libstdcxx-15.2.0-h934c35e_19.conda +https://conda.anaconda.org/conda-forge/linux-64/libuuid-2.42.1-h5347b49_0.conda +https://conda.anaconda.org/conda-forge/linux-64/libuv-1.52.1-h280c20c_0.conda +https://conda.anaconda.org/conda-forge/linux-64/libwebp-base-1.6.0-hd42ef1d_0.conda +https://conda.anaconda.org/conda-forge/linux-64/make-4.4.1-hb9d3cd8_2.conda +https://conda.anaconda.org/conda-forge/linux-64/ncurses-6.6-hdb14827_0.conda +https://conda.anaconda.org/conda-forge/linux-64/openssl-3.6.3-h35e630c_0.conda +https://conda.anaconda.org/conda-forge/linux-64/pthread-stubs-0.4-hb9d3cd8_1002.conda +https://conda.anaconda.org/conda-forge/linux-64/sed-4.10-h19d0853_0.conda +https://conda.anaconda.org/conda-forge/linux-64/tk-8.6.13-noxft_h366c992_103.conda +https://conda.anaconda.org/conda-forge/linux-64/xorg-libice-1.1.2-hb9d3cd8_0.conda +https://conda.anaconda.org/conda-forge/linux-64/xorg-libxau-1.0.12-hb03c661_1.conda +https://conda.anaconda.org/conda-forge/linux-64/xorg-libxdmcp-1.1.5-hb03c661_1.conda +https://conda.anaconda.org/conda-forge/linux-64/binutils_linux-64-2.45.1-default_h4852527_102.conda +https://conda.anaconda.org/conda-forge/linux-64/bwidget-1.10.1-ha770c72_1.conda +https://conda.anaconda.org/conda-forge/linux-64/graphite2-1.3.15-hecca717_0.conda +https://conda.anaconda.org/conda-forge/linux-64/icu-78.3-h33c6efd_0.conda +https://conda.anaconda.org/conda-forge/linux-64/lerc-4.1.0-hdb68285_0.conda +https://conda.anaconda.org/conda-forge/linux-64/libedit-3.1.20250104-pl5321h7949ede_0.conda +https://conda.anaconda.org/conda-forge/linux-64/libev-4.33-hd590300_2.conda +https://conda.anaconda.org/conda-forge/linux-64/libfreetype6-2.14.3-h73754d4_0.conda +https://conda.anaconda.org/conda-forge/linux-64/libgfortran-15.2.0-h69a702a_19.conda +https://conda.anaconda.org/conda-forge/linux-64/libsanitizer-15.2.0-h90f66d4_19.conda +https://conda.anaconda.org/conda-forge/linux-64/libssh2-1.11.1-hcf80075_0.conda +https://conda.anaconda.org/conda-forge/linux-64/libstdcxx-ng-15.2.0-hdf11a46_19.conda +https://conda.anaconda.org/conda-forge/linux-64/libxcb-1.17.0-h8a09558_0.conda +https://conda.anaconda.org/conda-forge/linux-64/pcre2-10.47-haa7fec5_0.conda +https://conda.anaconda.org/conda-forge/linux-64/pixman-0.46.4-h54a6638_1.conda +https://conda.anaconda.org/conda-forge/linux-64/readline-8.3-h853b02a_0.conda +https://conda.anaconda.org/conda-forge/linux-64/xorg-libsm-1.2.6-he73a12e_0.conda +https://conda.anaconda.org/conda-forge/linux-64/gcc_impl_linux-64-15.2.0-he0086c7_19.conda +https://conda.anaconda.org/conda-forge/linux-64/gmp-6.3.0-hac33072_2.conda +https://conda.anaconda.org/conda-forge/linux-64/krb5-1.22.2-ha1258a1_0.conda +https://conda.anaconda.org/conda-forge/linux-64/libfreetype-2.14.3-ha770c72_0.conda +https://conda.anaconda.org/conda-forge/linux-64/libglib-2.88.1-h0d30a3d_2.conda +https://conda.anaconda.org/conda-forge/linux-64/libnghttp2-1.68.1-h877daf1_0.conda +https://conda.anaconda.org/conda-forge/linux-64/libopenblas-0.3.33-pthreads_h94d23a6_0.conda +https://conda.anaconda.org/conda-forge/linux-64/libtiff-4.7.1-h9d88235_1.conda +https://conda.anaconda.org/conda-forge/linux-64/libxml2-16-2.15.3-hca6bf5a_0.conda +https://conda.anaconda.org/conda-forge/linux-64/xorg-libx11-1.8.13-he1eb515_0.conda +https://conda.anaconda.org/conda-forge/linux-64/fontconfig-2.18.1-h27c8c51_0.conda +https://conda.anaconda.org/conda-forge/linux-64/gcc_linux-64-15.2.0-h7be306e_26.conda +https://conda.anaconda.org/conda-forge/linux-64/gfortran_impl_linux-64-15.2.0-h281d09f_19.conda +https://conda.anaconda.org/conda-forge/linux-64/glpk-5.0-h445213a_0.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/gxx_impl_linux-64-15.2.0-hda75c37_19.conda +https://conda.anaconda.org/conda-forge/linux-64/libblas-3.11.0-8_h4a7cf45_openblas.conda +https://conda.anaconda.org/conda-forge/linux-64/libcurl-8.20.0-hcf29cc6_0.conda +https://conda.anaconda.org/conda-forge/linux-64/libxml2-2.15.3-h49c6c72_0.conda +https://conda.anaconda.org/conda-forge/linux-64/tktable-2.10-h5a7a40f_8.conda +https://conda.anaconda.org/conda-forge/linux-64/xorg-libxext-1.3.7-hb03c661_0.conda +https://conda.anaconda.org/conda-forge/linux-64/xorg-libxrender-0.9.12-hb9d3cd8_0.conda +https://conda.anaconda.org/conda-forge/linux-64/xorg-libxt-1.3.1-hb9d3cd8_0.conda +https://conda.anaconda.org/conda-forge/linux-64/cairo-1.18.4-he90730b_1.conda +https://conda.anaconda.org/conda-forge/linux-64/curl-8.20.0-hcf29cc6_0.conda +https://conda.anaconda.org/conda-forge/linux-64/gxx_linux-64-15.2.0-h92113df_26.conda +https://conda.anaconda.org/conda-forge/linux-64/libcblas-3.11.0-8_h0358290_openblas.conda +https://conda.anaconda.org/conda-forge/linux-64/liblapack-3.11.0-8_h47877c9_openblas.conda +https://conda.anaconda.org/conda-forge/linux-64/gsl-2.7-he838d99_0.tar.bz2 +https://conda.anaconda.org/conda-forge/linux-64/harfbuzz-14.2.1-h6083320_0.conda +https://conda.anaconda.org/conda-forge/linux-64/pango-1.56.4-hda50119_1.conda +https://conda.anaconda.org/conda-forge/linux-64/r-base-4.5.3-h502d0c9_2.conda +https://conda.anaconda.org/bioconda/noarch/bioconductor-assorthead-1.4.0-r45hdfd78af_0.conda +https://conda.anaconda.org/bioconda/linux-64/bioconductor-rhdf5lib-1.32.0-r45h01b2380_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-abind-1.4_8-r45hc72bb7e_1.conda +https://conda.anaconda.org/conda-forge/noarch/r-assertthat-0.2.1-r45hc72bb7e_6.conda +https://conda.anaconda.org/conda-forge/linux-64/r-backports-1.5.1-r45h54b55ab_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-base64enc-0.1_6-r45h54b55ab_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-bh-1.90.0_1-r45hc72bb7e_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-bit-4.6.0-r45h54b55ab_1.conda +https://conda.anaconda.org/conda-forge/linux-64/r-bitops-1.0_9-r45h54b55ab_1.conda +https://conda.anaconda.org/conda-forge/linux-64/r-cli-3.6.6-r45h3697838_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-clipr-0.8.1-r45hc72bb7e_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-cluster-2.1.8.2-r45heaba542_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-codetools-0.2_20-r45hc72bb7e_2.conda +https://conda.anaconda.org/conda-forge/linux-64/r-colorspace-2.1_2-r45h54b55ab_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-commonmark-2.0.0-r45h54b55ab_1.conda +https://conda.anaconda.org/conda-forge/noarch/r-cpp11-0.5.5-r45h785f33e_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-crayon-1.5.3-r45hc72bb7e_2.conda +https://conda.anaconda.org/conda-forge/linux-64/r-curl-7.1.0-r45h10955f1_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-data.table-1.17.8-r45h1c8cec4_1.conda +https://conda.anaconda.org/conda-forge/noarch/r-dbi-1.3.0-r45hc72bb7e_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-deldir-2.0_4-r45heaba542_2.conda +https://conda.anaconda.org/conda-forge/linux-64/r-digest-0.6.39-r45h3697838_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-dotcall64-1.2-r45heaba542_1.conda +https://conda.anaconda.org/conda-forge/noarch/r-evaluate-1.0.5-r45hc72bb7e_1.conda +https://conda.anaconda.org/conda-forge/linux-64/r-fansi-1.0.7-r45h54b55ab_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-farver-2.1.2-r45h3697838_2.conda +https://conda.anaconda.org/conda-forge/linux-64/r-fastmap-1.2.0-r45h3697838_2.conda +https://conda.anaconda.org/conda-forge/linux-64/r-fnn-1.1.4.1-r45h3697838_2.conda +https://conda.anaconda.org/conda-forge/linux-64/r-fs-2.1.0-r45h3697838_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-generics-0.1.4-r45hc72bb7e_1.conda +https://conda.anaconda.org/conda-forge/linux-64/r-glue-1.8.1-r45h54b55ab_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-goftest-1.2_3-r45h54b55ab_4.conda +https://conda.anaconda.org/conda-forge/linux-64/r-gtools-3.9.5-r45h54b55ab_2.conda +https://conda.anaconda.org/conda-forge/noarch/r-ica-1.0_3-r45hc72bb7e_4.conda +https://conda.anaconda.org/conda-forge/linux-64/r-jsonlite-2.0.0-r45h54b55ab_1.conda +https://conda.anaconda.org/conda-forge/linux-64/r-kernsmooth-2.23_26-r45ha0a88a1_1.conda +https://conda.anaconda.org/conda-forge/noarch/r-labeling-0.4.3-r45hc72bb7e_2.conda +https://conda.anaconda.org/conda-forge/linux-64/r-lattice-0.22_9-r45h54b55ab_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-listenv-0.10.1-r45hc72bb7e_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-magrittr-2.0.5-r45h54b55ab_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-mass-7.3_65-r45h54b55ab_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-matrixstats-1.5.0-r45h54b55ab_1.conda +https://conda.anaconda.org/conda-forge/linux-64/r-mime-0.13-r45h54b55ab_1.conda +https://conda.anaconda.org/conda-forge/noarch/r-otel-0.2.0-r45hc72bb7e_1.conda +https://conda.anaconda.org/conda-forge/linux-64/r-parallelly-1.47.0-r45h54b55ab_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-pbapply-1.7_4-r45hc72bb7e_1.conda +https://conda.anaconda.org/conda-forge/noarch/r-pkgconfig-2.0.3-r45hc72bb7e_5.conda +https://conda.anaconda.org/conda-forge/linux-64/r-png-0.1_9-r45haf2892b_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-polyclip-1.10_7-r45h3697838_1.conda +https://conda.anaconda.org/conda-forge/linux-64/r-ps-1.9.3-r45h54b55ab_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-r6-2.6.1-r45hc72bb7e_1.conda +https://conda.anaconda.org/conda-forge/linux-64/r-rann-2.6.2-r45h3697838_2.conda +https://conda.anaconda.org/conda-forge/linux-64/r-rappdirs-0.3.4-r45h54b55ab_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-rcolorbrewer-1.1_3-r45h785f33e_4.conda +https://conda.anaconda.org/conda-forge/linux-64/r-rcpp-1.1.1_1.1-r45h3697838_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-rcppprogress-0.4.2-r45hc72bb7e_5.conda +https://conda.anaconda.org/conda-forge/noarch/r-rematch-2.0.0-r45hc72bb7e_2.conda +https://conda.anaconda.org/conda-forge/linux-64/r-rlang-1.2.0-r45h3697838_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-rprojroot-2.1.1-r45hc72bb7e_1.conda +https://conda.anaconda.org/conda-forge/noarch/r-rstudioapi-0.19.0-r45hc72bb7e_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-s7-0.2.2-r45h54b55ab_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-sourcetools-0.1.7_2-r45h3697838_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-spatstat.utils-3.2_3-r45h54b55ab_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-stringi-1.8.7-r45h3d52c89_2.conda +https://conda.anaconda.org/conda-forge/linux-64/r-sys-3.4.3-r45h54b55ab_1.conda +https://conda.anaconda.org/conda-forge/noarch/r-tensor-1.5.1-r45hc72bb7e_1.conda +https://conda.anaconda.org/conda-forge/linux-64/r-utf8-1.2.6-r45h54b55ab_1.conda +https://conda.anaconda.org/conda-forge/linux-64/r-uuid-1.2_2-r45h54b55ab_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-viridislite-0.4.3-r45hc72bb7e_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-withr-3.0.2-r45hc72bb7e_1.conda +https://conda.anaconda.org/conda-forge/linux-64/r-xfun-0.57-r45h3697838_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-xtable-1.8_8-r45hc72bb7e_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-yaml-2.3.12-r45h54b55ab_0.conda +https://conda.anaconda.org/bioconda/noarch/bioconductor-biocgenerics-0.56.0-r45hdfd78af_2.conda +https://conda.anaconda.org/bioconda/noarch/bioconductor-matrixgenerics-1.22.0-r45hdfd78af_1.conda +https://conda.anaconda.org/bioconda/linux-64/bioconductor-rhdf5filters-1.22.0-r45ha27e39d_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-askpass-1.2.1-r45h54b55ab_1.conda +https://conda.anaconda.org/conda-forge/linux-64/r-bit64-4.8.2-r45h54b55ab_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-cachem-1.1.0-r45h54b55ab_2.conda +https://conda.anaconda.org/conda-forge/linux-64/r-catools-1.18.3-r45h3697838_1.conda +https://conda.anaconda.org/conda-forge/linux-64/r-ellipsis-0.3.3-r45h54b55ab_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-globals-0.19.1-r45hc72bb7e_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-here-1.0.2-r45hc72bb7e_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-hexbin-1.28.5-r45heaba542_1.conda +https://conda.anaconda.org/conda-forge/noarch/r-highr-0.12-r45hc72bb7e_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-isoband-0.3.0-r45h3697838_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-later-1.4.8-r45h3697838_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-lazyeval-0.2.3-r45h54b55ab_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-lifecycle-1.0.5-r45hc72bb7e_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-matrix-1.7_5-r45h0e4624f_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-munsell-0.5.1-r45hc72bb7e_2.conda +https://conda.anaconda.org/conda-forge/linux-64/r-nlme-3.1_169-r45heaba542_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-plyr-1.8.9-r45h3697838_3.conda +https://conda.anaconda.org/conda-forge/noarch/r-prettyunits-1.2.0-r45hc72bb7e_2.conda +https://conda.anaconda.org/conda-forge/linux-64/r-processx-3.9.0-r45h54b55ab_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-progressr-0.19.0-r45hc72bb7e_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-rcppannoy-0.0.23-r45h3697838_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-rcpparmadillo-15.2.7_1-r45h3704496_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-rcpphnsw-0.7.0-r45h3697838_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-rcpptoml-0.2.3-r45h3697838_1.conda +https://conda.anaconda.org/conda-forge/linux-64/r-rtsne-0.17-r45hf1899b2_3.conda +https://conda.anaconda.org/conda-forge/linux-64/r-sitmo-2.0.2-r45h3697838_4.conda +https://conda.anaconda.org/conda-forge/linux-64/r-sp-2.2_1-r45h54b55ab_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-spam-2.11_4-r45h2ddecb4_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-spatstat.univar-3.2_0-r45h54b55ab_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-timechange-0.4.0-r45h3697838_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-tinytex-0.59-r45hc72bb7e_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-tzdb-0.5.0-r45h3697838_2.conda +https://conda.anaconda.org/conda-forge/linux-64/r-xml2-1.5.2-r45he78afff_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-zoo-1.8_15-r45h54b55ab_0.conda +https://conda.anaconda.org/bioconda/linux-64/bioconductor-biobase-2.70.0-r45h01b2380_0.conda +https://conda.anaconda.org/bioconda/linux-64/bioconductor-rhdf5-2.54.1-r45h882589e_0.conda +https://conda.anaconda.org/bioconda/linux-64/bioconductor-s4vectors-0.48.0-r45h01b2380_1.conda +https://conda.anaconda.org/bioconda/linux-64/bioconductor-sparsematrixstats-1.22.0-r45ha27e39d_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-callr-3.7.6-r45hc72bb7e_2.conda +https://conda.anaconda.org/conda-forge/linux-64/r-dqrng-0.3.2-r45h3697838_2.conda +https://conda.anaconda.org/conda-forge/noarch/r-future-1.70.0-r45hc72bb7e_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-gplots-3.3.0-r45hc72bb7e_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-gtable-0.3.6-r45hc72bb7e_1.conda +https://conda.anaconda.org/conda-forge/linux-64/r-htmltools-0.5.9-r45h3697838_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-irlba-2.3.7-r45h0e4624f_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-knitr-1.51-r45hc72bb7e_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-lmtest-0.9_40-r45heaba542_4.conda +https://conda.anaconda.org/conda-forge/linux-64/r-lubridate-1.9.5-r45h54b55ab_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-memoise-2.0.1-r45hc72bb7e_4.conda +https://conda.anaconda.org/conda-forge/linux-64/r-openssl-2.4.2-r45h68c19f5_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-promises-1.5.0-r45hc72bb7e_1.conda +https://conda.anaconda.org/conda-forge/linux-64/r-rcppeigen-0.3.4.0.2-r45h3704496_1.conda +https://conda.anaconda.org/conda-forge/linux-64/r-reticulate-1.46.0-r45h3697838_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-scales-1.4.0-r45hc72bb7e_1.conda +https://conda.anaconda.org/conda-forge/noarch/r-spatstat.data-3.1_9-r45hc72bb7e_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-spatstat.sparse-3.2_0-r45h54b55ab_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-survival-3.8_6-r45h54b55ab_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-systemfonts-1.3.2-r45h74f4acd_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-vctrs-0.7.3-r45h3697838_0.conda +https://conda.anaconda.org/bioconda/linux-64/bioconductor-iranges-2.44.0-r45h01b2380_1.conda +https://conda.anaconda.org/conda-forge/noarch/r-blob-1.3.0-r45hc72bb7e_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-conflicted-1.2.0-r45h785f33e_3.conda +https://conda.anaconda.org/conda-forge/noarch/r-crosstalk-1.2.2-r45hc72bb7e_1.conda +https://conda.anaconda.org/conda-forge/noarch/r-fitdistrplus-1.2_6-r45hc72bb7e_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-fontawesome-0.5.3-r45hc72bb7e_1.conda +https://conda.anaconda.org/conda-forge/noarch/r-future.apply-1.20.2-r45hc72bb7e_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-ggplot2-4.0.3-r45h785f33e_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-gridextra-2.3-r45hc72bb7e_1007.conda +https://conda.anaconda.org/conda-forge/noarch/r-hms-1.1.4-r45hc72bb7e_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-httpuv-1.6.17-r45h6d565e7_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-httr-1.4.8-r45hc72bb7e_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-ids-1.0.1-r45hc72bb7e_5.conda +https://conda.anaconda.org/conda-forge/linux-64/r-igraph-2.3.2-r45hf411e2a_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-jquerylib-0.1.4-r45hc72bb7e_4.conda +https://conda.anaconda.org/conda-forge/noarch/r-pillar-1.11.1-r45hc72bb7e_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-purrr-1.2.2-r45h54b55ab_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-rocr-1.0_12-r45hc72bb7e_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-rspectra-0.16_2-r45h3704496_2.conda +https://conda.anaconda.org/conda-forge/linux-64/r-sass-0.4.10-r45h3697838_1.conda +https://conda.anaconda.org/conda-forge/linux-64/r-spatstat.geom-3.8_1-r45h54b55ab_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-stringr-1.6.0-r45h785f33e_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-textshaping-1.0.5-r45h74f4acd_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-tidyselect-1.2.1-r45hc72bb7e_2.conda +https://conda.anaconda.org/bioconda/linux-64/bioconductor-s4arrays-1.10.1-r45h01b2380_0.conda +https://conda.anaconda.org/bioconda/noarch/bioconductor-seqinfo-1.0.0-r45hdfd78af_0.conda +https://conda.anaconda.org/bioconda/linux-64/bioconductor-xvector-0.50.0-r45h01b2380_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-bslib-0.11.0-r45hc72bb7e_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-cowplot-1.2.0-r45hc72bb7e_2.conda +https://conda.anaconda.org/conda-forge/noarch/r-gargle-1.6.1-r45h785f33e_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-ggrepel-0.9.8-r45h3697838_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-ggridges-0.5.7-r45hc72bb7e_1.conda +https://conda.anaconda.org/conda-forge/linux-64/r-leidenbase-0.1.36-r45ha11a66c_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-patchwork-1.3.2-r45hc72bb7e_1.conda +https://conda.anaconda.org/conda-forge/noarch/r-progress-1.2.3-r45hc72bb7e_2.conda +https://conda.anaconda.org/conda-forge/linux-64/r-ragg-1.5.2-r45h9f1dc4d_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-reshape2-1.4.5-r45h3697838_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-scattermore-1.2-r45h3697838_5.conda +https://conda.anaconda.org/conda-forge/noarch/r-selectr-0.5_1-r45hc72bb7e_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-seuratobject-5.4.0-r45h3697838_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-spatstat.random-3.4_5-r45h3697838_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-tibble-3.3.1-r45h54b55ab_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-uwot-0.2.4-r45h3697838_0.conda +https://conda.anaconda.org/bioconda/linux-64/bioconductor-genomicranges-1.62.1-r45h01b2380_0.conda +https://conda.anaconda.org/bioconda/linux-64/bioconductor-sparsearray-1.10.8-r45h01b2380_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-cellranger-1.1.0-r45hc72bb7e_1008.conda +https://conda.anaconda.org/conda-forge/linux-64/r-dplyr-1.2.1-r45h3697838_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-fastdummies-1.7.6-r45hc72bb7e_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-forcats-1.0.1-r45hc72bb7e_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-googledrive-2.1.2-r45hc72bb7e_1.conda +https://conda.anaconda.org/conda-forge/noarch/r-rematch2-2.1.2-r45hc72bb7e_5.conda +https://conda.anaconda.org/conda-forge/noarch/r-rmarkdown-2.31-r45hc72bb7e_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-rvest-1.0.5-r45hc72bb7e_1.conda +https://conda.anaconda.org/conda-forge/noarch/r-shiny-1.13.0-r45h785f33e_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-spatstat.explore-3.8_0-r45h54b55ab_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-vroom-1.7.1-r45h3697838_0.conda +https://conda.anaconda.org/bioconda/linux-64/bioconductor-delayedarray-0.36.0-r45h01b2380_0.conda +https://conda.anaconda.org/bioconda/linux-64/bioconductor-h5mread-1.2.1-r45h01b2380_1.conda +https://conda.anaconda.org/conda-forge/noarch/r-dtplyr-1.3.3-r45hc72bb7e_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-googlesheets4-1.1.2-r45h785f33e_1.conda +https://conda.anaconda.org/conda-forge/noarch/r-htmlwidgets-1.6.4-r45h785f33e_4.conda +https://conda.anaconda.org/conda-forge/noarch/r-miniui-0.1.2-r45hc72bb7e_1.conda +https://conda.anaconda.org/conda-forge/linux-64/r-readr-2.2.0-r45h3697838_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-readxl-1.5.0-r45h10e25cc_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-reprex-2.1.1-r45hc72bb7e_2.conda +https://conda.anaconda.org/conda-forge/linux-64/r-sctransform-0.4.3-r45h3704496_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-tidyr-1.3.2-r45h3697838_0.conda +https://conda.anaconda.org/bioconda/linux-64/bioconductor-beachmat-2.26.0-r45ha27e39d_0.conda +https://conda.anaconda.org/bioconda/noarch/bioconductor-delayedmatrixstats-1.32.0-r45hdfd78af_0.conda +https://conda.anaconda.org/bioconda/noarch/bioconductor-hdf5array-1.38.0-r45hdfd78af_0.conda +https://conda.anaconda.org/bioconda/noarch/bioconductor-summarizedexperiment-1.40.0-r45hdfd78af_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-broom-1.0.13-r45hc72bb7e_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-dbplyr-2.5.2-r45hc72bb7e_0.conda +https://conda.anaconda.org/conda-forge/linux-64/r-haven-2.5.5-r45h6d565e7_1.conda +https://conda.anaconda.org/conda-forge/noarch/r-plotly-4.12.0-r45hc72bb7e_0.conda +https://conda.anaconda.org/bioconda/noarch/bioconductor-singlecellexperiment-1.32.0-r45hdfd78af_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-modelr-0.1.11-r45hc72bb7e_3.conda +https://conda.anaconda.org/conda-forge/linux-64/r-seurat-5.5.0-r45h3697838_0.conda +https://conda.anaconda.org/bioconda/linux-64/bioconductor-glmgampoi-1.22.0-r45ha27e39d_0.conda +https://conda.anaconda.org/conda-forge/noarch/r-tidyverse-2.0.0-r45h785f33e_3.conda diff --git a/workflow/envs/downsample_clusters.yml b/workflow/envs/downsample_clusters.yml new file mode 100644 index 0000000..dd402ca --- /dev/null +++ b/workflow/envs/downsample_clusters.yml @@ -0,0 +1,11 @@ +channels: + - conda-forge + - bioconda + +channel_priority: strict + +dependencies: + - r-base + - r-seurat + - r-tidyverse + - bioconductor-glmgampoi diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk index 023708c..0a2d73b 100644 --- a/workflow/rules/common.smk +++ b/workflow/rules/common.smk @@ -2,8 +2,10 @@ import os from pathlib import Path -REQUIRED_CONFIG_KEYS = { +BASE_REQUIRED_CONFIG_KEYS = { "conda-channel-priority", +} +PREPROCESS_REQUIRED_CONFIG_KEYS = { "sampleTable", "workflow_seed", "emptydrop_removal_methods", @@ -14,6 +16,7 @@ REQUIRED_CONFIG_KEYS = { "min_ncount", "max_mtdna", } +ALLOWED_WORKFLOW_MODES = {"preprocess", "preprocess_and_downsample", "downsample_only"} ALLOWED_EMPTYDROP_METHODS = {"tenx", "emptydrops"} ALLOWED_DECON_METHODS = {"soupx", "cellbender_fromraw"} ALLOWED_DOUBLET_METHODS = {"doubletfinder", "scdblfinder"} @@ -32,6 +35,11 @@ def require_positive_int(config_values, key, errors): errors.append(f"{key} must be a positive integer") +def require_optional_positive_int(config_values, key, errors): + if key in config_values: + require_positive_int(config_values, key, errors) + + def validate_method_list(config_values, key, allowed_values, errors): value = config_values.get(key) if not isinstance(value, list) or not value: @@ -61,30 +69,83 @@ def validate_method_list(config_values, key, allowed_values, errors): return normalized +def validate_downsample_targets(config_values, errors): + if "downsampleTargets" not in config_values: + return + + value = config_values["downsampleTargets"] + if value == "all": + return + if not isinstance(value, list) or not value: + errors.append("downsampleTargets must be 'all' or a non-empty list") + return + for item in value: + if not isinstance(item, str) or not item.strip(): + errors.append("downsampleTargets contains a non-string or empty value") + return + + duplicates = sorted({item for item in value if value.count(item) > 1}) + if duplicates: + errors.append("downsampleTargets contains duplicate value(s): " + ", ".join(duplicates)) + + def validate_workflow_config(config_values): errors = [] - missing_keys = sorted(REQUIRED_CONFIG_KEYS - set(config_values.keys())) + workflow_mode = config_values.get("workflow_mode", "preprocess") + if workflow_mode not in ALLOWED_WORKFLOW_MODES: + errors.append( + "workflow_mode must be one of: " + ", ".join(sorted(ALLOWED_WORKFLOW_MODES)) + ) + workflow_mode = "preprocess" + + required_keys = set(BASE_REQUIRED_CONFIG_KEYS) + if workflow_mode in {"preprocess", "preprocess_and_downsample"}: + required_keys.update(PREPROCESS_REQUIRED_CONFIG_KEYS) + + missing_keys = sorted(required_keys - set(config_values.keys())) if missing_keys: errors.append("missing required key(s): " + ", ".join(missing_keys)) - require_non_empty_string(config_values, "sampleTable", errors) if config_values.get("conda-channel-priority") != "strict": errors.append("conda-channel-priority must be 'strict'") - workflow_seed = config_values.get("workflow_seed") - if isinstance(workflow_seed, bool) or not isinstance(workflow_seed, int): - errors.append("workflow_seed must be an integer") - - emptydrop_methods = validate_method_list(config_values, "emptydrop_removal_methods", ALLOWED_EMPTYDROP_METHODS, errors) - decon_methods = validate_method_list(config_values, "ambient_decon_methods", ALLOWED_DECON_METHODS, errors) - doublet_methods = validate_method_list(config_values, "doublet_removal_methods", ALLOWED_DOUBLET_METHODS, errors) - posthoc_methods = validate_method_list(config_values, "posthoc_methods", ALLOWED_POSTHOC_METHODS, errors) - - require_positive_int(config_values, "min_nfeature", errors) - require_positive_int(config_values, "min_ncount", errors) - max_mtdna = config_values.get("max_mtdna") - if isinstance(max_mtdna, bool) or not isinstance(max_mtdna, (int, float)) or not 0 <= max_mtdna <= 100: - errors.append("max_mtdna must be a number between 0 and 100") + emptydrop_methods = [] + decon_methods = [] + doublet_methods = [] + posthoc_methods = [] + if workflow_mode in {"preprocess", "preprocess_and_downsample"}: + require_non_empty_string(config_values, "sampleTable", errors) + + workflow_seed = config_values.get("workflow_seed") + if isinstance(workflow_seed, bool) or not isinstance(workflow_seed, int): + errors.append("workflow_seed must be an integer") + + emptydrop_methods = validate_method_list(config_values, "emptydrop_removal_methods", ALLOWED_EMPTYDROP_METHODS, errors) + decon_methods = validate_method_list(config_values, "ambient_decon_methods", ALLOWED_DECON_METHODS, errors) + doublet_methods = validate_method_list(config_values, "doublet_removal_methods", ALLOWED_DOUBLET_METHODS, errors) + posthoc_methods = validate_method_list(config_values, "posthoc_methods", ALLOWED_POSTHOC_METHODS, errors) + + require_positive_int(config_values, "min_nfeature", errors) + require_positive_int(config_values, "min_ncount", errors) + max_mtdna = config_values.get("max_mtdna") + if isinstance(max_mtdna, bool) or not isinstance(max_mtdna, (int, float)) or not 0 <= max_mtdna <= 100: + errors.append("max_mtdna must be a number between 0 and 100") + + if workflow_mode in {"preprocess_and_downsample", "downsample_only"}: + if "downsampleSeuratObjectDir" in config_values: + require_non_empty_string(config_values, "downsampleSeuratObjectDir", errors) + if "downsampleResultsDir" in config_values: + require_non_empty_string(config_values, "downsampleResultsDir", errors) + require_optional_positive_int(config_values, "nDownsampleReplicates", errors) + if "workflowSeed" in config_values: + seed = config_values["workflowSeed"] + if isinstance(seed, bool) or not isinstance(seed, int): + errors.append("workflowSeed must be an integer") + if "downsampleRate" in config_values: + rate = config_values["downsampleRate"] + if isinstance(rate, bool) or not isinstance(rate, (int, float)) or not 0 < rate <= 1: + errors.append("downsampleRate must be > 0 and <= 1") + validate_downsample_targets(config_values, errors) if "resultsDir" in config_values: require_non_empty_string(config_values, "resultsDir", errors) @@ -93,6 +154,7 @@ def validate_workflow_config(config_values): raise ValueError("Invalid workflow config: " + "; ".join(errors)) return { + "workflow_mode": workflow_mode, "emptydrop_methods": emptydrop_methods, "decon_methods": decon_methods, "doublet_methods": doublet_methods, @@ -187,3 +249,47 @@ def marker_chunk_inputs(wildcards): prefix=wildcards.prefix, cluster=cluster_ids, ) + + +def marker_targets(prefixes): + return [f"{RESULTS_DIR}/{prefix}_markergenes.csv" for prefix in prefixes] + + +def rds_targets(prefixes): + return [f"{RESULTS_DIR}/{prefix}.rds" for prefix in prefixes] + + +def select_downsample_inputs(inputs_by_target): + requested = config.get("downsampleTargets", "all") + if requested == "all" or requested == ["all"]: + return dict(inputs_by_target) + + missing = sorted(set(requested) - set(inputs_by_target)) + if missing: + raise ValueError( + "downsampleTargets contains unknown target(s): " + + ", ".join(missing) + + "; available targets are: " + + ", ".join(sorted(inputs_by_target)) + ) + return {target: inputs_by_target[target] for target in requested} + + +def downsample_inputs_from_preprocess_outputs(): + inputs = {} + for path in PREPROCESS_SEURAT_TARGETS: + target = Path(path).stem + if target in inputs and inputs[target] != path: + raise ValueError(f"duplicate downsample target name: {target}") + inputs[target] = path + return inputs + + +def downsample_inputs_from_external_dir(): + targets = sorted( + glob_wildcards(f"{DOWNSAMPLE_SEURAT_OBJECT_DIR}/{{downsample_target}}.rds").downsample_target + ) + return { + target: f"{DOWNSAMPLE_SEURAT_OBJECT_DIR}/{target}.rds" + for target in targets + } diff --git a/workflow/rules/downsample_clusters.smk b/workflow/rules/downsample_clusters.smk new file mode 100644 index 0000000..d729dad --- /dev/null +++ b/workflow/rules/downsample_clusters.smk @@ -0,0 +1,26 @@ +N_DOWNSAMPLE_REPLICATES = int(config.get("nDownsampleReplicates", 100)) + + +rule downsample_clusters: + input: + seurat_object=lambda wildcards: DOWNSAMPLE_INPUTS_BY_TARGET[wildcards.downsample_target], + script="workflow/scripts/downsample_clusters.R" + output: + tsv=f"{DOWNSAMPLE_RESULTS_DIR}/{{downsample_target}}_clusterdownsampling.tsv" + log: + f"{DOWNSAMPLE_RESULTS_DIR}/logs/downsample_clusters/{{downsample_target}}.log" + params: + downsample_rate=lambda wildcards: float(config.get("downsampleRate", 0.8)), + seed=lambda wildcards: int(config.get("workflowSeed", config.get("workflow_seed", 12345))), + n_replicates=N_DOWNSAMPLE_REPLICATES + conda: + "../envs/downsample_clusters.yml" + wildcard_constraints: + downsample_target=DOWNSAMPLE_TARGET_REGEX + resources: + mem_mb=lambda wildcards, attempt: int(24000 * (2 ** (attempt - 1))), + runtime=lambda wildcards, attempt: int(270 * (2 ** (attempt - 1))) + shell: + """ + SCRNASEQ_DOWNSAMPLE_SEED={params.seed} Rscript {input.script} {input.seurat_object} {output.tsv} {params.n_replicates} {params.downsample_rate} > {log} 2>&1 + """ diff --git a/workflow/rules/markers.smk b/workflow/rules/markers.smk index 3897ace..d22b29c 100644 --- a/workflow/rules/markers.smk +++ b/workflow/rules/markers.smk @@ -52,6 +52,7 @@ rule combine_markers: shell: """ SCRNASEQ_PREPROCESS_SEED={params.seed} Rscript {input.script} {output} {input.markers} > {log} 2>&1 + rm -f {input.markers} rm -f {params.cluster_ids} rm -f {params.nclusters} rm -rf {params.manifest} diff --git a/workflow/scripts/downsample_clusters.R b/workflow/scripts/downsample_clusters.R new file mode 100644 index 0000000..8167624 --- /dev/null +++ b/workflow/scripts/downsample_clusters.R @@ -0,0 +1,137 @@ +args <- commandArgs(trailingOnly = TRUE) +if (length(args) < 3) { + stop("Usage: downsample_clusters.R [downsample_rate]", call. = FALSE) +} + +seurat_rds <- args[1] +output <- args[2] +n_replicates <- suppressWarnings(as.integer(args[3])) +downsample_rate <- if (length(args) >= 4) as.numeric(args[4]) else as.numeric(Sys.getenv("SCRNASEQ_DOWNSAMPLE_RATE", "0.8")) +workflow_seed <- suppressWarnings(as.integer(Sys.getenv("SCRNASEQ_DOWNSAMPLE_SEED", "12345"))) + +if (length(n_replicates) != 1 || is.na(n_replicates) || n_replicates < 1) { + stop("n_replicates must be a positive integer", call. = FALSE) +} +if (length(downsample_rate) != 1 || is.na(downsample_rate) || downsample_rate <= 0 || downsample_rate > 1) { + stop("downsample_rate must be > 0 and <= 1", call. = FALSE) +} +if (length(workflow_seed) != 1 || is.na(workflow_seed)) { + stop("SCRNASEQ_DOWNSAMPLE_SEED must be an integer", call. = FALSE) +} +set.seed(workflow_seed) + +suppressPackageStartupMessages({ + library("tidyverse") + library("Seurat") + library("glmGamPoi") +}) +options(future.globals.maxSize = 2 * 1024^3) + +JaccardSimilarity <- function(set1, set2) { + intersect_length <- length(intersect(set1, set2)) + union_length <- length(set1) + length(set2) - intersect_length + intersect_length / union_length +} + +RandomSubsetData <- function(object, rate, random.subset.seed = NULL, ...) { + ncells <- nrow(object@meta.data) + ncells.subsample <- round(ncells * rate) + + set.seed(random.subset.seed) + + selected.cells <- sample(colnames(object), ncells.subsample) + object <- subset(object, cells = selected.cells, ...) + return(object) +} + +SubSampleReSCTSeuratObject <- function(seurat_obj, subrate, replicate_seed) { + set.seed(replicate_seed) + subsampled_seurat <- RandomSubsetData( + seurat_obj, + rate = subrate, + random.subset.seed = replicate_seed + ) + subsampled_seurat$presub_clusters <- seurat_obj@meta.data[ + colnames(subsampled_seurat), + "seurat_clusters", + 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) + } + 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) + subsampled_seurat <- FindClusters( + subsampled_seurat, + random.seed = replicate_seed, + verbose = FALSE + ) + return(subsampled_seurat) +} + +GetJaccardMaxByCluster <- function(seurat_obj, bootstrap) { + jaccard_max_stats <- tibble::tibble( + clusterid = factor(), + max_jaccard = numeric(), + bootstrap_number = integer() + ) + + 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)) + ) %>% + dplyr::pull(jaccard) %>% + unlist() %>% + max() + + jaccard_max_stats <- jaccard_max_stats %>% + tibble::add_row( + clusterid = original_cluster, + max_jaccard = maxstat, + bootstrap_number = bootstrap + ) + } + + return(jaccard_max_stats) +} + +seurat_obj <- readRDS(seurat_rds) +if (!"seurat_clusters" %in% colnames(seurat_obj@meta.data)) { + stop("Input Seurat object is missing required metadata column: seurat_clusters", call. = FALSE) +} + +replicate_results <- vector("list", n_replicates) +for (replicate in seq_len(n_replicates)) { + t0 <- Sys.time() + replicate_seed <- workflow_seed + replicate + subsampled_obj <- SubSampleReSCTSeuratObject(seurat_obj, downsample_rate, replicate_seed) + replicate_results[[replicate]] <- GetJaccardMaxByCluster(subsampled_obj, replicate) + rm(subsampled_obj) + gc(verbose = FALSE) + message(sprintf( + "[%s] replicate %d/%d done: %.1f sec", + format(Sys.time(), "%H:%M:%S"), replicate, n_replicates, + as.numeric(Sys.time() - t0, units = "secs") + )) +} + +jaccard_max_stats <- dplyr::bind_rows(replicate_results) + +dir.create(dirname(output), showWarnings = FALSE, recursive = TRUE) +write_tsv(jaccard_max_stats, output)