Simplify application of ft fire

fiberseq · Dec 11, 2024 · 83ba801 · 83ba801
1 parent c3d8e1f
commit 83ba801
Show file tree

Hide file tree

Showing 6 changed files with 22 additions and 42 deletions.
diff --git a/workflow/Snakefile b/workflow/Snakefile
@@ -130,7 +130,7 @@ rule all:
             v=VERSION,
         ),
         # model results
-        expand(rules.merged_fire_bam.output.cram, sm=MANIFEST.index, v=VERSION),
+        expand(rules.fire.output.cram, sm=MANIFEST.index, v=VERSION),
         expand(rules.fire_sites.output, sm=MANIFEST.index, v=VERSION),
         # Stats and Tables
         expand(rules.fires_in_peaks.output.txt, sm=MANIFEST.index, v=VERSION),

diff --git a/workflow/rules/apply-model.smk b/workflow/rules/apply-model.smk
@@ -4,64 +4,44 @@
 rule fire:
     input:
         bam=ancient(get_input_bam),
+        ref=ancient(REF),
     output:
-        bam=temp("temp/{sm}/fire/{chrom}.fire.bam"),
-    threads: 4
+        cram="results/{sm}/{sm}-fire-{v}-filtered.cram",
+        crai="results/{sm}/{sm}-fire-{v}-filtered.cram.crai",
+    threads: 32
     resources:
-        mem_mb=8 * 1024,
+        mem_mb=32 * 1024,
+        runtime=600,
     params:
         min_msp=config.get("min_msp", 10),
         min_ave_msp_size=config.get("min_ave_msp_size", 10),
         use_ont=USE_ONT,
         flag=FILTER_FLAG,
+    benchmark:
+        "results/{sm}/additional-outputs-{v}/benchmarks/{sm}-fire-bam.txt"
     conda:
         DEFAULT_ENV
     shell:
         """
-        samtools view -F {params.flag} -u -@ {threads} {input.bam} {wildcards.chrom} \
-            | {FT_EXE} fire -t {threads} \
+        samtools view -@ {threads} -u -F {params.flag} {input.bam} \
+            | {FT_EXE} fire -F {params.flag} -t {threads} \
                 {params.use_ont} \
                 --min-msp {params.min_msp} \
                 --min-ave-msp-size {params.min_ave_msp_size} \
                 --skip-no-m6a \
-                - {output.bam}
-        """
-
-
-rule merged_fire_bam:
-    input:
-        ref=ancient(REF),
-        fai=ancient(FAI),
-        bams=expand(rules.fire.output.bam, chrom=get_chroms(), allow_missing=True),
-    output:
-        cram="results/{sm}/{sm}-fire-{v}-filtered.cram",
-        crai="results/{sm}/{sm}-fire-{v}-filtered.cram.crai",
-    threads: 16
-    resources:
-        mem_mb=16 * 1024,
-        runtime=300,
-    conda:
-        DEFAULT_ENV
-    benchmark:
-        "results/{sm}/additional-outputs-{v}/benchmarks/{sm}-merged-fire-bam.txt"
-    priority: 100
-    shell:
-        """
-        samtools merge -@ {threads} -u {input.bams} -o - \
+                - - \
             | samtools view -C -@ {threads} -T {input.ref} \
                 --output-fmt-option embed_ref=1 \
             | samtools view -C -@ {threads} -T {input.ref} \
                 --output-fmt-option embed_ref=1 \
                 --input-fmt-option required_fields=0x1bff \
                 --write-index -o {output.cram}
-        # the second samtools view of CRAM file is needed to drop the quality scores
-        # this halves the size of the CRAM file
         """
 
 
 rule fire_sites_chrom:
     input:
-        cram=rules.merged_fire_bam.output.cram,
+        cram=rules.fire.output.cram,
     output:
         bed=temp("temp/{sm}/chrom/{v}-{chrom}.sorted.bed.gz"),
     threads: 4

diff --git a/workflow/rules/coverages.smk b/workflow/rules/coverages.smk
@@ -5,8 +5,8 @@ rule genome_bedgraph:
     input:
         ref=ancient(REF),
         fai=ancient(FAI),
-        cram=rules.merged_fire_bam.output.cram,
-        crai=rules.merged_fire_bam.output.crai,
+        cram=rules.fire.output.cram,
+        crai=rules.fire.output.crai,
     output:
         bg=temp("temp/{sm}/coverage/{sm}-{v}.bed.gz"),
         tbi=temp("temp/{sm}/coverage/{sm}-{v}.bed.gz.tbi"),
@@ -53,8 +53,8 @@ rule coverage:
 #
 rule fiber_locations_chromosome:
     input:
-        cram=rules.merged_fire_bam.output.cram,
-        crai=rules.merged_fire_bam.output.crai,
+        cram=rules.fire.output.cram,
+        crai=rules.fire.output.crai,
     output:
         bed=temp("temp/{sm}/coverage/{v}-{chrom}.fiber-locations.bed.gz"),
     threads: 4

diff --git a/workflow/rules/decorated-reads.smk b/workflow/rules/decorated-reads.smk
@@ -6,8 +6,8 @@ ITEMS_PER_SLOT = 1024 * 8
 
 rule decorate_fibers_chromosome:
     input:
-        cram=rules.merged_fire_bam.output.cram,
-        crai=rules.merged_fire_bam.output.crai,
+        cram=rules.fire.output.cram,
+        crai=rules.fire.output.crai,
     output:
         bed=temp("temp/{sm}/decorate/{v}-{chrom}.bed.gz"),
         decorated=temp("temp/{sm}/decorate/{v}-{chrom}.dec.bed.gz"),

diff --git a/workflow/rules/fire-peaks.smk b/workflow/rules/fire-peaks.smk
@@ -25,7 +25,7 @@ rule filtered_and_shuffled_fiber_locations_chromosome:
 
 rule shuffled_pileup_chromosome:
     input:
-        cram=rules.merged_fire_bam.output.cram,
+        cram=rules.fire.output.cram,
         shuffled=rules.filtered_and_shuffled_fiber_locations_chromosome.output.shuffled,
     output:
         bed=temp("temp/{sm}/shuffle/{v}-{chrom}.pileup.bed.gz"),
@@ -93,7 +93,7 @@ rule fdr_table:
 # coverage_H2  fire_coverage_H2  score_H2   nuc_coverage_H2 msp_coverage_H2
 rule pileup_chromosome:
     input:
-        bam=rules.merged_fire_bam.output.cram,
+        bam=rules.fire.output.cram,
     output:
         bed=temp("temp/{sm}/{v}-{chrom}.pileup.bed.gz"),
     threads: 4

diff --git a/workflow/rules/stats.smk b/workflow/rules/stats.smk
@@ -52,7 +52,7 @@ rule fires_in_peaks:
 
 rule ft_qc:
     input:
-        cram=rules.merged_fire_bam.output.cram,
+        cram=rules.fire.output.cram,
     output:
         tbl="results/{sm}/{sm}-fire-{v}-qc.tbl.gz",
     conda: