Skip to content

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

Open
mrvollger wants to merge 4 commits into
mainfrom
feat-per-sample-refs
Open

feat: per-sample reference support via manifest ref/ref_name columns#71
mrvollger wants to merge 4 commits into
mainfrom
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.

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.

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>
Mitchell R. Vollger and others added 2 commits August 14, 2026 16:15
…enome file

Every multi-chromosome bedtools sort now uses -g with the per-sample
genome file (exclude_from_shuffle, clustering_vs_null, wide and
one-percent peaks, hap-differences track), and the two sorted intersects
in fires_in_peaks that ran without -g are anchored too — an unanchored
-sorted sweep only worked while both inputs happened to share an order.
exclude_from_shuffle also drops exclude-file contigs that are absent
from the sample's bam header, since bedtools sort -g rejects undeclared
chromosomes and those records were dead weight. Per-chromosome rules
keep plain position sorts; chromosome order cannot apply within one
chromosome.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
verify-multi now asserts full content equality between the test and
test-rev samples (same reads, reversed header) for peaks and pileup —
the one check that catches a silent ordering regression. The manifest
reader uses index_col=False with ParserWarning promoted to an error, so
a uniformly extra column fails loudly instead of column-shifting or
silently dropping a field; the malformed-row check covers NaN and empty
cells. config/README.md documents that outputs follow bam header order
(not the lexicographic order of earlier FIRE versions) and that old
results directories must be rerun from scratch.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
@mrvollger

Copy link
Copy Markdown
Member Author

Second multi-agent review (ordering-focused, with empirical bedtools/bigtools/tabix experiments using chromosome names where header and lexicographic order disagree): no ordering-correctness bugs found. Decisive check: test vs test-rev (same reads, reversed header) outputs are byte-identical modulo row order — 1,180 peak rows and 3.75M pileup rows — and bedtools provably fails loudly, not silently, on order inconsistencies.

Three confirmed findings, fixed in the latest commit: (1) verify-multi.sh now asserts full content equality between test and test-rev, the one assertion that catches a silent ordering regression; (2) the output-order change (bam header order instead of lexicographic) is now documented in config/README.md with migration guidance — anchor downstream bedtools intersect -sorted with -g, and rerun pre-existing results dirs from scratch; (3) a manifest row with a uniformly extra column previously column-shifted silently via pandas index inference — now a loud parse error, with a fixture.

Migration notes (for the release)

  • Output bed files follow the bam header order (hg38: chr1, chr2, …), not the lexicographic order (chr1, chr10, …) of earlier FIRE versions.
  • Downstream bedtools intersect -sorted between FIRE outputs and lexicographically sorted annotations needs -g (or re-sorted inputs).
  • Rerun old results directories from scratch; resuming a pre-0.3 results dir mixes sort orders and fails loudly mid-run.
  • Input bams must stay readable for the lifetime of the results (headers are read at every pipeline start).

The content-equality guard checks file existence up front, so missing
files fail loudly instead of comparing two empty streams. The manifest
parse error gains a plain-language hint, the README notes that manifest
paths cannot contain spaces or quotes (the whitespace-separated format
cannot express them), and the err-nan-cell fixture description matches
the validator that actually fires.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
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