Skip to content

feat: per-sample reference support via manifest ref/ref_name columns - #68

Closed
mrvollger wants to merge 4 commits into
remove-leviofrom
feat-per-sample-refs
Closed

feat: per-sample reference support via manifest ref/ref_name columns#68
mrvollger wants to merge 4 commits into
remove-leviofrom
feat-per-sample-refs

Conversation

@mrvollger

@mrvollger mrvollger commented Aug 14, 2026

Copy link
Copy Markdown
Member

Summary

  • The manifest can now have four columns: sample, bam, ref, ref_name. Each sample can use its own reference genome.
  • The old two-column manifest keeps working; the config then supplies ref and ref_name for all samples. A . cell in the new columns falls back to the config value for that row. Filled manifest cells override config values (with an INFO line).
  • Chromosome names and lengths now come only from each sample's BAM header, never from the fasta .fai. A BAM whose contigs are a subset of its fasta now works instead of producing silently empty per-chromosome output.
  • A new genome_file rule writes a per-sample chrom.sizes file from the BAM header; all bedtools/bigtools sizes arguments use it. The fasta/fai remain only as true references (CRAM encoding, mosdepth).
  • All reference and chromosome lookups are functions of the {sm} wildcard (get_ref(wc), get_ref_name(wc), get_chroms(wc), get_excludes(wc)), backed by an lru_cached pysam header reader. Manifest validation runs at parse time with FIRE:-prefixed errors (duplicate samples, pairing rules, malformed rows, missing ref/fai/bam).
  • The trackhub now writes each sample's own ref_name into its genomes.txt; the hg38/GRCh38 auto-excludes are keyed per sample.

Supersedes #44 (reimplemented from scratch per review feedback there: BAM header as the only chromosome source, {sm}-keyed lookups, no dead code).

Behavior change

Because chromosome lists are needed to build the DAG, every snakemake invocation (including dry-runs) now opens all manifest BAMs. Previously a completed run tolerated archived input BAMs. This mirrors the existing parse-time ref/fai checks and is documented in config/README.md.

Testing

  • pixi run test-dry (new, runs first in CI): 21-check dry-run matrix — back-compat, 4-column, sentinel, override-INFO, levio parse coverage, and one error fixture per validator branch, plus rendered-shell assertions covering both excludes branches.
  • pixi run test + pixi run test-verify (new in CI): the single-sample regression completes 49/49 rules with golden counts identical to main (1,181-line peaks file, 46,223 fire elements) — the chromosome-source swap is behaviorally invisible.
  • pixi run test-multi (new, local): full three-sample run (136/136 rules) with per-sample genomes.txt (hg38/testref/testref2 via sentinel), the empty-excludes branch, and a generated chr20-only BAM proving per-sample chromosome sets.

Notes for review

  • snakefmt's shell formatting became active with the current pixi env, so pixi run fmt reformatted shell blocks across the workflow — mechanical whitespace noise beyond the functional changes. Two pre-existing ((cat $(...) subshells in fire-peaks.smk were disambiguated to ( (cat ... (shfmt correctly rejects the former as arithmetic).
  • The unused leviosam2 chain-lift (DSA) rules are removed entirely: nothing consumed their outputs, chain mode has been broken since rule fire stopped producing a bam, and a single global chain: cannot coexist with per-sample references.
  • CI golden counts were measured locally on macOS; if the ubuntu conda builds differ, test-verify's numbers may need hardening from this PR's first CI run.

🤖 Generated with Claude Code

https://claude.ai/code/session_01NMLEvpsc6e62iPmUQKQz6E

@mrvollger

Copy link
Copy Markdown
Member Author

Update per review: chromosome lists now preserve BAM header order instead of sorting lexicographically — header order is the data's sort order. genome_bedgraph no longer re-sorts the mosdepth output, and the fiber_locations sorted intersect is anchored with -g. A fourth test-multi sample with a reversed header (chr21 before chr20) asserts the order survives into concatenated outputs. Regression counts unchanged (49/49, 1,181/46,223); multi run 184/184 with all assertions green.

snakefmt's shell formatting (shfmt) is active in the current pixi env
and reformats every shell block. Two pre-existing '((cat $(...)'
subshells are disambiguated to '( (cat ...' because shfmt correctly
rejects the former as arithmetic. No functional changes.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
@mrvollger
mrvollger force-pushed the feat-per-sample-refs branch from d543389 to 22c51a9 Compare August 14, 2026 18:54
@mrvollger

