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

.github/workflows/main.yml

@@ -4,7 +4,6 @@ on:
   push:
     branches: [main]
   pull_request:
-    branches: [main]
 
 jobs:
   # pixi test
@@ -16,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. 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.
 ```
 manifest: config/config.tbl
 ```

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,44 @@
+# 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). |
+| `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-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,153 @@
+#!/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"
+
+# sample-attributed rendered-shell checks for exclude_from_shuffle;
+# -R forces the rule so cached results do not hide the rendering. The
+# rendered shell contains the sample-scoped output path, so paragraphs
+# (awk RS='') can be attributed 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: test gets hg38 blacklists + the config exclude,
+# test2 gets only the config exclude (a leak of hg38 paths into test2 fails)
+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 references {params.exclude} twice (the [ -n ] guard
+# and gunzip), so each exclude path renders exactly twice; more means a
+# get_excludes 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: test2 exercises the truly-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
+
+# a leftover genome file that no longer matches the bam header must be
+# removed at parse time (the staleness guard for the ancient() bam input)
+mkdir -p temp/test
+printf 'chrBOGUS\t12345\n' >temp/test/test.genome
+out=$(snk test.yaml 2>&1) || {
+    echo "FAIL: stale-genome dry-run did not build"
+    FAILURES=$((FAILURES + 1))
+}
+if ! grep -qF "removing stale genome file" <<<"$out"; then
+    echo "FAIL: stale genome file was not reported"
+    FAILURES=$((FAILURES + 1))
+elif [ -f temp/test/test.genome ]; then
+    echo "FAIL: stale genome file was not removed"
+    FAILURES=$((FAILURES + 1))
+else
+    echo "ok (stale guard): mismatched genome file removed at parse time"
+fi
+
+# error cases, one per validator branch; substrings are specific enough
+# that 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-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"