Preflight cbendercheck dev

config/config.yaml

@@ -12,10 +12,29 @@ workflow_seed: 12345
 # ID in samplesheet.tsv, remove that ID from this list so the new data is processed.
 excluded_samples: []
 
+# Preflight quality gate (runs before any job): counts the called cells in each sample's
+# CellRanger filtered_feature_bc_matrix and handles samples with too few to cluster.
+#   preflight_mode: "skip" drops them from the DAG and reports them - re-checked every run, so a
+#     re-sequenced sample with more cells is picked up automatically; "error" hard-stops the run;
+#     "warn" only reports; "off" disables the check.
+# Only the CellRanger filtered count is known preflight; the runtime min-cells guard and the
+# results/low_quality_samples/ quarantine remain the backstop for the emptydrops/cellbender arms.
+preflight_min_cells: 100
+preflight_mode: "skip"
+
 emptydrop_removal_methods: ["tenx","emptydrops"]
 ambient_decon_methods: ["soupx","cellbender_fromraw"]
 doublet_removal_methods: ["doubletfinder", "scdblfinder"]
 posthoc_methods: ["threshold", "mad"]
+
+# CellBender (cellbender_fromraw arm) adaptive re-run. After the initial run, the HTML report's
+# automated assessment is parsed: if the learning curve looks normal, that run is kept; if the
+# report recommends halving the learning rate and cellbender_adaptive_rerun is true, CellBender is
+# re-run at half cellbender_learning_rate and the re-run is kept only if ITS curve looks normal
+# (otherwise the initial run is kept). results/cellbender/cellbender_adaptive_summary.tsv reports
+# which samples needed a re-run and whether it resolved the issue.
+cellbender_learning_rate: 0.0001      # CellBender's default; the re-run uses exactly half of this
+cellbender_adaptive_rerun: true       # set false to always keep the single initial run
 min_nfeature: 200
 min_ncount: 500
 max_mtdna: 5

tests/test_cellbender_adaptive.py

