diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 3236e25..53d31e8 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -43,4 +43,3 @@ jobs: run: snakemake -s workflow/Snakefile -d .test --lint --workflow-profile none - name: Dry-run (decoupled smoke test) run: snakemake -s workflow/Snakefile -d .test --dry-run --workflow-profile none - diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml new file mode 100644 index 0000000..0a8346d --- /dev/null +++ b/.pre-commit-config.yaml @@ -0,0 +1,14 @@ +exclude: ^(demos/|config/.*\.(tsv|csv)|.*\.svg$) +repos: + - repo: https://github.com/snakemake/snakefmt + rev: v2.0.3 # Replace by any tag/version ≥v0.6.0 : https://github.com/snakemake/snakefmt/releases + hooks: + - id: snakefmt + + - repo: https://github.com/pre-commit/pre-commit-hooks + rev: v4.6.0 + hooks: + - id: check-yaml + exclude: ^mkdocs\.yml$ + - id: end-of-file-fixer + - id: trailing-whitespace diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index 5a40831..9e82639 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -28,19 +28,38 @@ workflow/ └── schemas/*.yaml config + sample-sheet validation config/ config.yaml, README.md (config reference), sample sheets docs/ the MkDocs site -.test/ tiny synthetic decoupled fixture for CI +.test/ tiny synthetic decoupled fixture for CI reproduction/ public-data worked examples ``` ## Development setup +Snakemake >= 8 and Conda/Mamba are the only host requirements. We also use snakefmt +and pre-commit hooks to ensure clean diffs and comply with Snakemake linting and formatting +standards. + ```bash git clone https://github.com/cbib/SpaceBlocks && cd SpaceBlocks -# Snakemake >= 8 and Conda/Mamba are the only host requirements. +# Install the development environment +conda create workflow/envs/dev.yaml +conda activate spaceblocks_dev snakemake -n --sdm conda # dry-run: builds the DAG, validates the config, provisions envs snakemake -s workflow/Snakefile -d .test -n --workflow-profile none # decoupled smoke test ``` +Install the pre-commit hooks (`pre-commit` is included in `dev.yaml`) so snakefmt and the +formatting checks run automatically on every commit: + +```bash +pre-commit install +``` + +To run them on demand without committing (e.g. against the whole repo): + +```bash +pre-commit run --all-files +``` + Each rule group has its own Conda env, provisioned by `--sdm conda`. Keep the loose `envs/*.yaml` as the maintainable surface; the `*_linux-64.lock` files are the exact-reproducibility surface (see [Environments](https://cbib.github.io/SpaceBlocks/environments/)). @@ -55,7 +74,7 @@ These are load-bearing, and most past bugs we experienced during development cam `params:`; the script reads `snakemake.params`, **never `config` directly**. After any change, cross-check that every param name matches between the `.smk` and its script. - **Resources scale with retries.** Don't hardcode resources in a rule. Every compute rule draws `mem_mb`/`runtime`/`threads` from - `config["resources"]` (with a `default`), and `mem_mb` grows with the attempt number. + `config["resources"]` (with a `default`), and `mem_mb` grows with the attempt number. - **The contract convention.** `obs["cell_id"]` must equal `obs_names` (as strings); downstream joins key on it. Head-produced contracts live at `SAMPLES_DIR/{sample}/{sample}_unfiltered.h5ad` (nested); decoupled contracts live at `contract_dir/{sample}.h5ad` (flat). @@ -65,7 +84,7 @@ These are load-bearing, and most past bugs we experienced during development cam - **Guard plotting.** Wrap plot generation in `try/except` (Python) / `tryCatch` (R) so one failed figure doesn't crash the rule. - **Palettes.** SpaceBlocks allows color palette customization directly from the config. Levels not referenced in config fall back to grey. - When no palette is configured for a column, leave colours to scanpy and drop any stale `*_colors` from `uns`. + When no palette is configured for a column, leave colours to scanpy and drop any stale `*_colors` from `uns`. ## Adding things @@ -94,7 +113,9 @@ snakemake -s workflow/Snakefile -d .test -n --workflow-profile none python -c "import yaml,jsonschema; jsonschema.validate(yaml.safe_load(open('config/config.yaml')), yaml.safe_load(open('workflow/schemas/config.schema.yaml'))); print('Configuration schema validation passed')" # 4. Docs build cleanly (only if you touched docs/) +# If you did not install the development environment, run: pip install mkdocs-material pymdown-extensions +# Build the docs mkdocs build --strict ``` diff --git a/LICENSE.md b/LICENSE.md index d385772..b58ab93 100644 --- a/LICENSE.md +++ b/LICENSE.md @@ -18,4 +18,4 @@ FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -SOFTWARE. \ No newline at end of file +SOFTWARE. diff --git a/README.md b/README.md index 59baa8a..8bb1edd 100644 --- a/README.md +++ b/README.md @@ -34,7 +34,7 @@ Its modular architecture supports both public and in-house ST datasets while rem The workflow is divided into (optional) technology-specific **HeadBlocks** and common **CoreBlocks** that streamline pre-, post-processing, and informative exploration of results.

SpaceBlocks workflow overview

Workflow overview and partial output from Xenium 5K demo example.