Copy link
Copy Markdown
Member Author

Rebased into two commits for reviewability: 096999d is formatting only (snakefmt shell formatting, no functional changes), and 22c51a9 is the actual feature. Review the second commit alone for a clean diff, or append ?w=1 to the files view to hide whitespace. The tree is byte-identical to the previously CI-validated state.

@mrvollger

Copy link
Copy Markdown
Member Author

Ran a multi-agent code review over the feature commit (17 agents, 4 lenses, every finding adversarially verified — 13 confirmed, 0 refuted). All fixed in e61ccd4:

Validation robustness — numeric or NA-like sample names (2023, NA) previously passed validation then crashed with raw TypeErrors at the wildcard constraint; the manifest is now read as literal text (dtype=str, keep_default_na=False). Null/empty config ref:/ref_name: values (a half-filled template) previously produced an error blaming valid manifest rows; they now get a config-specific message, and ref/fai existence uses isfile so a directory can't pass.

levio (chain mode) — the index rule wrote to a hardcoded shared path instead of its declared per-sample output (pre-existing), and the lift rules were un-instantiable because their outputs lacked {v} while the new cram input carries it. Both fixed minimally; chain mode remains unconsumed by rule all and a fuller redesign (per-ref index dedup, per-sample chains) is deliberately out of scope.

genome_file — back to an ancient() bam input, matching rule fire's policy: a re-copied bam with a fresh mtime must not cascade reruns through every genome-file consumer; the narrow stale window (interrupted run + bam swap) is the accepted trade, documented in a comment.

Test hardening — rendered-shell excludes checks are now sample-attributed (a leak of hg38 paths into the non-hg38 sample now fails, as does exclude duplication from the list-mutation class of bug); a new no-excludes.yaml fixture covers the truly-empty excludes branch that was previously claimed but untested; error-fixture substrings are branch-specific so a wrong error can't satisfy them; the chr21-contamination check is SIGPIPE-proof; the bam generator handles unmapped reads; all test pixi tasks use PIXI_PROJECT_ROOT consistently.

Full dry-run matrix (now 26 checks) and the four-sample multi assertions pass locally; CI rerunning.

Stacked pull requests have a non-main base and previously got no checks.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Nothing consumes the leviosam2 outputs, chain mode has been broken since
the fire rule stopped producing a bam (its input referenced a bam output
that no longer exists), and the rules carried further latent bugs (index
written to a hardcoded shared path instead of the declared output).
Delete the rules and the chain/levio_exe config options.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
The manifest can now carry ref and ref_name columns so each sample can
use its own reference; the two-column manifest with config-level values
keeps working, and a '.' cell falls back to the config value. Chromosome
names and lengths come only from each sample's BAM header, in header
order (the sort order of the data), and a generated per-sample genome
file replaces the fai for all bedtools and bigtools sizes arguments
(with a parse-time staleness guard, so a leftover genome file can never
disagree with its bam header). All reference and chromosome lookups are
functions of the sm wildcard. Includes the fixes from the multi-agent
code review (manifest cells parsed as literal text, config-specific
errors for null ref values, sample-attributed test assertions).
Supersedes #44.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
@mrvollger
mrvollger force-pushed the feat-per-sample-refs branch from 7e8a41f to 32b9e5d Compare August 14, 2026 21:09
@mrvollger
mrvollger changed the base branch from main to remove-levio August 14, 2026 21:09
@mrvollger

Copy link
Copy Markdown
Member Author

Restructured into a 3-PR stack: #69 (style: shell formatting + CI trigger for stacked PRs) → #70 (refactor: remove unused leviosam2) → this PR (feat: per-sample references, now a single clean commit containing only the feature, review fixes, and staleness guard). Merge order: #69, then #70 (base auto-retargets), then this. Each squash-merge keeps its conventional title, so release-please records all three correctly. The combined tree is identical to the previously CI-validated state except the one-line CI trigger change.

@mrvollger

Copy link
Copy Markdown
Member Author

GitHub auto-closed this PR when its stack base branch (remove-levio) was deleted on #70's merge, and closed PRs with a deleted base cannot be reopened. Continued in #71 — same head branch, same single feature commit, rebased onto main.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant