diff --git a/README.md b/README.md
index f8054a8..e3fe904 100644
--- a/README.md
+++ b/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.
diff --git a/config/config.yaml b/config/config.yaml
index 741101a..6df43fb 100644
--- a/config/config.yaml
+++ b/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"]
diff --git a/profiles/slurm/config.yaml b/profiles/slurm/config.yaml
index 14957f9..2356b5f 100644
--- a/profiles/slurm/config.yaml
+++ b/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
diff --git a/scrnaseq_preprocess_slurmrunner.sh b/scrnaseq_preprocess_slurmrunner.sh
index e11056a..ca5458a 100755
--- a/scrnaseq_preprocess_slurmrunner.sh
+++ b/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
diff --git a/workflow/Snakefile b/workflow/Snakefile
index bf22523..30126d4 100644
--- a/workflow/Snakefile
+++ b/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))
diff --git a/workflow/rules/cellbender2seurat.smk b/workflow/rules/cellbender2seurat.smk
index 4a536cc..a130e95 100644
--- a/workflow/rules/cellbender2seurat.smk
+++ b/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
diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk
index 0a2d73b..7eee7b1 100644
--- a/workflow/rules/common.smk
+++ b/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)
diff --git a/workflow/rules/doubletfinder.smk b/workflow/rules/doubletfinder.smk
index 9eb8431..cdb46e0 100644
--- a/workflow/rules/doubletfinder.smk
+++ b/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
diff --git a/workflow/rules/downsample_clusters.smk b/workflow/rules/downsample_clusters.smk
index d729dad..20e336f 100644
--- a/workflow/rules/downsample_clusters.smk
+++ b/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
diff --git a/workflow/rules/emptydrops.smk b/workflow/rules/emptydrops.smk
index 10825bc..b1673a9 100644
--- a/workflow/rules/emptydrops.smk
+++ b/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
diff --git a/workflow/rules/markers.smk b/workflow/rules/markers.smk
index d22b29c..0deb2f9 100644
--- a/workflow/rules/markers.smk
+++ b/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
diff --git a/workflow/rules/posthocfilter_mad.smk b/workflow/rules/posthocfilter_mad.smk
index 53c84ae..e288ea3 100644
--- a/workflow/rules/posthocfilter_mad.smk
+++ b/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
diff --git a/workflow/rules/posthocfilter_threshold.smk b/workflow/rules/posthocfilter_threshold.smk
index 4d4c1e0..c1ef039 100644
--- a/workflow/rules/posthocfilter_threshold.smk
+++ b/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"],
diff --git a/workflow/rules/scdblfinder.smk b/workflow/rules/scdblfinder.smk
index e63e807..a3a2555 100644
--- a/workflow/rules/scdblfinder.smk
+++ b/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
diff --git a/workflow/rules/soupx.smk b/workflow/rules/soupx.smk
index 2a3a431..c1095a5 100644
--- a/workflow/rules/soupx.smk
+++ b/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
diff --git a/workflow/rules/soupx_emptydrops.smk b/workflow/rules/soupx_emptydrops.smk
index 3bc9869..541e269 100644
--- a/workflow/rules/soupx_emptydrops.smk
+++ b/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
diff --git a/workflow/rules/tenx2seuratrds.smk b/workflow/rules/tenx2seuratrds.smk
index fa5ab8a..dd8a82a 100644
--- a/workflow/rules/tenx2seuratrds.smk
+++ b/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
diff --git a/workflow/scripts/cellbender2seurat.R b/workflow/scripts/cellbender2seurat.R
index 3400be9..a800aab 100644
--- a/workflow/scripts/cellbender2seurat.R
+++ b/workflow/scripts/cellbender2seurat.R
@@ -22,9 +22,14 @@ write_cluster_metadata <- function(seurat_obj, nclusters_output, cluster_ids_out
 }
 
 mat <- Read_CellBender_h5_Mat(cellbender_h5)
-seurat <- CreateSeuratObject(mat)
+# CellBender's *_filtered.h5 can still contain a rare barcode whose counts were
+# all subtracted to zero during background removal. A zero-count cell yields
+# log_umi = log10(0) = -Inf, which makes SCTransform's make_cell_attr abort.
+# min.features = 1 drops only these empty barcodes and leaves all real cells intact.
+seurat <- CreateSeuratObject(mat, min.features = 1)
 seurat[["percent.mt"]] <- PercentageFeatureSet(seurat, pattern = "(?i)^mt-")
 seurat <- SCTransform(seurat, vars.to.regress = "percent.mt", seed.use = WORKFLOW_SEED, verbose = FALSE)
+require_min_cells_for_pca(seurat, context = "cellbender2seurat")
 seurat <- RunPCA(seurat, seed.use = WORKFLOW_SEED, verbose = FALSE)
 seurat <- RunUMAP(seurat, dims = 1:30, seed.use = WORKFLOW_SEED)
 seurat <- FindNeighbors(seurat, dims = 1:30)
diff --git a/workflow/scripts/combine_markers.R b/workflow/scripts/combine_markers.R
index bf8ecc6..f8d7c29 100644
--- a/workflow/scripts/combine_markers.R
+++ b/workflow/scripts/combine_markers.R
@@ -9,5 +9,14 @@ marker_tables <- lapply(
   function(marker_path) read_csv(marker_path, show_col_types = FALSE)
 )
 
-combined_markers <- bind_rows(marker_tables)
+# Clusters with < 3 cells produce header-only (0-row) files. read_csv infers their
+# empty columns as logical, which would clash with the character/double types from
+# populated files under bind_rows. Drop empties before binding; if every cluster was
+# skipped, fall back to a single header-only table so the output still has the schema.
+non_empty <- marker_tables[vapply(marker_tables, nrow, integer(1)) > 0]
+combined_markers <- if (length(non_empty) > 0) {
+  bind_rows(non_empty)
+} else {
+  marker_tables[[1]]
+}
 write_csv(combined_markers, file = output_path)
diff --git a/workflow/scripts/doubletfinder.R b/workflow/scripts/doubletfinder.R
index 98ad1de..b86e7d4 100644
--- a/workflow/scripts/doubletfinder.R
+++ b/workflow/scripts/doubletfinder.R
@@ -1,3 +1,12 @@
+# SCALABILITY NOTE: doubletFinder() augments the full dataset with ~25% synthetic doublets
+# and then builds a DENSE pairwise distance matrix (fields::rdist) over all cells. That matrix
+# is O(N^2): ~78 GB at 77k cells, ~132 GB at 100k, ~298 GB at 150k, ~530 GB at 200k. paramSweep
+# is not the driver (it caps at a 10k-cell subsample); the main classification step is. This is
+# inherent to the DoubletFinder algorithm and cannot be tuned away here, so for large datasets
+# (roughly >100k cells, e.g. emptyDrops outputs) prefer scDblFinder (the workflow's other
+# doublet_removal_method), which does not materialize a full distance matrix and scales far
+# better. Memory for this rule is provisioned by input size in workflow/rules/doubletfinder.smk.
+
 args <- commandArgs(trailingOnly = TRUE)
 seurat <- args[1]
 output <- args[2]
@@ -50,6 +59,7 @@ gc()
 
 seurat_singlets[["percent.mt"]] <- PercentageFeatureSet(seurat_singlets, pattern = "(?i)^mt-")
 seurat_singlets <- SCTransform(seurat_singlets, vars.to.regress = "percent.mt", seed.use = WORKFLOW_SEED, verbose = FALSE)
+require_min_cells_for_pca(seurat_singlets, context = "doubletfinder")
 seurat_singlets <- RunPCA(seurat_singlets, seed.use = WORKFLOW_SEED, verbose = FALSE)
 seurat_singlets <- RunUMAP(seurat_singlets, dims = 1:30, seed.use = WORKFLOW_SEED)
 seurat_singlets <- FindNeighbors(seurat_singlets, dims = 1:30)
diff --git a/workflow/scripts/downsample_clusters.R b/workflow/scripts/downsample_clusters.R
index 0d37790..3ae459a 100644
--- a/workflow/scripts/downsample_clusters.R
+++ b/workflow/scripts/downsample_clusters.R
@@ -25,7 +25,7 @@ suppressPackageStartupMessages({
   library("Seurat")
   library("glmGamPoi")
 })
-options(future.globals.maxSize = 2 * 1024^3)
+options(future.globals.maxSize = 16 * 1024^3)
 
 JaccardSimilarity <- function(set1, set2) {
   intersect_length <- length(intersect(set1, set2))
@@ -113,6 +113,20 @@ if (!"seurat_clusters" %in% colnames(seurat_obj@meta.data)) {
   stop("Input Seurat object is missing required metadata column: seurat_clusters", call. = FALSE)
 }
 
+# Each replicate re-runs SCTransform/PCA/clustering from raw counts and only needs the
+# original cluster labels for the Jaccard comparison. The input also carries a full SCT
+# assay (a dense scale.data), PCA/UMAP embeddings and neighbor graphs from upstream
+# clustering; those would be pinned for the whole loop and re-copied into every subset(),
+# yet are never used here. Strip to a minimal counts-only object (measured ~36 -> ~30 GB
+# peak per replicate on a 77k-cell dataset; the gap grows with cell count). Results are
+# unchanged: SCTransform operates on the RNA counts, which are preserved exactly.
+keep_meta <- intersect(c("seurat_clusters", "percent.mt"), colnames(seurat_obj@meta.data))
+seurat_obj <- CreateSeuratObject(
+  counts = GetAssayData(seurat_obj, assay = "RNA", layer = "counts"),
+  meta.data = seurat_obj@meta.data[, keep_meta, drop = FALSE]
+)
+gc(verbose = FALSE)
+
 replicate_results <- vector("list", n_replicates)
 for (replicate in seq_len(n_replicates)) {
   t0 <- Sys.time()
diff --git a/workflow/scripts/emptydrops.R b/workflow/scripts/emptydrops.R
index 5ac7d5d..fb9d331 100644
--- a/workflow/scripts/emptydrops.R
+++ b/workflow/scripts/emptydrops.R
@@ -53,6 +53,7 @@ gc()
 seurat_droputil_filtered <- RenameAssays(seurat_droputil_filtered,originalexp="RNA")
 seurat_droputil_filtered[["percent.mt"]] <- PercentageFeatureSet(seurat_droputil_filtered, pattern = "(?i)^mt-")
 seurat_droputil_filtered <- SCTransform(seurat_droputil_filtered, vars.to.regress = "percent.mt", seed.use = WORKFLOW_SEED, verbose = FALSE)
+require_min_cells_for_pca(seurat_droputil_filtered, context = "emptydrops")
 seurat_droputil_filtered <- RunPCA(seurat_droputil_filtered, seed.use = WORKFLOW_SEED, verbose = FALSE)
 seurat_droputil_filtered <- RunUMAP(seurat_droputil_filtered, dims = 1:30, seed.use = WORKFLOW_SEED)
 seurat_droputil_filtered <- FindNeighbors(seurat_droputil_filtered, dims = 1:30)
diff --git a/workflow/scripts/find_markers.R b/workflow/scripts/find_markers.R
index 573bafb..b2cb18f 100644
--- a/workflow/scripts/find_markers.R
+++ b/workflow/scripts/find_markers.R
@@ -15,18 +15,42 @@ WORKFLOW_SEED <- set_workflow_seed()
 seurat_obj <- readRDS(rds_input)
 Idents(seurat_obj) <- "seurat_clusters"
 
-set.seed(WORKFLOW_SEED)
-markers <- FindMarkers(seurat_obj, ident.1 = cluster_id)
-markers_tbl <- markers %>%
-  as.data.frame() %>%
-  rownames_to_column(var = "genesymbol") %>%
-  as_tibble()
-
 workflow_name <- file_path_sans_ext(basename(output_path))
 workflow_name <- sub("_markergenes_cluster.*$", "", workflow_name)
 
-sig_markers <- markers_tbl %>%
-  filter(p_val_adj <= 0.05) %>%
-  mutate(cluster = cluster_id, workflow = workflow_name)
-
-write_csv(sig_markers, file = output_path)
+# FindMarkers (via ValidateCellGroups) requires >= 3 cells in the cluster.
+# Clustering can legitimately emit micro-clusters (outliers / residual doublets)
+# with 1-2 cells; those cannot yield markers. Skip them gracefully by writing a
+# schema-correct empty file so combine_markers can still concatenate all clusters.
+n_cells_in_cluster <- sum(Idents(seurat_obj) == cluster_id)
+empty_markers <- tibble(
+  genesymbol = character(),
+  p_val = numeric(),
+  avg_log2FC = numeric(),
+  pct.1 = numeric(),
+  pct.2 = numeric(),
+  p_val_adj = numeric(),
+  cluster = character(),
+  workflow = character()
+)
+
+if (n_cells_in_cluster < 3) {
+  message(sprintf(
+    "Cluster %s has %d cell(s) (< 3); skipping marker detection and writing empty output.",
+    cluster_id, n_cells_in_cluster
+  ))
+  write_csv(empty_markers, file = output_path)
+} else {
+  set.seed(WORKFLOW_SEED)
+  markers <- FindMarkers(seurat_obj, ident.1 = cluster_id)
+  markers_tbl <- markers %>%
+    as.data.frame() %>%
+    rownames_to_column(var = "genesymbol") %>%
+    as_tibble()
+
+  sig_markers <- markers_tbl %>%
+    filter(p_val_adj <= 0.05) %>%
+    mutate(cluster = cluster_id, workflow = workflow_name)
+
+  write_csv(sig_markers, file = output_path)
+}
diff --git a/workflow/scripts/posthocfilter_mad.R b/workflow/scripts/posthocfilter_mad.R
index eed8c92..2175c81 100644
--- a/workflow/scripts/posthocfilter_mad.R
+++ b/workflow/scripts/posthocfilter_mad.R
@@ -32,9 +32,20 @@ low_umi     <- isOutlier(qc$sum, nmads = 3, type = "lower")
 low_feature <- isOutlier(qc$detected, nmads = 3, type = "lower")
 discard <- high_mito | low_umi | low_feature
 cells_to_keep <- colnames(seurat)[!discard]
+rm(sce)
+gc(verbose = FALSE)
 seurat_filtered<- subset(seurat, cells = cells_to_keep)
+# subset() carries the upstream SCT assay, PCA/UMAP embeddings and neighbor graphs into the
+# filtered object; this rule recomputes them from RNA counts below. DietSeurat drops that
+# baggage (SCT assay, reductions, graphs) while preserving the RNA counts and full metadata.
+# Results are unchanged.
+DefaultAssay(seurat_filtered) <- "RNA"
+seurat_filtered <- DietSeurat(seurat_filtered, assays = "RNA", dimreducs = NULL, graphs = NULL)
+rm(seurat)
+gc(verbose = FALSE)
 
 seurat_filtered[["percent.mt"]] <- PercentageFeatureSet(seurat_filtered, pattern = "(?i)^mt-")
+require_min_cells_for_pca(seurat_filtered, context = "posthocfilter_mad")
 seurat_filtered <- SCTransform(seurat_filtered, vars.to.regress = "percent.mt", seed.use = WORKFLOW_SEED, verbose = FALSE)
 seurat_filtered <- RunPCA(seurat_filtered, seed.use = WORKFLOW_SEED, verbose = FALSE)
 seurat_filtered <- RunUMAP(seurat_filtered, dims = 1:30, seed.use = WORKFLOW_SEED)
diff --git a/workflow/scripts/posthocfilter_threshold.R b/workflow/scripts/posthocfilter_threshold.R
index 1680d5f..4b3562e 100644
--- a/workflow/scripts/posthocfilter_threshold.R
+++ b/workflow/scripts/posthocfilter_threshold.R
@@ -30,6 +30,17 @@ seurat <- readRDS(seurat)
 
 seurat_filtered <- subset(seurat, subset = nFeature_RNA > min_nfeature & nCount_RNA > min_ncount & percent.mt < max_mtdna)
 
+# subset() carries the upstream SCT assay, PCA/UMAP embeddings and neighbor graphs into the
+# filtered object, but this rule recomputes all of them from RNA counts below. DietSeurat drops
+# that baggage (SCT assay, reductions, graphs) while preserving the RNA counts AND the full cell
+# metadata, so the stale results are not held alongside the fresh ones and the saved object
+# keeps its upstream metadata columns. Results are unchanged.
+DefaultAssay(seurat_filtered) <- "RNA"
+seurat_filtered <- DietSeurat(seurat_filtered, assays = "RNA", dimreducs = NULL, graphs = NULL)
+rm(seurat)
+gc(verbose = FALSE)
+
+require_min_cells_for_pca(seurat_filtered, context = "posthocfilter_threshold")
 seurat_filtered <- SCTransform(seurat_filtered, vars.to.regress = "percent.mt", seed.use = WORKFLOW_SEED, verbose = FALSE)
 seurat_filtered <- RunPCA(seurat_filtered, seed.use = WORKFLOW_SEED, verbose = FALSE)
 seurat_filtered <- RunUMAP(seurat_filtered, dims = 1:30, seed.use = WORKFLOW_SEED)
diff --git a/workflow/scripts/quarantine_low_quality_samples.py b/workflow/scripts/quarantine_low_quality_samples.py
new file mode 100644
index 0000000..d6c01a6
--- /dev/null
+++ b/workflow/scripts/quarantine_low_quality_samples.py
@@ -0,0 +1,189 @@
+#!/usr/bin/env python3
+"""Quarantine outputs of samples flagged as low-quality (too few cells to cluster).
+
+The require_min_cells_for_pca() guard in workflow/scripts/silhouette_utils.R writes a
+stable "[LOW_QUALITY_SAMPLE]" token into a job's log when QC/filtering leaves too few
+cells to run PCA/clustering for a sample. This script scans the run logs for that token,
+identifies the affected samples, and moves ALL of each flagged sample's output files from
+results/<subdir>/ into results/low_quality_samples/<subdir>/, preserving the subdirectory
+structure.
+
+Job logs are deliberately left in results/logs/ (they explain why a sample was flagged,
+and keeping them makes detection stable across re-runs). The move is per-sample: if a
+sample is flagged in any branch/stage, every output carrying that sample name is moved.
+
+The flagged sample IDs are also written to <results>/<dest>/flagged_samples.txt as an
+ADVISORY record. This file does NOT drive the workflow. To actually skip a sample on
+future runs, add it to `excluded_samples` in config/config.yaml (an explicit, visible
+choice) - this script never edits the config, so a re-sequenced library reusing a sample
+ID is never silently skipped.
+
+Intended to run on completion of the Snakemake workflow (wired into the runner script).
+Safe to run repeatedly; a no-op when nothing is flagged.
+"""
+from __future__ import annotations
+
+import argparse
+import re
+import shutil
+import subprocess
+import sys
+from pathlib import Path
+
+# Must match the token emitted by require_min_cells_for_pca() in silhouette_utils.R.
+TOKEN = "[LOW_QUALITY_SAMPLE]"
+
+# Top-level results subdirectories that are never scanned as outputs and never moved.
+SKIP_DIRS = {"logs", "low_quality_samples", "low_quality_flags"}
+
+# Log subdirectories skipped when scanning for the token, purely for speed. The guard runs
+# only in the RunPCA scripts, so find_markers logs (the overwhelming majority - one per
+# cluster per prefix) can never carry the token. Skipping them is safe: at worst, if these
+# logs ever move, the scan just gets slower - it can never miss a real flag.
+SKIP_LOG_DIRS = {"markers"}
+
+
+def read_sample_ids(samplesheet: Path) -> list[str]:
+    """Sample IDs from the first (tab-separated) column of the sample sheet, minus header."""
+    ids: list[str] = []
+    with samplesheet.open() as fh:
+        next(fh, None)  # skip header row
+        for line in fh:
+            if not line.strip():
+                continue
+            ids.append(line.split("\t")[0].strip())
+    return ids
+
+
+def name_matches_sample(name: str, sid: str) -> bool:
+    """True if a file/dir name carries the sample ID as a whole token.
+
+    Sample names appear bounded by start/underscore/slash on the left and by
+    underscore/dot/end on the right, e.g. "..._cteleta.rds", "cteleta_..._matrix",
+    "cellbender_cteleta.h5". Bounded matching prevents one sample ID from matching
+    inside another (e.g. sc108 vs sc109, obOLsample1 vs obOLsample2).
+    """
+    return re.search(rf"(?:^|[_/]){re.escape(sid)}(?:[_.]|$)", name) is not None
+
+
+def _token_log_names(logs_dir: Path) -> list[str]:
+    """Basenames of log files containing the token. Uses grep (fast over a large, networked
+    log tree - this workflow can produce tens of thousands of logs) and falls back to a
+    pure-Python scan if grep is unavailable."""
+    exclude_args = [f"--exclude-dir={d}" for d in SKIP_LOG_DIRS]
+    try:
+        result = subprocess.run(
+            ["grep", "-rlF", *exclude_args, TOKEN, str(logs_dir)],
+            capture_output=True, text=True, check=False,
+        )
+        # rc 0 = matches, 1 = no matches; anything else -> fall back.
+        if result.returncode in (0, 1):
+            return [Path(line).name for line in result.stdout.splitlines() if line]
+    except (FileNotFoundError, OSError):
+        pass
+    names: list[str] = []
+    for log_path in logs_dir.rglob("*.log"):
+        if any(part in SKIP_LOG_DIRS for part in log_path.relative_to(logs_dir).parts):
+            continue
+        try:
+            if TOKEN in log_path.read_text(errors="replace"):
+                names.append(log_path.name)
+        except OSError:
+            continue
+    return names
+
+
+def find_flagged_samples(results_dir: Path, sample_ids: list[str]) -> set[str]:
+    logs_dir = results_dir / "logs"
+    flagged: set[str] = set()
+    if not logs_dir.is_dir():
+        return flagged
+    # Longest IDs first so the most specific sample name wins the attribution.
+    ordered = sorted(sample_ids, key=len, reverse=True)
+    for name in _token_log_names(logs_dir):
+        for sid in ordered:
+            if name_matches_sample(name, sid):
+                flagged.add(sid)
+                break
+    return flagged
+
+
+def quarantine(results_dir: Path, flagged: set[str], dest_name: str,
+               dry_run: bool) -> list[tuple[Path, Path]]:
+    dest_root = results_dir / dest_name
+    moves: list[tuple[Path, Path]] = []
+    for sub in sorted(p for p in results_dir.iterdir() if p.is_dir()):
+        if sub.name in SKIP_DIRS:
+            continue
+        for item in sorted(sub.iterdir()):
+            if not any(name_matches_sample(item.name, sid) for sid in flagged):
+                continue
+            dest_dir = dest_root / sub.name
+            dest = dest_dir / item.name
+            moves.append((item, dest))
+            if dry_run:
+                continue
+            dest_dir.mkdir(parents=True, exist_ok=True)
+            if dest.is_dir():
+                shutil.rmtree(dest)
+            elif dest.exists():
+                dest.unlink()
+            shutil.move(str(item), str(dest))
+    return moves
+
+
+def main(argv=None) -> int:
+    ap = argparse.ArgumentParser(description="Quarantine low-quality-sample outputs.")
+    ap.add_argument("--results-dir", default="results", type=Path,
+                    help="Workflow results directory (default: results).")
+    ap.add_argument("--samplesheet", default="samplesheet.tsv", type=Path,
+                    help="Sample sheet whose first column lists sample IDs (default: samplesheet.tsv).")
+    ap.add_argument("--dest-name", default="low_quality_samples",
+                    help="Subdirectory of results/ to move flagged outputs into.")
+    ap.add_argument("--dry-run", action="store_true",
+                    help="Report what would move without moving anything.")
+    args = ap.parse_args(argv)
+
+    if not args.results_dir.is_dir():
+        print(f"[quarantine] results dir not found: {args.results_dir}; nothing to do.")
+        return 0
+    if not args.samplesheet.is_file():
+        print(f"[quarantine] sample sheet not found: {args.samplesheet}; skipping quarantine.",
+              file=sys.stderr)
+        return 0
+
+    sample_ids = read_sample_ids(args.samplesheet)
+    flagged = find_flagged_samples(args.results_dir, sample_ids)
+    if not flagged:
+        print("[quarantine] no low-quality samples flagged; nothing to move.")
+        return 0
+
+    flagged_sorted = sorted(flagged)
+    print(f"[quarantine] flagged low-quality sample(s): {', '.join(flagged_sorted)}")
+
+    # Write the flagged manifest BEFORE moving anything. It records what was detected and is
+    # also read by verify_run_complete.py to excuse these samples from the run's pass/fail
+    # decision; writing it first keeps that decision correct even if a later move errors.
+    dest_root = args.results_dir / args.dest_name
+    manifest = dest_root / "flagged_samples.txt"
+    if not args.dry_run:
+        dest_root.mkdir(parents=True, exist_ok=True)
+        manifest.write_text("\n".join(flagged_sorted) + "\n")
+
+    moves = quarantine(args.results_dir, flagged, args.dest_name, args.dry_run)
+    verb = "would move" if args.dry_run else "moved"
+    for src, dest in moves:
+        print(f"[quarantine] {verb}: {src} -> {dest}")
+    print(f"[quarantine] {verb} {len(moves)} item(s) for {len(flagged)} flagged sample(s) "
+          f"into {args.results_dir / args.dest_name}/")
+
+    # The manifest is advisory for the DAG: it does NOT exclude samples. To skip these on
+    # future runs, add them to excluded_samples in config/config.yaml (an explicit choice).
+    print(f"[quarantine] advisory list {'would be ' if args.dry_run else ''}written to {manifest}")
+    print("[quarantine] to skip these samples on future runs, add them to "
+          "'excluded_samples' in config/config.yaml")
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())
diff --git a/workflow/scripts/scdblfinder.R b/workflow/scripts/scdblfinder.R
index e4bd4fa..81e22ed 100644
--- a/workflow/scripts/scdblfinder.R
+++ b/workflow/scripts/scdblfinder.R
@@ -34,6 +34,7 @@ rm(seurat)
 gc()
 seurat_singlets[["percent.mt"]] <- PercentageFeatureSet(seurat_singlets, pattern = "(?i)^mt-")
 seurat_singlets <- SCTransform(seurat_singlets, vars.to.regress = "percent.mt", seed.use = WORKFLOW_SEED, verbose = FALSE)
+require_min_cells_for_pca(seurat_singlets, context = "scdblfinder")
 seurat_singlets <- RunPCA(seurat_singlets, seed.use = WORKFLOW_SEED, verbose = FALSE)
 seurat_singlets <- RunUMAP(seurat_singlets, dims = 1:30, seed.use = WORKFLOW_SEED)
 seurat_singlets <- FindNeighbors(seurat_singlets, dims = 1:30)
diff --git a/workflow/scripts/silhouette_utils.R b/workflow/scripts/silhouette_utils.R
index 9bd1a4a..da7527e 100644
--- a/workflow/scripts/silhouette_utils.R
+++ b/workflow/scripts/silhouette_utils.R
@@ -7,6 +7,32 @@ set_workflow_seed <- function(seed = Sys.getenv("SCRNASEQ_PREPROCESS_SEED", "123
   seed
 }
 
+# Guard against too-few-cells before RunPCA. PCA (npcs, default 50) and the
+# downstream UMAP/neighbor graph require more cells than principal components;
+# when QC/filtering removes nearly all cells (e.g. a very low-depth sample that
+# fails a fixed count/feature threshold) RunPCA aborts with a cryptic
+# "max(nu, nv) must be strictly less than min(nrow(A), ncol(A))" SVD error.
+# Call this immediately before RunPCA to fail with an actionable message instead.
+# The message embeds the stable token "[LOW_QUALITY_SAMPLE]" so the post-workflow
+# quarantine step (workflow/scripts/quarantine_low_quality_samples.py) can detect
+# flagged samples from their logs; keep the token in sync with that script.
+require_min_cells_for_pca <- function(seurat_obj, context = "", npcs = 50L) {
+  n_cells <- ncol(seurat_obj)
+  if (n_cells <= npcs) {
+    prefix <- if (nzchar(context)) paste0(context, ": ") else ""
+    stop(sprintf(
+      paste0(
+        "%s[LOW_QUALITY_SAMPLE] only %d cell(s) remain - too few to compute %d principal ",
+        "components (RunPCA and downstream UMAP/clustering require more cells than PCs). ",
+        "This usually means upstream QC/filtering removed nearly all cells for this ",
+        "sample; consider dataset-specific thresholds or excluding this sample."
+      ),
+      prefix, n_cells, npcs
+    ), call. = FALSE)
+  }
+  invisible(n_cells)
+}
+
 add_silhouette_to_metadata <- function(
     seurat_obj,
     cluster_col = "seurat_clusters",
diff --git a/workflow/scripts/soupx.R b/workflow/scripts/soupx.R
index 163fcdb..867faf4 100644
--- a/workflow/scripts/soupx.R
+++ b/workflow/scripts/soupx.R
@@ -33,18 +33,37 @@ seurat_base <- readRDS(seurat_base)
 soup_channel <- SoupX::SoupChannel(tod = raw_matrix,toc=filtered_matrix,
                 is10X = TRUE)
 soup_channel$tod <- raw
-soup_channel <- SoupX::setClusters(soup_channel, 
-                                   clusters = as.factor(Idents(seurat_base)))
-soup_channel <- setDR(soup_channel, 
-                DR=Seurat::Embeddings(seurat_base, "umap"))
+# seurat_base carries a full SCT assay and neighbor graphs, but only its cluster labels and
+# UMAP embedding are needed here. Extract those and drop the object before the memory-heavy
+# autoEstCont/adjustCounts steps so its baggage is not held throughout. Results are unchanged.
+soup_clusters <- as.factor(Idents(seurat_base))
+soup_umap <- Seurat::Embeddings(seurat_base, "umap")
+rm(seurat_base)
+gc()
+soup_channel <- SoupX::setClusters(soup_channel, clusters = soup_clusters)
+soup_channel <- setDR(soup_channel, DR = soup_umap)
 set.seed(WORKFLOW_SEED)
-soup_channel <- autoEstCont(soup_channel)
+# autoEstCont aborts when it estimates an extremely high contamination fraction
+# (> 0.8), treating it as a likely estimation failure. Across many datasets this
+# hard stop kills otherwise-recoverable samples. Fall back to forceAccept = TRUE so
+# the estimated fraction is used and the sample proceeds, with a clear warning.
+soup_channel <- tryCatch(
+  autoEstCont(soup_channel),
+  error = function(e) {
+    message(
+      "autoEstCont failed (", conditionMessage(e),
+      "); retrying with forceAccept = TRUE."
+    )
+    autoEstCont(soup_channel, forceAccept = TRUE)
+  }
+)
 corrected_counts <- adjustCounts(soup_channel,roundToInt=TRUE)
 seurat_soupx <- CreateSeuratObject(counts = corrected_counts)
-rm(filtered_matrix, raw_matrix, seurat_base, soup_channel, corrected_counts)
+rm(filtered_matrix, raw_matrix, soup_channel, corrected_counts)
 gc()
 seurat_soupx[["percent.mt"]] <- PercentageFeatureSet(seurat_soupx, pattern = "(?i)^mt-")
 seurat_soupx <- SCTransform(seurat_soupx, vars.to.regress = "percent.mt", seed.use = WORKFLOW_SEED, verbose = FALSE)
+require_min_cells_for_pca(seurat_soupx, context = "soupx")
 seurat_soupx <- RunPCA(seurat_soupx, seed.use = WORKFLOW_SEED, verbose = FALSE)
 seurat_soupx <- RunUMAP(seurat_soupx, dims = 1:30, seed.use = WORKFLOW_SEED)
 seurat_soupx <- FindNeighbors(seurat_soupx, dims = 1:30)
diff --git a/workflow/scripts/tenx2seuratrds.R b/workflow/scripts/tenx2seuratrds.R
index 747708b..fad2a25 100644
--- a/workflow/scripts/tenx2seuratrds.R
+++ b/workflow/scripts/tenx2seuratrds.R
@@ -26,6 +26,7 @@ filtered_loaded <- Seurat::Read10X(filtered)
 seurat_obj <- CreateSeuratObject(counts = filtered_loaded)
 seurat_obj[["percent.mt"]] <- PercentageFeatureSet(seurat_obj, pattern = "(?i)^mt-")
 seurat_obj <- SCTransform(seurat_obj, vars.to.regress = "percent.mt", seed.use = WORKFLOW_SEED, verbose = FALSE)
+require_min_cells_for_pca(seurat_obj, context = "tenx2seuratrds")
 seurat_obj <- RunPCA(seurat_obj, seed.use = WORKFLOW_SEED, verbose = FALSE)
 seurat_obj <- RunUMAP(seurat_obj, dims = 1:30, seed.use = WORKFLOW_SEED)
 seurat_obj <- FindNeighbors(seurat_obj, dims = 1:30)
diff --git a/workflow/scripts/verify_run_complete.py b/workflow/scripts/verify_run_complete.py
new file mode 100644
index 0000000..705b7ee
--- /dev/null
+++ b/workflow/scripts/verify_run_complete.py
@@ -0,0 +1,94 @@
+#!/usr/bin/env python3
+"""Verify a workflow run produced its required outputs, excusing low-quality samples.
+
+Snakemake's own process exit code can be unreliable with the SLURM executor: a terminal job
+failure (e.g. an OOM that exhausts retries) does not always propagate to a non-zero exit, so
+the runner batch job can show COMPLETED when real work failed. This script gives the runner
+job a correct final state independent of that.
+
+It checks whether every required target (the rule-all inputs, passed via --targets-file)
+exists, and fails ONLY when a required output is missing for a sample that was NOT flagged as
+low quality. Failures confined to flagged low-quality samples - which are legitimately
+quarantined and cannot be diagnosed without running the workflow - are excused, so a run whose
+only problems are low-quality samples still reports success. A missing output for any other
+sample (a real failure) fails the run.
+
+Exit status: 0 if the run is complete for every non-low-quality sample; 1 otherwise. Intended
+to be invoked from the Snakefile onsuccess/onerror handlers, whose sys.exit() propagates this
+status to the runner job.
+"""
+from __future__ import annotations
+
+import argparse
+import re
+import sys
+from pathlib import Path
+
+
+def read_lines(path: Path) -> list[str]:
+    if not path.is_file():
+        return []
+    return [line.strip() for line in path.read_text().splitlines() if line.strip()]
+
+
+def read_sample_ids(samplesheet: Path) -> list[str]:
+    if not samplesheet.is_file():
+        return []
+    ids: list[str] = []
+    with samplesheet.open() as fh:
+        next(fh, None)  # header
+        for line in fh:
+            if line.strip():
+                ids.append(line.split("\t")[0].strip())
+    return ids
+
+
+def sample_of(path: str, sample_ids: list[str]) -> str | None:
+    """Sample ID carried by a target path (bounded token match, longest ID first)."""
+    name = Path(path).name
+    for sid in sorted(sample_ids, key=len, reverse=True):
+        if re.search(rf"(?:^|[_/]){re.escape(sid)}(?:[_.]|$)", name):
+            return sid
+    return None
+
+
+def main(argv=None) -> int:
+    ap = argparse.ArgumentParser(description="Verify run completeness, excusing low-quality samples.")
+    ap.add_argument("--results-dir", default="results", type=Path)
+    ap.add_argument("--samplesheet", default="samplesheet.tsv", type=Path)
+    ap.add_argument("--targets-file", required=True, type=Path,
+                    help="File listing the run's required target paths, one per line.")
+    ap.add_argument("--dest-name", default="low_quality_samples",
+                    help="Subdirectory of results/ holding quarantined low-quality outputs.")
+    args = ap.parse_args(argv)
+
+    targets = read_lines(args.targets_file)
+    if not targets:
+        print("[verify] no required targets to check; treating run as complete.")
+        return 0
+
+    flagged = set(read_lines(args.results_dir / args.dest_name / "flagged_samples.txt"))
+    sample_ids = read_sample_ids(args.samplesheet)
+
+    missing = [t for t in targets if not Path(t).exists()]
+    real_failures = [(t, sample_of(t, sample_ids)) for t in missing
+                     if sample_of(t, sample_ids) not in flagged]
+    excused = len(missing) - len(real_failures)
+
+    if real_failures:
+        print(f"[verify] RUN INCOMPLETE: {len(real_failures)} required output(s) missing for "
+              f"sample(s) not flagged as low quality - marking the run FAILED:", file=sys.stderr)
+        for target, sample in real_failures[:25]:
+            print(f"  [{sample or 'unknown-sample'}] {target}", file=sys.stderr)
+        if len(real_failures) > 25:
+            print(f"  ... and {len(real_failures) - 25} more", file=sys.stderr)
+        return 1
+
+    print(f"[verify] run complete: all required outputs present for non-low-quality samples "
+          f"({len(targets)} target(s) checked; {excused} missing output(s) excused for "
+          f"{len(flagged)} quarantined low-quality sample(s)).")
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())