@@ -0,0 +1,225 @@
+"""Adaptive CellBender re-run wrapper (workflow/scripts/cellbender_adaptive_run.py) and the
+run-level summary (cellbender_adaptive_summary.py).
+
+`assess_report` is unit-tested directly. The orchestration is tested by driving `main()` against a
+fake `cellbender` on PATH whose report verdict is scripted per run (initial vs re-run) through
+environment variables, so every branch -- keep initial, re-run and resolve, re-run and fail to
+resolve, adaptive off, and the hard error on an unrecognized report -- is exercised without a GPU."""
+
+import importlib.util
+import os
+from pathlib import Path
+
+import pytest
+
+SCRIPTS = Path(__file__).resolve().parents[1] / "workflow" / "scripts"
+
+
+def _load(module_name):
+    path = SCRIPTS / f"{module_name}.py"
+    spec = importlib.util.spec_from_file_location(module_name, path)
+    module = importlib.util.module_from_spec(spec)
+    spec.loader.exec_module(module)
+    return module
+
+
+run_mod = _load("cellbender_adaptive_run")
+summary_mod = _load("cellbender_adaptive_summary")
+
+
+# A fake `cellbender`: advertises --seed, and on each real run writes the .h5/_filtered.h5 outputs
+# plus a _report.html whose verdict is chosen by a per-invocation counter -- the first real run
+# reads $FAKE_INITIAL_PHRASE, the second reads $FAKE_RERUN_PHRASE (values: normal|rerun|unknown).
+FAKE_CELLBENDER = r"""#!/usr/bin/env bash
+set -euo pipefail
+
+if [[ "${1:-}" == "remove-background" && "${2:-}" == "--help" ]]; then
+    printf '%s\n' "cellbender remove-background" "  --seed INTEGER"
+    exit 0
+fi
+if [[ "${1:-}" != "remove-background" ]]; then
+    echo "unexpected cellbender invocation: $*" >&2
+    exit 2
+fi
+shift
+
+output=""
+while [[ $# -gt 0 ]]; do
+    case "$1" in
+        --output) output="$2"; shift 2 ;;
+        --input|--learning-rate|--seed) shift 2 ;;
+        *) shift ;;
+    esac
+done
+
+n=1
+if [[ -f "${FAKE_COUNTER}" ]]; then n=$(( $(cat "${FAKE_COUNTER}") + 1 )); fi
+echo "${n}" > "${FAKE_COUNTER}"
+if [[ "${n}" == "1" ]]; then phrase="${FAKE_INITIAL_PHRASE:-normal}"; else phrase="${FAKE_RERUN_PHRASE:-normal}"; fi
+
+case "${phrase}" in
+    normal)  summary="This learning curve looks normal." ;;
+    rerun)   summary="Consider re-running with half the current learning rate to compare the results." ;;
+    *)       summary="The assessment is inconclusive here." ;;
+esac
+
+mkdir -p "$(dirname "${output}")"
+printf 'fake raw for run %s\n' "${n}" > "${output}"
+printf 'fake filtered for run %s\n' "${n}" > "${output%.h5}_filtered.h5"
+printf '<html><body><h2>Automated assessment</h2><h2>Summary</h2><p>%s</p></body></html>\n' \
+    "${summary}" > "${output%.h5}_report.html"
+"""
+
+
+def _install_fake(tmp_path, monkeypatch, initial_phrase, rerun_phrase="normal"):
+    fake_bin = tmp_path / "bin"
+    fake_bin.mkdir()
+    fake = fake_bin / "cellbender"
+    fake.write_text(FAKE_CELLBENDER)
+    fake.chmod(0o755)
+    # Prepend so the fake wins as `cellbender` while bash/env/coreutils still resolve.
+    monkeypatch.setenv("PATH", f"{fake_bin}:{os.environ['PATH']}")
+    monkeypatch.setenv("FAKE_COUNTER", str(tmp_path / "counter"))
+    monkeypatch.setenv("FAKE_INITIAL_PHRASE", initial_phrase)
+    monkeypatch.setenv("FAKE_RERUN_PHRASE", rerun_phrase)
+
+
+def _run(tmp_path, adaptive="true"):
+    """Drive main() and return (status_dict, results_dir Path)."""
+    input_h5 = tmp_path / "raw.h5"
+    input_h5.write_text("raw input")
+    results = tmp_path / "results" / "cellbender"
+    run_mod.main([
+        "--input", str(input_h5),
+        "--output-base", str(results / "cellbender_test.h5"),
+        "--output-filtered", str(results / "cellbender_test_filtered.h5"),
+        "--report", str(results / "cellbender_test_report.html"),
+        "--status", str(results / "cellbender_test_adaptive_status.txt"),
+        "--sample", "test",
+        "--learning-rate", "0.0001",
+        "--adaptive", adaptive,
+        "--seed", "12345",
+        "--workdir", str(tmp_path / "scratch"),
+    ])
+    status = {}
+    for line in (results / "cellbender_test_adaptive_status.txt").read_text().splitlines():
+        key, value = line.split("\t", 1)
+        status[key] = value
+    return status, results
+
+
+# --- assess_report unit tests --------------------------------------------------------------------
+
+def test_assess_report_normal_ignores_tags_and_case():
+    html = "<div><h2>Summary</h2><p>This <b>learning</b> curve <i>looks</i> NORMAL.</p></div>"
+    assert run_mod.assess_report(html) == "normal"
+
+
+def test_assess_report_detects_rerun_suggestion():
+    html = "<p>Consider re-running with half the current learning rate to compare.</p>"
+    assert run_mod.assess_report(html) == "rerun"
+
+
+def test_assess_report_unknown_when_neither_phrase_present():
+    assert run_mod.assess_report("<p>Training finished. Elbo converged.</p>") == "unknown"
+
+
+# --- orchestration -------------------------------------------------------------------------------
+
+def test_normal_first_try_keeps_initial(tmp_path, monkeypatch):
+    _install_fake(tmp_path, monkeypatch, initial_phrase="normal")
+    status, results = _run(tmp_path)
+
+    assert status["outcome"] == "NORMAL_FIRST_TRY"
+    assert status["reran"] == "false"
+    assert status["kept_run"] == "initial"
+    assert status["rerun_learning_rate"] == "NA"
+    assert (results / "cellbender_test.h5").read_text() == "fake raw for run 1\n"
+    assert (results / "cellbender_test_report_initial.html").exists()
+    assert not (results / "cellbender_test_report_rerun.html").exists()
+
+
+def test_rerun_resolves_keeps_rerun(tmp_path, monkeypatch):
+    _install_fake(tmp_path, monkeypatch, initial_phrase="rerun", rerun_phrase="normal")
+    status, results = _run(tmp_path)
+
+    assert status["outcome"] == "RERUN_RESOLVED"
+    assert status["reran"] == "true"
+    assert status["kept_run"] == "rerun"
+    assert status["initial_learning_rate"] == "0.0001"
+    assert status["rerun_learning_rate"] == "5e-05"  # exactly half
+    # The kept output is the second (re-run) result, and both reports are archived.
+    assert (results / "cellbender_test.h5").read_text() == "fake raw for run 2\n"
+    assert (results / "cellbender_test_report_initial.html").exists()
+    assert (results / "cellbender_test_report_rerun.html").exists()
+
+
+def test_rerun_does_not_resolve_keeps_initial(tmp_path, monkeypatch):
+    _install_fake(tmp_path, monkeypatch, initial_phrase="rerun", rerun_phrase="rerun")
+    status, results = _run(tmp_path)
+
+    assert status["outcome"] == "RERUN_DID_NOT_RESOLVE"
+    assert status["reran"] == "true"
+    assert status["kept_run"] == "initial"
+    assert status["rerun_verdict"] == "rerun"
+    # Kept output falls back to the initial run even though a re-run happened.
+    assert (results / "cellbender_test.h5").read_text() == "fake raw for run 1\n"
+    assert (results / "cellbender_test_report_rerun.html").exists()
+
+
+def test_adaptive_off_does_not_rerun(tmp_path, monkeypatch):
+    _install_fake(tmp_path, monkeypatch, initial_phrase="rerun", rerun_phrase="normal")
+    status, results = _run(tmp_path, adaptive="false")
+
+    assert status["outcome"] == "RERUN_SUGGESTED_BUT_ADAPTIVE_OFF"
+    assert status["reran"] == "false"
+    assert status["kept_run"] == "initial"
+    assert (results / "cellbender_test.h5").read_text() == "fake raw for run 1\n"
+    assert not (results / "cellbender_test_report_rerun.html").exists()
+
+
+def test_unknown_initial_report_hard_errors(tmp_path, monkeypatch):
+    _install_fake(tmp_path, monkeypatch, initial_phrase="unknown")
+    with pytest.raises(SystemExit) as excinfo:
+        _run(tmp_path)
+    assert "could not classify" in str(excinfo.value)
+
+
+def test_unknown_rerun_report_hard_errors(tmp_path, monkeypatch):
+    _install_fake(tmp_path, monkeypatch, initial_phrase="rerun", rerun_phrase="unknown")
+    with pytest.raises(SystemExit) as excinfo:
+        _run(tmp_path)
+    assert "re-run report" in str(excinfo.value)
+
+
+# --- run-level summary ---------------------------------------------------------------------------
+
+def _write_status(path, **fields):
+    path.write_text("".join(f"{k}\t{v}\n" for k, v in fields.items()))
+    return path
+
+
+def test_summary_aggregates_status_files(tmp_path, capsys):
+    a = _write_status(
+        tmp_path / "a.txt", sample="alpha", outcome="NORMAL_FIRST_TRY",
+        initial_learning_rate="0.0001", initial_verdict="normal", reran="false",
+        rerun_learning_rate="NA", rerun_verdict="NA", kept_run="initial",
+    )
+    b = _write_status(
+        tmp_path / "b.txt", sample="beta", outcome="RERUN_DID_NOT_RESOLVE",
+        initial_learning_rate="0.0001", initial_verdict="rerun", reran="true",
+        rerun_learning_rate="5e-05", rerun_verdict="rerun", kept_run="initial",
+    )
+    out = tmp_path / "summary.tsv"
+    # Pass b before a to confirm the summary sorts rows by sample id.
+    summary_mod.main(["--output", str(out), str(b), str(a)])
+
+    lines = out.read_text().splitlines()
+    assert lines[0].split("\t") == list(summary_mod.COLUMNS)
+    assert lines[1].startswith("alpha\t")
+    assert lines[2].startswith("beta\t")
+
+    printed = capsys.readouterr().out
+    assert "2 sample(s)" in printed
+    assert "RERUN_DID_NOT_RESOLVE: 1" in printed
+    assert "beta" in printed  # named as re-ran and as unresolved

