From 650051148680328686788547ba4bc8df982c1c91 Mon Sep 17 00:00:00 2001 From: Adam Freedman Date: Tue, 4 Aug 2026 12:58:16 -0400 Subject: [PATCH 1/7] Fix OOM failures and data-driven errors in multi-dataset runs A 26-dataset run failed with ~190 errors across 13 rules. Diagnosis: 179 were OUT_OF_MEMORY (mislabeled "FAILED" by the SLURM executor), hitting the largest datasets across nearly every rule; the remaining 3 were genuine data-driven errors. Data guards: - cellbender2seurat.R: CreateSeuratObject(min.features=1) drops the rare denoised-to-zero barcode in CellBender *_filtered.h5 that made log_umi=-Inf abort SCTransform. - find_markers.R: skip clusters with <3 cells (FindMarkers minimum), writing a schema-correct empty file instead of crashing. - combine_markers.R: drop empty marker tables before bind_rows to avoid a logical/character type clash; handle the all-empty case. - soupx.R: fall back to autoEstCont(forceAccept=TRUE) when contamination estimation aborts (e.g. estimated fraction >0.8). Memory right-sizing (from observed peak RSS; MaxRSS undercounts killed jobs): - doubletfinder/_cellbender: scale by input size, floor 48GB (paramSweep peaks ~43x the input rds; 92GB observed for the largest sample). - downsample_clusters: scale by input size, floor 32GB (all replicates run in one job; ~41GB observed). - emptydrops 24->64GB; tenx2seuratrds/soupx/soupx_emptydrops/scdblfinder/ posthocfilter_mad/posthocfilter_threshold 24->48GB; find_markers 12->24GB. Right-sizing attempt-1 baselines also sidesteps the unreliable retry path (observed: ~half of OOM restart intents were never resubmitted). Co-Authored-By: Claude Opus 4.8 --- workflow/rules/doubletfinder.smk | 12 +++++- workflow/rules/downsample_clusters.smk | 5 ++- workflow/rules/emptydrops.smk | 2 +- workflow/rules/markers.smk | 2 +- workflow/rules/posthocfilter_mad.smk | 4 +- workflow/rules/posthocfilter_threshold.smk | 4 +- workflow/rules/scdblfinder.smk | 4 +- workflow/rules/soupx.smk | 2 +- workflow/rules/soupx_emptydrops.smk | 2 +- workflow/rules/tenx2seuratrds.smk | 2 +- workflow/scripts/cellbender2seurat.R | 6 ++- workflow/scripts/combine_markers.R | 11 ++++- workflow/scripts/find_markers.R | 48 ++++++++++++++++------ workflow/scripts/soupx.R | 15 ++++++- 14 files changed, 90 insertions(+), 29 deletions(-) diff --git a/workflow/rules/doubletfinder.smk b/workflow/rules/doubletfinder.smk index 9eb8431..a55d12a 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, 48 * 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, 48 * 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..31f08ec 100644 --- a/workflow/rules/downsample_clusters.smk +++ b/workflow/rules/downsample_clusters.smk @@ -18,7 +18,10 @@ rule downsample_clusters: wildcard_constraints: downsample_target=DOWNSAMPLE_TARGET_REGEX resources: - mem_mb=lambda wildcards, attempt: int(24000 * (2 ** (attempt - 1))), + # One job runs all replicates (subset + SCTransform + recluster) for a target, + # so peak memory tracks the input object size (~19x observed: 2.1 GB rds -> 41 GB). + # Scale by input size (floor 32 GB) so large datasets clear attempt 1. + mem_mb=lambda wildcards, input, attempt: int(max(32000, 24 * input.size_mb) * (2 ** (attempt - 1))), runtime=lambda wildcards, attempt: int(270 * (2 ** (attempt - 1))) shell: """ 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..aa64c8f 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(48000 * (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(48000 * (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..7da40e2 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(48000 * (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(48000 * (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..8301243 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(48000 * (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..0f339df 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(48000 * (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..3efa379 100644 --- a/workflow/scripts/cellbender2seurat.R +++ b/workflow/scripts/cellbender2seurat.R @@ -22,7 +22,11 @@ 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) seurat <- RunPCA(seurat, seed.use = WORKFLOW_SEED, verbose = FALSE) 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/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/soupx.R b/workflow/scripts/soupx.R index 163fcdb..fbce566 100644 --- a/workflow/scripts/soupx.R +++ b/workflow/scripts/soupx.R @@ -38,7 +38,20 @@ soup_channel <- SoupX::setClusters(soup_channel, soup_channel <- setDR(soup_channel, DR=Seurat::Embeddings(seurat_base, "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) From cf698526ff565751872303fc24763bc0267d91b4 Mon Sep 17 00:00:00 2001 From: Adam Freedman Date: Fri, 7 Aug 2026 11:53:01 -0400 Subject: [PATCH 2/7] Harden downsample_clusters, right-size memory, and improve scalability Second debugging round on 26-dataset runs: the prior fixes cut failures from 193 to 82 and reached 99.7%, but downsample_clusters and a few heavy-dataset jobs still failed. Root causes were the serial single-job downsampling design, under-sized memory for the largest (emptyDrops) cell counts, a stray small future limit, and an unguarded PCA on over-filtered samples. downsample_clusters (all 100 replicates run serially in one job): - Resize to a flat 64 GB / 15 h: peak is a single replicate's SCTransform footprint (~42 GB) and wall-time is ~100x per-replicate (worst ~10.2 h); input-size scaling mis-modeled this (small cteleta was among the heaviest). - Fix future.globals.maxSize 2 GB -> 16 GB (every other script already used 16 GB); the 2 GB cap hard-failed SCTransform on large subsamples (5 jobs). - Strip the input to a counts-only object before the loop: the upstream SCT assay / PCA / UMAP / graphs were pinned across all replicates and re-copied into every subset() but never used (~36 -> ~30 GB peak; grows with cells). Memory right-sizing for the largest emptyDrops cell counts: - doubletfinder / _cellbender: input scaling 48x -> 64x (amellifera peaked ~105 GB from DoubletFinder's O(N^2) distance matrix). - cellbender2seurat 24 -> 48 GB; posthocfilter (mad+threshold), soupx, soupx_emptydrops 48 -> 64 GB. Robustness guard: - require_min_cells_for_pca() in silhouette_utils.R, called before RunPCA in all RunPCA scripts. Over-filtered samples (e.g. cteleta: 21 of 54k cells survive the fixed threshold) now stop with an actionable message instead of a cryptic SVD error. Memory efficiency (results-identical): - posthocfilter (mad+threshold): rebuild a counts-only object after subset() so the upstream SCT/PCA/graph baggage is dropped before recomputation; posthocfilter_mad also frees the SingleCellExperiment QC copy early. - soupx: extract cluster labels + UMAP from seurat_base and drop it before autoEstCont/adjustCounts. Scalability docs: - Document DoubletFinder's O(N^2) memory in doubletfinder.R and README, and recommend scDblFinder for >100k-cell datasets. Co-Authored-By: Claude Opus 4.8 --- README.md | 6 ++++++ workflow/rules/cellbender2seurat.smk | 2 +- workflow/rules/doubletfinder.smk | 4 ++-- workflow/rules/downsample_clusters.smk | 14 ++++++++----- workflow/rules/posthocfilter_mad.smk | 4 ++-- workflow/rules/posthocfilter_threshold.smk | 4 ++-- workflow/rules/soupx.smk | 2 +- workflow/rules/soupx_emptydrops.smk | 2 +- workflow/scripts/cellbender2seurat.R | 1 + workflow/scripts/doubletfinder.R | 10 ++++++++++ workflow/scripts/downsample_clusters.R | 16 ++++++++++++++- workflow/scripts/emptydrops.R | 1 + workflow/scripts/posthocfilter_mad.R | 9 +++++++++ workflow/scripts/posthocfilter_threshold.R | 12 +++++++++++ workflow/scripts/scdblfinder.R | 1 + workflow/scripts/silhouette_utils.R | 23 ++++++++++++++++++++++ workflow/scripts/soupx.R | 16 ++++++++++----- workflow/scripts/tenx2seuratrds.R | 1 + 18 files changed, 108 insertions(+), 20 deletions(-) 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/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/doubletfinder.smk b/workflow/rules/doubletfinder.smk index a55d12a..cdb46e0 100644 --- a/workflow/rules/doubletfinder.smk +++ b/workflow/rules/doubletfinder.smk @@ -20,7 +20,7 @@ rule doubletfinder: # 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, 48 * input.size_mb) * (2 ** (attempt - 1))), + 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 @@ -49,7 +49,7 @@ rule doubletfinder_cellbender: # 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, 48 * input.size_mb) * (2 ** (attempt - 1))), + 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 31f08ec..20e336f 100644 --- a/workflow/rules/downsample_clusters.smk +++ b/workflow/rules/downsample_clusters.smk @@ -18,11 +18,15 @@ rule downsample_clusters: wildcard_constraints: downsample_target=DOWNSAMPLE_TARGET_REGEX resources: - # One job runs all replicates (subset + SCTransform + recluster) for a target, - # so peak memory tracks the input object size (~19x observed: 2.1 GB rds -> 41 GB). - # Scale by input size (floor 32 GB) so large datasets clear attempt 1. - mem_mb=lambda wildcards, input, attempt: int(max(32000, 24 * input.size_mb) * (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/posthocfilter_mad.smk b/workflow/rules/posthocfilter_mad.smk index aa64c8f..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(48000 * (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(48000 * (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 7da40e2..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(48000 * (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(48000 * (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/soupx.smk b/workflow/rules/soupx.smk index 8301243..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(48000 * (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 0f339df..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(48000 * (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/scripts/cellbender2seurat.R b/workflow/scripts/cellbender2seurat.R index 3efa379..a800aab 100644 --- a/workflow/scripts/cellbender2seurat.R +++ b/workflow/scripts/cellbender2seurat.R @@ -29,6 +29,7 @@ mat <- Read_CellBender_h5_Mat(cellbender_h5) 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/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/posthocfilter_mad.R b/workflow/scripts/posthocfilter_mad.R index eed8c92..cfb4cc5 100644 --- a/workflow/scripts/posthocfilter_mad.R +++ b/workflow/scripts/posthocfilter_mad.R @@ -32,9 +32,18 @@ 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 (and re-derives percent.mt), +# so rebuild a minimal counts-only object to avoid holding that baggage. Results are unchanged. +seurat_filtered <- CreateSeuratObject(counts = GetAssayData(seurat_filtered, assay = "RNA", layer = "counts")) +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..9444ba3 100644 --- a/workflow/scripts/posthocfilter_threshold.R +++ b/workflow/scripts/posthocfilter_threshold.R @@ -30,6 +30,18 @@ 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. Rebuild +# a minimal counts-only object (preserving percent.mt for the SCTransform regression) so +# that baggage is not held alongside the freshly computed results. Results are unchanged. +seurat_filtered <- CreateSeuratObject( + counts = GetAssayData(seurat_filtered, assay = "RNA", layer = "counts"), + meta.data = seurat_filtered@meta.data[, "percent.mt", drop = FALSE] +) +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/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..60c07be 100644 --- a/workflow/scripts/silhouette_utils.R +++ b/workflow/scripts/silhouette_utils.R @@ -7,6 +7,29 @@ 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. +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( + "%sonly %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 fbce566..867faf4 100644 --- a/workflow/scripts/soupx.R +++ b/workflow/scripts/soupx.R @@ -33,10 +33,15 @@ 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) # autoEstCont aborts when it estimates an extremely high contamination fraction # (> 0.8), treating it as a likely estimation failure. Across many datasets this @@ -54,10 +59,11 @@ soup_channel <- tryCatch( ) 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) From abe0cf1e8988105bf4a6885f932e6f71e674ba21 Mon Sep 17 00:00:00 2001 From: Adam Freedman Date: Mon, 10 Aug 2026 11:52:08 -0400 Subject: [PATCH 3/7] Add low-quality-sample quarantine and config-driven excluded_samples When post-hoc filtering leaves too few cells to cluster, the require_min_cells_for_pca guard stops the job. This adds two decoupled mechanisms to handle such samples cleanly. Quarantine (automatic, advisory): - The guard now embeds a stable "[LOW_QUALITY_SAMPLE]" token in its message. - New workflow/scripts/quarantine_low_quality_samples.py scans run logs for the token, identifies the affected samples, and moves ALL of each flagged sample's outputs from results// into results/low_quality_samples//, preserving structure (job logs stay in results/logs/ for debugging and stable detection). It writes an advisory results/low_quality_samples/flagged_samples.txt and points the user at excluded_samples. It never edits config and never drives the DAG. - Wired into the runner after the snakemake call (runs on completion, preserving the workflow's exit status). Exclusion (explicit, user-controlled): - New excluded_samples list in config/config.yaml (validated in common.smk) drops the listed sample IDs from the DAG, with a visible "[excluded_samples] skipping ..." message at startup and a note for IDs not found in the sample sheet. Also filters downsample targets (for downsample_only mode). - This is the ONLY thing that changes the DAG, so nothing is skipped automatically: a re-sequenced library that reuses a sample ID is processed normally unless the user deliberately lists that ID. Verified by dry-run: excluding 2 samples removes exactly their jobs (2394 -> 2210) with no leaked targets; quarantine moves flagged outputs, writes the advisory list, and leaves other samples untouched; idempotent and a no-op when nothing is flagged. Co-Authored-By: Claude Opus 4.8 --- config/config.yaml | 8 + scrnaseq_preprocess_slurmrunner.sh | 9 + workflow/Snakefile | 32 ++++ workflow/rules/common.smk | 24 +++ .../scripts/quarantine_low_quality_samples.py | 155 ++++++++++++++++++ workflow/scripts/silhouette_utils.R | 7 +- 6 files changed, 233 insertions(+), 2 deletions(-) create mode 100644 workflow/scripts/quarantine_low_quality_samples.py 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/scrnaseq_preprocess_slurmrunner.sh b/scrnaseq_preprocess_slurmrunner.sh index e11056a..7cb49f4 100755 --- a/scrnaseq_preprocess_slurmrunner.sh +++ b/scrnaseq_preprocess_slurmrunner.sh @@ -36,3 +36,12 @@ 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 +snakemake_status=$? + +# On completion (whether or not every job succeeded), move outputs of any samples flagged +# low-quality by the require_min_cells_for_pca guard (too few cells to cluster after +# filtering) into results/low_quality_samples/, preserving the results subdirectory layout. +python workflow/scripts/quarantine_low_quality_samples.py --results-dir results --samplesheet samplesheet.tsv + +# Preserve the workflow's exit status (the quarantine step should not mask a failed run). +exit $snakemake_status diff --git a/workflow/Snakefile b/workflow/Snakefile index bf22523..561ca9a 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 = {} diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk index 0a2d73b..564cd56 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, } diff --git a/workflow/scripts/quarantine_low_quality_samples.py b/workflow/scripts/quarantine_low_quality_samples.py new file mode 100644 index 0000000..07c18e5 --- /dev/null +++ b/workflow/scripts/quarantine_low_quality_samples.py @@ -0,0 +1,155 @@ +#!/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// into results/low_quality_samples//, 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 //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 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"} + + +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 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 log_path in logs_dir.rglob("*.log"): + try: + if TOKEN not in log_path.read_text(errors="replace"): + continue + except OSError: + continue + for sid in ordered: + if name_matches_sample(log_path.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)}") + 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}/") + + # Advisory record only; does NOT drive the DAG. To skip these on future runs, add them + # to excluded_samples in config/config.yaml (an explicit, visible choice). + 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") + 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/silhouette_utils.R b/workflow/scripts/silhouette_utils.R index 60c07be..da7527e 100644 --- a/workflow/scripts/silhouette_utils.R +++ b/workflow/scripts/silhouette_utils.R @@ -13,14 +13,17 @@ set_workflow_seed <- function(seed = Sys.getenv("SCRNASEQ_PREPROCESS_SEED", "123 # 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( - "%sonly %d cell(s) remain - too few to compute %d principal components ", - "(RunPCA and downstream UMAP/clustering require more cells than PCs). ", + "%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." ), From 932e42b98a9689b53cb442fda6f83712b7ac95ba Mon Sep 17 00:00:00 2001 From: Adam Freedman Date: Wed, 12 Aug 2026 12:01:37 -0400 Subject: [PATCH 4/7] Move execution policy to profile and quarantine to Snakefile handlers Previously the SLURM submission throttling lived on the runner command line and the low-quality-sample quarantine ran as a step appended after the snakemake call. Both are workflow behavior, so a run launched with a different runner script silently skipped them (as happened with a personal runner that lacked the quarantine step). Move both into the workflow so any launcher gets them via --workflow-profile / the Snakefile. - profiles/slurm/config.yaml: add keep-going, jobs, max-jobs-per-timespan, max-status-checks-per-second, latency-wait, rerun-incomplete alongside the existing retries. Remove stale downsample_cluster_replicate set-threads/set-resources entries (that rule was consolidated into downsample_clusters). - Snakefile: run the quarantine from BOTH onsuccess and onerror handlers, so it fires on completion whether or not every job succeeded (flagged samples fail the guard, so such runs end via onerror). It is intentionally not a rule - it moves other rules' outputs, which would fight Snakemake's dependency tracking. - scrnaseq_preprocess_slurmrunner.sh: slim to the invocation-level flags only. - quarantine_low_quality_samples.py: scan logs with grep and skip the find_markers log subtree (which cannot carry the guard token), cutting an 18k-file scan that timed out to a few seconds. Falls back to a Python scan if grep is unavailable. Co-Authored-By: Claude Opus 4.8 --- profiles/slurm/config.yaml | 16 ++++--- scrnaseq_preprocess_slurmrunner.sh | 14 ++---- workflow/Snakefile | 21 +++++++++ .../scripts/quarantine_low_quality_samples.py | 43 ++++++++++++++++--- 4 files changed, 71 insertions(+), 23 deletions(-) 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 7cb49f4..ca5458a 100755 --- a/scrnaseq_preprocess_slurmrunner.sh +++ b/scrnaseq_preprocess_slurmrunner.sh @@ -35,13 +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 -snakemake_status=$? - -# On completion (whether or not every job succeeded), move outputs of any samples flagged -# low-quality by the require_min_cells_for_pca guard (too few cells to cluster after -# filtering) into results/low_quality_samples/, preserving the results subdirectory layout. -python workflow/scripts/quarantine_low_quality_samples.py --results-dir results --samplesheet samplesheet.tsv - -# Preserve the workflow's exit status (the quarantine step should not mask a failed run). -exit $snakemake_status +# 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 561ca9a..8c45c3f 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -199,3 +199,24 @@ include: "rules/posthocfilter_mad.smk" include: "rules/posthocfilter_threshold.smk" include: "rules/markers.smk" include: "rules/downsample_clusters.smk" + + +# Post-workflow quarantine of low-quality samples. The require_min_cells_for_pca guard flags +# samples with too few cells to cluster; this moves all of their outputs into +# results/low_quality_samples/. It runs from BOTH handlers so it fires on completion whether +# or not every job succeeded (flagged samples fail the guard, so such runs end via onerror). +# It is intentionally not a rule: it moves other rules' outputs, which a rule cannot do +# without fighting Snakemake's dependency tracking. +_QUARANTINE_CMD = ( + "python workflow/scripts/quarantine_low_quality_samples.py " + f'--results-dir "{RESULTS_DIR}" ' + f"--samplesheet \"{config.get('sampleTable', 'samplesheet.tsv')}\"" +) + + +onsuccess: + shell(_QUARANTINE_CMD) + + +onerror: + shell(_QUARANTINE_CMD) diff --git a/workflow/scripts/quarantine_low_quality_samples.py b/workflow/scripts/quarantine_low_quality_samples.py index 07c18e5..b8ae6eb 100644 --- a/workflow/scripts/quarantine_low_quality_samples.py +++ b/workflow/scripts/quarantine_low_quality_samples.py @@ -26,6 +26,7 @@ import argparse import re import shutil +import subprocess import sys from pathlib import Path @@ -35,6 +36,12 @@ # 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.""" @@ -59,6 +66,33 @@ def name_matches_sample(name: str, sid: str) -> bool: 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() @@ -66,14 +100,9 @@ def find_flagged_samples(results_dir: Path, sample_ids: list[str]) -> set[str]: return flagged # Longest IDs first so the most specific sample name wins the attribution. ordered = sorted(sample_ids, key=len, reverse=True) - for log_path in logs_dir.rglob("*.log"): - try: - if TOKEN not in log_path.read_text(errors="replace"): - continue - except OSError: - continue + for name in _token_log_names(logs_dir): for sid in ordered: - if name_matches_sample(log_path.name, sid): + if name_matches_sample(name, sid): flagged.add(sid) break return flagged From 31ed6722ca1ad3f699479a396f44cceb92061961 Mon Sep 17 00:00:00 2001 From: Adam Freedman Date: Thu, 13 Aug 2026 12:48:30 -0400 Subject: [PATCH 5/7] Give the runner job a correct pass/fail state via a post-run verification With the SLURM executor, Snakemake can exit 0 even when a job fails terminally (e.g. an OOM that exhausts retries), so the runner batch job shows COMPLETED and a real failure escapes detection. This adds a verification that drives the true exit state, independent of Snakemake's own exit code, and honors the rule that low-quality-sample failures are excused. - verify_run_complete.py (new): checks that every required target (the rule-all inputs, passed via a targets file) exists. A missing output is a real failure only if its sample was NOT flagged low quality; failures confined to flagged samples - which are quarantined and cannot be diagnosed without running - are excused. Exit 1 on real failure, else 0. - Snakefile: onsuccess/onerror now call _finalize_run(), which quarantines, verifies, then sys.exit(verify_code). sys.exit() in a handler deterministically sets Snakemake's process exit code in BOTH directions (unlike shell, which can only force failure), so the runner job's COMPLETED/FAILED status reflects reality even when Snakemake would exit 0 on a terminal failure. Being in the handlers, it applies to any launcher. - quarantine: write flagged_samples.txt before moving outputs, so the verification can excuse flagged samples even if a later move errors. Detection is by missing required outputs, so it is executor-exit-code-agnostic and recovered retries (output present) correctly do not count as failures. Verified end-to-end: a hidden missing target for a non-flagged sample forces the process to exit 1 (runner FAILED) while a complete run exits 0, and failures confined to flagged samples are excused. Co-Authored-By: Claude Opus 4.8 --- workflow/Snakefile | 47 +++++++--- .../scripts/quarantine_low_quality_samples.py | 19 ++-- workflow/scripts/verify_run_complete.py | 94 +++++++++++++++++++ 3 files changed, 140 insertions(+), 20 deletions(-) create mode 100644 workflow/scripts/verify_run_complete.py diff --git a/workflow/Snakefile b/workflow/Snakefile index 8c45c3f..e079f3d 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -201,22 +201,43 @@ include: "rules/markers.smk" include: "rules/downsample_clusters.smk" -# Post-workflow quarantine of low-quality samples. The require_min_cells_for_pca guard flags -# samples with too few cells to cluster; this moves all of their outputs into -# results/low_quality_samples/. It runs from BOTH handlers so it fires on completion whether -# or not every job succeeded (flagged samples fail the guard, so such runs end via onerror). -# It is intentionally not a rule: it moves other rules' outputs, which a rule cannot do -# without fighting Snakemake's dependency tracking. -_QUARANTINE_CMD = ( - "python workflow/scripts/quarantine_low_quality_samples.py " - f'--results-dir "{RESULTS_DIR}" ' - f"--samplesheet \"{config.get('sampleTable', 'samplesheet.tsv')}\"" -) +# Post-workflow finalization, run from BOTH handlers so it fires on completion regardless of +# Snakemake's own (SLURM-executor-dependent) exit code: +# 1. Quarantine low-quality samples - move all outputs of samples the require_min_cells_for_pca +# guard flagged (too few cells to cluster) into results/low_quality_samples/, and record +# them in flagged_samples.txt. This is intentionally not a rule: it moves other rules' +# outputs, which a rule cannot do without fighting Snakemake's dependency tracking. +# 2. Verify completeness and DRIVE THE EXIT CODE - fail the run only if a required output is +# missing for a sample that was NOT flagged low quality (a real failure, e.g. an OOM that +# exhausted retries). Failures confined to flagged low-quality samples are excused. +# sys.exit() in a handler deterministically sets Snakemake's process exit code in either +# direction, so the runner batch job's COMPLETED/FAILED status reflects the true state even +# when Snakemake would otherwise exit 0 on a terminal failure. +def _finalize_run(): + import subprocess + import sys as _sys + from pathlib import Path as _Path + + samplesheet = str(config.get("sampleTable", "samplesheet.tsv")) + _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) onsuccess: - shell(_QUARANTINE_CMD) + _finalize_run() onerror: - shell(_QUARANTINE_CMD) + _finalize_run() diff --git a/workflow/scripts/quarantine_low_quality_samples.py b/workflow/scripts/quarantine_low_quality_samples.py index b8ae6eb..d6c01a6 100644 --- a/workflow/scripts/quarantine_low_quality_samples.py +++ b/workflow/scripts/quarantine_low_quality_samples.py @@ -160,6 +160,16 @@ def main(argv=None) -> int: 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: @@ -167,13 +177,8 @@ def main(argv=None) -> int: print(f"[quarantine] {verb} {len(moves)} item(s) for {len(flagged)} flagged sample(s) " f"into {args.results_dir / args.dest_name}/") - # Advisory record only; does NOT drive the DAG. To skip these on future runs, add them - # to excluded_samples in config/config.yaml (an explicit, visible choice). - 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") + # 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") 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()) From a0a9b64d3077cf626a08a95fa4fd17e5bb58d0d9 Mon Sep 17 00:00:00 2001 From: Adam Freedman Date: Fri, 14 Aug 2026 14:30:53 -0400 Subject: [PATCH 6/7] Fix CI: move finalize function to common.smk and gate it to SLURM runs The push failed the pytest and downsample-rule CI jobs for two reasons introduced with the run-verification handlers: - snakemake --lint flagged "mixed rules and functions in same snakefile" because _finalize_run was defined in the Snakefile. Move it into rules/common.smk as finalize_run(). - The onsuccess/onerror finalization verified against ALL_TARGETS on every run, so tests that build a single target (e.g. the marker checkpoint and downsample-rule tests) saw the other targets reported missing and the process forced to exit 1. Gate finalization on workflow.non_local_exec so it runs only under a remote (SLURM) executor - where it is needed to correct that executor's unreliable exit code - and is skipped for local runs (CI, tests, ad-hoc builds), which have reliable exit codes and may build only a subset of targets. Verified locally: snakemake --lint passes, and test_snakemake_lint, test_marker_checkpoint_ expansion, test_cellbender_rule, test_r_rule_execution, and the downsampling rule test all pass. Co-Authored-By: Claude Opus 4.8 --- workflow/Snakefile | 48 ++++++++++----------------------------- workflow/rules/common.smk | 37 ++++++++++++++++++++++++++++++ 2 files changed, 49 insertions(+), 36 deletions(-) diff --git a/workflow/Snakefile b/workflow/Snakefile index e079f3d..30126d4 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -201,43 +201,19 @@ include: "rules/markers.smk" include: "rules/downsample_clusters.smk" -# Post-workflow finalization, run from BOTH handlers so it fires on completion regardless of -# Snakemake's own (SLURM-executor-dependent) exit code: -# 1. Quarantine low-quality samples - move all outputs of samples the require_min_cells_for_pca -# guard flagged (too few cells to cluster) into results/low_quality_samples/, and record -# them in flagged_samples.txt. This is intentionally not a rule: it moves other rules' -# outputs, which a rule cannot do without fighting Snakemake's dependency tracking. -# 2. Verify completeness and DRIVE THE EXIT CODE - fail the run only if a required output is -# missing for a sample that was NOT flagged low quality (a real failure, e.g. an OOM that -# exhausted retries). Failures confined to flagged low-quality samples are excused. -# sys.exit() in a handler deterministically sets Snakemake's process exit code in either -# direction, so the runner batch job's COMPLETED/FAILED status reflects the true state even -# when Snakemake would otherwise exit 0 on a terminal failure. -def _finalize_run(): - import subprocess - import sys as _sys - from pathlib import Path as _Path - - samplesheet = str(config.get("sampleTable", "samplesheet.tsv")) - _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) - - +# 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() + finalize_run(RESULTS_DIR, ALL_TARGETS, str(config.get("sampleTable", "samplesheet.tsv")), + getattr(workflow, "non_local_exec", False)) onerror: - _finalize_run() + finalize_run(RESULTS_DIR, ALL_TARGETS, str(config.get("sampleTable", "samplesheet.tsv")), + getattr(workflow, "non_local_exec", False)) diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk index 564cd56..7eee7b1 100644 --- a/workflow/rules/common.smk +++ b/workflow/rules/common.smk @@ -317,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) From b65102acdff48dc961b8d92bd06b5b3b9008ffe5 Mon Sep 17 00:00:00 2001 From: Adam Freedman Date: Sun, 16 Aug 2026 15:23:28 -0400 Subject: [PATCH 7/7] Preserve cell metadata in posthocfilter memory strip (fix full-workflow CI) The posthocfilter memory optimization rebuilt the filtered object with CreateSeuratObject keeping only percent.mt, which dropped the upstream metadata columns (orig.ident extras, scDblFinder.class, upstream silhouette/purity, etc.). The full-workflow reference comparison checks metadata column names, so the saved posthocfilter objects no longer matched and the pytest --run-workflow CI job failed. Switch both posthocfilter scripts to DietSeurat(assays="RNA", dimreducs=NULL, graphs=NULL) (with DefaultAssay set to RNA first). This still drops the heavy SCT assay, PCA/UMAP embeddings and neighbor graphs before recomputation - the intended memory savings - but keeps the RNA counts and the full cell metadata intact, so the outputs match the reference exactly. Verified locally: pytest tests --run-workflow passes (37 passed, 12 skipped) reusing the prebuilt conda envs. Co-Authored-By: Claude Opus 4.8 --- workflow/scripts/posthocfilter_mad.R | 8 +++++--- workflow/scripts/posthocfilter_threshold.R | 15 +++++++-------- 2 files changed, 12 insertions(+), 11 deletions(-) diff --git a/workflow/scripts/posthocfilter_mad.R b/workflow/scripts/posthocfilter_mad.R index cfb4cc5..2175c81 100644 --- a/workflow/scripts/posthocfilter_mad.R +++ b/workflow/scripts/posthocfilter_mad.R @@ -36,9 +36,11 @@ 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 (and re-derives percent.mt), -# so rebuild a minimal counts-only object to avoid holding that baggage. Results are unchanged. -seurat_filtered <- CreateSeuratObject(counts = GetAssayData(seurat_filtered, assay = "RNA", layer = "counts")) +# 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) diff --git a/workflow/scripts/posthocfilter_threshold.R b/workflow/scripts/posthocfilter_threshold.R index 9444ba3..4b3562e 100644 --- a/workflow/scripts/posthocfilter_threshold.R +++ b/workflow/scripts/posthocfilter_threshold.R @@ -30,14 +30,13 @@ 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. Rebuild -# a minimal counts-only object (preserving percent.mt for the SCTransform regression) so -# that baggage is not held alongside the freshly computed results. Results are unchanged. -seurat_filtered <- CreateSeuratObject( - counts = GetAssayData(seurat_filtered, assay = "RNA", layer = "counts"), - meta.data = seurat_filtered@meta.data[, "percent.mt", drop = FALSE] -) +# 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)