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

.github/workflows/main.yml

@@ -15,4 +15,6 @@ jobs:
         with:
           pixi-version: v0.70.2
           cache: true
+      - run: pixi run test-dry
       - run: pixi run test
+      - run: pixi run test-verify

CONTRIBUTING.md

@@ -1,13 +1,22 @@
 # Contributing to FIRE
 
-PRs are welcome! Please run the test suite and the formatters before
-submitting:
+PRs are welcome! Please run the validation checks, the test suite, and the
+formatters before submitting:
 
 ```bash
+pixi run test-dry
 pixi run test
+pixi run test-verify
 pixi run fmt
 ```
 
+If your change touches manifest or reference handling, also run the
+multi-sample test:
+
+```bash
+pixi run test-multi
+```
+
 ## Conventional commits and releases
 
 Releases are automated with

config/README.md

@@ -2,15 +2,27 @@
 see `config.yaml` for an example.
 
 ## Required input options
-Reference name, this is the name that will be used as the genome name in the UCSC track hub, so be sure to use a valid UCSC genome name when possible. 
+Reference name, this is the name that will be used as the genome name in the UCSC track hub, so be sure to use a valid UCSC genome name when possible. A reference name of `hg38` or `GRCh38` also turns on the default `excludes` files (see below).
 ```
 ref_name: hg38
 ```
-Reference `fasta` file:
+Reference `fasta` file (a `.fai` index must exist beside it):
 ```
 ref: /path/to/hg38.fa
 ```
-Manifest of input sample(s), must have two white-space separated columns: sample name (`sample`) and input bam file path (`bam`). See `config.tbl` for an example. The `bam` file must be indexed and aligned to the reference genome in the `ref` option.
+Manifest of input sample(s), white-space separated with a header row. See `config.tbl` for an example. The two-column form gives a sample name (`sample`) and an input bam path (`bam`); every sample then uses the `ref` and `ref_name` from this config file:
+```
+sample	bam
+sample1	/path/to/sample1.bam
+```
+The manifest can instead carry per-sample references with two more columns, `ref` and `ref_name`. Add both columns together. Every cell must be non-empty. A cell containing `.` uses the config value for that row. Filled cells override the config values:
+```
+sample	bam	ref	ref_name
+sample1	/path/to/sample1.bam	/path/to/hg38.fa	hg38
+sample2	/path/to/sample2.bam	/path/to/chm13.fa	GCA_009914755.4
+sample3	/path/to/sample3.bam	.	.
+```
+Each `bam` file must be indexed and aligned to its reference genome. Manifest paths cannot contain spaces or quotes (the manifest is whitespace-separated). FIRE reads the chromosome names and lengths from the bam header, in header order, not from the fasta, so the fasta can contain extra contigs that the bam does not use. Because of this, FIRE opens every manifest bam when it starts, for every command including dry-runs — keep the input bams readable for the lifetime of the results. Output bed files follow the bam header order (for hg38: chr1, chr2, ...), not the lexicographic order of earlier FIRE versions; anchor downstream `bedtools intersect -sorted` calls with `-g`, and rerun old results directories from scratch rather than resuming them.
 ```
 manifest: config/config.tbl
 ```
