From 6db2985e596b1b94050bf1b01fa37155387fd776 Mon Sep 17 00:00:00 2001 From: "Mitchell R. Vollger" Date: Fri, 14 Aug 2026 14:08:52 -0700 Subject: [PATCH 1/4] feat: per-sample reference support via manifest ref/ref_name columns 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 --- .github/workflows/main.yml | 2 + CONTRIBUTING.md | 13 +- config/README.md | 18 +- pixi.lock | 259 ++++++++++++++--------- pixi.toml | 14 +- tests/README.md | 44 ++++ tests/config/err-config-ref-only.yaml | 2 + tests/config/err-dup-sample.tbl | 3 + tests/config/err-dup-sample.yaml | 3 + tests/config/err-keepchrs.yaml | 4 + tests/config/err-missing-bam.tbl | 2 + tests/config/err-missing-bam.yaml | 1 + tests/config/err-missing-fai.tbl | 2 + tests/config/err-missing-fai.yaml | 1 + tests/config/err-missing-refpath.tbl | 2 + tests/config/err-missing-refpath.yaml | 1 + tests/config/err-nan-cell.tbl | 3 + tests/config/err-nan-cell.yaml | 1 + tests/config/err-no-ref-anywhere.yaml | 1 + tests/config/err-ref-col-only.tbl | 2 + tests/config/err-ref-col-only.yaml | 1 + tests/config/err-sentinel-no-config.tbl | 2 + tests/config/err-sentinel-no-config.yaml | 1 + tests/config/extra-exclude.bed | 1 + tests/config/multi.tbl | 5 + tests/config/multi.yaml | 5 + tests/config/no-excludes.yaml | 1 + tests/config/no-fai.fa | 2 + tests/config/override-info.yaml | 3 + tests/config/sentinel.tbl | 2 + tests/config/sentinel.yaml | 3 + tests/config/single-4col.tbl | 2 + tests/config/single-4col.yaml | 1 + tests/config/two-sample.tbl | 3 + tests/config/two-sample.yaml | 3 + tests/dry-run-checks.sh | 153 +++++++++++++ tests/generate-test-bams.py | 69 ++++++ tests/verify-multi.sh | 52 +++++ tests/verify-results.sh | 28 +++ workflow/Snakefile | 17 +- workflow/envs/runner.yaml | 5 +- workflow/rules/apply-model.smk | 4 +- workflow/rules/common.smk | 253 +++++++++++++++++----- workflow/rules/coverages.smk | 40 +++- workflow/rules/decorated-reads.smk | 14 +- workflow/rules/fire-peaks.smk | 14 +- workflow/rules/stats.smk | 8 +- workflow/rules/track-hub.smk | 26 +-- 48 files changed, 883 insertions(+), 213 deletions(-) create mode 100644 tests/README.md create mode 100644 tests/config/err-config-ref-only.yaml create mode 100644 tests/config/err-dup-sample.tbl create mode 100644 tests/config/err-dup-sample.yaml create mode 100644 tests/config/err-keepchrs.yaml create mode 100644 tests/config/err-missing-bam.tbl create mode 100644 tests/config/err-missing-bam.yaml create mode 100644 tests/config/err-missing-fai.tbl create mode 100644 tests/config/err-missing-fai.yaml create mode 100644 tests/config/err-missing-refpath.tbl create mode 100644 tests/config/err-missing-refpath.yaml create mode 100644 tests/config/err-nan-cell.tbl create mode 100644 tests/config/err-nan-cell.yaml create mode 100644 tests/config/err-no-ref-anywhere.yaml create mode 100644 tests/config/err-ref-col-only.tbl create mode 100644 tests/config/err-ref-col-only.yaml create mode 100644 tests/config/err-sentinel-no-config.tbl create mode 100644 tests/config/err-sentinel-no-config.yaml create mode 100644 tests/config/extra-exclude.bed create mode 100644 tests/config/multi.tbl create mode 100644 tests/config/multi.yaml create mode 100644 tests/config/no-excludes.yaml create mode 100644 tests/config/no-fai.fa create mode 100644 tests/config/override-info.yaml create mode 100644 tests/config/sentinel.tbl create mode 100644 tests/config/sentinel.yaml create mode 100644 tests/config/single-4col.tbl create mode 100644 tests/config/single-4col.yaml create mode 100644 tests/config/two-sample.tbl create mode 100644 tests/config/two-sample.yaml create mode 100644 tests/dry-run-checks.sh create mode 100644 tests/generate-test-bams.py create mode 100644 tests/verify-multi.sh create mode 100644 tests/verify-results.sh diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index f5ba8c199..abf9d5d36 100644 --- a/.github/workflows/main.yml +++ b/.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 diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index 964441543..5008f9972 100644 --- a/CONTRIBUTING.md +++ b/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 diff --git a/config/README.md b/config/README.md index 01fe3a763..3a7de2e3d 100644 --- a/config/README.md +++ b/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 ``` diff --git a/pixi.lock b/pixi.lock index 060ac2f1c..d74bcf4a6 100644 --- a/pixi.lock +++ b/pixi.lock @@ -18,7 +18,7 @@ environments: - conda: https://conda.anaconda.org/conda-forge/noarch/boltons-25.0.0-pyhd8ed1ab_0.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/brotli-python-1.2.0-py312h67db365_0.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/bzip2-1.0.8-hda65f42_8.conda - - conda: https://conda.anaconda.org/conda-forge/linux-64/c-ares-1.34.5-hb9d3cd8_0.conda + - conda: https://conda.anaconda.org/conda-forge/linux-64/c-ares-1.34.8-h280c20c_1.conda - conda: https://conda.anaconda.org/conda-forge/noarch/ca-certificates-2025.11.12-hbd8a1cb_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/certifi-2025.11.12-pyhd8ed1ab_0.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/cffi-2.0.0-py312h460c074_1.conda @@ -53,7 +53,7 @@ environments: - conda: https://conda.anaconda.org/conda-forge/noarch/hpack-4.1.0-pyhd8ed1ab_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/humanfriendly-10.0-pyh707e725_8.conda - conda: https://conda.anaconda.org/conda-forge/noarch/hyperframe-6.1.0-pyhd8ed1ab_0.conda - - conda: https://conda.anaconda.org/conda-forge/linux-64/icu-75.1-he02047a_0.conda + - conda: https://conda.anaconda.org/conda-forge/linux-64/icu-78.3-py310h44b86e0_2.conda - conda: https://conda.anaconda.org/conda-forge/noarch/idna-3.11-pyhd8ed1ab_0.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/immutables-0.21-py312h4c3975b_2.conda - conda: https://conda.anaconda.org/conda-forge/noarch/iniconfig-2.3.0-pyhd8ed1ab_0.conda @@ -64,12 +64,13 @@ environments: - conda: https://conda.anaconda.org/conda-forge/noarch/jsonschema-specifications-2025.9.1-pyhcf101f3_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/jupyter_core-5.9.1-pyhc90fa1f_0.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/keyutils-1.6.3-hb9d3cd8_0.conda - - conda: https://conda.anaconda.org/conda-forge/linux-64/krb5-1.21.3-h659f571_0.conda + - conda: https://conda.anaconda.org/conda-forge/linux-64/krb5-1.22.2-hbde042b_1.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/ld_impl_linux-64-2.44-h1aa0949_5.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/libarchive-3.8.2-gpl_h7be2006_100.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/libblas-3.9.0-38_h4a7cf45_openblas.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/libcblas-3.9.0-38_h0358290_openblas.conda - - conda: https://conda.anaconda.org/conda-forge/linux-64/libcurl-8.17.0-h4e3cde8_0.conda + - conda: https://conda.anaconda.org/conda-forge/linux-64/libcurl-8.21.0-heca4667_4.conda + - conda: https://conda.anaconda.org/conda-forge/linux-64/libdeflate-1.25-hd45a770_1.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/libedit-3.1.20250104-pl5321h7949ede_0.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/libev-4.33-hd590300_2.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/libexpat-2.7.1-hecca717_0.conda @@ -82,22 +83,22 @@ environments: - conda: https://conda.anaconda.org/conda-forge/linux-64/libiconv-1.18-h3b78370_2.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/liblapack-3.9.0-38_h47877c9_openblas.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/liblapacke-3.9.0-38_h6ae95b6_openblas.conda - - conda: https://conda.anaconda.org/conda-forge/linux-64/liblzma-5.8.1-hb9d3cd8_2.conda + - conda: https://conda.anaconda.org/conda-forge/linux-64/liblzma-5.8.3-hb03c661_1.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/libmamba-2.3.3-h09abcb8_1.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/libmambapy-2.3.3-py312h9b917ee_1.conda - - conda: https://conda.anaconda.org/conda-forge/linux-64/libnghttp2-1.67.0-had1ee68_0.conda + - conda: https://conda.anaconda.org/conda-forge/linux-64/libnghttp2-1.68.1-h877daf1_0.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/libnsl-2.0.1-hb9d3cd8_1.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/libopenblas-0.3.30-pthreads_h94d23a6_3.conda + - conda: https://conda.anaconda.org/conda-forge/linux-64/libpsl-0.23.1-hf670292_0.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/libsolv-0.7.35-h9463b59_0.conda - - conda: https://conda.anaconda.org/conda-forge/linux-64/libsqlite-3.51.0-hee844dc_0.conda + - conda: https://conda.anaconda.org/conda-forge/linux-64/libsqlite-3.53.4-hf4e2dac_0.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/libssh2-1.11.1-hcf80075_0.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/libstdcxx-15.2.0-h8f9b012_7.conda - - conda: https://conda.anaconda.org/conda-forge/linux-64/libstdcxx-ng-15.2.0-h4852527_7.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/libuuid-2.41.2-he9a06e4_0.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/libxcrypt-4.4.36-hd590300_1.conda - - conda: https://conda.anaconda.org/conda-forge/linux-64/libxml2-16-2.15.1-ha9997c6_0.conda - - conda: https://conda.anaconda.org/conda-forge/linux-64/libxml2-2.15.1-h26afc86_0.conda - - conda: https://conda.anaconda.org/conda-forge/linux-64/libzlib-1.3.1-hb9d3cd8_2.conda + - conda: https://conda.anaconda.org/conda-forge/linux-64/libxml2-16-2.15.3-hca6bf5a_0.conda + - conda: https://conda.anaconda.org/conda-forge/linux-64/libxml2-2.15.3-h49c6c72_0.conda + - conda: https://conda.anaconda.org/conda-forge/linux-64/libzlib-1.3.2-h25fd6f3_3.conda - conda: https://conda.anaconda.org/conda-forge/noarch/logmuse-0.2.8-pyhd8ed1ab_1.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/lz4-c-1.10.0-h5888daf_1.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/lzo-2.10-h280c20c_1002.conda @@ -128,6 +129,7 @@ environments: - conda: https://conda.anaconda.org/conda-forge/linux-64/pydantic-core-2.41.5-py312h868fb18_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/pygments-2.19.2-pyhd8ed1ab_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/pyparsing-3.2.5-pyhcf101f3_0.conda + - conda: https://conda.anaconda.org/bioconda/linux-64/pysam-0.24.0-py312hf5ad864_1.conda - conda: https://conda.anaconda.org/conda-forge/noarch/pysocks-1.7.1-pyha55dd90_7.conda - conda: https://conda.anaconda.org/conda-forge/noarch/pytest-9.0.0-pyhcf101f3_0.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/python-3.12.12-hd63d673_1_cpython.conda @@ -271,6 +273,7 @@ environments: - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libcblas-3.11.0-9_hb0561ab_openblas.conda - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libcurl-8.21.0-hf618e03_4.conda - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libcxx-22.1.8-h55c6f16_0.conda + - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libdeflate-1.25-he7e0567_1.conda - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libedit-3.1.20250104-pl5321h26f1114_1.conda - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libev-4.33-h1a92334_3.conda - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libexpat-2.8.1-hf6b4638_1.conda @@ -329,6 +332,7 @@ environments: - conda: https://conda.anaconda.org/conda-forge/noarch/pygments-2.20.0-pyhd8ed1ab_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/pyparsing-3.3.2-pyhcf101f3_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/pyproject_hooks-1.2.0-pyhd8ed1ab_1.conda + - conda: https://conda.anaconda.org/bioconda/osx-arm64/pysam-0.24.0-py312h12d0683_1.conda - conda: https://conda.anaconda.org/conda-forge/noarch/pysocks-1.7.1-pyha55dd90_7.conda - conda: https://conda.anaconda.org/conda-forge/osx-arm64/python-3.12.13-hd1323d7_1_cpython.conda - conda: https://conda.anaconda.org/conda-forge/noarch/python-build-1.5.0-pyhc364b38_0.conda @@ -667,16 +671,18 @@ packages: license_family: BSD size: 124965 timestamp: 1785906749812 -- conda: https://conda.anaconda.org/conda-forge/linux-64/c-ares-1.34.5-hb9d3cd8_0.conda - sha256: f8003bef369f57396593ccd03d08a8e21966157269426f71e943f96e4b579aeb - md5: f7f0d6cc2dc986d42ac2689ec88192be +- conda: https://conda.anaconda.org/conda-forge/linux-64/c-ares-1.34.8-h280c20c_1.conda + sha256: 5139b6afbfaca91c47104ba9a6a40f81211e1c9200e96b73ce3a56f0ca1902f4 + md5: 2cef891b791040aab83e218c7d137679 depends: - __glibc >=2.17,<3.0.a0 - - libgcc >=13 + - libgcc >=14 + constrains: + - c-ares-static <0a0 license: MIT license_family: MIT - size: 206884 - timestamp: 1744127994291 + size: 226755 + timestamp: 1786116641939 - conda: https://conda.anaconda.org/conda-forge/osx-arm64/c-ares-1.34.8-h1a92334_1.conda sha256: 104b41473845649101ba8ebf8221c7431256465d34ce380b10e9a90558ed33ac md5: 7d9390a4d4b43f91652823d870f68065 @@ -1621,17 +1627,17 @@ packages: license_family: MIT size: 17397 timestamp: 1737618427549 -- conda: https://conda.anaconda.org/conda-forge/linux-64/icu-75.1-he02047a_0.conda - sha256: 71e750d509f5fa3421087ba88ef9a7b9be11c53174af3aa4d06aff4c18b38e8e - md5: 8b189310083baabfb622af68fd9d3ae3 +- conda: https://conda.anaconda.org/conda-forge/linux-64/icu-78.3-py310h44b86e0_2.conda + sha256: 9f07834f0c546ab14d885ce0366285f61f44e326c0edd1fc63b8294e113ae432 + md5: 72a381cbad04f24b1c2a43ef707f45b4 depends: - __glibc >=2.17,<3.0.a0 - - libgcc-ng >=12 - - libstdcxx-ng >=12 + - libstdcxx >=14 + - libgcc >=14 license: MIT license_family: MIT - size: 12129203 - timestamp: 1720853576813 + size: 14459115 + timestamp: 1786545741408 - conda: https://conda.anaconda.org/conda-forge/osx-arm64/icu-78.3-py310h579977c_2.conda sha256: 6cdb5dee54c72e56ab189fb3ad33cb28533553d42590e7e831160248f4416a43 md5: a5efc0b42bb8b42e97d0a29ae3e3c187 @@ -1819,20 +1825,21 @@ packages: license: LGPL-2.1-or-later size: 134088 timestamp: 1754905959823 -- conda: https://conda.anaconda.org/conda-forge/linux-64/krb5-1.21.3-h659f571_0.conda - sha256: 99df692f7a8a5c27cd14b5fb1374ee55e756631b9c3d659ed3ee60830249b238 - md5: 3f43953b7d3fb3aaa1d0d0723d91e368 +- conda: https://conda.anaconda.org/conda-forge/linux-64/krb5-1.22.2-hbde042b_1.conda + sha256: 9b07046870772f28740e3f6149f09ff222843733087a33c5540b169c6289652d + md5: 54157a1c8c0bb70f62dd0b17fba7e7f2 depends: - - keyutils >=1.6.1,<2.0a0 - - libedit >=3.1.20191231,<3.2.0a0 - - libedit >=3.1.20191231,<4.0a0 - - libgcc-ng >=12 - - libstdcxx-ng >=12 - - openssl >=3.3.1,<4.0a0 + - __glibc >=2.17,<3.0.a0 + - keyutils >=1.6.3,<2.0a0 + - libedit >=3.1.20250104,<3.2.0a0 + - libedit >=3.1.20250104,<4.0a0 + - libgcc >=14 + - libstdcxx >=14 + - openssl >=3.5.7,<4.0a0 license: MIT license_family: MIT - size: 1370023 - timestamp: 1719463201255 + size: 1388990 + timestamp: 1781859420533 - conda: https://conda.anaconda.org/conda-forge/osx-arm64/krb5-1.22.2-hfd3d5f3_1.conda sha256: c740e4a2e7247776a9883158fdab50ae0732c8f67f96d8f1db8ad9da5e0b5222 md5: 8780f41b013d19219faef9c82260744b @@ -1957,22 +1964,23 @@ packages: license_family: BSD size: 18110 timestamp: 1786058893756 -- conda: https://conda.anaconda.org/conda-forge/linux-64/libcurl-8.17.0-h4e3cde8_0.conda - sha256: 100e29ca864c32af15a5cc354f502d07b2600218740fdf2439fa7d66b50b3529 - md5: 01e149d4a53185622dc2e788281961f2 +- conda: https://conda.anaconda.org/conda-forge/linux-64/libcurl-8.21.0-heca4667_4.conda + sha256: aff8ef75636d0825ce34911e0bef2f26816e7afd090a9dcb4b9bacab75cbf584 + md5: 3f2fd5617cfacac49c85f6dc63842ea1 depends: - __glibc >=2.17,<3.0.a0 - - krb5 >=1.21.3,<1.22.0a0 + - krb5 >=1.22.2,<1.23.0a0 - libgcc >=14 - - libnghttp2 >=1.67.0,<2.0a0 + - libnghttp2 >=1.68.1,<2.0a0 + - libpsl >=0.23.0,<0.24.0a0 - libssh2 >=1.11.1,<2.0a0 - - libzlib >=1.3.1,<2.0a0 - - openssl >=3.5.4,<4.0a0 + - libzlib >=1.3.2,<2.0a0 + - openssl >=3.5.7,<4.0a0 - zstd >=1.5.7,<1.6.0a0 license: curl license_family: MIT - size: 460366 - timestamp: 1762333743748 + size: 480565 + timestamp: 1785500108494 - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libcurl-8.21.0-hf618e03_4.conda sha256: d25be36712d7f854a5f93513b9fbf1b0ee3219975b7fef4a3ee071b021b10167 md5: 66cf9c5003ee81ecdb0dc9f9df17bbe5 @@ -1998,6 +2006,25 @@ packages: license_family: Apache size: 569349 timestamp: 1781670209146 +- conda: https://conda.anaconda.org/conda-forge/linux-64/libdeflate-1.25-hd45a770_1.conda + sha256: 82e134c8a08b1eed9a2ed8ab578b89aa1730dcde3dea8dd87645ed0637878e54 + md5: 40f9b31aa9cf007789867df0decd0492 + depends: + - __glibc >=2.17,<3.0.a0 + - libgcc >=14 + license: MIT + license_family: MIT + size: 73710 + timestamp: 1785908694612 +- conda: https://conda.anaconda.org/conda-forge/osx-arm64/libdeflate-1.25-he7e0567_1.conda + sha256: d896f4aa4ce4c590c2838678cb1917356fdb461d2a189991c0280c818c362172 + md5: 78650d671cb56909bb3e5c13bce310f9 + depends: + - __osx >=11.0 + license: MIT + license_family: MIT + size: 55727 + timestamp: 1785909153744 - conda: https://conda.anaconda.org/conda-forge/linux-64/libedit-3.1.20250104-pl5321h7949ede_0.conda sha256: d789471216e7aba3c184cd054ed61ce3f6dac6f87a50ec69291b9297f8c18724 md5: c277e0a4d549b03ac1e9d6cbbe3d017b @@ -2242,17 +2269,17 @@ packages: license_family: BSD size: 18174 timestamp: 1786058909286 -- conda: https://conda.anaconda.org/conda-forge/linux-64/liblzma-5.8.1-hb9d3cd8_2.conda - sha256: f2591c0069447bbe28d4d696b7fcb0c5bd0b4ac582769b89addbcf26fb3430d8 - md5: 1a580f7796c7bf6393fddb8bbbde58dc +- conda: https://conda.anaconda.org/conda-forge/linux-64/liblzma-5.8.3-hb03c661_1.conda + sha256: 9787df8c22a59c9a70d3e5a10db9ad663485e75e9ccc3f09bd092cb7b95e0dab + md5: 1390b7c5ac0b1d8e447bc5efa6d3c8c2 depends: - __glibc >=2.17,<3.0.a0 - - libgcc >=13 + - libgcc >=14 constrains: - - xz 5.8.1.* + - xz 5.8.3.* license: 0BSD - size: 112894 - timestamp: 1749230047870 + size: 112995 + timestamp: 1786348617826 - conda: https://conda.anaconda.org/conda-forge/osx-arm64/liblzma-5.8.3-h8088a28_1.conda sha256: 23d0630046a3e8b164d8f80f2b74ed2605af2e7050ab9913018056402fae4311 md5: 8ab10323068b107661a4b9a4af84f3b5 @@ -2376,22 +2403,22 @@ packages: license: BSL-1.0 size: 39302 timestamp: 1786189396450 -- conda: https://conda.anaconda.org/conda-forge/linux-64/libnghttp2-1.67.0-had1ee68_0.conda - sha256: a4a7dab8db4dc81c736e9a9b42bdfd97b087816e029e221380511960ac46c690 - md5: b499ce4b026493a13774bcf0f4c33849 +- conda: https://conda.anaconda.org/conda-forge/linux-64/libnghttp2-1.68.1-h877daf1_0.conda + sha256: 663444d77a42f2265f54fb8b48c5450bfff4388d9c0f8253dd7855f0d993153f + md5: 2a45e7f8af083626f009645a6481f12d depends: - __glibc >=2.17,<3.0.a0 - - c-ares >=1.34.5,<2.0a0 + - c-ares >=1.34.6,<2.0a0 - libev >=4.33,<4.34.0a0 - libev >=4.33,<5.0a0 - libgcc >=14 - libstdcxx >=14 - libzlib >=1.3.1,<2.0a0 - - openssl >=3.5.2,<4.0a0 + - openssl >=3.5.5,<4.0a0 license: MIT license_family: MIT - size: 666600 - timestamp: 1756834976695 + size: 663344 + timestamp: 1773854035739 - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libnghttp2-1.68.1-h8f3e76b_0.conda sha256: 2bc7bc3978066f2c274ebcbf711850cc9ab92e023e433b9631958a098d11e10a md5: 6ea18834adbc3b33df9bd9fb45eaf95b @@ -2445,6 +2472,18 @@ packages: license_family: BSD size: 4318474 timestamp: 1784288246205 +- conda: https://conda.anaconda.org/conda-forge/linux-64/libpsl-0.23.1-hf670292_0.conda + sha256: 2e6405feb59e3f40a5a4641dfa73afda158046893aa73d478e23415e18997ece + md5: c8217f5bdd5b018087bdea081912cda5 + depends: + - libgcc >=14 + - libstdcxx >=14 + - __glibc >=2.17,<3.0.a0 + - icu >=78.3,<79.0a0 + license: MIT + license_family: MIT + size: 72505 + timestamp: 1786443494718 - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libpsl-0.23.1-h7a62e17_0.conda sha256: 0140cb0d059ac531e476b85a543668b69b7dffec16f0d26269ab6c9b70546921 md5: 75a4cdf128d016141db61732ea462276 @@ -2479,17 +2518,17 @@ packages: license_family: BSD size: 430365 timestamp: 1780057267477 -- conda: https://conda.anaconda.org/conda-forge/linux-64/libsqlite-3.51.0-hee844dc_0.conda - sha256: 4c992dcd0e34b68f843e75406f7f303b1b97c248d18f3c7c330bdc0bc26ae0b3 - md5: 729a572a3ebb8c43933b30edcc628ceb +- conda: https://conda.anaconda.org/conda-forge/linux-64/libsqlite-3.53.4-hf4e2dac_0.conda + sha256: 72023efc207fe681e26b65fc9d668062cf0b4f0eacf3431e6eb099b95c1f2efd + md5: df088a279cd5e6fd2790b4c196434da1 depends: - __glibc >=2.17,<3.0.a0 - - icu >=75.1,<76.0a0 + - icu >=78.3,<79.0a0 - libgcc >=14 - - libzlib >=1.3.1,<2.0a0 + - libzlib >=1.3.2,<2.0a0 license: blessing - size: 945576 - timestamp: 1762299687230 + size: 964200 + timestamp: 1785016112246 - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libsqlite-3.53.4-h1ae2325_0.conda sha256: 745662565e103f290e9dc4263bbd88285082f8cf699854fe2d5f1e35a4a0d326 md5: 0e3477c0c3e718dcf2eb74ccc8f68570 @@ -2534,15 +2573,6 @@ packages: license_family: GPL size: 3898269 timestamp: 1759968103436 -- conda: https://conda.anaconda.org/conda-forge/linux-64/libstdcxx-ng-15.2.0-h4852527_7.conda - sha256: 024fd46ac3ea8032a5ec3ea7b91c4c235701a8bf0e6520fe5e6539992a6bd05f - md5: f627678cf829bd70bccf141a19c3ad3e - depends: - - libstdcxx 15.2.0 h8f9b012_7 - license: GPL-3.0-only WITH GCC-exception-3.1 - license_family: GPL - size: 29343 - timestamp: 1759968157195 - conda: https://conda.anaconda.org/conda-forge/linux-64/libuuid-2.41.2-he9a06e4_0.conda sha256: e5ec6d2ad7eef538ddcb9ea62ad4346fde70a4736342c4ad87bd713641eb9808 md5: 80c07c68d2f6870250959dcc95b209d1 @@ -2561,21 +2591,21 @@ packages: license: LGPL-2.1-or-later size: 100393 timestamp: 1702724383534 -- conda: https://conda.anaconda.org/conda-forge/linux-64/libxml2-2.15.1-h26afc86_0.conda - sha256: ec0735ae56c3549149eebd7dc22c0bed91fd50c02eaa77ff418613ddda190aa8 - md5: e512be7dc1f84966d50959e900ca121f +- conda: https://conda.anaconda.org/conda-forge/linux-64/libxml2-2.15.3-h49c6c72_0.conda + sha256: 3bc5551720c58591f6ea1146f7d1539c734ed1c40e7b9f5cb8cb7e900c509aba + md5: 995d8c8bad2a3cc8db14675a153dec2b depends: - __glibc >=2.17,<3.0.a0 - - icu >=75.1,<76.0a0 + - icu >=78.3,<79.0a0 - libgcc >=14 - libiconv >=1.18,<2.0a0 - - liblzma >=5.8.1,<6.0a0 - - libxml2-16 2.15.1 ha9997c6_0 - - libzlib >=1.3.1,<2.0a0 + - liblzma >=5.8.3,<6.0a0 + - libxml2-16 2.15.3 hca6bf5a_0 + - libzlib >=1.3.2,<2.0a0 license: MIT license_family: MIT - size: 45283 - timestamp: 1761015644057 + size: 46810 + timestamp: 1776376751152 - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libxml2-2.15.3-h5654f7c_0.conda sha256: 2fe1d8de0854342ae9cabe408b476935f82f5636e153b3b497456264dc8ff3a1 md5: 8e037d73747d6fe34e12d7bcac10cf21 @@ -2590,22 +2620,22 @@ packages: license_family: MIT size: 41102 timestamp: 1776377119495 -- conda: https://conda.anaconda.org/conda-forge/linux-64/libxml2-16-2.15.1-ha9997c6_0.conda - sha256: 71436e72a286ef8b57d6f4287626ff91991eb03c7bdbe835280521791efd1434 - md5: e7733bc6785ec009e47a224a71917e84 +- conda: https://conda.anaconda.org/conda-forge/linux-64/libxml2-16-2.15.3-hca6bf5a_0.conda + sha256: 3d44f737c5ae52d5af32682cc1530df433f401f8e58a7533926536244127572a + md5: e79d2c2f24b027aa8d5ab1b1ba3061e7 depends: - __glibc >=2.17,<3.0.a0 - - icu >=75.1,<76.0a0 + - icu >=78.3,<79.0a0 - libgcc >=14 - libiconv >=1.18,<2.0a0 - - liblzma >=5.8.1,<6.0a0 - - libzlib >=1.3.1,<2.0a0 + - liblzma >=5.8.3,<6.0a0 + - libzlib >=1.3.2,<2.0a0 constrains: - - libxml2 2.15.1 + - libxml2 2.15.3 license: MIT license_family: MIT - size: 556302 - timestamp: 1761015637262 + size: 559775 + timestamp: 1776376739004 - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libxml2-16-2.15.3-h5ef1a60_0.conda sha256: ff75b84cdb9e8d123db2fa694a8ac2c2059516b6cbc98ac21fb68e235d0fd354 md5: 19edaa53885fc8205614b03da2482282 @@ -2621,18 +2651,17 @@ packages: license_family: MIT size: 466360 timestamp: 1776377102261 -- conda: https://conda.anaconda.org/conda-forge/linux-64/libzlib-1.3.1-hb9d3cd8_2.conda - sha256: d4bfe88d7cb447768e31650f06257995601f89076080e76df55e3112d4e47dc4 - md5: edb0dca6bc32e4f4789199455a1dbeb8 +- conda: https://conda.anaconda.org/conda-forge/linux-64/libzlib-1.3.2-h25fd6f3_3.conda + sha256: eb8a0db0aa570124f7d2a93d7c7f596e3390df5e047818d873baad32985fc736 + md5: 0de0122d9570a8ab637c6b73db268389 depends: - __glibc >=2.17,<3.0.a0 - - libgcc >=13 constrains: - - zlib 1.3.1 *_2 + - zlib 1.3.2 *_3 license: Zlib license_family: Other - size: 60963 - timestamp: 1727963148474 + size: 63713 + timestamp: 1785362952714 - conda: https://conda.anaconda.org/conda-forge/osx-arm64/libzlib-1.3.2-h8088a28_3.conda sha256: a18fa5d5bac452401459f966cf0d872224e8080c4ff93c77e168d43ab42ef9d7 md5: f39288f0ea63ae962e1a2e4f355a0d75 @@ -3400,6 +3429,38 @@ packages: license_family: MIT size: 15528 timestamp: 1733710122949 +- conda: https://conda.anaconda.org/bioconda/linux-64/pysam-0.24.0-py312hf5ad864_1.conda + sha256: cea331ad0156a060b4ee9f534daee29da451439679e5b4d2d1558ddaa4102f61 + md5: f94309496d6273c25f0f0a499b75b2fb + depends: + - __glibc >=2.17,<3.0.a0 + - bzip2 >=1.0.8,<2.0a0 + - libcurl >=8.20.0,<9.0a0 + - libdeflate >=1.25,<1.26.0a0 + - libgcc >=14 + - liblzma >=5.8.3,<6.0a0 + - libzlib >=1.3.2,<2.0a0 + - openssl >=3.5.7,<4.0a0 + - python >=3.12,<3.13.0a0 + - python_abi 3.12.* *_cp312 + license: MIT + size: 3908681 + timestamp: 1781873141976 +- conda: https://conda.anaconda.org/bioconda/osx-arm64/pysam-0.24.0-py312h12d0683_1.conda + sha256: bab067cb5ad9c4603783533b55c0ae8c999e2ae23183a4d150c3c71c9036840f + md5: 6959148824e8c370a1578914aa2f7d0c + depends: + - __osx >=11.0 + - bzip2 >=1.0.8,<2.0a0 + - libcurl >=8.20.0,<9.0a0 + - libdeflate >=1.25,<1.26.0a0 + - liblzma >=5.8.3,<6.0a0 + - libzlib >=1.3.2,<2.0a0 + - python >=3.12,<3.13.0a0 + - python_abi 3.12.* *_cp312 + license: MIT + size: 3660664 + timestamp: 1781872373639 - conda: https://conda.anaconda.org/conda-forge/noarch/pysocks-1.7.1-pyha55dd90_7.conda sha256: ba3b032fa52709ce0d9fd388f63d330a026754587a2f461117cac9ab73d8d0d8 md5: 461219d1a5bd61342293efa2c0c90eac diff --git a/pixi.toml b/pixi.toml index ff7499ff9..9e1d148a5 100644 --- a/pixi.toml +++ b/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 = "*" diff --git a/tests/README.md b/tests/README.md new file mode 100644 index 000000000..36e475cf6 --- /dev/null +++ b/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. | diff --git a/tests/config/err-config-ref-only.yaml b/tests/config/err-config-ref-only.yaml new file mode 100644 index 000000000..9e84fb19f --- /dev/null +++ b/tests/config/err-config-ref-only.yaml @@ -0,0 +1,2 @@ +manifest: test.tbl +ref: test.fa.gz diff --git a/tests/config/err-dup-sample.tbl b/tests/config/err-dup-sample.tbl new file mode 100644 index 000000000..50cc1fa84 --- /dev/null +++ b/tests/config/err-dup-sample.tbl @@ -0,0 +1,3 @@ +sample bam +test test.cram +test test.cram diff --git a/tests/config/err-dup-sample.yaml b/tests/config/err-dup-sample.yaml new file mode 100644 index 000000000..487ec70bd --- /dev/null +++ b/tests/config/err-dup-sample.yaml @@ -0,0 +1,3 @@ +manifest: ../tests/config/err-dup-sample.tbl +ref: test.fa.gz +ref_name: hg38 diff --git a/tests/config/err-keepchrs.yaml b/tests/config/err-keepchrs.yaml new file mode 100644 index 000000000..97f254ab3 --- /dev/null +++ b/tests/config/err-keepchrs.yaml @@ -0,0 +1,4 @@ +ref: test.fa.gz +ref_name: hg38 +manifest: test.tbl +keep_chromosomes: "chrX" diff --git a/tests/config/err-missing-bam.tbl b/tests/config/err-missing-bam.tbl new file mode 100644 index 000000000..ef0be9158 --- /dev/null +++ b/tests/config/err-missing-bam.tbl @@ -0,0 +1,2 @@ +sample bam ref ref_name +test nope.cram test.fa.gz hg38 diff --git a/tests/config/err-missing-bam.yaml b/tests/config/err-missing-bam.yaml new file mode 100644 index 000000000..912850091 --- /dev/null +++ b/tests/config/err-missing-bam.yaml @@ -0,0 +1 @@ +manifest: ../tests/config/err-missing-bam.tbl diff --git a/tests/config/err-missing-fai.tbl b/tests/config/err-missing-fai.tbl new file mode 100644 index 000000000..08caed378 --- /dev/null +++ b/tests/config/err-missing-fai.tbl @@ -0,0 +1,2 @@ +sample bam ref ref_name +test test.cram ../tests/config/no-fai.fa hg38 diff --git a/tests/config/err-missing-fai.yaml b/tests/config/err-missing-fai.yaml new file mode 100644 index 000000000..54d2f4923 --- /dev/null +++ b/tests/config/err-missing-fai.yaml @@ -0,0 +1 @@ +manifest: ../tests/config/err-missing-fai.tbl diff --git a/tests/config/err-missing-refpath.tbl b/tests/config/err-missing-refpath.tbl new file mode 100644 index 000000000..83d42911e --- /dev/null +++ b/tests/config/err-missing-refpath.tbl @@ -0,0 +1,2 @@ +sample bam ref ref_name +test test.cram nope.fa hg38 diff --git a/tests/config/err-missing-refpath.yaml b/tests/config/err-missing-refpath.yaml new file mode 100644 index 000000000..4bbcfa376 --- /dev/null +++ b/tests/config/err-missing-refpath.yaml @@ -0,0 +1 @@ +manifest: ../tests/config/err-missing-refpath.tbl diff --git a/tests/config/err-nan-cell.tbl b/tests/config/err-nan-cell.tbl new file mode 100644 index 000000000..7f2924577 --- /dev/null +++ b/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 diff --git a/tests/config/err-nan-cell.yaml b/tests/config/err-nan-cell.yaml new file mode 100644 index 000000000..906d3274a --- /dev/null +++ b/tests/config/err-nan-cell.yaml @@ -0,0 +1 @@ +manifest: ../tests/config/err-nan-cell.tbl diff --git a/tests/config/err-no-ref-anywhere.yaml b/tests/config/err-no-ref-anywhere.yaml new file mode 100644 index 000000000..0cb6d126b --- /dev/null +++ b/tests/config/err-no-ref-anywhere.yaml @@ -0,0 +1 @@ +manifest: test.tbl diff --git a/tests/config/err-ref-col-only.tbl b/tests/config/err-ref-col-only.tbl new file mode 100644 index 000000000..dec20c0a8 --- /dev/null +++ b/tests/config/err-ref-col-only.tbl @@ -0,0 +1,2 @@ +sample bam ref +test test.cram test.fa.gz diff --git a/tests/config/err-ref-col-only.yaml b/tests/config/err-ref-col-only.yaml new file mode 100644 index 000000000..50f5a8801 --- /dev/null +++ b/tests/config/err-ref-col-only.yaml @@ -0,0 +1 @@ +manifest: ../tests/config/err-ref-col-only.tbl diff --git a/tests/config/err-sentinel-no-config.tbl b/tests/config/err-sentinel-no-config.tbl new file mode 100644 index 000000000..b5a0612b2 --- /dev/null +++ b/tests/config/err-sentinel-no-config.tbl @@ -0,0 +1,2 @@ +sample bam ref ref_name +test test.cram . . diff --git a/tests/config/err-sentinel-no-config.yaml b/tests/config/err-sentinel-no-config.yaml new file mode 100644 index 000000000..d8cb5ece0 --- /dev/null +++ b/tests/config/err-sentinel-no-config.yaml @@ -0,0 +1 @@ +manifest: ../tests/config/err-sentinel-no-config.tbl diff --git a/tests/config/extra-exclude.bed b/tests/config/extra-exclude.bed new file mode 100644 index 000000000..44655c5ff --- /dev/null +++ b/tests/config/extra-exclude.bed @@ -0,0 +1 @@ +chr20 0 10000 diff --git a/tests/config/multi.tbl b/tests/config/multi.tbl new file mode 100644 index 000000000..cd6f818fa --- /dev/null +++ b/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 diff --git a/tests/config/multi.yaml b/tests/config/multi.yaml new file mode 100644 index 000000000..24ff66189 --- /dev/null +++ b/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 diff --git a/tests/config/no-excludes.yaml b/tests/config/no-excludes.yaml new file mode 100644 index 000000000..cc1265e71 --- /dev/null +++ b/tests/config/no-excludes.yaml @@ -0,0 +1 @@ +manifest: ../tests/config/two-sample.tbl diff --git a/tests/config/no-fai.fa b/tests/config/no-fai.fa new file mode 100644 index 000000000..f5fc23851 --- /dev/null +++ b/tests/config/no-fai.fa @@ -0,0 +1,2 @@ +>tiny +ACGT diff --git a/tests/config/override-info.yaml b/tests/config/override-info.yaml new file mode 100644 index 000000000..e77abe1d5 --- /dev/null +++ b/tests/config/override-info.yaml @@ -0,0 +1,3 @@ +manifest: ../tests/config/single-4col.tbl +ref: test.fa.gz +ref_name: hg38 diff --git a/tests/config/sentinel.tbl b/tests/config/sentinel.tbl new file mode 100644 index 000000000..b5a0612b2 --- /dev/null +++ b/tests/config/sentinel.tbl @@ -0,0 +1,2 @@ +sample bam ref ref_name +test test.cram . . diff --git a/tests/config/sentinel.yaml b/tests/config/sentinel.yaml new file mode 100644 index 000000000..13f516a20 --- /dev/null +++ b/tests/config/sentinel.yaml @@ -0,0 +1,3 @@ +manifest: ../tests/config/sentinel.tbl +ref: test.fa.gz +ref_name: hg38 diff --git a/tests/config/single-4col.tbl b/tests/config/single-4col.tbl new file mode 100644 index 000000000..306a70b75 --- /dev/null +++ b/tests/config/single-4col.tbl @@ -0,0 +1,2 @@ +sample bam ref ref_name +test test.cram test.fa.gz hg38 diff --git a/tests/config/single-4col.yaml b/tests/config/single-4col.yaml new file mode 100644 index 000000000..0c7cde333 --- /dev/null +++ b/tests/config/single-4col.yaml @@ -0,0 +1 @@ +manifest: ../tests/config/single-4col.tbl diff --git a/tests/config/two-sample.tbl b/tests/config/two-sample.tbl new file mode 100644 index 000000000..040266d51 --- /dev/null +++ b/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 diff --git a/tests/config/two-sample.yaml b/tests/config/two-sample.yaml new file mode 100644 index 000000000..31cc2a867 --- /dev/null +++ b/tests/config/two-sample.yaml @@ -0,0 +1,3 @@ +manifest: ../tests/config/two-sample.tbl +excludes: + - ../tests/config/extra-exclude.bed diff --git a/tests/dry-run-checks.sh b/tests/dry-run-checks.sh new file mode 100644 index 000000000..d97ad04cc --- /dev/null +++ b/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" diff --git a/tests/generate-test-bams.py b/tests/generate-test-bams.py new file mode 100644 index 000000000..d734c2f67 --- /dev/null +++ b/tests/generate-test-bams.py @@ -0,0 +1,69 @@ +"""Make the generated test BAMs from the test CRAM. + +Two BAMs, each skipped when already present: + +- generated/test-chr20.bam: header lists only chr20, so the sample + exercises the BAM-contigs-are-a-subset-of-the-fasta case. +- generated/test-rev.bam: header lists chr21 before chr20, so the sample + exercises header-order preservation (header order differs from + lexicographic order). + +Run with the working directory set to fire-test-data (see the test-multi +pixi task). +""" + +from pathlib import Path + +import pysam + +CRAM = "test.cram" +REF = "test.fa.gz" +OUT_DIR = Path("generated") + + +def remap(read: pysam.AlignedSegment, out: pysam.AlignmentFile) -> pysam.AlignedSegment: + """Point the read's tids at `out`'s header by contig name.""" + if read.reference_name is not None: + read.reference_id = out.get_tid(read.reference_name) + if read.next_reference_id >= 0 and read.next_reference_name is not None: + read.next_reference_id = out.get_tid(read.next_reference_name) + return read + + +def make_chr20() -> None: + out = OUT_DIR / "test-chr20.bam" + if out.exists(): + print(f"{out} already present, skipping generation") + return + with pysam.AlignmentFile(CRAM, "rc", reference_filename=REF) as cram: + header = cram.header.to_dict() + header["SQ"] = [sq for sq in header["SQ"] if sq["SN"] == "chr20"] + with pysam.AlignmentFile(out, "wb", header=header) as bam: + for read in cram.fetch("chr20"): + bam.write(remap(read, bam)) + pysam.index(str(out)) + print(f"wrote {out}") + + +def make_reversed() -> None: + out = OUT_DIR / "test-rev.bam" + if out.exists(): + print(f"{out} already present, skipping generation") + return + unsorted = OUT_DIR / "test-rev.unsorted.bam" + with pysam.AlignmentFile(CRAM, "rc", reference_filename=REF) as cram: + header = cram.header.to_dict() + header["SQ"] = list(reversed(header["SQ"])) + with pysam.AlignmentFile(unsorted, "wb", header=header) as bam: + for read in cram.fetch(until_eof=True): + bam.write(remap(read, bam)) + pysam.sort("-o", str(out), str(unsorted)) + unsorted.unlink() + pysam.index(str(out)) + print(f"wrote {out}") + + +if __name__ == "__main__": + OUT_DIR.mkdir(exist_ok=True) + make_chr20() + make_reversed() diff --git a/tests/verify-multi.sh b/tests/verify-multi.sh new file mode 100644 index 000000000..e4222e707 --- /dev/null +++ b/tests/verify-multi.sh @@ -0,0 +1,52 @@ +#!/usr/bin/env bash +# Checks for the multi-sample run (four samples, per-sample references). +# Run with the working directory set to fire-test-data after `pixi run test-multi`. +set -euo pipefail + +V="v$(echo "$PIXI_PROJECT_VERSION" | cut -d. -f1-2)" +FAILURES=0 + +assert_genome() { + local sm=$1 name=$2 + local genomes="results/$sm/trackHub-$V/genomes.txt" + if ! grep -qx "genome $name" "$genomes"; then + echo "FAIL: $genomes does not contain 'genome $name'" + FAILURES=$((FAILURES + 1)) + else + echo "ok: $sm uses genome $name" + fi +} + +assert_genome test hg38 +assert_genome test2 testref +assert_genome test20 testref2 +assert_genome test-rev hg38 + +# the chr20-only sample must have no chr21 peaks; grep -c reads the whole +# stream, so no SIGPIPE can mask a hit under pipefail +peaks20="results/test20/test20-fire-$V-peaks.bed.gz" +chr21_count=$(gunzip -c "$peaks20" | cut -f 1 | { grep -cx "chr21" || true; }) +if [ "$chr21_count" -gt 0 ]; then + echo "FAIL: $peaks20 contains $chr21_count chr21 records" + FAILURES=$((FAILURES + 1)) +else + echo "ok: test20 output is chr20 only" +fi + +# the reversed-header sample must keep bam header order (chr21 first) in +# concatenated per-chromosome outputs +pileup_rev="results/test-rev/test-rev-fire-$V-pileup.bed.gz" +# || true absorbs the SIGPIPE that head sends up the pipeline +first_chrom=$(gunzip -c "$pileup_rev" | grep -v "^#" | head -n 1 | cut -f 1 || true) +if [ "$first_chrom" != "chr21" ]; then + echo "FAIL: $pileup_rev starts with $first_chrom, expected chr21 (header order)" + FAILURES=$((FAILURES + 1)) +else + echo "ok: test-rev output preserves bam header order" +fi + +if [ "$FAILURES" -gt 0 ]; then + echo "verify-multi: $FAILURES failure(s)" + exit 1 +fi +echo "verify-multi: all checks passed" diff --git a/tests/verify-results.sh b/tests/verify-results.sh new file mode 100644 index 000000000..1566c2d29 --- /dev/null +++ b/tests/verify-results.sh @@ -0,0 +1,28 @@ +#!/usr/bin/env bash +# Golden-count checks for the single-sample regression test. +# Run with the working directory set to fire-test-data after `pixi run test`. +set -euo pipefail + +V="v$(echo "$PIXI_PROJECT_VERSION" | cut -d. -f1-2)" +PEAKS="results/test/test-fire-$V-peaks.bed.gz" +ELEMENTS="results/test/additional-outputs-$V/fire-peaks/test-$V-fire-elements.bed.gz" +GENOMES="results/test/trackHub-$V/genomes.txt" + +peaks=$(gunzip -c "$PEAKS" | wc -l | tr -d ' ') +elements=$(gunzip -c "$ELEMENTS" | wc -l | tr -d ' ') +echo "peaks file lines: $peaks (expected 1181)" +echo "fire elements: $elements (expected 46223)" + +FAILURES=0 +[ "$peaks" -eq 1181 ] || FAILURES=$((FAILURES + 1)) +[ "$elements" -eq 46223 ] || FAILURES=$((FAILURES + 1)) +if ! grep -qx "genome hg38" "$GENOMES"; then + echo "FAIL: $GENOMES does not contain 'genome hg38'" + FAILURES=$((FAILURES + 1)) +fi + +if [ "$FAILURES" -gt 0 ]; then + echo "verify-results: $FAILURES failure(s)" + exit 1 +fi +echo "verify-results: all checks passed" diff --git a/workflow/Snakefile b/workflow/Snakefile index 4cb4ca372..1cf3304eb 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -1,5 +1,6 @@ -# pandas is a dependency of snakemake so a runner environment is not needed -# if I ever add this back because of new depedancies this must be the first line of the Snakefile +# pandas and pysam come from the pixi runner environment (pixi.toml [dependencies]); +# envs/runner.yaml mirrors them for users who run snakemake without pixi. +# If a conda runner env is ever re-enabled, this must be the first line of the Snakefile: # conda: "envs/runner.yaml" # TODO get the container working @@ -32,11 +33,7 @@ if config.get("full-version", False): MAX_THREADS = config.get("max_threads", 4) SORT_THREADS = config.get("sort_threads", 8) -# reference genome and reference regions -REF = get_ref() -FAI = get_fai() -REF_NAME = config["ref_name"] -EXCLUDES = get_excludes() +# chromosome filtering; references are per-sample (see the manifest columns) KEEP_CHRS = config.get("keep_chromosomes", ".*") # coverage requirements @@ -44,9 +41,9 @@ MIN_COVERAGE = config.get("min_coverage", 4) COVERAGE_WITHIN_N_SD = config.get("coverage_within_n_sd", 5) # sample, haplotype, and chromosome wildcard building -FAI_DF = get_fai_df() DEFAULT_ENV = config.get("env", "../envs/env.yaml") MANIFEST = get_manifest() +remove_stale_genome_files() MIN_FIRE_FDR = config.get("min_fire_fdr", 0.10) # FDR / peak calling thresholds @@ -100,9 +97,9 @@ include: "rules/track-hub.smk" wildcard_constraints: - chrom="|".join(get_chroms()), + chrom="|".join(re.escape(chrom) for chrom in all_chroms()), call="|".join(["msp", "m6a"]), - sm="|".join(MANIFEST.index), + sm="|".join(re.escape(sm) for sm in MANIFEST.index), types="|".join(types), fdr=r"\d+", hp="|".join(haps), diff --git a/workflow/envs/runner.yaml b/workflow/envs/runner.yaml index dddf584ec..b13f49efe 100644 --- a/workflow/envs/runner.yaml +++ b/workflow/envs/runner.yaml @@ -5,7 +5,4 @@ channels: dependencies: - numpy - pandas==2.2.3 -#- tqdm -#- pip -#- pip: -#- pysam==0.22.1 + - pysam diff --git a/workflow/rules/apply-model.smk b/workflow/rules/apply-model.smk index af721bfd2..18ec9ad50 100644 --- a/workflow/rules/apply-model.smk +++ b/workflow/rules/apply-model.smk @@ -4,7 +4,7 @@ rule fire: input: bam=ancient(get_input_bam), - ref=ancient(REF), + ref=lambda wc: ancient(get_ref(wc)), output: cram="results/{sm}/{sm}-fire-{v}-filtered.cram", crai="results/{sm}/{sm}-fire-{v}-filtered.cram.crai", @@ -75,7 +75,7 @@ rule fire_sites_chrom: rule fire_sites: input: beds=expand( - rules.fire_sites_chrom.output.bed, chrom=get_chroms(), allow_missing=True + rules.fire_sites_chrom.output.bed, chrom=get_chroms, allow_missing=True ), output: bed="results/{sm}/additional-outputs-{v}/fire-peaks/{sm}-{v}-fire-elements.bed.gz", diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk index 2dec22b79..00211edf4 100644 --- a/workflow/rules/common.smk +++ b/workflow/rules/common.smk @@ -1,29 +1,95 @@ import re import logging import sys +from functools import lru_cache + +import pysam + +# marks a manifest cell that falls back to the config value +MANIFEST_NA = "." + + +@lru_cache(maxsize=None) +def _bam_contigs(bam): + """Read (name, length) pairs from the header of a BAM/CRAM, once per file.""" + verbosity = pysam.set_verbosity(0) + try: + with pysam.AlignmentFile(bam, check_sq=False, require_index=False) as f: + contigs = tuple(zip(f.references, f.lengths)) + except (OSError, ValueError) as e: + raise ValueError(f"FIRE: cannot read input bam {bam}: {e}") from e + finally: + pysam.set_verbosity(verbosity) + if not contigs: + raise ValueError( + f"FIRE: input bam {bam} has no reference sequences (@SQ) in its header; " + "FIRE requires an aligned bam" + ) + return contigs + -FIRST_REPORT = True +@lru_cache(maxsize=None) +def _sample_chroms(sm): + """Filtered chromosome names for one sample, in BAM header order. + Header order is the sort order of the data (mosdepth output, coordinate + sorted reads), so it must be preserved for the sorted bedtools + operations downstream. + """ + min_contig_length = config.get("min_contig_length", 0) + try: + contigs = _bam_contigs(MANIFEST.loc[sm, "bam"]) + except ValueError as e: + raise ValueError(f"{e} (sample '{sm}')") from e + skipped = [name for name, length in contigs if length < min_contig_length] + if skipped: + print( + f"WARNING: {sm}: skipping contigs with length < {min_contig_length:,}: {skipped}", + file=sys.stderr, + ) + chroms = tuple( + name + for name, length in contigs + if length >= min_contig_length + and "chrUn_" not in name + and "_random" not in name + and re.fullmatch(KEEP_CHRS, name) + ) + print(f"INFO: {sm}: using {len(chroms)} chromosomes", file=sys.stderr) + if not chroms: + raise ValueError( + f"FIRE: no chromosomes left for sample '{sm}' after filtering. " + "Check the keep_chromosomes and min_contig_length options in config.yaml. " + f"The bam header contains: {[name for name, _ in contigs]}" + ) + return chroms -def get_ref(): - if "ref" not in config: - raise ValueError("FIRE: ref parameter is missing in config.yaml") - ref = config["ref"] - if not os.path.exists(ref): - raise ValueError(f"FIRE: reference file {ref} does not exist") - return os.path.abspath(ref) +def get_ref(wc): + return MANIFEST.loc[wc.sm, "ref"] -def get_fai(): - fai = f"{get_ref()}.fai" - if not os.path.exists(fai): - raise ValueError(f"FIRE: reference index file {fai} does not exist") - return fai +def get_fai(wc): + return f"{get_ref(wc)}.fai" -def get_excludes(): - excludes = config.get("excludes", []) - if REF_NAME == "hg38" or REF_NAME == "GRCh38": + +def get_ref_name(wc): + return MANIFEST.loc[wc.sm, "ref_name"] + + +def get_chroms(wc): + return list(_sample_chroms(wc.sm)) + + +def all_chroms(): + # sorted only for a deterministic wildcard-constraint regex; the + # alternation order has no effect on matching + return sorted({chrom for sm in MANIFEST.index for chrom in _sample_chroms(sm)}) + + +def get_excludes(wc): + excludes = list(config.get("excludes", [])) + if get_ref_name(wc) in ["hg38", "GRCh38"]: files = [ "../annotations/hg38.gap.bed.gz", "../annotations/hg38.blacklist.ENCFF356LFX.bed.gz", @@ -33,47 +99,110 @@ def get_excludes(): return excludes -def get_fai_df(): - fai = get_fai() - return pd.read_csv(fai, sep="\t", names=["chr", "length", "x", "y", "z"]) +def _config_ref_value(col): + """The config value for ref/ref_name, or None when absent or empty.""" + value = config.get(col) + if value is None or str(value) == "": + return None + return str(value) -def get_chroms(): - global FIRST_REPORT - min_contig_length = config.get("min_contig_length", 0) - skipped_contigs = FAI_DF["chr"][FAI_DF["length"] < min_contig_length] - if len(skipped_contigs) > 0 and FIRST_REPORT: +def _fill_manifest_refs(manifest): + """Fill and validate the ref and ref_name manifest columns.""" + empty_in_config = [ + col + for col in ["ref", "ref_name"] + if col in config and _config_ref_value(col) is None + ] + if empty_in_config: + raise ValueError( + f"FIRE: config options {empty_in_config} are empty in config.yaml; " + "set a value or remove the key" + ) + in_manifest = [col for col in ["ref", "ref_name"] if col in manifest.columns] + in_config = [ + col for col in ["ref", "ref_name"] if _config_ref_value(col) is not None + ] + if len(in_manifest) == 1: + raise ValueError( + "FIRE: manifest columns 'ref' and 'ref_name' must be provided together " + f"(found only '{in_manifest[0]}')" + ) + if len(in_config) == 1: + raise ValueError( + "FIRE: config options 'ref' and 'ref_name' must be provided together " + f"(found only '{in_config[0]}')" + ) + if not in_manifest and not in_config: + raise ValueError( + "FIRE: no reference specified: add 'ref' and 'ref_name' columns to the " + "manifest, or set 'ref' and 'ref_name' in config.yaml" + ) + if in_manifest and in_config: print( - f"WARNING: Skipping contigs with length < {min_contig_length:,}: {skipped_contigs}", + "INFO: manifest ref/ref_name columns override config-level values", file=sys.stderr, ) + for col in ["ref", "ref_name"]: + if col not in manifest.columns: + manifest[col] = _config_ref_value(col) + continue + sentinel = manifest[col] == MANIFEST_NA + if sentinel.any(): + if _config_ref_value(col) is None: + missing = manifest.index[sentinel].tolist() + raise ValueError( + f"FIRE: samples {missing} use '{MANIFEST_NA}' for '{col}' in the " + f"manifest, but '{col}' is not set in config.yaml" + ) + manifest.loc[sentinel, col] = _config_ref_value(col) + return manifest - chroms = FAI_DF["chr"][FAI_DF["length"] >= min_contig_length] - chroms = sorted([chrom for chrom in chroms if "chrUn_" not in chrom]) - chroms = [chrom for chrom in chroms if "_random" not in chrom] - chroms = [chrom for chrom in chroms if re.fullmatch(KEEP_CHRS, chrom)] - - if FIRST_REPORT: - FIRST_REPORT = False - print(f"INFO: Using N chromosomes: {len(chroms)}", file=sys.stderr) - if len(chroms) == 0: +def get_manifest(): + manifest_path = config.get("manifest") + if manifest_path is None: + raise ValueError("FIRE: manifest parameter is missing in config.yaml") + if not os.path.exists(manifest_path): + raise ValueError(f"FIRE: manifest file {manifest_path} does not exist") + try: + # dtype=str + keep_default_na=False keep every cell as literal text: + # numeric sample names stay strings, a sample named NA stays "NA", + # and missing trailing cells parse as "" instead of NaN + manifest = pd.read_csv( + manifest_path, sep=r"\s+", comment="#", dtype=str, keep_default_na=False + ) + except (pd.errors.ParserError, pd.errors.EmptyDataError) as e: + raise ValueError(f"FIRE: cannot parse manifest {manifest_path}: {e}") from e + for col in ["sample", "bam"]: + if col not in manifest.columns: + raise ValueError( + f"FIRE: manifest must have 'sample' and 'bam' columns; " + f"found: {list(manifest.columns)}" + ) + if len(manifest) == 0: + raise ValueError(f"FIRE: manifest {manifest_path} has no samples") + dups = manifest["sample"][manifest["sample"].duplicated()].tolist() + if dups: + raise ValueError(f"FIRE: duplicate sample names in manifest: {dups}") + for sm in manifest["sample"]: + if not re.fullmatch(r"[A-Za-z0-9_.-]+", sm): + raise ValueError(f"FIRE: sample name '{sm}' must match [A-Za-z0-9_.-]+") + manifest = manifest.set_index("sample") + manifest = _fill_manifest_refs(manifest) + malformed = manifest.index[manifest[["bam", "ref", "ref_name"]].eq("").any(axis=1)] + if len(malformed) > 0: raise ValueError( - f"No chromosomes left after filtering. Check your keep_chromosomes parameter in config.yaml. " - f"Your fai file contains the following chromosomes: {FAI_DF['chr']}" + f"FIRE: samples {malformed.tolist()} have missing or malformed manifest " + "fields; every row must fill all manifest columns " + f"(use '{MANIFEST_NA}' in ref/ref_name to fall back to the config value)" ) - return chroms - - -def get_manifest(): - manifest = config.get("manifest") - if manifest is None: - raise ValueError("manifest parameter is missing in config.yaml") - if not os.path.exists(manifest): - raise ValueError(f"Manifest file {manifest} does not exist") - manifest = pd.read_csv(config["manifest"], sep=r"\s+", comment="#").set_index( - "sample" - ) + manifest["ref"] = manifest["ref"].map(os.path.abspath) + for ref in manifest["ref"].unique(): + if not os.path.isfile(ref): + raise ValueError(f"FIRE: reference file {ref} does not exist") + if not os.path.isfile(f"{ref}.fai"): + raise ValueError(f"FIRE: reference index file {ref}.fai does not exist") return manifest @@ -81,6 +210,34 @@ def get_input_bam(wc): return MANIFEST.loc[wc.sm, "bam"] +def genome_file_content(sm): + """The genome (chrom sizes) file text for one sample: all bam header + contigs, in header order.""" + return "".join( + f"{name}\t{length}\n" for name, length in _bam_contigs(MANIFEST.loc[sm, "bam"]) + ) + + +def remove_stale_genome_files(): + """Delete leftover genome files that no longer match the bam header. + + The genome_file rule has an ancient() input, so a leftover temp file + from an interrupted run is not regenerated when the manifest bam is + swapped. Every parse reads the live bam headers anyway, so a + mismatched file is removed here, and the missing output then forces + regeneration before any consumer runs. + """ + for sm in MANIFEST.index: + path = f"temp/{sm}/{sm}.genome" + if not os.path.exists(path): + continue + with open(path) as f: + existing = f.read() + if existing != genome_file_content(sm): + print(f"INFO: {sm}: removing stale genome file {path}", file=sys.stderr) + os.remove(path) + + def get_mem_mb(wildcards, attempt): if attempt < 3: return attempt * 1024 * 32 diff --git a/workflow/rules/coverages.smk b/workflow/rules/coverages.smk index ec3b6f955..9fc949259 100644 --- a/workflow/rules/coverages.smk +++ b/workflow/rules/coverages.smk @@ -1,10 +1,26 @@ +rule genome_file: + input: + # ancient for the same reason as rule fire: a re-copied bam with a + # fresh mtime must not cascade reruns through every genome-file + # consumer. Staleness is impossible: remove_stale_genome_files() + # deletes any leftover file that no longer matches the bam header + # at parse time, which forces regeneration. + bam=lambda wc: ancient(get_input_bam(wc)), + output: + genome=temp("temp/{sm}/{sm}.genome"), + localrule: True + run: + with open(output.genome, "w") as out: + out.write(genome_file_content(wildcards.sm)) + + # # Coverage calculations # rule genome_bedgraph: input: - ref=ancient(REF), - fai=ancient(FAI), + ref=lambda wc: ancient(get_ref(wc)), + fai=lambda wc: ancient(get_fai(wc)), cram=rules.fire.output.cram, crai=rules.fire.output.crai, output: @@ -17,9 +33,10 @@ rule genome_bedgraph: threads: 16 shell: """ + # mosdepth output is position sorted in bam header order; keep that + # order, it is the order of every other per-sample file mosdepth -F 4 -f {input.ref} -t {threads} tmp {input.cram} bgzip -cd tmp.per-base.bed.gz \ - | LC_ALL=C sort --parallel={threads} -k1,1 -k2,2n -k3,3n -k4,4 \ | bgzip -@ {threads} \ >{output.bg} tabix -f -p bed {output.bg} @@ -43,7 +60,7 @@ rule coverage: params: coverage_within_n_sd=COVERAGE_WITHIN_N_SD, min_coverage=MIN_COVERAGE, - chroms=get_chroms(), + chroms=get_chroms, script: "../scripts/cov.py" @@ -76,12 +93,13 @@ rule fiber_locations: input: fibers=expand( rules.fiber_locations_chromosome.output.bed, - chrom=get_chroms(), + chrom=get_chroms, allow_missing=True, ), bg=rules.genome_bedgraph.output.bg, minimum=rules.coverage.output.minimum, maximum=rules.coverage.output.maximum, + genome=rules.genome_file.output.genome, output: bed=temp("temp/{sm}/coverage/{v}-fiber-locations.bed.gz"), bed_tbi=temp("temp/{sm}/coverage/{v}-fiber-locations.bed.gz.tbi"), @@ -104,7 +122,7 @@ rule fiber_locations: # get filtered fiber locations MIN=$(cat {input.minimum}) MAX=$(cat {input.maximum}) - bedtools intersect -header -sorted -v -f {params.max_frac_overlap} \ + bedtools intersect -header -sorted -g {input.genome} -v -f {params.max_frac_overlap} \ -a {output.bed} \ -b <(bgzip -cd {input.bg} | awk -v MAX="$MAX" -v MIN="$MIN" '$4 <= MIN || $4 >= MAX') \ | bgzip -@ {threads} \ @@ -119,18 +137,18 @@ rule fiber_locations: rule exclude_from_shuffle: input: filtered=rules.fiber_locations.output.filtered, - fai=ancient(FAI), + genome=rules.genome_file.output.genome, output: bed="results/{sm}/additional-outputs-{v}/coverage/exclude-from-shuffles.bed.gz", conda: DEFAULT_ENV threads: 4 params: - exclude=lambda wc: " ".join(EXCLUDES) if EXCLUDES else "", + exclude=lambda wc: " ".join(get_excludes(wc)), shell: """ ( - bedtools genomecov -bga -i {input.filtered} -g {input.fai} | awk '$4 == 0' + bedtools genomecov -bga -i {input.filtered} -g {input.genome} | awk '$4 == 0' if [ -n "{params.exclude}" ]; then gunzip -cf {params.exclude} fi @@ -148,7 +166,7 @@ rule unreliable_coverage_regions: bg=rules.genome_bedgraph.output.bg, minimum=rules.coverage.output.minimum, maximum=rules.coverage.output.maximum, - fai=ancient(FAI), + genome=rules.genome_file.output.genome, output: bed="results/{sm}/additional-outputs-{v}/coverage/unreliable-coverage-regions.bed.gz", bed_tbi="results/{sm}/additional-outputs-{v}/coverage/unreliable-coverage-regions.bed.gz.tbi", @@ -177,7 +195,7 @@ rule unreliable_coverage_regions: bedtools merge -i {output.bed} >{output.tmp} bigtools bedtobigbed \ -s start -a {params.bed3_as} \ - {output.tmp} {input.fai} {output.bb} + {output.tmp} {input.genome} {output.bb} # index tabix -f -p bed {output.bed} diff --git a/workflow/rules/decorated-reads.smk b/workflow/rules/decorated-reads.smk index 6a92ee80d..d985f2de7 100644 --- a/workflow/rules/decorated-reads.smk +++ b/workflow/rules/decorated-reads.smk @@ -32,10 +32,10 @@ rule decorate_fibers_1: input: bed=expand( rules.decorate_fibers_chromosome.output.bed, - chrom=get_chroms(), + chrom=get_chroms, allow_missing=True, ), - fai=ancient(FAI), + genome=rules.genome_file.output.genome, output: #bed=temp("temp/{sm}/fiber-calls/fire-fibers.bed.gz"), bb="results/{sm}/trackHub-{v}/bb/fire-fibers.bb", @@ -62,7 +62,7 @@ rule decorate_fibers_1: --block-size {params.block_size} --items-per-slot {params.items_per_slot} \ --nzooms {params.nzooms} \ -s start -a {params.bed_as} \ - - {input.fai} {output.bb} + - {input.genome} {output.bb} """ @@ -70,10 +70,10 @@ rule decorate_fibers_2: input: decorated=expand( rules.decorate_fibers_chromosome.output.decorated, - chrom=get_chroms(), + chrom=get_chroms, allow_missing=True, ), - fai=ancient(FAI), + genome=rules.genome_file.output.genome, output: bb="results/{sm}/trackHub-{v}/bb/fire-fiber-decorators.bb", #bed=temp("temp/{sm}/trackHub-{v}/bb/fire-fiber-decorators.bed.gz"), @@ -104,7 +104,7 @@ rule decorate_fibers_2: --block-size {params.block_size} --items-per-slot {params.items_per_slot} \ --nzooms {params.nzooms} \ -s start -a {params.dec_as} \ - - {input.fai} {output.bb} + - {input.genome} {output.bb} """ @@ -114,5 +114,5 @@ if False: cat {input.decorated} > {output.bed} bedToBigBed \ -allow1bpOverlap -type=bed12+ -as={params.dec_as} \ - {output.bed} {input.fai} {output.bb} + {output.bed} {input.genome} {output.bb} """ diff --git a/workflow/rules/fire-peaks.smk b/workflow/rules/fire-peaks.smk index fbccff0e5..a0e2d2dbd 100644 --- a/workflow/rules/fire-peaks.smk +++ b/workflow/rules/fire-peaks.smk @@ -3,7 +3,7 @@ rule filtered_and_shuffled_fiber_locations_chromosome: filtered=rules.fiber_locations.output.filtered, filtered_tbi=rules.fiber_locations.output.filtered_tbi, exclude=rules.exclude_from_shuffle.output.bed, - fai=ancient(FAI), + genome=rules.genome_file.output.genome, output: shuffled=temp("temp/{sm}/shuffle/{v}-{chrom}.fiber-locations-shuffled.bed.gz"), conda: @@ -16,7 +16,7 @@ rule filtered_and_shuffled_fiber_locations_chromosome: | bedtools shuffle -chrom -seed 42 \ -excl {input.exclude} \ -i - \ - -g {input.fai} \ + -g {input.genome} \ | sort -k1,1 -k2,2n -k3,3n -k4,4 \ | bgzip -@ {threads} \ >{output.shuffled} @@ -46,7 +46,7 @@ rule shuffled_pileup: input: beds=expand( rules.shuffled_pileup_chromosome.output.bed, - chrom=get_chroms(), + chrom=get_chroms, allow_missing=True, ), output: @@ -137,7 +137,7 @@ rule pileup: input: beds=expand( rules.fdr_track_chromosome.output.bed, - chrom=get_chroms(), + chrom=get_chroms, allow_missing=True, ), output: @@ -250,7 +250,7 @@ rule fire_peaks: input: beds=expand( rules.fdr_peaks_by_fire_elements_chromosome.output.bed, - chrom=get_chroms(), + chrom=get_chroms, allow_missing=True, ), output: @@ -282,7 +282,7 @@ rule wide_fire_peaks: input: bed=rules.fire_peaks.output.bed, track=rules.pileup.output.bed, - fai=ancient(FAI), + genome=rules.genome_file.output.genome, output: bed="results/{sm}/additional-outputs-{v}/fire-peaks/{sm}-fire-{v}-wide-peaks.bed.gz", tbi="results/{sm}/additional-outputs-{v}/fire-peaks/{sm}-fire-{v}-wide-peaks.bed.gz.tbi", @@ -311,7 +311,7 @@ rule wide_fire_peaks: bgzip -cd -@ 16 {output.bed} \ | bigtools bedtobigbed \ -s start -a {params.bed3_as} \ - - {input.fai} {output.bb} + - {input.genome} {output.bb} tabix -p bed {output.bed} """ diff --git a/workflow/rules/stats.smk b/workflow/rules/stats.smk index f462c8592..56cdc0764 100644 --- a/workflow/rules/stats.smk +++ b/workflow/rules/stats.smk @@ -4,7 +4,7 @@ rule clustering_vs_null: input: bed=rules.fire_sites.output.bed, - fai=ancient(FAI), + genome=rules.genome_file.output.genome, output: tmp=temp("temp/{sm}/tmp.pre.calls.bed"), null=temp("temp/{sm}/null.calls.bed"), @@ -15,11 +15,11 @@ rule clustering_vs_null: shell: """ bgzip -cd -@{threads} {input.bed} | cut -f 1-3 >{output.tmp} - bedtools shuffle -chrom -i {output.tmp} -g {input.fai} >{output.null} + bedtools shuffle -chrom -i {output.tmp} -g {input.genome} >{output.null} ( - bedtools genomecov -bg -i {output.tmp} -g {input.fai} | sed 's/$/\\tReal/g' - bedtools genomecov -bg -i {output.null} -g {input.fai} | sed 's/$/\\tNull/g' + bedtools genomecov -bg -i {output.tmp} -g {input.genome} | sed 's/$/\\tReal/g' + bedtools genomecov -bg -i {output.null} -g {input.genome} | sed 's/$/\\tNull/g' ) \ | bedtools sort \ | bgzip -@ {threads} \ diff --git a/workflow/rules/track-hub.smk b/workflow/rules/track-hub.smk index 1be0fcb55..8dd036a6b 100644 --- a/workflow/rules/track-hub.smk +++ b/workflow/rules/track-hub.smk @@ -1,7 +1,7 @@ rule percent_accessible: input: bed=rules.pileup.output.bed, - fai=ancient(FAI), + genome=rules.genome_file.output.genome, output: tmp=temp("temp/{sm}/{hp}/{v}-percent.accessible.bed"), bw="results/{sm}/trackHub-{v}/bw/{hp}.percent.accessible.bw", @@ -13,7 +13,7 @@ rule percent_accessible: params: suffix=get_hap_col_suffix, nzooms=NZOOMS, - chrom=get_chroms()[0], + chrom=lambda wc: get_chroms(wc)[0], shell: """ bgzip -cd {input.bed} \ @@ -32,14 +32,14 @@ rule percent_accessible: bigtools bedgraphtobigwig \ --nzooms {params.nzooms} -s start \ - {output.tmp} {input.fai} {output.bw} + {output.tmp} {input.genome} {output.bw} """ rule element_coverages_bw: input: bed=rules.pileup.output.bed, - fai=ancient(FAI), + genome=rules.genome_file.output.genome, output: bw="results/{sm}/trackHub-{v}/bw/{hp}.{el_type}.coverage.bw", conda: @@ -54,14 +54,14 @@ rule element_coverages_bw: | grep -v "^#" \ | bigtools bedgraphtobigwig \ -s start --nzooms {params.nzooms} \ - - {input.fai} {output.bw} + - {input.genome} {output.bw} """ rule fdr_track_to_bw: input: bed=rules.pileup.output.bed, - fai=ancient(FAI), + genome=rules.genome_file.output.genome, output: bw="results/{sm}/trackHub-{v}/bw/{col}.bw", conda: @@ -75,14 +75,14 @@ rule fdr_track_to_bw: | grep -v "^#" \ | bigtools bedgraphtobigwig \ -s start --nzooms {params.nzooms} \ - - {input.fai} {output.bw} + - {input.genome} {output.bw} """ rule fire_peaks_bb: input: bed=rules.fire_peaks.output.bed, - fai=ancient(FAI), + genome=rules.genome_file.output.genome, output: bb="results/{sm}/trackHub-{v}/bb/fire-peaks.bb", conda: @@ -98,14 +98,14 @@ rule fire_peaks_bb: | rg -v '^#' \ | bigtools bedtobigbed \ -a {params.bedfmt} -s start \ - - {input.fai} {output.bb} + - {input.genome} {output.bb} """ rule hap_differences_track: input: bed9=rules.hap_differences.output.bed9, - fai=ancient(FAI), + genome=rules.genome_file.output.genome, output: bb="results/{sm}/trackHub-{v}/bb/hap_differences.bb", conda: @@ -114,7 +114,7 @@ rule hap_differences_track: resources: mem_mb=get_mem_mb, params: - chrom=get_chroms()[0], + chrom=lambda wc: get_chroms(wc)[0], bed9_as=workflow.source_path("../templates/bed9.as"), shell: """ @@ -124,7 +124,7 @@ rule hap_differences_track: ) \ | bigtools bedtobigbed \ -s start -a {params.bed9_as} \ - - {input.fai} {output.bb} + - {input.genome} {output.bb} """ @@ -140,7 +140,7 @@ rule trackhub: resources: load=get_load, params: - ref=REF_NAME, + ref=get_ref_name, script=workflow.source_path("../scripts/trackhub.py"), description=workflow.source_path("../templates/fire-description.html"), shell: From 0c05d742d7a1eb51dc3d8a81f37440b603e45944 Mon Sep 17 00:00:00 2001 From: "Mitchell R. Vollger" Date: Fri, 14 Aug 2026 16:15:04 -0700 Subject: [PATCH 2/4] fix: anchor all whole-genome bed sorts and sorted intersects to the genome file MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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 --- workflow/rules/coverages.smk | 3 ++- workflow/rules/fire-peaks.smk | 5 +++-- workflow/rules/stats.smk | 7 ++++--- workflow/rules/track-hub.smk | 2 +- 4 files changed, 10 insertions(+), 7 deletions(-) diff --git a/workflow/rules/coverages.smk b/workflow/rules/coverages.smk index 9fc949259..c6ce51383 100644 --- a/workflow/rules/coverages.smk +++ b/workflow/rules/coverages.smk @@ -154,7 +154,8 @@ rule exclude_from_shuffle: fi ) \ | cut -f 1-3 \ - | bedtools sort \ + | awk 'NR==FNR {{keep[$1] = 1; next}} keep[$1]' {input.genome} - \ + | bedtools sort -g {input.genome} \ | bedtools merge \ | bgzip -@ {threads} \ >{output.bed} diff --git a/workflow/rules/fire-peaks.smk b/workflow/rules/fire-peaks.smk index a0e2d2dbd..ad2146861 100644 --- a/workflow/rules/fire-peaks.smk +++ b/workflow/rules/fire-peaks.smk @@ -303,7 +303,7 @@ rule wide_fire_peaks: | bioawk -tc hdr 'NR==1 || (NF>0 && $coverage>0 && $fire_coverage/$coverage>={params.min_frac_acc})' ) \ | cut -f 1-3 \ - | bedtools sort \ + | bedtools sort -g {input.genome} \ | bedtools merge -d {params.nuc_size} \ | bgzip -@ {threads} \ >{output.bed} @@ -321,6 +321,7 @@ rule one_percent_fire_peaks: input: bed=rules.fire_peaks.output.bed, track=rules.pileup.output.bed, + genome=rules.genome_file.output.genome, output: bed="results/{sm}/additional-outputs-{v}/fire-peaks/one-percent-FDR/{sm}-fire-{v}-01-fire-peaks.bed.gz", tbi="results/{sm}/additional-outputs-{v}/fire-peaks/one-percent-FDR/{sm}-fire-{v}-01-fire-peaks.bed.gz.tbi", @@ -344,7 +345,7 @@ rule one_percent_fire_peaks: bioawk -tc hdr '$FDR<=0.01' {input.track} ) \ | cut -f 1-3 \ - | bedtools sort \ + | bedtools sort -g {input.genome} \ | bedtools merge -d {params.nuc_size} \ | bgzip -@ {threads} \ >{output.wide} diff --git a/workflow/rules/stats.smk b/workflow/rules/stats.smk index 56cdc0764..ac797e5b1 100644 --- a/workflow/rules/stats.smk +++ b/workflow/rules/stats.smk @@ -21,7 +21,7 @@ rule clustering_vs_null: bedtools genomecov -bg -i {output.tmp} -g {input.genome} | sed 's/$/\\tReal/g' bedtools genomecov -bg -i {output.null} -g {input.genome} | sed 's/$/\\tNull/g' ) \ - | bedtools sort \ + | bedtools sort -g {input.genome} \ | bgzip -@ {threads} \ >{output.bed} """ @@ -32,6 +32,7 @@ rule fires_in_peaks: fire=rules.fire_sites.output.bed, exclude=rules.unreliable_coverage_regions.output.bed, peaks=rules.fire_peaks.output.bed, + genome=rules.genome_file.output.genome, output: tmp=temp("temp/{sm}/tmp.FIREs-{v}-in-peaks.bed"), txt="results/{sm}/additional-outputs-{v}/fire-peaks/{sm}-{v}-fires-in-peaks.txt", @@ -42,13 +43,13 @@ rule fires_in_peaks: script=workflow.source_path("../scripts/percent-in-clusters.sh"), shell: """ - bedtools intersect -sorted -a {input.fire} -b {input.exclude} -v >{output.tmp} + bedtools intersect -sorted -g {input.genome} -a {input.fire} -b {input.exclude} -v >{output.tmp} echo "Total # of FIREs within normal coverage regions" >>{output.txt} wc -l {output.tmp} >>{output.txt} echo "# of FIREs within peaks" >>{output.txt} - bedtools intersect -sorted -u -a {input.fire} -b {input.peaks} | wc -l >>{output.txt} + bedtools intersect -sorted -g {input.genome} -u -a {input.fire} -b {input.peaks} | wc -l >>{output.txt} """ diff --git a/workflow/rules/track-hub.smk b/workflow/rules/track-hub.smk index 8dd036a6b..2e995fce5 100644 --- a/workflow/rules/track-hub.smk +++ b/workflow/rules/track-hub.smk @@ -120,7 +120,7 @@ rule hap_differences_track: """ ( printf "{params.chrom}\t0\t1\tfake\t100\t+\t0\t1\t230,230,230\\n" - bedtools sort -i {input.bed9} + bedtools sort -g {input.genome} -i {input.bed9} ) \ | bigtools bedtobigbed \ -s start -a {params.bed9_as} \ From 0d7193730ef27d7fa6d7403785df20b522972e82 Mon Sep 17 00:00:00 2001 From: "Mitchell R. Vollger" Date: Fri, 14 Aug 2026 17:10:30 -0700 Subject: [PATCH 3/4] fix: address ordering re-review findings MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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 --- config/README.md | 4 ++-- tests/README.md | 3 ++- tests/config/err-extra-column.tbl | 2 ++ tests/config/err-extra-column.yaml | 1 + tests/dry-run-checks.sh | 1 + tests/verify-multi.sh | 20 +++++++++++++++++++ workflow/rules/common.smk | 32 +++++++++++++++++++++++------- 7 files changed, 53 insertions(+), 10 deletions(-) create mode 100644 tests/config/err-extra-column.tbl create mode 100644 tests/config/err-extra-column.yaml diff --git a/config/README.md b/config/README.md index 3a7de2e3d..6e4769f0f 100644 --- a/config/README.md +++ b/config/README.md @@ -22,7 +22,7 @@ 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. +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. 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 ``` @@ -54,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 ``` diff --git a/tests/README.md b/tests/README.md index 36e475cf6..ce9e697d0 100644 --- a/tests/README.md +++ b/tests/README.md @@ -35,7 +35,8 @@ Error cases (`err-*`, the dry run must fail with a specific message): | `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-nan-cell.yaml` + `.tbl` | Short manifest row (missing cells); rejected at parse time. | +| `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. | diff --git a/tests/config/err-extra-column.tbl b/tests/config/err-extra-column.tbl new file mode 100644 index 000000000..029893491 --- /dev/null +++ b/tests/config/err-extra-column.tbl @@ -0,0 +1,2 @@ +sample bam ref ref_name +test test.cram test.fa.gz hg38 extra diff --git a/tests/config/err-extra-column.yaml b/tests/config/err-extra-column.yaml new file mode 100644 index 000000000..31932f4f3 --- /dev/null +++ b/tests/config/err-extra-column.yaml @@ -0,0 +1 @@ +manifest: ../tests/config/err-extra-column.tbl diff --git a/tests/dry-run-checks.sh b/tests/dry-run-checks.sh index d97ad04cc..4ab5902d9 100644 --- a/tests/dry-run-checks.sh +++ b/tests/dry-run-checks.sh @@ -139,6 +139,7 @@ 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" diff --git a/tests/verify-multi.sh b/tests/verify-multi.sh index e4222e707..b884e7f77 100644 --- a/tests/verify-multi.sh +++ b/tests/verify-multi.sh @@ -45,6 +45,26 @@ else echo "ok: test-rev output preserves bam header order" fi +# test-rev holds the same reads as test with a reversed header, so output +# CONTENT must be identical modulo row order — the only assertion that can +# catch a silent (non-crashing) ordering bug +assert_same_content() { + local a=$1 b=$2 label=$3 + if diff <(gunzip -c "$a" | grep -v "^#" | LC_ALL=C sort) \ + <(gunzip -c "$b" | grep -v "^#" | LC_ALL=C sort) >/dev/null; then + echo "ok: $label content identical between test and test-rev" + else + echo "FAIL: $label content differs between test and test-rev" + FAILURES=$((FAILURES + 1)) + fi +} +assert_same_content \ + "results/test/test-fire-$V-peaks.bed.gz" \ + "results/test-rev/test-rev-fire-$V-peaks.bed.gz" "peaks" +assert_same_content \ + "results/test/test-fire-$V-pileup.bed.gz" \ + "results/test-rev/test-rev-fire-$V-pileup.bed.gz" "pileup" + if [ "$FAILURES" -gt 0 ]; then echo "verify-multi: $FAILURES failure(s)" exit 1 diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk index 00211edf4..d8959b497 100644 --- a/workflow/rules/common.smk +++ b/workflow/rules/common.smk @@ -1,6 +1,7 @@ import re import logging import sys +import warnings from functools import lru_cache import pysam @@ -167,12 +168,28 @@ def get_manifest(): raise ValueError(f"FIRE: manifest file {manifest_path} does not exist") try: # dtype=str + keep_default_na=False keep every cell as literal text: - # numeric sample names stay strings, a sample named NA stays "NA", - # and missing trailing cells parse as "" instead of NaN - manifest = pd.read_csv( - manifest_path, sep=r"\s+", comment="#", dtype=str, keep_default_na=False - ) - except (pd.errors.ParserError, pd.errors.EmptyDataError) as e: + # numeric sample names stay strings and a sample named NA stays "NA" + # (missing trailing cells still parse as NaN; the malformed-row check + # below catches both NaN and ""). index_col=False stops pandas from + # silently treating the first field as an index when every data row + # has one extra column; promoting ParserWarning to an error turns + # the resulting silent field drop into a loud failure + with warnings.catch_warnings(): + warnings.simplefilter("error", pd.errors.ParserWarning) + manifest = pd.read_csv( + manifest_path, + sep=r"\s+", + comment="#", + dtype=str, + keep_default_na=False, + index_col=False, + engine="python", + ) + except ( + pd.errors.ParserError, + pd.errors.EmptyDataError, + pd.errors.ParserWarning, + ) as e: raise ValueError(f"FIRE: cannot parse manifest {manifest_path}: {e}") from e for col in ["sample", "bam"]: if col not in manifest.columns: @@ -190,7 +207,8 @@ def get_manifest(): raise ValueError(f"FIRE: sample name '{sm}' must match [A-Za-z0-9_.-]+") manifest = manifest.set_index("sample") manifest = _fill_manifest_refs(manifest) - malformed = manifest.index[manifest[["bam", "ref", "ref_name"]].eq("").any(axis=1)] + ref_cols = manifest[["bam", "ref", "ref_name"]] + malformed = manifest.index[(ref_cols.isna() | ref_cols.eq("")).any(axis=1)] if len(malformed) > 0: raise ValueError( f"FIRE: samples {malformed.tolist()} have missing or malformed manifest " From 9d78a07636096372492e588823ad963d82b11dca Mon Sep 17 00:00:00 2001 From: "Mitchell R. Vollger" Date: Fri, 14 Aug 2026 17:39:05 -0700 Subject: [PATCH 4/4] fix: close targeted-review findings on the final fix commit 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 --- config/README.md | 2 +- psub-err.txt | 5 +++++ tests/README.md | 2 +- tests/verify-multi.sh | 8 ++++++++ workflow/rules/common.smk | 5 ++++- 5 files changed, 19 insertions(+), 3 deletions(-) create mode 100644 psub-err.txt diff --git a/config/README.md b/config/README.md index 6e4769f0f..1650630d9 100644 --- a/config/README.md +++ b/config/README.md @@ -22,7 +22,7 @@ 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. 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. +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 ``` diff --git a/psub-err.txt b/psub-err.txt new file mode 100644 index 000000000..44a1c9dc1 --- /dev/null +++ b/psub-err.txt @@ -0,0 +1,5 @@ +gunzip: can't stat: missing-a.bed.gz (missing-a.bed.gz.gz): No such file or directory +gunzip: can't stat: missing-b.bed.gz (missing-b.bed.gz.gz): No such file or directory +gunzip: can't stat: missing-b.bed.gz (missing-b.bed.gz.gz): No such file or directory +gunzip: truncated.bed.gz: unexpected end of file +gunzip: truncated.bed.gz: uncompress failed diff --git a/tests/README.md b/tests/README.md index ce9e697d0..3b10bc798 100644 --- a/tests/README.md +++ b/tests/README.md @@ -35,7 +35,7 @@ Error cases (`err-*`, the dry run must fail with a specific message): | `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); rejected at parse time. | +| `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. | diff --git a/tests/verify-multi.sh b/tests/verify-multi.sh index b884e7f77..df856e860 100644 --- a/tests/verify-multi.sh +++ b/tests/verify-multi.sh @@ -50,6 +50,14 @@ fi # catch a silent (non-crashing) ordering bug assert_same_content() { local a=$1 b=$2 label=$3 + # missing inputs must fail loudly, not compare two empty streams + for f in "$a" "$b"; do + if [ ! -f "$f" ]; then + echo "FAIL: $label missing file $f" + FAILURES=$((FAILURES + 1)) + return + fi + done if diff <(gunzip -c "$a" | grep -v "^#" | LC_ALL=C sort) \ <(gunzip -c "$b" | grep -v "^#" | LC_ALL=C sort) >/dev/null; then echo "ok: $label content identical between test and test-rev" diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk index d8959b497..42c668290 100644 --- a/workflow/rules/common.smk +++ b/workflow/rules/common.smk @@ -190,7 +190,10 @@ def get_manifest(): pd.errors.EmptyDataError, pd.errors.ParserWarning, ) as e: - raise ValueError(f"FIRE: cannot parse manifest {manifest_path}: {e}") from e + raise ValueError( + f"FIRE: cannot parse manifest {manifest_path}: {e} (check that every " + "row has the same number of whitespace-separated fields as the header)" + ) from e for col in ["sample", "bam"]: if col not in manifest.columns: raise ValueError(