- + A full SpaceBlocks run takes three inputs: 1. **ST formatted AnnData objects**. Either generated from the HeadBlocks, or manually formatted as a standardized h5ad AnnData object. For brevity, **we refer to each of these objects as THE CONTRACT**. Their structure (count matrix, spatial coordinates and optional region annotations) is validated (`validate_input`) before downstream analyses. diff --git a/config/config.yaml b/config/config.yaml index e8c641e..6bc2735 100644 --- a/config/config.yaml +++ b/config/config.yaml @@ -9,7 +9,7 @@ post_processing_outdir: /scratch/CBiB/soterino/smk_visiumhd_results logdir: /scratch/CBiB/soterino/smk_logs spaceranger_processing_outdir: /scratch/CBiB/soterino/smk_spaceranger_results # dir with region annotation GeoJSONs (optional, recommended) -geojson_path: /scratch/CBiB/soterino/geojson +geojson_path: /scratch/CBiB/soterino/geojson # only when using "mode: decoupled" — dir of pre-built contract h5ads contract_dir: '' @@ -128,9 +128,9 @@ analysis: # ─── Annotation ────────────────────────────────────────────────────────────── # (optional, recommended) predefined cell type markers, screened during leiden to ease annotation -snakemake_cell_markers: config/snakemake_cell_markers.tsv +snakemake_cell_markers: config/snakemake_cell_markers.tsv # (default) cluster → cell-type TSV (fill template from generate_annotation_template) -cluster_annotations: config/cluster_annotations.tsv +cluster_annotations: config/cluster_annotations.tsv # (optional) use external tsv with annotations external_annotation: enabled: false diff --git a/docs/getting-started.md b/docs/getting-started.md index a7f3590..7311c42 100644 --- a/docs/getting-started.md +++ b/docs/getting-started.md @@ -87,7 +87,7 @@ The output (see [Outputs](outputs.md) for what every file is, and [the output tr Annotation is the most limiting step in every high-throughput single cell pipeline. It conditions downstream analysis and, thus, it is a breakpoint between CoreBlocks within the SpaceBlocks workflow. -Through the development of SpaceBlocks, we have tested several automatic annotation tools in Visium HD data. The results in our own data had very limited accuracy and, thus, manual annotation is the default option. +Through the development of SpaceBlocks, we have tested several automatic annotation tools in Visium HD data. The results in our own data had very limited accuracy and, thus, manual annotation is the default option. ### 5.1 Manual annotation diff --git a/workflow/Snakefile b/workflow/Snakefile index 73c4e29..9df8b03 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -14,7 +14,10 @@ from snakemake.utils import validate, min_version min_version("8.0") + configfile: "config/config.yaml" + + # --- Public-data examples (mode: decoupled) ------------------------------------- # To run a worked example instead of your own data, comment the line above and # uncomment ONE of these (see docs -> Reproduce & examples): @@ -25,32 +28,34 @@ configfile: "config/config.yaml" validate(config, schema="schemas/config.schema.yaml") + # Helper functions include: "rules/common.smk" + # ── Globals ────────────────────────────────────────────────────────────────── -ANALYSIS = config.get("analysis", {}) +ANALYSIS = config.get("analysis", {}) _region_colors = ANALYSIS.get("region_colors", {}) -REGION_LEVELS = ANALYSIS.get("region_levels", list(_region_colors.keys())) -LOGDIR = config["logdir"] -OUTDIR_SR = config["spaceranger_processing_outdir"] -OUTDIR_PP = config["post_processing_outdir"] -GEOJ_DIR = config["geojson_path"] +REGION_LEVELS = ANALYSIS.get("region_levels", list(_region_colors.keys())) +LOGDIR = config["logdir"] +OUTDIR_SR = config["spaceranger_processing_outdir"] +OUTDIR_PP = config["post_processing_outdir"] +GEOJ_DIR = config["geojson_path"] SAMPLES_DIR = f"{OUTDIR_PP}/Samples" RANDOM_SEED = config.get("random_seed", 42) # Xenium 5K head specifics -XENIUM5K = config.get("xenium5k", {}) or {} +XENIUM5K = config.get("xenium5k", {}) or {} XENIUM_ZARR_DIR = XENIUM5K.get("zarr_dir", "") or OUTDIR_SR # Atera head specifics (ALPHA — the platform ships H2 2026 and its output format is # still a preview; see docs/rules.md). -ATERA = config.get("atera", {}) or {} +ATERA = config.get("atera", {}) or {} ATERA_ZARR_DIR = ATERA.get("zarr_dir", "") or OUTDIR_SR # The registered H&E is OPTIONAL and ships separately from the outs/ bundle. Both the # image and its alignment matrix are needed for the transform to be defined, so the # H&E QuPath rule is included only when both are configured. -ATERA_HAS_HE = bool(ATERA.get("he_image")) and bool(ATERA.get("he_alignment")) +ATERA_HAS_HE = bool(ATERA.get("he_image")) and bool(ATERA.get("he_alignment")) # MERSCOPE (Vizgen MERFISH) head specifics MERSCOPE = config.get("merscope", {}) or {} @@ -61,15 +66,19 @@ DEFAULT_MER_CHANNELS = ["DAPI", "PolyT", "Cellbound1", "Cellbound2", "Cellbound3 MODE = str(config.get("mode", "visiumhd")).lower() KNOWN_HEADS = ("visiumhd", "xenium5k", "atera", "merscope") if MODE != "decoupled" and MODE not in KNOWN_HEADS: - sys.exit(f"[config error] mode must be 'decoupled' or a head technology " - f"{list(KNOWN_HEADS)}, got '{MODE}'.") -IS_DECOUPLED = (MODE == "decoupled") + sys.exit( + f"[config error] mode must be 'decoupled' or a head technology " + f"{list(KNOWN_HEADS)}, got '{MODE}'." + ) +IS_DECOUPLED = MODE == "decoupled" # A half-configured H&E is a silent no-op otherwise: the rule would be skipped and the # user would annotate an image the pipeline never reads. if MODE == "atera" and bool(ATERA.get("he_image")) != bool(ATERA.get("he_alignment")): - sys.exit("[config error] atera.he_image and atera.he_alignment must be set together " - "(or both left empty to annotate the morphology composite instead).") + sys.exit( + "[config error] atera.he_image and atera.he_alignment must be set together " + "(or both left empty to annotate the morphology composite instead)." + ) # ── Input contract (modular head → core hand-off) ──────────────────────────── config.setdefault("contract", {}) @@ -78,8 +87,12 @@ if IS_DECOUPLED: config["contract"]["unfiltered_h5ad"] = f"{CONTRACT_DIR}/{{sample}}.h5ad" else: CONTRACT_DIR = config.get("contract_dir", "") or SAMPLES_DIR - config["contract"]["unfiltered_h5ad"] = f"{CONTRACT_DIR}/{{sample}}/{{sample}}_unfiltered.h5ad" -config["outdir"] = SAMPLES_DIR # validate_input / qc_sweep write under the per-sample tree + config["contract"][ + "unfiltered_h5ad" + ] = f"{CONTRACT_DIR}/{{sample}}/{{sample}}_unfiltered.h5ad" +config["outdir"] = ( + SAMPLES_DIR # validate_input / qc_sweep write under the per-sample tree +) # ── Sample sheets ──────────────────────────────────────────────────────────── CORE_SAMPLES = None @@ -89,7 +102,9 @@ if _core_sheet: sys.exit(f"[config error] core_samples file not found: {_core_sheet}") CORE_SAMPLES = pd.read_csv(_core_sheet, sep="\t", dtype=str, comment="#").fillna("") if "sample" not in CORE_SAMPLES.columns: - sys.exit(f"[config error] core_samples ({_core_sheet}) needs a 'sample' column.") + sys.exit( + f"[config error] core_samples ({_core_sheet}) needs a 'sample' column." + ) validate(CORE_SAMPLES, schema="schemas/core_samples.schema.yaml") CORE_SAMPLES = CORE_SAMPLES.set_index("sample", drop=False) SAMPLE_IDS = list(CORE_SAMPLES["sample"]) @@ -101,7 +116,7 @@ SAMPLES = {} _head_sheet = config.get("samples", "") or "" if os.path.isfile(_head_sheet): samples_df = ( - pd.read_csv(_head_sheet, dtype=str) # keep numeric sample_ids as strings + pd.read_csv(_head_sheet, dtype=str) # keep numeric sample_ids as strings .fillna("") .set_index("sample_id", drop=False) ) @@ -111,15 +126,17 @@ if os.path.isfile(_head_sheet): # Derive the sample list from the head sheet only when no core sheet is set if CORE_SAMPLES is None: if not SAMPLES: - sys.exit("[config error] No core_samples sheet configured and no readable " - "visiumhd_samples.csv found — cannot determine the sample list.") + sys.exit( + "[config error] No core_samples sheet configured and no readable " + "visiumhd_samples.csv found — cannot determine the sample list." + ) SAMPLE_IDS = list(SAMPLES.keys()) # Optional per-column colour palettes (sample + design columns) for downstream visualization SAMPLE_COLORS = config.get("sample_colors", {}) or {} # Integration key on integrate_samples -INTEGRATION = config.get("integration", {}) or {} +INTEGRATION = config.get("integration", {}) or {} INTEGRATE_KEY = str(INTEGRATION.get("integrate_key", "sample") or "sample") # Extra-annotation columns (core-sheet columns beyond `sample`) surfaced to downstream plots/annotation tracks @@ -136,8 +153,11 @@ EXTRA_ANNO_COLS, EXTRA_ANNO_VALS, EXTRA_ANNO_COLORS = _extra_annot_color_triples # ── Contract presence gate (decoupled mode) ────────────────────────────────── if IS_DECOUPLED: _pat = config["contract"]["unfiltered_h5ad"] - _missing = [(s, _pat.format(sample=s)) for s in SAMPLE_IDS - if not os.path.isfile(_pat.format(sample=s))] + _missing = [ + (s, _pat.format(sample=s)) + for s in SAMPLE_IDS + if not os.path.isfile(_pat.format(sample=s)) + ] if _missing: _lst = "\n".join(f" - {s}: {p}" for s, p in _missing) sys.exit( @@ -145,7 +165,8 @@ if IS_DECOUPLED: f"sample(s) in the input directory:\n{_lst}\n" "Every sample in the core sheet must already have a contract present as " "/.h5ad (mode: decoupled does not build them). Fix the " - "sheet or contract_dir, or set mode to a head technology (e.g. visiumhd).") + "sheet or contract_dir, or set mode to a head technology (e.g. visiumhd)." + ) RESOLUTIONS = [ str(r) @@ -157,13 +178,15 @@ RESOLUTIONS = [ ] # Primary/external annotation -EXTERNAL_ENABLED = bool((config.get("external_annotation", {}) or {}).get("enabled", False)) +EXTERNAL_ENABLED = bool( + (config.get("external_annotation", {}) or {}).get("enabled", False) +) PRIMARY_ANNOT_TYPE = "external_annotation" if EXTERNAL_ENABLED else "tsv_annotation" _ANNOT_TYPE_TO_COL = { - "tsv_annotation": "cell_type_tsv", + "tsv_annotation": "cell_type_tsv", "external_annotation": "cell_type_external", - "ingest_annotation": "cell_type_ingest", - "refined_annotation": "cell_type_refined", + "ingest_annotation": "cell_type_ingest", + "refined_annotation": "cell_type_refined", } DEFAULT_ANNOT_COL = _ANNOT_TYPE_TO_COL.get(PRIMARY_ANNOT_TYPE, "cell_type_tsv") ANNOT_TYPES = list(config.get("annotation_types") or [PRIMARY_ANNOT_TYPE]) @@ -180,9 +203,11 @@ check_external_annotation() # Consistency: analysing the external annotation downstream requires it to be produced. if "external_annotation" in ANNOT_TYPES and not EXTERNAL_ENABLED: - sys.exit("[config error] annotation_types includes 'external_annotation' but " - "external_annotation.enabled is false — cell_type_external would be absent " - "downstream. Enable external_annotation or drop it from annotation_types.") + sys.exit( + "[config error] annotation_types includes 'external_annotation' but " + "external_annotation.enabled is false — cell_type_external would be absent " + "downstream. Enable external_annotation or drop it from annotation_types." + ) HAS_CLUSTER_ANNOT = _cluster_annotations_available() HAS_INGEST_REF = bool(config.get("ingest_ref", "")) @@ -200,23 +225,33 @@ SUBCOMPARTMENTS = list(config.get("subcompartments", {}).keys()) # ── Rules ──────────────────────────────────────────────────────────────────── # Headblock (technology-specific): builds the contract h5ad. if MODE == "visiumhd": + include: "rules/spaceranger_count_vhd.smk" include: "rules/generate_qupath_vhd.smk" include: "rules/prepare_input_vhd.smk" + elif MODE == "xenium5k": + include: "rules/convert_zarr_x5k.smk" include: "rules/generate_qupath_x5k.smk" include: "rules/prepare_input_x5k.smk" + elif MODE == "atera": + include: "rules/convert_zarr_ate.smk" include: "rules/generate_qupath_ate.smk" + if ATERA_HAS_HE: + include: "rules/generate_qupath_he_ate.smk" include: "rules/prepare_input_ate.smk" + elif MODE == "merscope": + include: "rules/generate_qupath_mer.smk" include: "rules/prepare_input_mer.smk" + # Coreblock — step 1: preprocessing. include: "rules/validate_input.smk" include: "rules/qc_sweep.smk" @@ -225,7 +260,6 @@ include: "rules/leiden_analysis.smk" include: "rules/generate_annotation_template.smk" include: "rules/ingest_ref.smk" include: "rules/spatial_niches.smk" - # Coreblock — step 2: postprocessing. include: "rules/annotate_cells.smk" include: "rules/integrate_samples.smk" @@ -234,25 +268,33 @@ include: "rules/pseudobulk_de.smk" include: "rules/neighbourhood_analysis.smk" include: "rules/subcluster.smk" include: "rules/sample_report.smk" - # Coreblock — step 3: exploration. include: "rules/explore_genes.smk" # ── Targets ────────────────────────────────────────────────────────────────── if MODE == "visiumhd": - QUPATH_IMAGES = expand(rules.generate_qupath_vhd.output.qupath_image, sample=SAMPLE_IDS) + QUPATH_IMAGES = expand( + rules.generate_qupath_vhd.output.qupath_image, sample=SAMPLE_IDS + ) elif MODE == "xenium5k": - QUPATH_IMAGES = expand(rules.generate_qupath_x5k.output.qupath_image, sample=SAMPLE_IDS) + QUPATH_IMAGES = expand( + rules.generate_qupath_x5k.output.qupath_image, sample=SAMPLE_IDS + ) elif MODE == "atera": - QUPATH_IMAGES = expand(rules.generate_qupath_ate.output.qupath_image, sample=SAMPLE_IDS) + QUPATH_IMAGES = expand( + rules.generate_qupath_ate.output.qupath_image, sample=SAMPLE_IDS + ) if ATERA_HAS_HE: # The H&E-on-morphology-grid background is the annotation image (aligned with # the cells; the raw _he.tiff was dropped as a distracting, error-prone output). - QUPATH_IMAGES += expand(rules.generate_qupath_he_ate.output.he_background, - sample=SAMPLE_IDS) + QUPATH_IMAGES += expand( + rules.generate_qupath_he_ate.output.he_background, sample=SAMPLE_IDS + ) elif MODE == "merscope": - QUPATH_IMAGES = expand(rules.generate_qupath_mer.output.qupath_image, sample=SAMPLE_IDS) + QUPATH_IMAGES = expand( + rules.generate_qupath_mer.output.qupath_image, sample=SAMPLE_IDS + ) else: QUPATH_IMAGES = [] @@ -262,6 +304,7 @@ rule all: get_all_targets, rules.generate_annotation_template.output.template, + rule run_preprocessing: """STEP 1 — validate, preprocess, leiden, annotation template, ingest, niches.""" input: @@ -269,50 +312,85 @@ rule run_preprocessing: expand(rules.preprocess_umap.output.metadata, sample=SAMPLE_IDS), expand(rules.preprocess_umap.output.report, sample=SAMPLE_IDS), *QUPATH_IMAGES, - *(expand(rules.leiden_analysis.output.res_dir, sample=SAMPLE_IDS, resolution=RESOLUTIONS) - if RUN_LEIDEN_ANALYSIS else []), + *( + expand( + rules.leiden_analysis.output.res_dir, + sample=SAMPLE_IDS, + resolution=RESOLUTIONS, + ) + if RUN_LEIDEN_ANALYSIS + else [] + ), rules.generate_annotation_template.output.template, - *(expand(rules.ingest_ref.output.adata_ingested, sample=SAMPLE_IDS) - if HAS_INGEST_REF else []), + *( + expand(rules.ingest_ref.output.adata_ingested, sample=SAMPLE_IDS) + if HAS_INGEST_REF + else [] + ), *([rules.spatial_niches.output.concatenated] if SPATIAL_NICHES_ENABLED else []), + rule run_postprocessing: """STEP 2 — annotation + DE + neighbourhood + integration + report.""" input: expand(rules.annotate_cells.output.adata_annot, sample=SAMPLE_IDS), - *(expand(rules.pseudobulk_de.output.results_dir, - annot_type=ANNOT_TYPES, analysis_level=ANALYSIS_LEVELS) if RUN_DE else []), - expand(rules.neighbourhood_analysis.output.results_dir, - sample=SAMPLE_IDS, annot_type=ANNOT_TYPES), + *( + expand( + rules.pseudobulk_de.output.results_dir, + annot_type=ANNOT_TYPES, + analysis_level=ANALYSIS_LEVELS, + ) + if RUN_DE + else [] + ), + expand( + rules.neighbourhood_analysis.output.results_dir, + sample=SAMPLE_IDS, + annot_type=ANNOT_TYPES, + ), rules.integrate_samples.output.concatenated, rules.integrate_samples.output.harmony, rules.integrate_samples.output.sketched, rules.sample_report.output.report, expand(rules.subcluster.output.sub_dir, subcompartment=SUBCOMPARTMENTS), + rule run_exploration: """STEP 3 — subcompartment re-analysis + gene/signature exploration.""" input: - *([rules.explore_genes_integrated.output.ranges] if HAS_GENE_EXPLORATION else []), - *(expand(rules.explore_genes_sample.output.done, sample=SAMPLE_IDS) - if HAS_GENE_EXPLORATION else []), + *( + [rules.explore_genes_integrated.output.ranges] + if HAS_GENE_EXPLORATION + else [] + ), + *( + expand(rules.explore_genes_sample.output.done, sample=SAMPLE_IDS) + if HAS_GENE_EXPLORATION + else [] + ), + rule subcluster_all: """Run all subcompartment analyses.""" input: expand(rules.subcluster.output.sub_dir, subcompartment=SUBCOMPARTMENTS), + rule qupath_images: """HEAD: generate the hires images to annotate in QuPath (run FIRST, annotate, - export each {sample}_tissue_hires_image.geojson into geojson_path, then continue).""" + export each {sample}_tissue_hires_image.geojson into geojson_path, then continue). + + """ input: QUPATH_IMAGES, + rule qc_sweep_all: """Optional pre-filtering QC diagnostic for every sample (no filtering).""" input: expand(rules.qc_sweep.output.violins_png, sample=SAMPLE_IDS), + rule explore_genes: """Generate gene/signature exploration PNGs (integrated + per sample).""" input: diff --git a/workflow/envs/dev.yaml b/workflow/envs/dev.yaml new file mode 100644 index 0000000..77d8837 --- /dev/null +++ b/workflow/envs/dev.yaml @@ -0,0 +1,17 @@ +name: spaceblocks_dev +channels: + - conda-forge + - bioconda + +dependencies: + # Base snakemake environment + - snakemake>=9.13 + + # Formatting and linting tools + - pre_commit + - bioconda::snakefmt + + # Documentation tools + - mkdocs + - pymdown-extensions + - mkdocs-material diff --git a/workflow/envs/visiumhd.yaml b/workflow/envs/visiumhd.yaml index 89646c0..69bb6e8 100644 --- a/workflow/envs/visiumhd.yaml +++ b/workflow/envs/visiumhd.yaml @@ -40,4 +40,4 @@ dependencies: - pip: - pyclustree>=0.4 - pydeseq2>=0.4 - - geosketch>=1.3 \ No newline at end of file + - geosketch>=1.3 diff --git a/workflow/rules/annotate_cells.smk b/workflow/rules/annotate_cells.smk index a12617b..22be00f 100644 --- a/workflow/rules/annotate_cells.smk +++ b/workflow/rules/annotate_cells.smk @@ -11,6 +11,17 @@ rule annotate_cells: output: adata_annot=f"{SAMPLES_DIR}/{{sample}}/adata_{{sample}}_annotated.h5ad", plots_dir=directory(f"{SAMPLES_DIR}/{{sample}}/annotation"), + log: + out=f"{LOGDIR}/annotate_cells/{{sample}}.out", + err=f"{LOGDIR}/annotate_cells/{{sample}}.err", + benchmark: + f"{LOGDIR}/benchmarks/annotate_cells/{{sample}}.tsv" + conda: + "../envs/visiumhd.yaml" + threads: get_resource("annotate_cells", "threads") + resources: + mem_mb=mem_mb_attempt("annotate_cells"), + runtime=get_resource("annotate_cells", "runtime"), params: sample_id=lambda wc: wc.sample, min_cells_per_type=ANALYSIS.get("min_cells_per_type", 15), @@ -22,17 +33,5 @@ rule annotate_cells: region_colors=ANALYSIS.get("region_colors", {}), dpi=ANALYSIS.get("plot_dpi", 300), niche_column=GENE_EXPLORATION.get("niche_column", ""), - log: - out=f"{LOGDIR}/annotate_cells/{{sample}}.out", - err=f"{LOGDIR}/annotate_cells/{{sample}}.err", - benchmark: - f"{LOGDIR}/benchmarks/annotate_cells/{{sample}}.tsv" - conda: - "../envs/visiumhd.yaml" - threads: - get_resource("annotate_cells", "threads") - resources: - mem_mb=mem_mb_attempt("annotate_cells"), - runtime=get_resource("annotate_cells", "runtime"), script: "../scripts/annotate_cells.py" diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk index 0f8e1a4..2eb1183 100644 --- a/workflow/rules/common.smk +++ b/workflow/rules/common.smk @@ -25,9 +25,11 @@ def mem_mb_attempt(rule_name, cap_factor=None): killed for OOM is automatically resubmitted with more memory. Pairs with the global `retries` set in the execution profile. attempt starts at 1 (= base).""" base = get_resource(rule_name, "mem_mb") + def _mem(wildcards, attempt): factor = attempt if cap_factor is None else min(attempt, cap_factor) return int(base * factor) + return _mem @@ -42,12 +44,16 @@ def check_external_annotation(): return col = cfg.get("column", "") if not col: - sys.exit("[config error] external_annotation.enabled is true but 'column' is empty.") + sys.exit( + "[config error] external_annotation.enabled is true but 'column' is empty." + ) meta_dir = config.get("precomputed_metadata_dir", "") or "" if not meta_dir: - sys.exit("[config error] external_annotation.enabled requires 'precomputed_metadata_dir' " - f"to point at a directory of metadata_{{sample}}.tsv files carrying the '{col}' " - "column (the pipeline's own metadata does not contain external labels).") + sys.exit( + "[config error] external_annotation.enabled requires 'precomputed_metadata_dir' " + f"to point at a directory of metadata_{{sample}}.tsv files carrying the '{col}' " + "column (the pipeline's own metadata does not contain external labels)." + ) missing_file, missing_col = [], [] for s in SAMPLE_IDS: f = os.path.join(meta_dir, f"metadata_{s}.tsv") @@ -61,11 +67,17 @@ def check_external_annotation(): if col not in header: missing_col.append(f) if missing_file or missing_col: - msg = [f"[config error] external_annotation enabled (column '{col}') but not in place:"] + msg = [ + f"[config error] external_annotation enabled (column '{col}') but not in place:" + ] if missing_file: - msg.append(" missing metadata file(s):\n " + "\n ".join(missing_file)) + msg.append( + " missing metadata file(s):\n " + "\n ".join(missing_file) + ) if missing_col: - msg.append(f" column '{col}' absent in:\n " + "\n ".join(missing_col)) + msg.append( + f" column '{col}' absent in:\n " + "\n ".join(missing_col) + ) sys.exit("\n".join(msg)) @@ -94,7 +106,9 @@ def _extra_annot_color_triples(): sc_ = SAMPLE_COLORS if isinstance(SAMPLE_COLORS, dict) else {} for c in EXTRA_ANNOT_COLUMNS: for v, hexc in (sc_.get(c, {}) or {}).items(): - cols.append(c); vals.append(str(v)); colors.append(str(hexc)) + cols.append(c) + vals.append(str(v)) + colors.append(str(hexc)) return cols, vals, colors @@ -126,18 +140,23 @@ def xenium_dir_for(sample): {sample} pattern (e.g. 'data/xenium/{sample}').""" pat = (config.get("xenium5k", {}) or {}).get("xenium_dir", "") if not pat: - sys.exit("[config error] mode 'xenium5k' requires xenium5k.xenium_dir " - "(a {sample} pattern to each Xenium output bundle).") + sys.exit( + "[config error] mode 'xenium5k' requires xenium5k.xenium_dir " + "(a {sample} pattern to each Xenium output bundle)." + ) return pat.format(sample=sample) + # ── HEAD (Atera) ───────────────────────────────────────────────────────────── def atera_dir_for(sample): """Atera output bundle directory for a sample, from the atera.atera_dir {sample} pattern (e.g. 'data/atera/{sample}/outs').""" pat = (config.get("atera", {}) or {}).get("atera_dir", "") if not pat: - sys.exit("[config error] mode 'atera' requires atera.atera_dir " - "(a {sample} pattern to each Atera outs/ bundle).") + sys.exit( + "[config error] mode 'atera' requires atera.atera_dir " + "(a {sample} pattern to each Atera outs/ bundle)." + ) return pat.format(sample=sample) @@ -150,8 +169,10 @@ def he_file_for(sample, key, required=True): pat = (config.get("atera", {}) or {}).get(key, "") if not pat: if required: - sys.exit(f"[config error] atera.{key} is required when the optional H&E " - "annotation image is enabled (atera.he_image + atera.he_alignment).") + sys.exit( + f"[config error] atera.{key} is required when the optional H&E " + "annotation image is enabled (atera.he_image + atera.he_alignment)." + ) return "" return pat.format(sample=sample) @@ -163,13 +184,17 @@ def _ate_prepare_inputs(wildcards): inputs = { "done": rules.convert_zarr_ate.output.done.format(sample=wildcards.sample), "qupath_meta": rules.generate_qupath_ate.output.qupath_meta.format( - sample=wildcards.sample), + sample=wildcards.sample + ), } if ATERA_HAS_HE: # Only the embedded background is consumed here; the raw-H&E QuPath image and its # polygon affine were dropped, so he_meta is no longer an input. - inputs["he_background"] = rules.generate_qupath_he_ate.output.he_background.format( - sample=wildcards.sample) + inputs["he_background"] = ( + rules.generate_qupath_he_ate.output.he_background.format( + sample=wildcards.sample + ) + ) # Track the region GeoJSON when present, so editing/renaming it retriggers the # contract build. It is read by filename (regions are optional), not hard-required — # same "declare only when present" pattern as _preprocess_inputs' precomputed_meta. @@ -180,6 +205,7 @@ def _ate_prepare_inputs(wildcards): break return inputs + # ── HEAD (MERSCOPE) input helper ───────────────────────────────────────────── def merscope_dir_for(sample): """MERSCOPE region directory for a sample, from the merscope.merscope_dir @@ -188,8 +214,10 @@ def merscope_dir_for(sample): micron_to_mosaic_pixel_transform.csv).""" pat = (config.get("merscope", {}) or {}).get("merscope_dir", "") if not pat: - sys.exit("[config error] mode 'merscope' requires merscope.merscope_dir " - "(a {sample} pattern to each MERSCOPE region directory).") + sys.exit( + "[config error] mode 'merscope' requires merscope.merscope_dir " + "(a {sample} pattern to each MERSCOPE region directory)." + ) return pat.format(sample=sample) @@ -226,7 +254,9 @@ def _preprocess_inputs(wc): if _ext.get("enabled") and not _ext.get("keep_unannotated", True): _meta_dir = config.get("precomputed_metadata_dir", "") or "" if _meta_dir: - inputs["external_meta"] = os.path.join(_meta_dir, f"metadata_{wc.sample}.tsv") + inputs["external_meta"] = os.path.join( + _meta_dir, f"metadata_{wc.sample}.tsv" + ) return inputs @@ -238,8 +268,9 @@ def _spatial_niches_inputs(wc): rule reloads from its own output dir, so nothing is added here (no circular dependency).""" inputs = { - "adatas": expand(f"{SAMPLES_DIR}/{{sample}}/adata_{{sample}}.h5ad", - sample=SAMPLE_IDS), + "adatas": expand( + f"{SAMPLES_DIR}/{{sample}}/adata_{{sample}}.h5ad", sample=SAMPLE_IDS + ), } sn = config.get("spatial_niches", {}) if sn.get("use_precomputed", False): @@ -279,10 +310,13 @@ def get_all_targets(wildcards): targets += expand(rules.preprocess_umap.output.adata, sample=SAMPLE_IDS) targets += expand(rules.preprocess_umap.output.metadata, sample=SAMPLE_IDS) targets += expand(rules.preprocess_umap.output.report, sample=SAMPLE_IDS) - targets += list(QUPATH_IMAGES) # head QuPath image(s) for the active mode (or []) + targets += list(QUPATH_IMAGES) # head QuPath image(s) for the active mode (or []) if RUN_LEIDEN_ANALYSIS: - targets += expand(rules.leiden_analysis.output.res_dir, - sample=SAMPLE_IDS, resolution=RESOLUTIONS) + targets += expand( + rules.leiden_analysis.output.res_dir, + sample=SAMPLE_IDS, + resolution=RESOLUTIONS, + ) if SPATIAL_NICHES_ENABLED: targets.append(rules.spatial_niches.output.concatenated) targets.append(rules.spatial_niches.output.plots_dir) @@ -291,13 +325,22 @@ def get_all_targets(wildcards): if HAS_INGEST_REF: targets += expand(rules.ingest_ref.output.adata_ingested, sample=SAMPLE_IDS) targets += expand(rules.annotate_cells.output.adata_annot, sample=SAMPLE_IDS) - targets += expand(rules.neighbourhood_analysis.output.results_dir, - sample=SAMPLE_IDS, annot_type=ANNOT_TYPES) - targets += expand(rules.pseudobulk_aggregate.output.agg_dir, - annot_type=ANNOT_TYPES, analysis_level=ANALYSIS_LEVELS) + targets += expand( + rules.neighbourhood_analysis.output.results_dir, + sample=SAMPLE_IDS, + annot_type=ANNOT_TYPES, + ) + targets += expand( + rules.pseudobulk_aggregate.output.agg_dir, + annot_type=ANNOT_TYPES, + analysis_level=ANALYSIS_LEVELS, + ) if RUN_DE: - targets += expand(rules.pseudobulk_de.output.results_dir, - annot_type=ANNOT_TYPES, analysis_level=ANALYSIS_LEVELS) + targets += expand( + rules.pseudobulk_de.output.results_dir, + annot_type=ANNOT_TYPES, + analysis_level=ANALYSIS_LEVELS, + ) targets.append(rules.integrate_samples.output.concatenated) targets.append(rules.integrate_samples.output.harmony) targets.append(rules.integrate_samples.output.sketched) @@ -306,5 +349,7 @@ def get_all_targets(wildcards): targets.append(rules.explore_genes_integrated.output.ranges) targets += expand(rules.explore_genes_sample.output.done, sample=SAMPLE_IDS) if SUBCOMPARTMENTS: - targets += expand(rules.subcluster.output.sub_dir, subcompartment=SUBCOMPARTMENTS) + targets += expand( + rules.subcluster.output.sub_dir, subcompartment=SUBCOMPARTMENTS + ) return targets diff --git a/workflow/rules/convert_zarr_ate.smk b/workflow/rules/convert_zarr_ate.smk index 7813503..da619a7 100644 --- a/workflow/rules/convert_zarr_ate.smk +++ b/workflow/rules/convert_zarr_ate.smk @@ -8,9 +8,6 @@ rule convert_zarr_ate: atera_dir=lambda wc: atera_dir_for(wc.sample), output: done=f"{ATERA_ZARR_DIR}/{{sample}}/{{sample}}.zarr.done", - params: - sample_id=lambda wc: wc.sample, - zarr_path=lambda wc, output: str(output.done)[:-len(".done")], log: out=f"{LOGDIR}/convert_zarr_ate/{{sample}}.out", err=f"{LOGDIR}/convert_zarr_ate/{{sample}}.err", @@ -18,10 +15,12 @@ rule convert_zarr_ate: f"{LOGDIR}/benchmarks/convert_zarr_ate/{{sample}}.tsv" conda: "../envs/atera.yaml" - threads: - get_resource("convert_zarr_ate", "threads") + threads: get_resource("convert_zarr_ate", "threads") resources: mem_mb=mem_mb_attempt("convert_zarr_ate"), runtime=get_resource("convert_zarr_ate", "runtime"), + params: + sample_id=lambda wc: wc.sample, + zarr_path=lambda wc, output: str(output.done)[: -len(".done")], script: "../scripts/convert_zarr_ate.py" diff --git a/workflow/rules/convert_zarr_x5k.smk b/workflow/rules/convert_zarr_x5k.smk index 8e4e2b6..8d7c5c6 100644 --- a/workflow/rules/convert_zarr_x5k.smk +++ b/workflow/rules/convert_zarr_x5k.smk @@ -7,9 +7,6 @@ rule convert_zarr_x5k: xenium_dir=lambda wc: xenium_dir_for(wc.sample), output: done=f"{XENIUM_ZARR_DIR}/{{sample}}/{{sample}}.zarr.done", - params: - sample_id=lambda wc: wc.sample, - zarr_path=lambda wc, output: str(output.done)[:-len(".done")], log: out=f"{LOGDIR}/convert_zarr_x5k/{{sample}}.out", err=f"{LOGDIR}/convert_zarr_x5k/{{sample}}.err", @@ -17,10 +14,12 @@ rule convert_zarr_x5k: f"{LOGDIR}/benchmarks/convert_zarr_x5k/{{sample}}.tsv" conda: "../envs/xenium5k.yaml" - threads: - get_resource("convert_zarr_x5k", "threads") + threads: get_resource("convert_zarr_x5k", "threads") resources: mem_mb=mem_mb_attempt("convert_zarr_x5k"), runtime=get_resource("convert_zarr_x5k", "runtime"), + params: + sample_id=lambda wc: wc.sample, + zarr_path=lambda wc, output: str(output.done)[: -len(".done")], script: "../scripts/convert_zarr_x5k.py" diff --git a/workflow/rules/explore_genes.smk b/workflow/rules/explore_genes.smk index cdd7f68..8915d73 100644 --- a/workflow/rules/explore_genes.smk +++ b/workflow/rules/explore_genes.smk @@ -12,16 +12,6 @@ rule explore_genes_integrated: output: ranges=f"{OUTDIR_PP}/gene_exploration/expression_ranges.tsv", done=touch(f"{OUTDIR_PP}/gene_exploration/.integrated_done"), - params: - outdir=lambda wc, output: os.path.dirname(output.ranges), - annot_key=GENE_EXPLORATION.get("annot_key") or DEFAULT_ANNOT_COL, - aucell_fraction=GENE_EXPLORATION.get("aucell_max_rank_fraction", 0.05), - niche_column=GENE_EXPLORATION.get("niche_column", ""), - dpi=GENE_EXPLORATION.get("dpi", ANALYSIS.get("plot_dpi", 300)), - annotation_colors=config.get("annotation_colors", {}), - region_colors=ANALYSIS.get("region_colors", {}), - extra_annot_columns=EXTRA_ANNOT_COLUMNS, - sample_colors=SAMPLE_COLORS, log: out=f"{LOGDIR}/explore_genes/integrated.out", err=f"{LOGDIR}/explore_genes/integrated.err", @@ -29,11 +19,20 @@ rule explore_genes_integrated: f"{LOGDIR}/benchmarks/explore_genes/integrated.tsv" conda: "../envs/pseudobulk_aggregate.yaml" - threads: - get_resource("explore_genes_integrated", "threads") + threads: get_resource("explore_genes_integrated", "threads") resources: mem_mb=mem_mb_attempt("explore_genes_integrated"), runtime=get_resource("explore_genes_integrated", "runtime"), + params: + outdir=lambda wc, output: os.path.dirname(output.ranges), + annot_key=GENE_EXPLORATION.get("annot_key") or DEFAULT_ANNOT_COL, + aucell_fraction=GENE_EXPLORATION.get("aucell_max_rank_fraction", 0.05), + niche_column=GENE_EXPLORATION.get("niche_column", ""), + dpi=GENE_EXPLORATION.get("dpi", ANALYSIS.get("plot_dpi", 300)), + annotation_colors=config.get("annotation_colors", {}), + region_colors=ANALYSIS.get("region_colors", {}), + extra_annot_columns=EXTRA_ANNOT_COLUMNS, + sample_colors=SAMPLE_COLORS, script: "../scripts/explore_genes_integrated.py" @@ -54,15 +53,6 @@ rule explore_genes_sample: queries=GENE_EXPLORATION.get("queries", "") or [], output: done=touch(f"{OUTDIR_PP}/gene_exploration/.sentinels/{{sample}}.done"), - params: - outdir=lambda wc, output: os.path.dirname(os.path.dirname(output.done)), - sample_id=lambda wc: wc.sample, - annot_key=GENE_EXPLORATION.get("annot_key") or DEFAULT_ANNOT_COL, - niche_column=GENE_EXPLORATION.get("niche_column", ""), - aucell_fraction=GENE_EXPLORATION.get("aucell_max_rank_fraction", 0.05), - dpi=GENE_EXPLORATION.get("dpi", ANALYSIS.get("plot_dpi", 300)), - annotation_colors=config.get("annotation_colors", {}), - region_colors=ANALYSIS.get("region_colors", {}), log: out=f"{LOGDIR}/explore_genes/{{sample}}.out", err=f"{LOGDIR}/explore_genes/{{sample}}.err", @@ -70,10 +60,18 @@ rule explore_genes_sample: f"{LOGDIR}/benchmarks/explore_genes/{{sample}}.tsv" conda: "../envs/pseudobulk_aggregate.yaml" - threads: - get_resource("explore_genes_sample", "threads") + threads: get_resource("explore_genes_sample", "threads") resources: mem_mb=mem_mb_attempt("explore_genes_sample"), runtime=get_resource("explore_genes_sample", "runtime"), + params: + outdir=lambda wc, output: os.path.dirname(os.path.dirname(output.done)), + sample_id=lambda wc: wc.sample, + annot_key=GENE_EXPLORATION.get("annot_key") or DEFAULT_ANNOT_COL, + niche_column=GENE_EXPLORATION.get("niche_column", ""), + aucell_fraction=GENE_EXPLORATION.get("aucell_max_rank_fraction", 0.05), + dpi=GENE_EXPLORATION.get("dpi", ANALYSIS.get("plot_dpi", 300)), + annotation_colors=config.get("annotation_colors", {}), + region_colors=ANALYSIS.get("region_colors", {}), script: "../scripts/explore_genes_sample.py" diff --git a/workflow/rules/generate_annotation_template.smk b/workflow/rules/generate_annotation_template.smk index 7646cf3..3ca3d75 100644 --- a/workflow/rules/generate_annotation_template.smk +++ b/workflow/rules/generate_annotation_template.smk @@ -12,14 +12,14 @@ rule generate_annotation_template: """ output: template=f"{OUTDIR_PP}/cluster_annotations_template.tsv", - params: - sample_ids=SAMPLE_IDS, - resolutions=RESOLUTIONS, - max_clusters=50, log: out=f"{LOGDIR}/generate_annotation_template.out", err=f"{LOGDIR}/generate_annotation_template.err", conda: "../envs/visiumhd.yaml" + params: + sample_ids=SAMPLE_IDS, + resolutions=RESOLUTIONS, + max_clusters=50, script: "../scripts/generate_annotation_template.py" diff --git a/workflow/rules/generate_qupath_ate.smk b/workflow/rules/generate_qupath_ate.smk index 902c9ff..be5b74f 100644 --- a/workflow/rules/generate_qupath_ate.smk +++ b/workflow/rules/generate_qupath_ate.smk @@ -10,11 +10,6 @@ rule generate_qupath_ate: output: qupath_image=f"{SAMPLES_DIR}/{{sample}}/QuPath_image/{{sample}}_morphology.tiff", qupath_meta=f"{SAMPLES_DIR}/{{sample}}/QuPath_image/{{sample}}_morphology_scalefactors.json", - params: - sample_id=lambda wc: wc.sample, - zarr_path=lambda wc, input: str(input.done)[:-len(".done")], - qupath_pyramid_level=ATERA.get("qupath_pyramid_level", 3), - pixel_size_um=ATERA.get("pixel_size_um", 0.2125), log: out=f"{LOGDIR}/generate_qupath_ate/{{sample}}.out", err=f"{LOGDIR}/generate_qupath_ate/{{sample}}.err", @@ -22,10 +17,14 @@ rule generate_qupath_ate: f"{LOGDIR}/benchmarks/generate_qupath_ate/{{sample}}.tsv" conda: "../envs/atera.yaml" - threads: - get_resource("generate_qupath_ate", "threads") + threads: get_resource("generate_qupath_ate", "threads") resources: mem_mb=mem_mb_attempt("generate_qupath_ate"), runtime=get_resource("generate_qupath_ate", "runtime"), + params: + sample_id=lambda wc: wc.sample, + zarr_path=lambda wc, input: str(input.done)[: -len(".done")], + qupath_pyramid_level=ATERA.get("qupath_pyramid_level", 3), + pixel_size_um=ATERA.get("pixel_size_um", 0.2125), script: "../scripts/generate_qupath_ate.py" diff --git a/workflow/rules/generate_qupath_he_ate.smk b/workflow/rules/generate_qupath_he_ate.smk index e303ff8..c843036 100644 --- a/workflow/rules/generate_qupath_he_ate.smk +++ b/workflow/rules/generate_qupath_he_ate.smk @@ -21,16 +21,6 @@ rule generate_qupath_he_ate: output: qupath_meta=f"{SAMPLES_DIR}/{{sample}}/QuPath_image/{{sample}}_he_scalefactors.json", he_background=f"{SAMPLES_DIR}/{{sample}}/QuPath_image/{{sample}}_he_background.tiff", - params: - sample_id=lambda wc: wc.sample, - zarr_path=lambda wc, input: str(input.done)[:-len(".done")], - # Optional: not an input: — absence must not block the rule, it only disables QA. - he_keypoints=lambda wc: he_file_for(wc.sample, "he_keypoints", required=False), - he_pyramid_level=ATERA.get("he_pyramid_level", 4), - # Must match prepare_input_ate's level: the background is warped onto that grid. - hires_pyramid_level=ATERA.get("hires_pyramid_level", 3), - pixel_size_um=ATERA.get("pixel_size_um", 0.2125), - residual_warn_px=ATERA.get("he_residual_warn_px", 50), log: out=f"{LOGDIR}/generate_qupath_he_ate/{{sample}}.out", err=f"{LOGDIR}/generate_qupath_he_ate/{{sample}}.err", @@ -38,10 +28,19 @@ rule generate_qupath_he_ate: f"{LOGDIR}/benchmarks/generate_qupath_he_ate/{{sample}}.tsv" conda: "../envs/atera.yaml" - threads: - get_resource("generate_qupath_he_ate", "threads") + threads: get_resource("generate_qupath_he_ate", "threads") resources: mem_mb=mem_mb_attempt("generate_qupath_he_ate"), runtime=get_resource("generate_qupath_he_ate", "runtime"), + params: + sample_id=lambda wc: wc.sample, + zarr_path=lambda wc, input: str(input.done)[: -len(".done")], + # Optional: not an input: — absence must not block the rule, it only disables QA. + he_keypoints=lambda wc: he_file_for(wc.sample, "he_keypoints", required=False), + he_pyramid_level=ATERA.get("he_pyramid_level", 4), + # Must match prepare_input_ate's level: the background is warped onto that grid. + hires_pyramid_level=ATERA.get("hires_pyramid_level", 3), + pixel_size_um=ATERA.get("pixel_size_um", 0.2125), + residual_warn_px=ATERA.get("he_residual_warn_px", 50), script: "../scripts/generate_qupath_he_ate.py" diff --git a/workflow/rules/generate_qupath_mer.smk b/workflow/rules/generate_qupath_mer.smk index c4e7912..c3162fa 100644 --- a/workflow/rules/generate_qupath_mer.smk +++ b/workflow/rules/generate_qupath_mer.smk @@ -3,21 +3,19 @@ rule generate_qupath_mer: MERSCOPE Headblock — composite a morphology RGB TIFF (+ pixel<->micron JSON and a grey background reused by the contract) from the Vizgen mosaic OME-TIFFs for QuPath region annotation. Analogous to generate_qupath_x5k, but reads the mosaic directly (no zarr). + """ input: # anchor on the always-present transform CSV; the script discovers the mosaics transform=lambda wc: os.path.join( - merscope_dir_for(wc.sample), "images", "micron_to_mosaic_pixel_transform.csv"), + merscope_dir_for(wc.sample), + "images", + "micron_to_mosaic_pixel_transform.csv", + ), output: qupath_image=f"{SAMPLES_DIR}/{{sample}}/QuPath_image/{{sample}}_morphology.tiff", qupath_meta=f"{SAMPLES_DIR}/{{sample}}/QuPath_image/{{sample}}_morphology_scalefactors.json", background=f"{SAMPLES_DIR}/{{sample}}/QuPath_image/{{sample}}_background.npy", - params: - sample_id=lambda wc: wc.sample, - merscope_dir=lambda wc: merscope_dir_for(wc.sample), - z_index=MERSCOPE.get("z_index", 3), - hires_pixel_size_um=MERSCOPE.get("hires_pixel_size_um", 1.0), - channels=MERSCOPE.get("channels", DEFAULT_MER_CHANNELS), log: out=f"{LOGDIR}/generate_qupath_mer/{{sample}}.out", err=f"{LOGDIR}/generate_qupath_mer/{{sample}}.err", @@ -25,10 +23,15 @@ rule generate_qupath_mer: f"{LOGDIR}/benchmarks/generate_qupath_mer/{{sample}}.tsv" conda: "../envs/xenium5k.yaml" - threads: - get_resource("generate_qupath_mer", "threads") + threads: get_resource("generate_qupath_mer", "threads") resources: mem_mb=mem_mb_attempt("generate_qupath_mer"), runtime=get_resource("generate_qupath_mer", "runtime"), + params: + sample_id=lambda wc: wc.sample, + merscope_dir=lambda wc: merscope_dir_for(wc.sample), + z_index=MERSCOPE.get("z_index", 3), + hires_pixel_size_um=MERSCOPE.get("hires_pixel_size_um", 1.0), + channels=MERSCOPE.get("channels", DEFAULT_MER_CHANNELS), script: "../scripts/generate_qupath_mer.py" diff --git a/workflow/rules/generate_qupath_vhd.smk b/workflow/rules/generate_qupath_vhd.smk index b825271..7b08307 100644 --- a/workflow/rules/generate_qupath_vhd.smk +++ b/workflow/rules/generate_qupath_vhd.smk @@ -4,18 +4,17 @@ rule generate_qupath_vhd: sr_done=f"{OUTDIR_SR}/{{sample}}/.done", output: qupath_image=f"{SAMPLES_DIR}/{{sample}}/QuPath_image/{{sample}}_tissue_hires_image.png", - params: - sample_id=lambda wc: wc.sample, - sr_outdir=lambda wc, input: os.path.dirname(input.sr_done), log: out=f"{LOGDIR}/generate_qupath_vhd/{{sample}}.out", err=f"{LOGDIR}/generate_qupath_vhd/{{sample}}.err", conda: "../envs/visiumhd.yaml" - threads: - get_resource("generate_qupath_vhd", "threads") + threads: get_resource("generate_qupath_vhd", "threads") resources: mem_mb=mem_mb_attempt("generate_qupath_vhd"), runtime=get_resource("generate_qupath_vhd", "runtime"), + params: + sample_id=lambda wc: wc.sample, + sr_outdir=lambda wc, input: os.path.dirname(input.sr_done), script: "../scripts/generate_qupath_vhd.py" diff --git a/workflow/rules/generate_qupath_x5k.smk b/workflow/rules/generate_qupath_x5k.smk index c7a023d..5d833ff 100644 --- a/workflow/rules/generate_qupath_x5k.smk +++ b/workflow/rules/generate_qupath_x5k.smk @@ -8,10 +8,6 @@ rule generate_qupath_x5k: output: qupath_image=f"{SAMPLES_DIR}/{{sample}}/QuPath_image/{{sample}}_morphology.tiff", qupath_meta=f"{SAMPLES_DIR}/{{sample}}/QuPath_image/{{sample}}_morphology_scalefactors.json", - params: - sample_id=lambda wc: wc.sample, - zarr_path=lambda wc, input: str(input.done)[:-len(".done")], - qupath_pyramid_level=XENIUM5K.get("qupath_pyramid_level", 3), log: out=f"{LOGDIR}/generate_qupath_x5k/{{sample}}.out", err=f"{LOGDIR}/generate_qupath_x5k/{{sample}}.err", @@ -19,10 +15,13 @@ rule generate_qupath_x5k: f"{LOGDIR}/benchmarks/generate_qupath_x5k/{{sample}}.tsv" conda: "../envs/xenium5k.yaml" - threads: - get_resource("generate_qupath_x5k", "threads") + threads: get_resource("generate_qupath_x5k", "threads") resources: mem_mb=mem_mb_attempt("generate_qupath_x5k"), runtime=get_resource("generate_qupath_x5k", "runtime"), + params: + sample_id=lambda wc: wc.sample, + zarr_path=lambda wc, input: str(input.done)[: -len(".done")], + qupath_pyramid_level=XENIUM5K.get("qupath_pyramid_level", 3), script: "../scripts/generate_qupath_x5k.py" diff --git a/workflow/rules/ingest_ref.smk b/workflow/rules/ingest_ref.smk index 736be98..2f7b569 100644 --- a/workflow/rules/ingest_ref.smk +++ b/workflow/rules/ingest_ref.smk @@ -17,11 +17,6 @@ rule ingest_ref: output: adata_ingested=f"{SAMPLES_DIR}/{{sample}}/adata_{{sample}}_ingested.h5ad", plots_dir=directory(f"{SAMPLES_DIR}/{{sample}}/ingest"), - params: - sample_id=lambda wc: wc.sample, - ref_label_key=config.get("ingest_ref_label_key", "cell_type"), - de_n_genes=ANALYSIS.get("de_n_genes", 10), - annotation_colors=config.get("annotation_colors", {}), log: out=f"{LOGDIR}/ingest_ref/{{sample}}.out", err=f"{LOGDIR}/ingest_ref/{{sample}}.err", @@ -29,10 +24,14 @@ rule ingest_ref: f"{LOGDIR}/benchmarks/ingest_ref/{{sample}}.tsv" conda: "../envs/visiumhd.yaml" - threads: - get_resource("ingest_ref", "threads") + threads: get_resource("ingest_ref", "threads") resources: mem_mb=mem_mb_attempt("ingest_ref"), runtime=get_resource("ingest_ref", "runtime"), + params: + sample_id=lambda wc: wc.sample, + ref_label_key=config.get("ingest_ref_label_key", "cell_type"), + de_n_genes=ANALYSIS.get("de_n_genes", 10), + annotation_colors=config.get("annotation_colors", {}), script: "../scripts/ingest_ref.py" diff --git a/workflow/rules/integrate_samples.smk b/workflow/rules/integrate_samples.smk index 43260bd..120240e 100644 --- a/workflow/rules/integrate_samples.smk +++ b/workflow/rules/integrate_samples.smk @@ -13,6 +13,17 @@ rule integrate_samples: concatenated=f"{OUTDIR_PP}/integrated_samples/concatenated.h5ad", harmony=f"{OUTDIR_PP}/integrated_samples/harmony_integrated.h5ad", sketched=f"{OUTDIR_PP}/integrated_samples/sketched.h5ad", + log: + out=f"{LOGDIR}/integrate_samples/integrate.out", + err=f"{LOGDIR}/integrate_samples/integrate.err", + benchmark: + f"{LOGDIR}/benchmarks/integrate_samples/integrate.tsv" + conda: + "../envs/visiumhd.yaml" + threads: get_resource("integrate_samples", "threads") + resources: + mem_mb=mem_mb_attempt("integrate_samples"), + runtime=get_resource("integrate_samples", "runtime"), params: sample_ids=SAMPLE_IDS, n_neighbors=ANALYSIS.get("n_neighbors", 10), @@ -25,17 +36,5 @@ rule integrate_samples: extra_annot_columns=EXTRA_ANNOT_COLUMNS, sample_colors=SAMPLE_COLORS, integrate_key=INTEGRATE_KEY, - log: - out=f"{LOGDIR}/integrate_samples/integrate.out", - err=f"{LOGDIR}/integrate_samples/integrate.err", - benchmark: - f"{LOGDIR}/benchmarks/integrate_samples/integrate.tsv" - conda: - "../envs/visiumhd.yaml" - threads: - get_resource("integrate_samples", "threads") - resources: - mem_mb=mem_mb_attempt("integrate_samples"), - runtime=get_resource("integrate_samples", "runtime"), script: "../scripts/integrate_samples.py" diff --git a/workflow/rules/leiden_analysis.smk b/workflow/rules/leiden_analysis.smk index 5483ae4..8081868 100644 --- a/workflow/rules/leiden_analysis.smk +++ b/workflow/rules/leiden_analysis.smk @@ -11,15 +11,6 @@ rule leiden_analysis: cell_markers=config["snakemake_cell_markers"], output: res_dir=directory(f"{SAMPLES_DIR}/{{sample}}/leiden_resolution_{{resolution}}"), - params: - sample_id=lambda wc: wc.sample, - resolution=lambda wc: wc.resolution, - de_n_genes=ANALYSIS.get("de_n_genes", 10), - resolution_scan_min=ANALYSIS.get("resolution_scan_min", 0.2), - resolution_scan_max=ANALYSIS.get("resolution_scan_max", 0.8), - resolution_scan_step=ANALYSIS.get("resolution_scan_step", 0.1), - annotation_colors=config.get("annotation_colors", {}), - region_colors=ANALYSIS.get("region_colors", {}), log: out=f"{LOGDIR}/leiden_analysis/{{sample}}_res{{resolution}}.out", err=f"{LOGDIR}/leiden_analysis/{{sample}}_res{{resolution}}.err", @@ -27,10 +18,18 @@ rule leiden_analysis: f"{LOGDIR}/benchmarks/leiden_analysis/{{sample}}_res{{resolution}}.tsv" conda: "../envs/visiumhd.yaml" - threads: - get_resource("leiden_analysis", "threads") + threads: get_resource("leiden_analysis", "threads") resources: mem_mb=mem_mb_attempt("leiden_analysis"), runtime=get_resource("leiden_analysis", "runtime"), + params: + sample_id=lambda wc: wc.sample, + resolution=lambda wc: wc.resolution, + de_n_genes=ANALYSIS.get("de_n_genes", 10), + resolution_scan_min=ANALYSIS.get("resolution_scan_min", 0.2), + resolution_scan_max=ANALYSIS.get("resolution_scan_max", 0.8), + resolution_scan_step=ANALYSIS.get("resolution_scan_step", 0.1), + annotation_colors=config.get("annotation_colors", {}), + region_colors=ANALYSIS.get("region_colors", {}), script: "../scripts/leiden_analysis.py" diff --git a/workflow/rules/neighbourhood_analysis.smk b/workflow/rules/neighbourhood_analysis.smk index 4b4e172..9845add 100644 --- a/workflow/rules/neighbourhood_analysis.smk +++ b/workflow/rules/neighbourhood_analysis.smk @@ -7,11 +7,9 @@ rule neighbourhood_analysis: input: adata=rules.annotate_cells.output.adata_annot, output: - results_dir=directory(f"{SAMPLES_DIR}/{{sample}}/neighbourhood_analysis/{{annot_type}}"), - params: - sample_id=lambda wc: wc.sample, - annot_type=lambda wc: wc.annot_type, - annotation_colors=config.get("annotation_colors", {}), + results_dir=directory( + f"{SAMPLES_DIR}/{{sample}}/neighbourhood_analysis/{{annot_type}}" + ), log: out=f"{LOGDIR}/neighbourhood_analysis/{{sample}}_{{annot_type}}.out", err=f"{LOGDIR}/neighbourhood_analysis/{{sample}}_{{annot_type}}.err", @@ -19,10 +17,13 @@ rule neighbourhood_analysis: f"{LOGDIR}/benchmarks/neighbourhood_analysis/{{sample}}_{{annot_type}}.tsv" conda: "../envs/visiumhd.yaml" - threads: - get_resource("neighbourhood_analysis", "threads") + threads: get_resource("neighbourhood_analysis", "threads") resources: mem_mb=mem_mb_attempt("neighbourhood_analysis"), runtime=get_resource("neighbourhood_analysis", "runtime"), + params: + sample_id=lambda wc: wc.sample, + annot_type=lambda wc: wc.annot_type, + annotation_colors=config.get("annotation_colors", {}), script: "../scripts/neighbourhood_analysis.py" diff --git a/workflow/rules/prepare_input_ate.smk b/workflow/rules/prepare_input_ate.smk index 3cd0c91..c3b1e60 100644 --- a/workflow/rules/prepare_input_ate.smk +++ b/workflow/rules/prepare_input_ate.smk @@ -11,12 +11,6 @@ rule prepare_input_ate: unpack(_ate_prepare_inputs), output: h5ad=config["contract"]["unfiltered_h5ad"], - params: - sample_id=lambda wc: wc.sample, - zarr_path=lambda wc, input: str(input.done)[:-len(".done")], - geojson_dir=GEOJ_DIR, - hires_pyramid_level=ATERA.get("hires_pyramid_level", 3), - pixel_size_um=ATERA.get("pixel_size_um", 0.2125), log: out=f"{LOGDIR}/prepare_input_ate/{{sample}}.out", err=f"{LOGDIR}/prepare_input_ate/{{sample}}.err", @@ -24,10 +18,15 @@ rule prepare_input_ate: f"{LOGDIR}/benchmarks/prepare_input_ate/{{sample}}.tsv" conda: "../envs/atera.yaml" - threads: - get_resource("prepare_input_ate", "threads") + threads: get_resource("prepare_input_ate", "threads") resources: mem_mb=mem_mb_attempt("prepare_input_ate"), runtime=get_resource("prepare_input_ate", "runtime"), + params: + sample_id=lambda wc: wc.sample, + zarr_path=lambda wc, input: str(input.done)[: -len(".done")], + geojson_dir=GEOJ_DIR, + hires_pyramid_level=ATERA.get("hires_pyramid_level", 3), + pixel_size_um=ATERA.get("pixel_size_um", 0.2125), script: "../scripts/prepare_input_ate.py" diff --git a/workflow/rules/prepare_input_mer.smk b/workflow/rules/prepare_input_mer.smk index dabc0fd..d3829df 100644 --- a/workflow/rules/prepare_input_mer.smk +++ b/workflow/rules/prepare_input_mer.smk @@ -6,14 +6,10 @@ rule prepare_input_mer: (so the mosaic is read once and the contract shares the QuPath grid). """ input: - qupath_meta=rules.generate_qupath_mer.output.qupath_meta, # px<->µm mapping + p0 - background=rules.generate_qupath_mer.output.background, # grey contract image + qupath_meta=rules.generate_qupath_mer.output.qupath_meta, # px<->µm mapping + p0 + background=rules.generate_qupath_mer.output.background, # grey contract image output: h5ad=config["contract"]["unfiltered_h5ad"], - params: - sample_id=lambda wc: wc.sample, - merscope_dir=lambda wc: merscope_dir_for(wc.sample), - geojson_dir=GEOJ_DIR, log: out=f"{LOGDIR}/prepare_input_mer/{{sample}}.out", err=f"{LOGDIR}/prepare_input_mer/{{sample}}.err", @@ -21,10 +17,13 @@ rule prepare_input_mer: f"{LOGDIR}/benchmarks/prepare_input_mer/{{sample}}.tsv" conda: "../envs/xenium5k.yaml" - threads: - get_resource("prepare_input_mer", "threads") + threads: get_resource("prepare_input_mer", "threads") resources: mem_mb=mem_mb_attempt("prepare_input_mer"), runtime=get_resource("prepare_input_mer", "runtime"), + params: + sample_id=lambda wc: wc.sample, + merscope_dir=lambda wc: merscope_dir_for(wc.sample), + geojson_dir=GEOJ_DIR, script: "../scripts/prepare_input_mer.py" diff --git a/workflow/rules/prepare_input_vhd.smk b/workflow/rules/prepare_input_vhd.smk index dab1b6e..b3950ac 100644 --- a/workflow/rules/prepare_input_vhd.smk +++ b/workflow/rules/prepare_input_vhd.smk @@ -5,10 +5,6 @@ rule prepare_input_vhd: hires_png=rules.generate_qupath_vhd.output.qupath_image, output: h5ad=config["contract"]["unfiltered_h5ad"], - params: - sample_id=lambda wc: wc.sample, - sr_outdir=lambda wc, input: os.path.dirname(input.sr_done), - geojson_path=GEOJ_DIR, log: out=f"{LOGDIR}/prepare_input_vhd/{{sample}}.out", err=f"{LOGDIR}/prepare_input_vhd/{{sample}}.err", @@ -16,10 +12,13 @@ rule prepare_input_vhd: f"{LOGDIR}/benchmarks/prepare_input_vhd/{{sample}}.tsv" conda: "../envs/visiumhd.yaml" - threads: - get_resource("prepare_input_vhd", "threads") + threads: get_resource("prepare_input_vhd", "threads") resources: mem_mb=mem_mb_attempt("prepare_input_vhd"), runtime=get_resource("prepare_input_vhd", "runtime"), + params: + sample_id=lambda wc: wc.sample, + sr_outdir=lambda wc, input: os.path.dirname(input.sr_done), + geojson_path=GEOJ_DIR, script: "../scripts/prepare_input_vhd.py" diff --git a/workflow/rules/prepare_input_x5k.smk b/workflow/rules/prepare_input_x5k.smk index 2f744e0..f9bffc5 100644 --- a/workflow/rules/prepare_input_x5k.smk +++ b/workflow/rules/prepare_input_x5k.smk @@ -2,18 +2,13 @@ rule prepare_input_x5k: """ Xenium5k Headblock — zarr + QuPath GeoJSON → the standardized UNFILTERED CONTRACT h5ad. Analogous to prepare_input_vhd for Visium HD; writes an AnnData.h5ad contract file per sample. + """ input: - done=rules.convert_zarr_x5k.output.done, # zarr completion marker - qupath_meta=rules.generate_qupath_x5k.output.qupath_meta, # geojson px→µm scale + done=rules.convert_zarr_x5k.output.done, # zarr completion marker + qupath_meta=rules.generate_qupath_x5k.output.qupath_meta, # geojson px→µm scale output: h5ad=config["contract"]["unfiltered_h5ad"], - params: - sample_id=lambda wc: wc.sample, - zarr_path=lambda wc, input: str(input.done)[:-len(".done")], - geojson_dir=GEOJ_DIR, - hires_pyramid_level=XENIUM5K.get("hires_pyramid_level", 3), - pixel_size_um=XENIUM5K.get("pixel_size_um", 0.2125), log: out=f"{LOGDIR}/prepare_input_x5k/{{sample}}.out", err=f"{LOGDIR}/prepare_input_x5k/{{sample}}.err", @@ -21,10 +16,15 @@ rule prepare_input_x5k: f"{LOGDIR}/benchmarks/prepare_input_x5k/{{sample}}.tsv" conda: "../envs/xenium5k.yaml" - threads: - get_resource("prepare_input_x5k", "threads") + threads: get_resource("prepare_input_x5k", "threads") resources: mem_mb=mem_mb_attempt("prepare_input_x5k"), runtime=get_resource("prepare_input_x5k", "runtime"), + params: + sample_id=lambda wc: wc.sample, + zarr_path=lambda wc, input: str(input.done)[: -len(".done")], + geojson_dir=GEOJ_DIR, + hires_pyramid_level=XENIUM5K.get("hires_pyramid_level", 3), + pixel_size_um=XENIUM5K.get("pixel_size_um", 0.2125), script: "../scripts/prepare_input_x5k.py" diff --git a/workflow/rules/preprocess_umap.smk b/workflow/rules/preprocess_umap.smk index b11ac12..cf057f2 100644 --- a/workflow/rules/preprocess_umap.smk +++ b/workflow/rules/preprocess_umap.smk @@ -6,16 +6,27 @@ rule preprocess_umap: adata=f"{SAMPLES_DIR}/{{sample}}/adata_{{sample}}.h5ad", metadata=f"{SAMPLES_DIR}/{{sample}}/metadata_{{sample}}.tsv", report=f"{SAMPLES_DIR}/{{sample}}/{{sample}}_report.tsv", + log: + out=f"{LOGDIR}/preprocess_umap/{{sample}}.out", + err=f"{LOGDIR}/preprocess_umap/{{sample}}.err", + benchmark: + f"{LOGDIR}/benchmarks/preprocess_umap/{{sample}}.tsv" + conda: + "../envs/visiumhd.yaml" + threads: get_resource("preprocess_umap", "threads") + resources: + mem_mb=mem_mb_attempt("preprocess_umap"), + runtime=get_resource("preprocess_umap", "runtime"), params: sample_id=lambda wc: wc.sample, - sample_meta=lambda wc: core_sample_meta(wc.sample), # design columns → obs + report + sample_meta=lambda wc: core_sample_meta(wc.sample), # design columns → obs + report # analysis.* filtering defaults; the script may override any of these # per sample from the optional thresholds_tsv (absent → these defaults). min_counts=ANALYSIS.get("min_counts", 1), min_cells=ANALYSIS.get("min_cells", 3), min_genes=ANALYSIS.get("min_genes", 100), - max_counts=ANALYSIS.get("max_counts", None), # optional upper bound; None = off - max_pct_mt=ANALYSIS.get("max_pct_mt", None), # optional upper bound; None = off + max_counts=ANALYSIS.get("max_counts", None), # optional upper bound; None = off + max_pct_mt=ANALYSIS.get("max_pct_mt", None), # optional upper bound; None = off n_neighbors=ANALYSIS.get("n_neighbors", 10), n_pcs=ANALYSIS.get("n_pcs", 30), resolution_scan_min=ANALYSIS.get("resolution_scan_min", 0.2), @@ -27,19 +38,9 @@ rule preprocess_umap: precomputed_metadata_dir=config.get("precomputed_metadata_dir", ""), external_enabled=EXTERNAL_ENABLED, external_column=(config.get("external_annotation", {}) or {}).get("column", ""), - keep_unannotated=(config.get("external_annotation", {}) or {}).get("keep_unannotated", True), + keep_unannotated=(config.get("external_annotation", {}) or {}).get( + "keep_unannotated", True + ), region_colors=ANALYSIS.get("region_colors", {}), - log: - out=f"{LOGDIR}/preprocess_umap/{{sample}}.out", - err=f"{LOGDIR}/preprocess_umap/{{sample}}.err", - benchmark: - f"{LOGDIR}/benchmarks/preprocess_umap/{{sample}}.tsv" - conda: - "../envs/visiumhd.yaml" - threads: - get_resource("preprocess_umap", "threads") - resources: - mem_mb=mem_mb_attempt("preprocess_umap"), - runtime=get_resource("preprocess_umap", "runtime"), script: "../scripts/preprocess_umap.py" diff --git a/workflow/rules/pseudobulk_aggregate.smk b/workflow/rules/pseudobulk_aggregate.smk index 64494ae..312f2d0 100644 --- a/workflow/rules/pseudobulk_aggregate.smk +++ b/workflow/rules/pseudobulk_aggregate.smk @@ -15,9 +15,20 @@ rule pseudobulk_aggregate: agg_dir=directory( f"{OUTDIR_PP}/pseudobulk/{{annot_type}}/{{analysis_level}}/aggregated" ), + log: + out=f"{LOGDIR}/pseudobulk_aggregate/{{annot_type}}_{{analysis_level}}.out", + err=f"{LOGDIR}/pseudobulk_aggregate/{{annot_type}}_{{analysis_level}}.err", + benchmark: + f"{LOGDIR}/benchmarks/pseudobulk_aggregate/{{annot_type}}_{{analysis_level}}.tsv" wildcard_constraints: annot_type="tsv_annotation|ingest_annotation", analysis_level="by_region|by_celltype_region|by_niche_region", + conda: + "../envs/pseudobulk_aggregate.yaml" + threads: get_resource("pseudobulk_aggregate", "threads") + resources: + mem_mb=mem_mb_attempt("pseudobulk_aggregate"), + runtime=get_resource("pseudobulk_aggregate", "runtime"), params: annot_type=lambda wc: wc.annot_type, analysis_level=lambda wc: wc.analysis_level, @@ -28,17 +39,5 @@ rule pseudobulk_aggregate: dpi=ANALYSIS.get("plot_dpi", 300), extra_annot_columns=EXTRA_ANNOT_COLUMNS, sample_colors=SAMPLE_COLORS, - log: - out=f"{LOGDIR}/pseudobulk_aggregate/{{annot_type}}_{{analysis_level}}.out", - err=f"{LOGDIR}/pseudobulk_aggregate/{{annot_type}}_{{analysis_level}}.err", - benchmark: - f"{LOGDIR}/benchmarks/pseudobulk_aggregate/{{annot_type}}_{{analysis_level}}.tsv" - conda: - "../envs/pseudobulk_aggregate.yaml" - threads: - get_resource("pseudobulk_aggregate", "threads") - resources: - mem_mb=mem_mb_attempt("pseudobulk_aggregate"), - runtime=get_resource("pseudobulk_aggregate", "runtime"), script: "../scripts/pseudobulk_aggregate.py" diff --git a/workflow/rules/pseudobulk_de.smk b/workflow/rules/pseudobulk_de.smk index 347a5f4..15c5536 100644 --- a/workflow/rules/pseudobulk_de.smk +++ b/workflow/rules/pseudobulk_de.smk @@ -9,9 +9,20 @@ rule pseudobulk_de: results_dir=directory( f"{OUTDIR_PP}/pseudobulk/{{annot_type}}/{{analysis_level}}/de_results" ), + log: + out=f"{LOGDIR}/pseudobulk_de/{{annot_type}}_{{analysis_level}}.out", + err=f"{LOGDIR}/pseudobulk_de/{{annot_type}}_{{analysis_level}}.err", + benchmark: + f"{LOGDIR}/benchmarks/pseudobulk_de/{{annot_type}}_{{analysis_level}}.tsv" wildcard_constraints: annot_type="tsv_annotation|refined_annotation|ingest_annotation", analysis_level="by_region|by_celltype_region|by_niche_region", + conda: + "../envs/pseudobulk_de.yaml" + threads: get_resource("pseudobulk_de", "threads") + resources: + mem_mb=mem_mb_attempt("pseudobulk_de"), + runtime=get_resource("pseudobulk_de", "runtime"), params: annot_type=lambda wc: wc.annot_type, analysis_level=lambda wc: wc.analysis_level, @@ -29,17 +40,5 @@ rule pseudobulk_de: extra_anno_col_names=EXTRA_ANNO_COLS, extra_anno_values=EXTRA_ANNO_VALS, extra_anno_colors=EXTRA_ANNO_COLORS, - log: - out=f"{LOGDIR}/pseudobulk_de/{{annot_type}}_{{analysis_level}}.out", - err=f"{LOGDIR}/pseudobulk_de/{{annot_type}}_{{analysis_level}}.err", - benchmark: - f"{LOGDIR}/benchmarks/pseudobulk_de/{{annot_type}}_{{analysis_level}}.tsv" - conda: - "../envs/pseudobulk_de.yaml" - threads: - get_resource("pseudobulk_de", "threads") - resources: - mem_mb=mem_mb_attempt("pseudobulk_de"), - runtime=get_resource("pseudobulk_de", "runtime"), script: "../scripts/pseudobulk_de.R" diff --git a/workflow/rules/qc_sweep.smk b/workflow/rules/qc_sweep.smk index 7b7a2e2..8e77c07 100644 --- a/workflow/rules/qc_sweep.smk +++ b/workflow/rules/qc_sweep.smk @@ -9,9 +9,11 @@ rule qc_sweep: input: h5ad=lambda wc: _CONTRACT["unfiltered_h5ad"].format(sample=wc.sample), validation=f"{_OUT}/{{sample}}/validation/input_validation.json", - ingest_ref=lambda wc: (config.get("ingest_ref", "") - if _QC.get("ingest_enabled") and config.get("ingest_ref") - else []), + ingest_ref=lambda wc: ( + config.get("ingest_ref", "") + if _QC.get("ingest_enabled") and config.get("ingest_ref") + else [] + ), output: violins_png=f"{_OUT}/{{sample}}/qc_sweep/qc_violins.png", low_genes_png=f"{_OUT}/{{sample}}/qc_sweep/qc_remove_low_genes.png", @@ -22,6 +24,15 @@ rule qc_sweep: joint_png=f"{_OUT}/{{sample}}/qc_sweep/qc_joint_scatter.png", summary_tsv=f"{_OUT}/{{sample}}/qc_sweep/qc_thresholds_summary.tsv", ingest_tsv=f"{_OUT}/{{sample}}/qc_sweep/qc_ingest_removed.tsv", + log: + out=f"{_OUT}/logs/{{sample}}/qc_sweep.log", + err=f"{_OUT}/logs/{{sample}}/qc_sweep.err", + conda: + "../envs/visiumhd.yaml" + threads: get_resource("qc_sweep", "threads") + resources: + mem_mb=mem_mb_attempt("qc_sweep"), + runtime=get_resource("qc_sweep", "runtime"), params: sample_id=lambda wc: wc.sample, thresholds=lambda wc: _thresholds_for(wc.sample), @@ -31,15 +42,5 @@ rule qc_sweep: ref_label_key=config.get("ingest_ref_label_key", "cell_type"), region_colors=config.get("analysis", {}).get("region_colors", {}), region_levels=config.get("analysis", {}).get("region_levels", []), - log: - out=f"{_OUT}/logs/{{sample}}/qc_sweep.log", - err=f"{_OUT}/logs/{{sample}}/qc_sweep.err", - conda: - "../envs/visiumhd.yaml" - threads: - get_resource("qc_sweep", "threads") - resources: - mem_mb=mem_mb_attempt("qc_sweep"), - runtime=get_resource("qc_sweep", "runtime"), script: "../scripts/qc_sweep.py" diff --git a/workflow/rules/sample_report.smk b/workflow/rules/sample_report.smk index eada07e..d8868c0 100644 --- a/workflow/rules/sample_report.smk +++ b/workflow/rules/sample_report.smk @@ -10,12 +10,6 @@ rule sample_report: annotated=expand(rules.annotate_cells.output.adata_annot, sample=SAMPLE_IDS), output: report=f"{OUTDIR_PP}/integrated_samples/samples_report.pdf", - params: - sample_ids=SAMPLE_IDS, - annotation_colors=config.get("annotation_colors", {}), - region_colors=ANALYSIS.get("region_colors", {}), - dpi=ANALYSIS.get("plot_dpi", 300), - niche_column=GENE_EXPLORATION.get("niche_column", ""), log: out=f"{LOGDIR}/sample_report/sample_report.out", err=f"{LOGDIR}/sample_report/sample_report.err", @@ -23,10 +17,15 @@ rule sample_report: f"{LOGDIR}/benchmarks/sample_report/sample_report.tsv" conda: "../envs/visiumhd.yaml" - threads: - get_resource("sample_report", "threads") + threads: get_resource("sample_report", "threads") resources: mem_mb=mem_mb_attempt("sample_report"), runtime=get_resource("sample_report", "runtime"), + params: + sample_ids=SAMPLE_IDS, + annotation_colors=config.get("annotation_colors", {}), + region_colors=ANALYSIS.get("region_colors", {}), + dpi=ANALYSIS.get("plot_dpi", 300), + niche_column=GENE_EXPLORATION.get("niche_column", ""), script: "../scripts/sample_report.py" diff --git a/workflow/rules/spaceranger_count_vhd.smk b/workflow/rules/spaceranger_count_vhd.smk index 38e743d..72559ce 100644 --- a/workflow/rules/spaceranger_count_vhd.smk +++ b/workflow/rules/spaceranger_count_vhd.smk @@ -10,6 +10,7 @@ rule spaceranger_count_vhd: as the output of this rule due to Martian restrictions. These files are instead tracked in downstream rules, giving Snakemake proper provenance. """ + # Space Ranger is proprietary and user-installed (path in config["spaceranger"]); # the visiumhd env below is declared only to satisfy `snakemake --lint` and to # provide the wrapper's Python — it does NOT contain Space Ranger itself. @@ -18,9 +19,7 @@ rule spaceranger_count_vhd: cytaimage=lambda wc: ( f"{SAMPLES[wc.sample]['fastq_dir']}/images/{wc.sample}_cyta.tiff" ), - image=lambda wc: ( - f"{SAMPLES[wc.sample]['fastq_dir']}/images/{wc.sample}.tiff" - ), + image=lambda wc: (f"{SAMPLES[wc.sample]['fastq_dir']}/images/{wc.sample}.tiff"), loupe_alignment=lambda wc: ( f"{SAMPLES[wc.sample]['fastq_dir']}/images/{wc.sample}.json" ), @@ -28,46 +27,45 @@ rule spaceranger_count_vhd: probe_set=config["probe_set"], output: done_flag=f"{OUTDIR_SR}/{{sample}}/.done", - params: - # outdir derived from the rule output (not a hardcoded prefix) so it - # cannot drift and is correct on non-shared filesystems. - outdir=lambda wc, output: os.path.dirname(output.done_flag), - slide=lambda wc: SAMPLES[wc.sample]["slide"], - area=lambda wc: SAMPLES[wc.sample]["area"], - spaceranger=config["spaceranger"], - fastqs_formatted=fastq_dirs_comma_separated, - conda: - "../envs/visiumhd.yaml" log: out=f"{LOGDIR}/spaceranger_count_vhd/{{sample}}.out", err=f"{LOGDIR}/spaceranger_count_vhd/{{sample}}.err", benchmark: f"{LOGDIR}/benchmarks/spaceranger_count_vhd/{{sample}}.tsv" - threads: - get_resource("spaceranger_count_vhd", "threads") + conda: + "../envs/visiumhd.yaml" + threads: get_resource("spaceranger_count_vhd", "threads") resources: mem_mb=mem_mb_attempt("spaceranger_count_vhd"), runtime=get_resource("spaceranger_count_vhd", "runtime"), + params: + # outdir derived from the rule output (not a hardcoded prefix) so it + # cannot drift and is correct on non-shared filesystems. + outdir=lambda wc, output: os.path.dirname(output.done_flag), + slide=lambda wc: SAMPLES[wc.sample]["slide"], + area=lambda wc: SAMPLES[wc.sample]["area"], + spaceranger=config["spaceranger"], + fastqs_formatted=fastq_dirs_comma_separated, shell: """ ( - echo "=== Space Ranger version ===" - {params.spaceranger} --version + echo "=== Space Ranger version ===" + {params.spaceranger} --version - {params.spaceranger} count \ - --id={wildcards.sample} \ - --transcriptome={input.transcriptome} \ - --fastqs={params.fastqs_formatted} \ - --sample={wildcards.sample} \ - --cytaimage={input.cytaimage} \ - --image={input.image} \ - --slide={params.slide} \ - --area={params.area} \ - --loupe-alignment={input.loupe_alignment} \ - --probe-set={input.probe_set} \ - --output-dir={params.outdir} \ - --create-bam=false - ) > {log.out} 2> {log.err} + {params.spaceranger} count \ + --id={wildcards.sample} \ + --transcriptome={input.transcriptome} \ + --fastqs={params.fastqs_formatted} \ + --sample={wildcards.sample} \ + --cytaimage={input.cytaimage} \ + --image={input.image} \ + --slide={params.slide} \ + --area={params.area} \ + --loupe-alignment={input.loupe_alignment} \ + --probe-set={input.probe_set} \ + --output-dir={params.outdir} \ + --create-bam=false + ) >{log.out} 2>{log.err} touch {output.done_flag} """ diff --git a/workflow/rules/spatial_niches.smk b/workflow/rules/spatial_niches.smk index 5100d22..60258c6 100644 --- a/workflow/rules/spatial_niches.smk +++ b/workflow/rules/spatial_niches.smk @@ -16,9 +16,21 @@ rule spatial_niches: unpack(_spatial_niches_inputs), output: concatenated=f"{OUTDIR_PP}/spatial_niches/spatial_niches_concatenated.h5ad", - niche_tsvs=expand(f"{OUTDIR_PP}/spatial_niches/tsv/niche_{{sample}}.tsv", - sample=SAMPLE_IDS), + niche_tsvs=expand( + f"{OUTDIR_PP}/spatial_niches/tsv/niche_{{sample}}.tsv", sample=SAMPLE_IDS + ), plots_dir=directory(f"{OUTDIR_PP}/spatial_niches/plots"), + log: + out=f"{LOGDIR}/spatial_niches/spatial_niches.out", + err=f"{LOGDIR}/spatial_niches/spatial_niches.err", + benchmark: + f"{LOGDIR}/benchmarks/spatial_niches/spatial_niches.tsv" + conda: + "../envs/spatial_niches.yaml" + threads: get_resource("spatial_niches", "threads") + resources: + mem_mb=mem_mb_attempt("spatial_niches"), + runtime=get_resource("spatial_niches", "runtime"), params: sample_ids=SAMPLE_IDS, random_seed=RANDOM_SEED, @@ -42,17 +54,5 @@ rule spatial_niches: annotation_colors=config.get("annotation_colors", {}), region_colors=ANALYSIS.get("region_colors", {}), dpi=ANALYSIS.get("plot_dpi", 300), - log: - out=f"{LOGDIR}/spatial_niches/spatial_niches.out", - err=f"{LOGDIR}/spatial_niches/spatial_niches.err", - benchmark: - f"{LOGDIR}/benchmarks/spatial_niches/spatial_niches.tsv" - conda: - "../envs/spatial_niches.yaml" - threads: - get_resource("spatial_niches", "threads") - resources: - mem_mb=mem_mb_attempt("spatial_niches"), - runtime=get_resource("spatial_niches", "runtime"), script: "../scripts/spatial_niches.py" diff --git a/workflow/rules/subcluster.smk b/workflow/rules/subcluster.smk index 5bd968e..d0f0e01 100644 --- a/workflow/rules/subcluster.smk +++ b/workflow/rules/subcluster.smk @@ -1,5 +1,6 @@ SUBCOMPARTMENTS = list(config.get("subcompartments", {}).keys()) + rule subcluster: """ Subset and subcluster a cell compartment from the concatenated dataset. @@ -11,15 +12,32 @@ rule subcluster: adata=rules.integrate_samples.output.concatenated, output: sub_dir=directory(f"{OUTDIR_PP}/Subcompartments/{{subcompartment}}"), + log: + out=f"{LOGDIR}/subcluster/{{subcompartment}}.out", + err=f"{LOGDIR}/subcluster/{{subcompartment}}.err", + benchmark: + f"{LOGDIR}/benchmarks/subcluster/{{subcompartment}}.tsv" wildcard_constraints: subcompartment="|".join(SUBCOMPARTMENTS) if SUBCOMPARTMENTS else "NONE", + conda: + "../envs/visiumhd.yaml" + threads: get_resource("subcluster", "threads") + resources: + mem_mb=mem_mb_attempt("subcluster"), + runtime=get_resource("subcluster", "runtime"), params: subcompartment=lambda wc: wc.subcompartment, strings=lambda wc: config["subcompartments"][wc.subcompartment]["strings"], annot_col=config.get("subcompartment_annot_col") or DEFAULT_ANNOT_COL, - resolution_min=lambda wc: config["subcompartments"][wc.subcompartment].get("resolution_min", 0.2), - resolution_max=lambda wc: config["subcompartments"][wc.subcompartment].get("resolution_max", 1.0), - resolution_step=lambda wc: config["subcompartments"][wc.subcompartment].get("resolution_step", 0.2), + resolution_min=lambda wc: config["subcompartments"][wc.subcompartment].get( + "resolution_min", 0.2 + ), + resolution_max=lambda wc: config["subcompartments"][wc.subcompartment].get( + "resolution_max", 1.0 + ), + resolution_step=lambda wc: config["subcompartments"][wc.subcompartment].get( + "resolution_step", 0.2 + ), n_neighbors=ANALYSIS.get("n_neighbors", 10), n_pcs=ANALYSIS.get("n_pcs", 30), de_n_genes=ANALYSIS.get("de_n_genes", 10), @@ -31,17 +49,5 @@ rule subcluster: niche_column=GENE_EXPLORATION.get("niche_column", ""), extra_annot_columns=EXTRA_ANNOT_COLUMNS, sample_colors=SAMPLE_COLORS, - log: - out=f"{LOGDIR}/subcluster/{{subcompartment}}.out", - err=f"{LOGDIR}/subcluster/{{subcompartment}}.err", - benchmark: - f"{LOGDIR}/benchmarks/subcluster/{{subcompartment}}.tsv" - conda: - "../envs/visiumhd.yaml" - threads: - get_resource("subcluster", "threads") - resources: - mem_mb=mem_mb_attempt("subcluster"), - runtime=get_resource("subcluster", "runtime"), script: "../scripts/subcluster.py" diff --git a/workflow/rules/validate_input.smk b/workflow/rules/validate_input.smk index ef8a856..f16b614 100644 --- a/workflow/rules/validate_input.smk +++ b/workflow/rules/validate_input.smk @@ -33,6 +33,15 @@ rule validate_input: h5ad=lambda wc: _CONTRACT["unfiltered_h5ad"].format(sample=wc.sample), output: report=f"{_OUT}/{{sample}}/validation/input_validation.json", + log: + out=f"{_OUT}/logs/{{sample}}/validate_input.log", + err=f"{_OUT}/logs/{{sample}}/validate_input.err", + conda: + "../envs/visiumhd.yaml" + threads: get_resource("validate_input", "threads") + resources: + mem_mb=mem_mb_attempt("validate_input"), + runtime=get_resource("validate_input", "runtime"), params: sample_id=lambda wc: wc.sample, sample_key=_CONTRACT.get("sample_key", "sample"), @@ -43,15 +52,5 @@ rule validate_input: external_enabled=EXTERNAL_ENABLED, external_column=(config.get("external_annotation", {}) or {}).get("column", ""), external_meta_dir=config.get("precomputed_metadata_dir", ""), - log: - out=f"{_OUT}/logs/{{sample}}/validate_input.log", - err=f"{_OUT}/logs/{{sample}}/validate_input.err", - conda: - "../envs/visiumhd.yaml" - threads: - get_resource("validate_input", "threads") - resources: - mem_mb=mem_mb_attempt("validate_input"), - runtime=get_resource("validate_input", "runtime"), script: "../scripts/validate_input.py"