diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 01e869b86..037f5d4e8 100644 --- a/.github/workflows/main.yml +++ b/.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 diff --git a/CHANGELOG.md b/CHANGELOG.md index 5769cd383..f87e1cb6b 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/config/README.md b/config/README.md index 2dda19989..01fe3a763 100644 --- a/config/README.md +++ b/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: ``` diff --git a/pixi.toml b/pixi.toml index cd59a75eb..dc98c8410 100644 --- a/pixi.toml +++ b/pixi.toml @@ -1,37 +1,16 @@ [workspace] authors = ["Mitchell Robert Vollger "] 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 = [ diff --git a/workflow/Snakefile b/workflow/Snakefile index dc1e50fdb..967e69371 100644 --- a/workflow/Snakefile +++ b/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 diff --git a/workflow/envs/env.yaml b/workflow/envs/env.yaml index 07984df4c..bc938e0cc 100644 --- a/workflow/envs/env.yaml +++ b/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 diff --git a/workflow/rules/coverages.smk b/workflow/rules/coverages.smk index 9a0aa1394..e7eb5ea89 100644 --- a/workflow/rules/coverages.smk +++ b/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 \ diff --git a/workflow/rules/fire-peaks.smk b/workflow/rules/fire-peaks.smk index 8093a0858..e1be2144a 100644 --- a/workflow/rules/fire-peaks.smk +++ b/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} """ diff --git a/workflow/scripts/fdr-table.py b/workflow/scripts/fdr-table.py index 3c185ead3..aa4f3505a 100644 --- a/workflow/scripts/fdr-table.py +++ b/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}" diff --git a/workflow/scripts/merge_fire_peaks.py b/workflow/scripts/merge_fire_peaks.py index 8867cec37..71d84ce29 100755 --- a/workflow/scripts/merge_fire_peaks.py +++ b/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