@@ -42,7 +54,7 @@ excludes:
     - annotations/cnvs.bed
 ```
 
-Reference contigs smaller than this length are skipped by the FIRE pipeline. Default is `0`.
+Contigs in the bam header smaller than this length are skipped by the FIRE pipeline. Default is `0`.
 ```
 min_contig_length: 0
 ```

pixi.toml

@@ -8,9 +8,16 @@ version = "0.2.0"
 
 [tasks]
 fmt = "ruff format . && taplo format pixi.toml && snakefmt workflow/"
-test-data = { cmd = '''bash -c 'if [ -f "$INIT_CWD/fire-test-data/test.cram" ]; then echo "test data already present, skipping download"; else mkdir -p "$INIT_CWD/fire-test-data" && rclone sync ":s3,env_auth=false,provider=Other,endpoint=s3.kopah.orci.washington.edu:stergachis/public/FIRE/test-data" "$INIT_CWD/fire-test-data/"; fi' ''' }
-test-clean = { cmd = '''bash -c 'cd "$INIT_CWD/fire-test-data" && rm -rf results temp .snakemake' ''' }
-test = { cmd = '''bash -c 'cd "$INIT_CWD/fire-test-data" && rm -rf results temp && snakemake -s "$PIXI_PROJECT_ROOT/workflow/Snakefile" --configfile test.yaml -k' ''', depends-on = [
+test-data = { cmd = '''bash -c 'if [ -f "$PIXI_PROJECT_ROOT/fire-test-data/test.cram" ]; then echo "test data already present, skipping download"; else mkdir -p "$PIXI_PROJECT_ROOT/fire-test-data" && rclone sync ":s3,env_auth=false,provider=Other,endpoint=s3.kopah.orci.washington.edu:stergachis/public/FIRE/test-data" "$PIXI_PROJECT_ROOT/fire-test-data/"; fi' ''' }
+test-clean = { cmd = '''bash -c 'cd "$PIXI_PROJECT_ROOT/fire-test-data" && rm -rf results temp .snakemake' ''' }
+test = { cmd = '''bash -c 'cd "$PIXI_PROJECT_ROOT/fire-test-data" && rm -rf results temp && snakemake -s "$PIXI_PROJECT_ROOT/workflow/Snakefile" --configfile test.yaml -k' ''', depends-on = [
+  "test-data",
+], clean-env = true }
+test-dry = { cmd = '''bash -c 'cd "$PIXI_PROJECT_ROOT/fire-test-data" && bash "$PIXI_PROJECT_ROOT/tests/dry-run-checks.sh"' ''', depends-on = [
+  "test-data",
+], clean-env = true }
+test-verify = { cmd = '''bash -c 'cd "$PIXI_PROJECT_ROOT/fire-test-data" && bash "$PIXI_PROJECT_ROOT/tests/verify-results.sh"' ''' }
+test-multi = { cmd = '''bash -c 'cd "$PIXI_PROJECT_ROOT/fire-test-data" && rm -rf results temp && python "$PIXI_PROJECT_ROOT/tests/generate-test-bams.py" && snakemake -s "$PIXI_PROJECT_ROOT/workflow/Snakefile" --configfile "$PIXI_PROJECT_ROOT/tests/config/multi.yaml" -k && bash "$PIXI_PROJECT_ROOT/tests/verify-multi.sh"' ''', depends-on = [
   "test-data",
 ], clean-env = true }
 fire = { cmd = [
@@ -34,6 +41,7 @@ slurm = { cmd = [
 
 [dependencies]
 conda = "*"
+pysam = "*"
 snakemake = "==8.21"
 snakemake-executor-plugin-slurm = ">=0.11.2"
 snakefmt = "*"

tests/README.md

@@ -0,0 +1,45 @@
+# FIRE test fixtures and scripts
+
+Everything in this directory is test material. Nothing here is an example
+configuration for real use — see `config/` for that. All configs reference
+the downloaded `fire-test-data/` directory and run with that directory as
+the working directory.
+
+## Entry points (pixi tasks)
+
+| Task | What it does |
+|---|---|
+| `pixi run test-dry` | Runs `dry-run-checks.sh`: the full validation matrix below as fast `snakemake -n` checks. Runs first in CI. |
+| `pixi run test` | Full single-sample regression run (unchanged, uses `fire-test-data/test.yaml`). |
+| `pixi run test-verify` | Runs `verify-results.sh`: golden output counts for the regression run. Runs in CI after `test`. |
+| `pixi run test-multi` | Local only. Generates a chr20-only bam and a reversed-header bam (`generate-test-bams.py`), runs the four-sample `config/multi.yaml`, then checks it with `verify-multi.sh`. |
+
+## Fixtures in `config/`
+
+Positive cases (the dry run must pass):
+
+| Fixture | Exercises |
+|---|---|
+| `single-4col.yaml` + `.tbl` | Four-column manifest, no config-level reference. |
+| `two-sample.yaml` + `.tbl` | Two samples, one with `ref_name: hg38` (auto-excludes branch) and one without (empty-excludes branch). Sets `excludes:` to guard the list-copy bug in `get_excludes`. |
+| `sentinel.yaml` + `.tbl` | A `.` cell in `ref`/`ref_name` falls back to the config value. |
+| `override-info.yaml` | Manifest columns and config values both set; asserts the override INFO line. |
+| `multi.yaml` + `.tbl` | Four samples for `test-multi`: `generated/test-chr20.bam` (bam contigs are a subset of the fasta), `generated/test-rev.bam` (header order chr21,chr20 — checks that outputs keep bam header order), and a sentinel row. |
+| `no-excludes.yaml` | Two samples with no `excludes:` key; the rendered shell must show the truly-empty excludes branch. |
+| `extra-exclude.bed` | Small exclude file used by `two-sample.yaml` and `multi.yaml`. |
+
+Error cases (`err-*`, the dry run must fail with a specific message):
+
+| Fixture | Validator branch |
+|---|---|
+| `err-no-ref-anywhere.yaml` | No reference in the manifest or the config. |
+| `err-ref-col-only.yaml` + `.tbl` | Manifest has `ref` without `ref_name`. |
+| `err-config-ref-only.yaml` | Config has `ref` without `ref_name`. |
+| `err-nan-cell.yaml` + `.tbl` | Short manifest row (missing cells); caught by the malformed-row check. |
+| `err-extra-column.yaml` + `.tbl` | Manifest row with one extra column (would silently column-shift without index_col=False). |
+| `err-sentinel-no-config.yaml` + `.tbl` | `.` cell with no config value to fall back to. |
+| `err-dup-sample.yaml` + `.tbl` | Duplicate sample names. |
+| `err-missing-refpath.yaml` + `.tbl` | Reference fasta does not exist. |
+| `err-missing-fai.yaml` + `.tbl` + `no-fai.fa` | Fasta exists but has no `.fai`. |
+| `err-missing-bam.yaml` + `.tbl` | Input bam does not exist. |
+| `err-keepchrs.yaml` | `keep_chromosomes` filters out every chromosome. |

tests/config/err-config-ref-only.yaml

@@ -0,0 +1,2 @@
+manifest: test.tbl
+ref: test.fa.gz

tests/config/err-dup-sample.tbl

@@ -0,0 +1,3 @@
+sample	bam
+test	test.cram
+test	test.cram

tests/config/err-dup-sample.yaml

@@ -0,0 +1,3 @@
+manifest: ../tests/config/err-dup-sample.tbl
+ref: test.fa.gz
+ref_name: hg38

tests/config/err-extra-column.tbl

@@ -0,0 +1,2 @@
+sample	bam	ref	ref_name
+test	test.cram	test.fa.gz	hg38	extra

tests/config/err-extra-column.yaml

@@ -0,0 +1 @@
+manifest: ../tests/config/err-extra-column.tbl

tests/config/err-keepchrs.yaml

@@ -0,0 +1,4 @@
+ref: test.fa.gz
+ref_name: hg38
+manifest: test.tbl
+keep_chromosomes: "chrX"

tests/config/err-missing-bam.tbl

@@ -0,0 +1,2 @@
+sample	bam	ref	ref_name
+test	nope.cram	test.fa.gz	hg38

tests/config/err-missing-bam.yaml

@@ -0,0 +1 @@
+manifest: ../tests/config/err-missing-bam.tbl

tests/config/err-missing-fai.tbl

@@ -0,0 +1,2 @@
+sample	bam	ref	ref_name
+test	test.cram	../tests/config/no-fai.fa	hg38

tests/config/err-missing-fai.yaml

@@ -0,0 +1 @@
+manifest: ../tests/config/err-missing-fai.tbl

tests/config/err-missing-refpath.tbl

@@ -0,0 +1,2 @@
+sample	bam	ref	ref_name
+test	test.cram	nope.fa	hg38

tests/config/err-missing-refpath.yaml

@@ -0,0 +1 @@
+manifest: ../tests/config/err-missing-refpath.tbl

tests/config/err-nan-cell.tbl

@@ -0,0 +1,3 @@
+sample	bam	ref	ref_name
+test	test.cram	test.fa.gz	hg38
+test2	test.cram	test.fa.gz

tests/config/err-nan-cell.yaml

@@ -0,0 +1 @@
+manifest: ../tests/config/err-nan-cell.tbl

tests/config/err-no-ref-anywhere.yaml

@@ -0,0 +1 @@
+manifest: test.tbl

tests/config/err-ref-col-only.tbl

@@ -0,0 +1,2 @@
+sample	bam	ref
+test	test.cram	test.fa.gz

tests/config/err-ref-col-only.yaml

@@ -0,0 +1 @@
+manifest: ../tests/config/err-ref-col-only.tbl

tests/config/err-sentinel-no-config.tbl

@@ -0,0 +1,2 @@
+sample	bam	ref	ref_name
+test	test.cram	.	.

tests/config/err-sentinel-no-config.yaml

@@ -0,0 +1 @@
+manifest: ../tests/config/err-sentinel-no-config.tbl

tests/config/extra-exclude.bed

@@ -0,0 +1 @@
+chr20	0	10000

tests/config/multi.tbl

@@ -0,0 +1,5 @@
+sample	bam	ref	ref_name
+test	test.cram	test.fa.gz	hg38
+test2	test.cram	test.fa.gz	testref
+test20	generated/test-chr20.bam	test.fa.gz	.
+test-rev	generated/test-rev.bam	test.fa.gz	hg38

tests/config/multi.yaml

@@ -0,0 +1,5 @@
+manifest: ../tests/config/multi.tbl
+ref: test.fa.gz
+ref_name: testref2
+excludes:
+  - ../tests/config/extra-exclude.bed

tests/config/no-excludes.yaml

@@ -0,0 +1 @@
+manifest: ../tests/config/two-sample.tbl

tests/config/no-fai.fa

@@ -0,0 +1,2 @@
+>tiny
+ACGT

tests/config/override-info.yaml

@@ -0,0 +1,3 @@
+manifest: ../tests/config/single-4col.tbl
+ref: test.fa.gz
+ref_name: hg38

tests/config/sentinel.tbl

@@ -0,0 +1,2 @@
+sample	bam	ref	ref_name
+test	test.cram	.	.

tests/config/sentinel.yaml

@@ -0,0 +1,3 @@
+manifest: ../tests/config/sentinel.tbl
+ref: test.fa.gz
+ref_name: hg38

tests/config/single-4col.tbl

@@ -0,0 +1,2 @@
+sample	bam	ref	ref_name
+test	test.cram	test.fa.gz	hg38

tests/config/single-4col.yaml

@@ -0,0 +1 @@
+manifest: ../tests/config/single-4col.tbl

tests/config/two-sample.tbl

@@ -0,0 +1,3 @@
+sample	bam	ref	ref_name
+test	test.cram	test.fa.gz	hg38
+test2	test.cram	test.fa.gz	testref

tests/config/two-sample.yaml

@@ -0,0 +1,3 @@
+manifest: ../tests/config/two-sample.tbl
+excludes:
+  - ../tests/config/extra-exclude.bed

tests/dry-run-checks.sh

@@ -0,0 +1,138 @@
+#!/usr/bin/env bash
+# Dry-run validation matrix for the FIRE manifest and reference handling.
+# Run with the working directory set to fire-test-data (see the test-dry pixi task).
+set -uo pipefail
+
+SNAKEFILE="$PIXI_PROJECT_ROOT/workflow/Snakefile"
+CFG="$PIXI_PROJECT_ROOT/tests/config"
+FAILURES=0
+
+snk() {
+    snakemake -s "$SNAKEFILE" -n -q rules --configfile "$@"
+}
+
+expect_pass() {
+    local config=$1
+    if ! out=$(snk "$config" 2>&1); then
+        echo "FAIL (expected pass): $config"
+        echo "$out" | tail -5
+        FAILURES=$((FAILURES + 1))
+    else
+        echo "ok (pass): $config"
+    fi
+}
+
+expect_fail() {
+    local config=$1 substring=$2
+    if out=$(snk "$config" 2>&1); then
+        echo "FAIL (expected failure): $config"
+        FAILURES=$((FAILURES + 1))
+    elif ! grep -qF "$substring" <<<"$out"; then
+        echo "FAIL (wrong message): $config, wanted: $substring"
+        echo "$out" | tail -5
+        FAILURES=$((FAILURES + 1))
+    else
+        echo "ok (fail): $config"
+    fi
+}
+
+expect_stderr() {
+    local config=$1 substring=$2
+    if ! out=$(snk "$config" 2>&1); then
+        echo "FAIL (expected pass): $config"
+        FAILURES=$((FAILURES + 1))
+    elif ! grep -qF "$substring" <<<"$out"; then
+        echo "FAIL (missing stderr line): $config, wanted: $substring"
+        FAILURES=$((FAILURES + 1))
+    else
+        echo "ok (stderr): $config"
+    fi
+}
+
+# positive cases
+expect_pass test.yaml
+expect_pass "$CFG/single-4col.yaml"
+expect_pass "$CFG/two-sample.yaml"
+expect_pass "$CFG/sentinel.yaml"
+expect_stderr "$CFG/override-info.yaml" "manifest ref/ref_name columns override config-level values"
+
+# rendered-shell checks for exclude_from_shuffle, attributed per sample.
+# -R forces the rule, so cached results cannot hide the rendering. The
+# rendered shell contains the sample-scoped output path. Paragraph mode
+# in awk (RS='') attributes each block to one sample.
+shell_block() {
+    # config, sample -> the rendered exclude_from_shuffle shell for sample
+    snakemake -s "$SNAKEFILE" -n -p -R exclude_from_shuffle --configfile "$1" 2>&1 \
+        | awk -v RS='' -v sm="results/$2/" '/bedtools genomecov/ && index($0, sm)'
+}
+
+check_block() {
+    local label=$1 block=$2 must=$3 must_not=$4
+    if [ -z "$block" ]; then
+        echo "FAIL: no rendered exclude_from_shuffle shell for $label"
+        FAILURES=$((FAILURES + 1))
+        return
+    fi
+    if [ -n "$must" ] && ! grep -qF "$must" <<<"$block"; then
+        echo "FAIL: $label rendered shell is missing: $must"
+        FAILURES=$((FAILURES + 1))
+        return
+    fi
+    if [ -n "$must_not" ] && grep -qF "$must_not" <<<"$block"; then
+        echo "FAIL: $label rendered shell must not contain: $must_not"
+        FAILURES=$((FAILURES + 1))
+        return
+    fi
+    echo "ok (rendered shell): $label"
+}
+
+# with config excludes set, sample test gets the hg38 blacklists and the
+# config exclude. Sample test2 gets only the config exclude. A leak of
+# hg38 paths into test2 fails the check.
+block_test=$(shell_block "$CFG/two-sample.yaml" test)
+block_test2=$(shell_block "$CFG/two-sample.yaml" test2)
+check_block "test (hg38 excludes)" "$block_test" "hg38.blacklist" ""
+check_block "test (config exclude)" "$block_test" "extra-exclude.bed" ""
+check_block "test2 (no hg38 leak)" "$block_test2" "extra-exclude.bed" "hg38.blacklist"
+# the shell template uses {params.exclude} two times: in the [ -n ]
+# guard and in gunzip. Each exclude path must render exactly two times.
+# More than two means a mutation bug duplicated the list.
+gap_count=$(grep -oF "hg38.gap.bed.gz" <<<"$block_test" | wc -l | tr -d ' ')
+if [ "$gap_count" -ne 2 ]; then
+    echo "FAIL: hg38.gap.bed.gz appears $gap_count times for test, expected 2"
+    FAILURES=$((FAILURES + 1))
+else
+    echo "ok (rendered shell): no exclude duplication for test"
+fi
+
+# without config excludes, sample test2 uses the empty excludes branch.
+# The [ -n ] guard must render with an empty parameter.
+block_empty=$(shell_block "$CFG/no-excludes.yaml" test2)
+check_block "test2 (empty excludes)" "$block_empty" "" "gunzip -cf ."
+if grep -qF 'if [ -n "" ]' <<<"$block_empty"; then
+    echo "ok (rendered shell): empty excludes guard renders for test2"
+else
+    echo "FAIL: empty excludes guard did not render empty for test2"
+    FAILURES=$((FAILURES + 1))
+fi
+
+
+# error cases, one per validator branch. Each substring is specific,
+# so a wrong error cannot satisfy the assertion.
+expect_fail "$CFG/err-no-ref-anywhere.yaml" "no reference specified"
+expect_fail "$CFG/err-ref-col-only.yaml" "manifest columns 'ref' and 'ref_name' must be provided together"
+expect_fail "$CFG/err-config-ref-only.yaml" "config options 'ref' and 'ref_name' must be provided together"
+expect_fail "$CFG/err-nan-cell.yaml" "missing or malformed manifest fields"
+expect_fail "$CFG/err-extra-column.yaml" "cannot parse manifest"
+expect_fail "$CFG/err-sentinel-no-config.yaml" "is not set in config.yaml"
+expect_fail "$CFG/err-dup-sample.yaml" "duplicate sample names"
+expect_fail "$CFG/err-missing-refpath.yaml" "reference file"
+expect_fail "$CFG/err-missing-fai.yaml" "reference index file"
+expect_fail "$CFG/err-missing-bam.yaml" "cannot read input bam"
+expect_fail "$CFG/err-keepchrs.yaml" "no chromosomes left"
+
+if [ "$FAILURES" -gt 0 ]; then
+    echo "dry-run-checks: $FAILURES failure(s)"
+    exit 1
+fi
+echo "dry-run-checks: all checks passed"