Debug

README.md

@@ -66,5 +66,11 @@ The main Snakemake entrypoint supports three `workflow_mode` values in `config/c
 
 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.
 
+## Scalability notes
+
+Most rules scale roughly linearly in cell count (dominated by the `SCTransform` working set), and per-rule memory requests in `workflow/rules/*.smk` are sized accordingly. One exception matters for large datasets:
+
+- **DoubletFinder memory grows as O(N²).** `doubletfinder` augments the full dataset with ~25% synthetic doublets and builds a dense pairwise distance matrix over all cells: roughly 78 GB at 77k cells, 132 GB at 100k, 298 GB at 150k, and 530 GB at 200k. This is inherent to the algorithm and cannot be tuned away. For large inputs (roughly >100k cells — e.g. emptyDrops cell calls, which are often much larger than the CellRanger filtered set) prefer **scDblFinder** via the `doublet_removal_methods` config, as it does not materialize a full distance matrix and scales far better.
+
 ## Tests
 For information on how to run the test suite, or run the workflow in test mode, see tests/README.md.

config/config.yaml

@@ -4,6 +4,14 @@ sampleTable: "samplesheet.tsv"
 workflow_mode: "preprocess_and_downsample" # {"preprocess", "preprocess_and_downsample", "downsample_only"}
 workflow_seed: 12345
 
+# Sample IDs to skip when building the workflow (they are dropped from the DAG, so their
+# outputs are neither generated nor regenerated). Use this to exclude samples flagged as
+# low quality after a run (see the run log and results/low_quality_samples/flagged_samples.txt)
+# or any sample you want to hold out. Exclusion is explicit and logged at startup - nothing
+# is skipped automatically. IMPORTANT: if you re-sequence a bad library and reuse its sample
+# ID in samplesheet.tsv, remove that ID from this list so the new data is processed.
+excluded_samples: []
+
 emptydrop_removal_methods: ["tenx","emptydrops"]
 ambient_decon_methods: ["soupx","cellbender_fromraw"]
 doublet_removal_methods: ["doubletfinder", "scdblfinder"]

profiles/slurm/config.yaml

@@ -1,13 +1,22 @@
 use-singularity: True
 singularity-args: "\"--nv\""
+
+# Execution policy (workflow-level, so it applies to every invocation regardless of the
+# launcher script). Submission/status-check rates are throttled to avoid overwhelming the
+# SLURM controller when thousands of jobs are queued.
 retries: 2
+keep-going: True
+jobs: 1500
+max-jobs-per-timespan: "10/1s"
+max-status-checks-per-second: 5
+latency-wait: 120
+rerun-incomplete: True
 
 set-threads:
     tenx2seuratrds: 1
     soupx: 1
     scdblfinder: 1
     emptydrops: 1
-    downsample_cluster_replicate: 1
     downsample_clusters: 1
 
 set-resources:
@@ -46,11 +55,6 @@ set-resources:
         tasks: 1
         cpus_per_task: 1
 
-    downsample_cluster_replicate:
-        nodes: 1
-        tasks: 1
-        cpus_per_task: 1
-
     downsample_clusters:
         nodes: 1
         tasks: 1

scrnaseq_preprocess_slurmrunner.sh

@@ -35,4 +35,7 @@ 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 --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
+# Execution policy (retries, keep-going, job/submission throttling, latency-wait,
+# rerun-incomplete) lives in profiles/slurm/config.yaml, and low-quality-sample quarantine
+# runs from the Snakefile's onsuccess/onerror handlers - so both apply to any launcher.
+snakemake --conda-prefix $PATH_TO_MY_CONDA_ENVS --snakefile workflow/Snakefile --configfile config/config.yaml --use-conda --workflow-profile profiles/slurm --profile cannon

workflow/Snakefile

@@ -1,4 +1,5 @@
 import re
+import sys
 from pathlib import Path
 
 import pandas as pd
@@ -16,10 +17,29 @@ config.setdefault("min_nfeature", 1)
 config.setdefault("min_ncount", 1)
 config.setdefault("max_mtdna", 100)
 
+EXCLUDED_SAMPLES = validated_config["excluded_samples"]
 if RUN_PREPROCESS:
     sampleinfo = pd.read_table(config["sampleTable"], dtype={"sampleid": str})
     sampleinfo = validate_sample_sheet(sampleinfo, config["sampleTable"])
     SAMPLES = list(sampleinfo.loc[:, "sampleid"])
+    if EXCLUDED_SAMPLES:
+        present = [s for s in SAMPLES if s in EXCLUDED_SAMPLES]
+        absent = [s for s in EXCLUDED_SAMPLES if s not in SAMPLES]
+        if present:
+            print(
+                "[excluded_samples] skipping "
+                f"{len(present)} sample(s) from the DAG: {', '.join(present)}",
+                file=sys.stderr,
+            )
+        if absent:
+            print(
+                "[excluded_samples] note: "
+                f"{len(absent)} excluded id(s) are not in the sample sheet and were ignored: "
+                f"{', '.join(absent)}",
+                file=sys.stderr,
+            )
+        SAMPLES = [s for s in SAMPLES if s not in EXCLUDED_SAMPLES]
+        sampleinfo = sampleinfo[~sampleinfo["sampleid"].isin(EXCLUDED_SAMPLES)].reset_index(drop=True)
 else:
     sampleinfo = pd.DataFrame({"sampleid": [], "tenx_datadir": []})
     SAMPLES = []
@@ -126,6 +146,18 @@ if RUN_DOWNSAMPLE:
             "No Seurat RDS inputs were found for downsampling. "
             "For downsample_only mode, set downsampleSeuratObjectDir to a directory containing *.rds files."
         )
+    # Drop downsample targets belonging to excluded samples. In preprocess_and_downsample
+    # this is already handled via SAMPLES; it matters for downsample_only, where targets come
+    # from external .rds files rather than the (filtered) sample sheet.
+    if EXCLUDED_SAMPLES:
+        DOWNSAMPLE_INPUTS_BY_TARGET = {
+            target: path
+            for target, path in DOWNSAMPLE_INPUTS_BY_TARGET.items()
+            if not any(
+                re.search(rf"(?:^|[_/]){re.escape(s)}(?:[_.]|$)", target)
+                for s in EXCLUDED_SAMPLES
+            )
+        }
 else:
     DOWNSAMPLE_INPUTS_BY_TARGET = {}
 
@@ -167,3 +199,21 @@ include: "rules/posthocfilter_mad.smk"
 include: "rules/posthocfilter_threshold.smk"
 include: "rules/markers.smk"
 include: "rules/downsample_clusters.smk"
+
+
+# Post-workflow finalization (defined in rules/common.smk), run from BOTH handlers so it fires
+# on completion regardless of Snakemake's own exit code. It quarantines low-quality samples
+# and verifies completeness, then drives the process exit code so the runner batch job's
+# COMPLETED/FAILED status reflects the true outcome (a real failure, e.g. an OOM that exhausted
+# retries, fails the run; failures confined to flagged low-quality samples are excused).
+# It runs only for non-local (SLURM) executor runs, where the verification is needed to correct
+# that executor's unreliable exit code; local runs (tests, ad-hoc builds) have reliable exit
+# codes and may intentionally build only a subset of targets, so they are skipped.
+onsuccess:
+    finalize_run(RESULTS_DIR, ALL_TARGETS, str(config.get("sampleTable", "samplesheet.tsv")),
+                 getattr(workflow, "non_local_exec", False))
+
+
+onerror:
+    finalize_run(RESULTS_DIR, ALL_TARGETS, str(config.get("sampleTable", "samplesheet.tsv")),
+                 getattr(workflow, "non_local_exec", False))

workflow/rules/cellbender2seurat.smk

@@ -12,7 +12,7 @@ rule cellbender2seurat:
     conda:
         "../envs/tenx2seuratrds.yml"
     resources:
-        mem_mb = lambda wildcards, attempt: int(24000 * (2 ** (attempt - 1))),
+        mem_mb = lambda wildcards, attempt: int(48000 * (2 ** (attempt - 1))),
         runtime = lambda wildcards, attempt: int(480* (2 ** (attempt - 1)))
     params:
         seed=WORKFLOW_SEED

workflow/rules/common.smk

@@ -89,6 +89,27 @@ def validate_downsample_targets(config_values, errors):
         errors.append("downsampleTargets contains duplicate value(s): " + ", ".join(duplicates))
 
 
+def validate_excluded_samples(config_values, errors):
+    """Normalize the optional user-controlled excluded_samples list (default empty).
+
+    Sample IDs listed here are dropped from the DAG. Exclusion is explicit and visible so
+    nothing is skipped automatically (e.g. a re-sequenced library reusing a sample ID is
+    processed normally unless the user deliberately lists that ID here)."""
+    value = config_values.get("excluded_samples", [])
+    if value is None:
+        return []
+    if not isinstance(value, list):
+        errors.append("excluded_samples must be a list of sample id strings")
+        return []
+    normalized = []
+    for item in value:
+        if not isinstance(item, str) or not item.strip():
+            errors.append("excluded_samples contains a non-string or empty value")
+            continue
+        normalized.append(item.strip())
+    return sorted(set(normalized))
+
+
 def validate_workflow_config(config_values):
     errors = []
     workflow_mode = config_values.get("workflow_mode", "preprocess")
@@ -150,6 +171,8 @@ def validate_workflow_config(config_values):
     if "resultsDir" in config_values:
         require_non_empty_string(config_values, "resultsDir", errors)
 
+    excluded_samples = validate_excluded_samples(config_values, errors)
+
     if errors:
         raise ValueError("Invalid workflow config: " + "; ".join(errors))
 
@@ -159,6 +182,7 @@ def validate_workflow_config(config_values):
         "decon_methods": decon_methods,
         "doublet_methods": doublet_methods,
         "posthoc_methods": posthoc_methods,
+        "excluded_samples": excluded_samples,
     }
 
 
@@ -293,3 +317,40 @@ def downsample_inputs_from_external_dir():
         target: f"{DOWNSAMPLE_SEURAT_OBJECT_DIR}/{target}.rds"
         for target in targets
     }
+
+
+def finalize_run(results_dir, all_targets, samplesheet, enabled):
+    """Post-run quarantine + completeness verification that drives the process exit code.
+
+    Runs from the Snakefile's onsuccess/onerror handlers. It (1) quarantines low-quality
+    samples and (2) verifies that every required target exists, failing the run only when a
+    required output is missing for a sample that was NOT flagged low quality. sys.exit() in a
+    handler deterministically sets Snakemake's exit code either way, so the runner batch job's
+    COMPLETED/FAILED state reflects the true outcome even when Snakemake would exit 0 on a
+    terminal failure.
+
+    `enabled` should be True only for non-local (SLURM) runs: the verification exists to
+    correct that executor's unreliable exit code, whereas local runs (tests, ad-hoc builds)
+    have reliable exit codes and may intentionally build only a subset of targets.
+    """
+    if not enabled:
+        return
+
+    import subprocess
+    import sys as _sys
+    from pathlib import Path as _Path
+
+    _Path(results_dir).mkdir(parents=True, exist_ok=True)
+    targets_file = _Path(results_dir) / ".run_targets.txt"
+    targets_file.write_text("\n".join(all_targets) + "\n")
+
+    subprocess.run([
+        _sys.executable, "workflow/scripts/quarantine_low_quality_samples.py",
+        "--results-dir", results_dir, "--samplesheet", samplesheet,
+    ])
+    result = subprocess.run([
+        _sys.executable, "workflow/scripts/verify_run_complete.py",
+        "--results-dir", results_dir, "--samplesheet", samplesheet,
+        "--targets-file", str(targets_file),
+    ])
+    _sys.exit(result.returncode)

workflow/rules/doubletfinder.smk

@@ -16,7 +16,11 @@ rule doubletfinder:
         decon_method="soupx",
         empty_method="tenx|emptydrops"
     resources:
-        mem_mb = lambda wildcards, attempt: int(24000 * (2 ** (attempt - 1))),
+        # DoubletFinder's paramSweep (sct=TRUE) synthesizes ~25% artificial doublets and
+        # re-runs SCTransform across a pK sweep, so peak memory is a large multiple of the
+        # input object (~43x observed: 2.1 GB rds -> 92 GB). A flat baseline is either
+        # wasteful for small samples or fatal for large ones; scale by input size (floor 48 GB).
+        mem_mb = lambda wildcards, input, attempt: int(max(48000, 64 * input.size_mb) * (2 ** (attempt - 1))),
         runtime = lambda wildcards, attempt: int(480* (2 ** (attempt - 1)))
     params:
         seed=WORKFLOW_SEED
@@ -41,7 +45,11 @@ rule doubletfinder_cellbender:
     conda:
         "../envs/doubletfinder.yml"
     resources:
-        mem_mb = lambda wildcards, attempt: int(24000 * (2 ** (attempt - 1))),
+        # DoubletFinder's paramSweep (sct=TRUE) synthesizes ~25% artificial doublets and
+        # re-runs SCTransform across a pK sweep, so peak memory is a large multiple of the
+        # input object (~43x observed: 2.1 GB rds -> 92 GB). A flat baseline is either
+        # wasteful for small samples or fatal for large ones; scale by input size (floor 48 GB).
+        mem_mb = lambda wildcards, input, attempt: int(max(48000, 64 * input.size_mb) * (2 ** (attempt - 1))),
         runtime = lambda wildcards, attempt: int(480* (2 ** (attempt - 1)))
     params:
         seed=WORKFLOW_SEED

workflow/rules/downsample_clusters.smk

@@ -18,8 +18,15 @@ rule downsample_clusters:
     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)))
+        # All 100 replicates run serially in one job (rm + gc between iterations), so the
+        # job's peak memory is a SINGLE replicate's SCTransform footprint. Observed per-
+        # replicate peak is ~42 GB and does NOT track input rds size (small datasets such as
+        # cteleta are among the heaviest), so a flat baseline beats input-scaling here.
+        # Wall-time is ~100x the per-replicate cost (worst observed ~370 s/rep -> ~10.2 h).
+        # 64 GB / 15 h cover the worst case on attempt 1; OOM/TIMEOUT restarts are not
+        # reliably resubmitted, so baselines must not depend on the retry escalation.
+        mem_mb=lambda wildcards, attempt: int(64000 * (2 ** (attempt - 1))),
+        runtime=lambda wildcards, attempt: int(900 * (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

workflow/rules/emptydrops.smk

@@ -13,7 +13,7 @@ rule emptydrops:
     conda:
         "../envs/emptydrops.yml"
     resources:
-        mem_mb = lambda wildcards, attempt: int(24000 * (2 ** (attempt - 1))),
+        mem_mb = lambda wildcards, attempt: int(64000 * (2 ** (attempt - 1))),
         runtime = lambda wildcards, attempt: int(480* (2 ** (attempt - 1)))
     params:
         seed=WORKFLOW_SEED

workflow/rules/markers.smk

@@ -23,7 +23,7 @@ rule find_markers:
     conda:
         "../envs/tenx2seuratrds.yml"
     resources:
-        mem_mb=lambda wildcards, attempt: int(12000 * (2 ** (attempt - 1))),
+        mem_mb=lambda wildcards, attempt: int(24000 * (2 ** (attempt - 1))),
         runtime=lambda wildcards, attempt: int(240 * (2 ** (attempt - 1)))
     params:
         seed=WORKFLOW_SEED

workflow/rules/posthocfilter_mad.smk

@@ -16,7 +16,7 @@ rule posthocfilter_mad:
         decon_method="soupx",
         empty_method="tenx|emptydrops"
     resources:
-        mem_mb = lambda wildcards, attempt: int(24000 * (2 ** (attempt - 1))),
+        mem_mb = lambda wildcards, attempt: int(64000 * (2 ** (attempt - 1))),
         runtime = lambda wildcards, attempt: int(480* (2 ** (attempt - 1)))
     params:
         seed=WORKFLOW_SEED
@@ -42,7 +42,7 @@ rule posthocfilter_mad_cellbender:
     wildcard_constraints:
         doublet_method="doubletfinder|scdblfinder"
     resources:
-        mem_mb = lambda wildcards, attempt: int(24000 * (2 ** (attempt - 1))),
+        mem_mb = lambda wildcards, attempt: int(64000 * (2 ** (attempt - 1))),
         runtime = lambda wildcards, attempt: int(480* (2 ** (attempt - 1)))
     params:
         seed=WORKFLOW_SEED

workflow/rules/posthocfilter_threshold.smk

@@ -16,7 +16,7 @@ rule posthocfilter_threshold:
         decon_method="soupx",
         empty_method="tenx|emptydrops"
     resources:
-        mem_mb = lambda wildcards, attempt: int(24000 * (2 ** (attempt - 1))),
+        mem_mb = lambda wildcards, attempt: int(64000 * (2 ** (attempt - 1))),
         runtime = lambda wildcards, attempt: int(480* (2 ** (attempt - 1)))
     params:
         min_numfeatures = config["min_nfeature"],
@@ -45,7 +45,7 @@ rule posthocfilter_threshold_cellbender:
     wildcard_constraints:
         doublet_method="doubletfinder|scdblfinder"
     resources:
-        mem_mb = lambda wildcards, attempt: int(24000 * (2 ** (attempt - 1))),
+        mem_mb = lambda wildcards, attempt: int(64000 * (2 ** (attempt - 1))),
         runtime = lambda wildcards, attempt: int(480* (2 ** (attempt - 1)))
     params:
         min_numfeatures = config["min_nfeature"],

workflow/rules/scdblfinder.smk

@@ -15,7 +15,7 @@ rule scdblfinder:
         decon_method="soupx",
         empty_method="tenx|emptydrops"
     resources:
-        mem_mb = lambda wildcards, attempt: int(24000 * (2 ** (attempt - 1))),
+        mem_mb = lambda wildcards, attempt: int(48000 * (2 ** (attempt - 1))),
         runtime = lambda wildcards, attempt: int(480* (2 ** (attempt - 1)))
     params:
         seed=WORKFLOW_SEED
@@ -39,7 +39,7 @@ rule scdblfinder_cellbender:
     conda:
         "../envs/scdblfinder.yml"
     resources:
-        mem_mb = lambda wildcards, attempt: int(24000 * (2 ** (attempt - 1))),
+        mem_mb = lambda wildcards, attempt: int(48000 * (2 ** (attempt - 1))),
         runtime = lambda wildcards, attempt: int(480* (2 ** (attempt - 1)))
     params:
         seed=WORKFLOW_SEED

workflow/rules/soupx.smk

@@ -14,7 +14,7 @@ rule soupx:
     conda:
         "../envs/soupx.yml"
     resources:
-        mem_mb = lambda wildcards, attempt: int(24000 * (2 ** (attempt - 1))),
+        mem_mb = lambda wildcards, attempt: int(64000 * (2 ** (attempt - 1))),
         runtime = lambda wildcards, attempt: int(480* (2 ** (attempt - 1)))
     params:
         seed=WORKFLOW_SEED

workflow/rules/soupx_emptydrops.smk

@@ -14,7 +14,7 @@ rule soupx_emptydrops:
     conda:
         "../envs/soupx.yml"
     resources:
-        mem_mb = lambda wildcards, attempt: int(24000 * (2 ** (attempt - 1))),
+        mem_mb = lambda wildcards, attempt: int(64000 * (2 ** (attempt - 1))),
         runtime = lambda wildcards, attempt: int(480* (2 ** (attempt - 1)))
     params:
         seed=WORKFLOW_SEED

workflow/rules/tenx2seuratrds.smk

@@ -12,7 +12,7 @@ rule tenx2seuratrds:
     conda:
         "../envs/tenx2seuratrds.yml"
     resources:
-        mem_mb = lambda wildcards, attempt: int(24000 * (2 ** (attempt - 1))),
+        mem_mb = lambda wildcards, attempt: int(48000 * (2 ** (attempt - 1))),
         runtime = lambda wildcards, attempt: int(480* (2 ** (attempt - 1)))
     params:
         seed=WORKFLOW_SEED