tests/test_cellbender_rule.py

@@ -62,6 +62,8 @@
 printf 'fake raw cellbender output for %s\\n' "${input}" > "${output}"
 filtered="${output%.h5}_filtered.h5"
 printf 'fake filtered cellbender output for %s\\n' "${input}" > "${filtered}"
+report="${output%.h5}_report.html"
+printf '<html><body><h2>Automated assessment</h2><h2>Summary</h2><p>This learning curve looks normal.</p></body></html>\\n' > "${report}"
 """
 
 
@@ -86,6 +88,8 @@ def test_cellbender_rule_uses_cellbender_filtered_output_convention(tmp_path):
     results_dir = tmp_path / "results"
     base_output = results_dir / "cellbender" / "cellbender_test.h5"
     filtered_output = results_dir / "cellbender" / "cellbender_test_filtered.h5"
+    report_output = results_dir / "cellbender" / "cellbender_test_report.html"
+    status_output = results_dir / "cellbender" / "cellbender_test_adaptive_status.txt"
 
     env = os.environ.copy()
     env["PATH"] = f"{fake_bin}{os.pathsep}{env.get('PATH', '')}"
@@ -125,3 +129,9 @@ def test_cellbender_rule_uses_cellbender_filtered_output_convention(tmp_path):
     assert result.returncode == 0, result.stdout + result.stderr
     assert base_output.read_text().startswith("fake raw cellbender output")
     assert filtered_output.read_text().startswith("fake filtered cellbender output")
+    # The adaptive wrapper also produces the report and per-sample status; a normal-looking
+    # learning curve means the single initial run is kept, with no re-run.
+    assert "This learning curve looks normal" in report_output.read_text()
+    status = status_output.read_text()
+    assert "outcome\tNORMAL_FIRST_TRY" in status
+    assert "reran\tfalse" in status

tests/test_conda_container_validation.py

@@ -136,7 +136,14 @@ def test_workflow_container_declarations_are_explicit_and_recognized():
             containers.append((rule_path, uri))
 
     assert containers, "no workflow container declarations found"
-    assert containers == [(root / "workflow/rules/cellbender.smk", CELLBENDER_CONTAINER_URI)]
+    # The `cellbender` rule and the `cellbender_adaptive_summary` rule (which reuses the same
+    # image for a stdlib aggregation) are the only container declarations, and both pin the
+    # identical CellBender digest.
+    cellbender_smk = root / "workflow/rules/cellbender.smk"
+    assert containers == [
+        (cellbender_smk, CELLBENDER_CONTAINER_URI),
+        (cellbender_smk, CELLBENDER_CONTAINER_URI),
+    ]
     for _, uri in containers:
         assert uri.startswith("docker://")
         assert ":" in uri.removeprefix("docker://"), f"container URI is missing a tag: {uri}"

tests/test_config_validation.py

@@ -114,6 +114,10 @@ def test_default_config_has_required_keys_and_valid_values():
     assert isinstance(config["min_ncount"], int) and config["min_ncount"] > 0
     assert isinstance(config["max_mtdna"], (int, float)) and 0 <= config["max_mtdna"] <= 100
 
+    assert config.get("preflight_mode", "skip") in {"off", "warn", "skip", "error"}
+    if "preflight_min_cells" in config:
+        assert isinstance(config["preflight_min_cells"], int) and config["preflight_min_cells"] > 0
+
 
 @pytest.mark.parametrize(
     "mutate, expected_message, override_results_dir",
@@ -128,6 +132,11 @@ def test_default_config_has_required_keys_and_valid_values():
         (lambda cfg: cfg.update({"resultsDir": ""}), "resultsDir must be a non-empty string", False),
         (lambda cfg: cfg.update({"workflow_mode": "bad_mode"}), "workflow_mode must be one of", True),
         (lambda cfg: cfg.update({"workflow_mode": "downsample_only", "downsampleRate": 1.5}), "downsampleRate must be > 0 and <= 1", True),
+        (lambda cfg: cfg.update({"preflight_min_cells": 0}), "preflight_min_cells must be a positive integer", True),
+        (lambda cfg: cfg.update({"preflight_mode": "bogus"}), "preflight_mode must be one of", True),
+        (lambda cfg: cfg.update({"cellbender_learning_rate": 0}), "cellbender_learning_rate must be a positive number", True),
+        (lambda cfg: cfg.update({"cellbender_learning_rate": "fast"}), "cellbender_learning_rate must be a positive number", True),
+        (lambda cfg: cfg.update({"cellbender_adaptive_rerun": "yes"}), "cellbender_adaptive_rerun must be a boolean", True),
     ],
 )
 def test_invalid_config_fails_early_with_clear_message(tmp_path, mutate, expected_message, override_results_dir):