Update fibertools-rs to 0.13.x and v0.2.0 housekeeping

.github/workflows/main.yml

@@ -12,7 +12,7 @@ jobs:
     runs-on: ubuntu-latest
     steps:
       - uses: actions/checkout@v4
-      - uses: prefix-dev/setup-pixi@v0.8.1
+      - uses: prefix-dev/setup-pixi@v0.9.6
         with:
           pixi-version: v0.70.2
           cache: true

CHANGELOG.md

@@ -2,6 +2,18 @@
 
 All notable changes to this project will be documented in this file.
 
+## v0.2.0
+
+- `min_frac_accessible` now defaults to 0.10: FIRE peaks must have at least 10%
+  of reads accessible in addition to passing FDR peak calling (set to 0.0 to
+  restore the old behavior)
+- Update fibertools-rs to 0.13.x (requires `ft pileup --rgn`; also updates samtools,
+  htslib, and bedtools pins in `workflow/envs/env.yaml`)
+- Polars fixes: explicit schemas for chrom/score columns (#48, #52, #53, #56)
+- `exclude_from_shuffle` now works when no `excludes` are configured (non-hg38 refs)
+- Test data now downloaded via rclone; test runs start from a clean state
+  (`test-clean` task added)
+
 ## v0.1.2
 
 - fix #34

config/README.md

@@ -57,9 +57,9 @@ Forgo the use of FDR peak calling and instead call peaks for regions with at lea
 ```
 min_per_acc_peak = 0.25
 ```
-Apply a percent actuation filter on top of the FDR peak calling. Default is `0.0` for no filter.
+Apply a percent actuation filter on top of the FDR peak calling. Default is `0.10`; set to `0.0` for no filter.
 ```
-min_frac_accessible: 0.0
+min_frac_accessible: 0.10
 ```
 Process only chromosomes matching this regular expression:
 ```

pixi.toml

@@ -1,37 +1,16 @@
 [workspace]
 authors = ["Mitchell Robert Vollger <mrvollger@gmail.com>"]
 channels = ["conda-forge", "bioconda"]
-description = "Add a short description here"
+description = "A Snakemake pipeline for calling FIRE peaks using fibertools-rs."
 name = "FIRE"
 platforms = ["osx-64", "linux-64"]
-version = "0.1.2"
+version = "0.2.0"
 
 [tasks]
 fmt = "ruff format . && taplo format pixi.toml && snakefmt workflow/"
-test-data = { cmd = [
-  "cd",
-  "$INIT_CWD",
-  "&&",
-  "mkdir",
-  "-p",
-  "fire-test-data",
-  "&&",
-  "rclone",
-  "sync",
-  ":s3,env_auth=false,provider=Other,endpoint=s3.kopah.orci.washington.edu:stergachis/public/FIRE/test-data",
-  "fire-test-data/",
-] }
-test = { cmd = [
-  "cd",
-  "$INIT_CWD/fire-test-data",
-  "&&",
-  "snakemake",
-  "-s",
-  "$PIXI_PROJECT_ROOT/workflow/Snakefile",
-  "--configfile",
-  "test.yaml",
-  "-k",
-], depends-on = [
+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",
 ], clean-env = true }
 fire = { cmd = [

workflow/Snakefile

@@ -56,7 +56,7 @@ if MIN_PER_ACC_PEAK is None:
     MIN_PER_ACC_PEAK = 0.0
 else:
     MAX_PEAK_FDR = 1.0
-MIN_FRAC_ACCESSIBLE = config.get("min_frac_accessible", 0)
+MIN_FRAC_ACCESSIBLE = config.get("min_frac_accessible", 0.10)
 
 # data filtering
 FILTER_FLAG = config.get("samtools-filter-flag", "260")  # 2308

workflow/envs/env.yaml

@@ -4,10 +4,10 @@ channels:
   - bioconda
   - defaults
 dependencies:
-  - samtools==1.19.1
-  - htslib==1.19.1
-  - bedtools==2.31
-  - bioconda::fibertools-rs==0.6
+  - samtools==1.21
+  - htslib==1.21
+  - bedtools==2.31.1
+  - bioconda::fibertools-rs>=0.13,<0.14
   - hck>=0.9.2
   - bioawk
   - ripgrep

workflow/rules/coverages.smk

@@ -126,13 +126,14 @@ rule exclude_from_shuffle:
     conda:
         DEFAULT_ENV
     params:
-        exclude=EXCLUDES,
+        exclude=lambda wc: " ".join(EXCLUDES) if EXCLUDES else "",
     shell:
         """
-
-        ( \
-            bedtools genomecov -bga -i {input.filtered} -g {input.fai} | awk '$4 == 0'; \
-            less {params.exclude} \
+        (
+            bedtools genomecov -bga -i {input.filtered} -g {input.fai} | awk '$4 == 0'
+            if [ -n "{params.exclude}" ]; then
+                gunzip -cf {params.exclude}
+            fi
         ) \
             | cut -f 1-3 \
             | bedtools sort \

workflow/rules/fire-peaks.smk

@@ -34,7 +34,7 @@ rule shuffled_pileup_chromosome:
         DEFAULT_ENV
     shell:
         """
-        {FT_EXE} pileup {input.cram} {wildcards.chrom} -t {threads} \
+        {FT_EXE} pileup {input.cram} --rgn {wildcards.chrom} -t {threads} \
             --fiber-coverage --shuffle {input.shuffled} \
             --no-msp --no-nuc \
             | bgzip -@ {threads} \
@@ -105,7 +105,7 @@ rule pileup_chromosome:
         """
         {FT_EXE} pileup -t {threads} \
             --haps --fiber-coverage \
-            {input.bam} {wildcards.chrom} \
+            {input.bam} --rgn {wildcards.chrom} \
             | bgzip -@ {threads} \
             > {output.bed}
         """

workflow/scripts/fdr-table.py

@@ -59,7 +59,7 @@ def read_pileup_file(infile, nrows):
     # add scema overrides for the score columns
     # Build schema overrides keyed by positional column names (column_1, column_2, ...)
     # because polars infers schema BEFORE new_columns is applied when has_header=False.
-    # Keying on '#chrom' / 'score' here would be silently ignored. 
+    # Keying on '#chrom' / 'score' here would be silently ignored.
     schema_overrides = {}
     for col_idx, col_name in enumerate(header, start=1):
         positional = f"column_{col_idx}"

workflow/scripts/merge_fire_peaks.py

@@ -124,7 +124,12 @@ def main(
     logger.setLevel(log_level)
 
     inf = io.StringIO(sys.stdin.read())
-    df = pl.read_csv(inf, separator="\t", null_values=".", schema_overrides={"#chrom": pl.Utf8},)
+    df = pl.read_csv(
+        inf,
+        separator="\t",
+        null_values=".",
+        schema_overrides={"#chrom": pl.Utf8},
+    )
     if df.shape[0] == 0:
         logging.info("No peaks to merge")
         return 0