diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 037f5d4e8..f5ba8c199 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -4,7 +4,6 @@ on: push: branches: [main] pull_request: - branches: [main] jobs: # pixi test diff --git a/workflow/rules/apply-model.smk b/workflow/rules/apply-model.smk index 8a2a6bf7d..af721bfd2 100644 --- a/workflow/rules/apply-model.smk +++ b/workflow/rules/apply-model.smk @@ -8,6 +8,10 @@ rule fire: output: cram="results/{sm}/{sm}-fire-{v}-filtered.cram", crai="results/{sm}/{sm}-fire-{v}-filtered.cram.crai", + benchmark: + "results/{sm}/additional-outputs-{v}/benchmarks/{sm}-fire-bam.txt" + conda: + DEFAULT_ENV threads: 32 resources: mem_mb=32 * 1024, @@ -17,10 +21,6 @@ rule fire: 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 -@ {threads} -u -F {params.flag} {input.bam} \ @@ -37,7 +37,7 @@ rule fire: --input-fmt-option required_fields=0x1bff \ --write-index -o {output.cram} - # check if the cram file has zero reads + # check if the cram file has zero reads reads_in_header=$(samtools view {output.cram} | head | wc -l || true) if [ $reads_in_header -eq 0 ]; then printf "\nNo reads passed filters exiting...\n\nPlease review https://fiberseq.github.io/quick-start.html to make sure the input BAM has been correctly processed.\n\n" @@ -51,9 +51,9 @@ rule fire_sites_chrom: cram=rules.fire.output.cram, output: bed=temp("temp/{sm}/chrom/{v}-{chrom}.sorted.bed.gz"), - threads: 4 conda: DEFAULT_ENV + threads: 4 resources: mem_mb=16 * 1024, params: @@ -62,13 +62,13 @@ rule fire_sites_chrom: """ samtools view -@ {threads} -u {input.cram} {wildcards.chrom} \ | {FT_EXE} fire -t {threads} --extract - \ - | LC_ALL=C sort --parallel={threads} \ - -k1,1 -k2,2n -k3,3n -k4,4 \ - | bioawk -tc hdr '$10<={params.min_fdr}' \ - | (grep '\\S' || true) \ - | (grep -v '^#' || true) \ - | bgzip -@ {threads} \ - > {output.bed} + | LC_ALL=C sort --parallel={threads} \ + -k1,1 -k2,2n -k3,3n -k4,4 \ + | bioawk -tc hdr '$10<={params.min_fdr}' \ + | (grep '\\S' || true) \ + | (grep -v '^#' || true) \ + | bgzip -@ {threads} \ + >{output.bed} """ @@ -79,12 +79,12 @@ rule fire_sites: ), output: bed="results/{sm}/additional-outputs-{v}/fire-peaks/{sm}-{v}-fire-elements.bed.gz", - threads: 1 conda: DEFAULT_ENV + threads: 1 shell: """ - cat {input.beds} > {output.bed} + cat {input.beds} >{output.bed} """ @@ -93,9 +93,9 @@ rule fire_sites_index: bed=rules.fire_sites.output.bed, output: tbi=rules.fire_sites.output.bed + ".tbi", - threads: 1 conda: DEFAULT_ENV + threads: 1 shell: """ tabix -p bed {input.bed} diff --git a/workflow/rules/coverages.smk b/workflow/rules/coverages.smk index e7eb5ea89..ec3b6f955 100644 --- a/workflow/rules/coverages.smk +++ b/workflow/rules/coverages.smk @@ -10,18 +10,18 @@ rule genome_bedgraph: output: bg=temp("temp/{sm}/coverage/{sm}-{v}.bed.gz"), tbi=temp("temp/{sm}/coverage/{sm}-{v}.bed.gz.tbi"), - threads: 16 shadow: "minimal" conda: DEFAULT_ENV + threads: 16 shell: - """ + """ 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 \ + | LC_ALL=C sort --parallel={threads} -k1,1 -k2,2n -k3,3n -k4,4 \ | bgzip -@ {threads} \ - > {output.bg} + >{output.bg} tabix -f -p bed {output.bg} """ @@ -33,13 +33,13 @@ rule coverage: cov="results/{sm}/additional-outputs-{v}/coverage/{sm}-{v}-median-coverage.txt", minimum="results/{sm}/additional-outputs-{v}/coverage/{sm}-{v}-minimum-coverage.txt", maximum="results/{sm}/additional-outputs-{v}/coverage/{sm}-{v}-maximum-coverage.txt", + benchmark: + "results/{sm}/additional-outputs-{v}/benchmarks/coverage/{sm}.txt" conda: "../envs/python.yaml" threads: 1 resources: mem_mb=64 * 1024, - benchmark: - "results/{sm}/additional-outputs-{v}/benchmarks/coverage/{sm}.txt" params: coverage_within_n_sd=COVERAGE_WITHIN_N_SD, min_coverage=MIN_COVERAGE, @@ -57,18 +57,18 @@ rule fiber_locations_chromosome: crai=rules.fire.output.crai, output: bed=temp("temp/{sm}/coverage/{v}-{chrom}.fiber-locations.bed.gz"), - threads: 4 conda: DEFAULT_ENV + threads: 4 shell: """ # get fiber locations (samtools view -@ {threads} -u {input.cram} {wildcards.chrom} \ | {FT_EXE} extract -t {threads} -s --all - \ - | hck -F '#ct' -F st -F en -F fiber -F strand -F HP ) \ + | hck -F '#ct' -F st -F en -F fiber -F strand -F HP) \ | (grep -v "^#" || true) \ | bgzip -@ {threads} \ - > {output.bed} + >{output.bed} """ @@ -91,24 +91,24 @@ rule fiber_locations: filtered_tbi=temp( "temp/{sm}/coverage/filtered-for-coverage/{v}-fiber-locations.bed.gz.tbi" ), - threads: 4 conda: DEFAULT_ENV + threads: 4 params: max_frac_overlap=0.2, shell: """ - cat {input.fibers} > {output.bed} + cat {input.fibers} >{output.bed} tabix -f -p bed {output.bed} - + # get filtered fiber locations MIN=$(cat {input.minimum}) MAX=$(cat {input.maximum}) bedtools intersect -header -sorted -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} \ - > {output.filtered} + | bgzip -@ {threads} \ + >{output.filtered} tabix -f -p bed {output.filtered} """ @@ -122,9 +122,9 @@ rule exclude_from_shuffle: fai=ancient(FAI), output: bed="results/{sm}/additional-outputs-{v}/coverage/exclude-from-shuffles.bed.gz", - threads: 4 conda: DEFAULT_ENV + threads: 4 params: exclude=lambda wc: " ".join(EXCLUDES) if EXCLUDES else "", shell: @@ -139,7 +139,7 @@ rule exclude_from_shuffle: | bedtools sort \ | bedtools merge \ | bgzip -@ {threads} \ - > {output.bed} + >{output.bed} """ @@ -154,12 +154,12 @@ rule unreliable_coverage_regions: bed_tbi="results/{sm}/additional-outputs-{v}/coverage/unreliable-coverage-regions.bed.gz.tbi", tmp=temp("temp/{sm}/additional-outputs-{v}/unreliable-coverage-regions.bed"), bb="results/{sm}/trackHub-{v}/bb/unreliable-coverage-regions.bb", + conda: + DEFAULT_ENV threads: 4 params: min_len=MIN_UNRELIABLE_COVERAGE_LEN, bed3_as=workflow.source_path("../templates/bed3.as"), - conda: - DEFAULT_ENV shell: """ MIN=$(cat {input.minimum}) @@ -170,15 +170,15 @@ rule unreliable_coverage_regions: | bedtools merge -i - \ | awk '$3-$2 >= {params.min_len}' \ | bgzip -@ {threads} \ - > {output.bed} + >{output.bed} # bigbed # for some reason bigtools gives a too many files open error when reading from stdin - bedtools merge -i {output.bed} > {output.tmp} + bedtools merge -i {output.bed} >{output.tmp} bigtools bedtobigbed \ - -s start -a {params.bed3_as} \ - {output.tmp} {input.fai} {output.bb} + -s start -a {params.bed3_as} \ + {output.tmp} {input.fai} {output.bb} - # index + # index tabix -f -p bed {output.bed} """ diff --git a/workflow/rules/decorated-reads.smk b/workflow/rules/decorated-reads.smk index b39eb5f6f..6a92ee80d 100644 --- a/workflow/rules/decorated-reads.smk +++ b/workflow/rules/decorated-reads.smk @@ -11,20 +11,20 @@ rule decorate_fibers_chromosome: output: bed=temp("temp/{sm}/decorate/{v}-{chrom}.bed.gz"), decorated=temp("temp/{sm}/decorate/{v}-{chrom}.dec.bed.gz"), - threads: 4 - resources: - mem_mb=get_mem_mb, # the following steps can take a while so this helps the pipeline start this earlier. priority: 10 conda: DEFAULT_ENV + threads: 4 + resources: + mem_mb=get_mem_mb, shell: """ samtools view -@ {threads} -u {input.cram} {wildcards.chrom} \ | {FT_EXE} track-decorators -t {threads} --bed12 {output.bed} \ | sort -k1,1 -k2,2n -k3,3n -k4,4 \ | bgzip -@ {threads} \ - > {output.decorated} + >{output.decorated} """ @@ -41,12 +41,12 @@ rule decorate_fibers_1: bb="results/{sm}/trackHub-{v}/bb/fire-fibers.bb", benchmark: "results/{sm}/additional-outputs-{v}/benchmarks/decorate_fibers_1/{sm}.txt" - threads: 8 - resources: - runtime=240, priority: 10 conda: DEFAULT_ENV + threads: 8 + resources: + runtime=240, params: bed_as=workflow.source_path("../templates/bed12_filter.as"), nzooms=NZOOMS, @@ -79,12 +79,12 @@ rule decorate_fibers_2: #bed=temp("temp/{sm}/trackHub-{v}/bb/fire-fiber-decorators.bed.gz"), benchmark: "results/{sm}/additional-outputs-{v}/benchmarks/decorate_fibers_2/{sm}.txt" - threads: 8 - resources: - runtime=60 * 16, priority: 10 conda: DEFAULT_ENV + threads: 8 + resources: + runtime=60 * 16, params: dec_as=workflow.source_path("../templates/decoration.as"), nzooms=NZOOMS, diff --git a/workflow/rules/fire-peaks.smk b/workflow/rules/fire-peaks.smk index e1be2144a..fbccff0e5 100644 --- a/workflow/rules/fire-peaks.smk +++ b/workflow/rules/fire-peaks.smk @@ -6,9 +6,9 @@ rule filtered_and_shuffled_fiber_locations_chromosome: fai=ancient(FAI), output: shuffled=temp("temp/{sm}/shuffle/{v}-{chrom}.fiber-locations-shuffled.bed.gz"), - threads: 4 conda: DEFAULT_ENV + threads: 4 shell: """ tabix {input.filtered} {wildcards.chrom} \ @@ -17,9 +17,9 @@ rule filtered_and_shuffled_fiber_locations_chromosome: -excl {input.exclude} \ -i - \ -g {input.fai} \ - | sort -k1,1 -k2,2n -k3,3n -k4,4 \ + | sort -k1,1 -k2,2n -k3,3n -k4,4 \ | bgzip -@ {threads} \ - > {output.shuffled} + >{output.shuffled} """ @@ -29,16 +29,16 @@ rule shuffled_pileup_chromosome: shuffled=rules.filtered_and_shuffled_fiber_locations_chromosome.output.shuffled, output: bed=temp("temp/{sm}/shuffle/{v}-{chrom}.pileup.bed.gz"), - threads: 4 conda: DEFAULT_ENV + threads: 4 shell: """ {FT_EXE} pileup {input.cram} --rgn {wildcards.chrom} -t {threads} \ --fiber-coverage --shuffle {input.shuffled} \ --no-msp --no-nuc \ | bgzip -@ {threads} \ - > {output.bed} + >{output.bed} """ @@ -52,12 +52,12 @@ rule shuffled_pileup: output: bed=temp("temp/{sm}/shuffle/{v}-shuffled-pileup.bed.gz"), tbi=temp("temp/{sm}/shuffle/{v}-shuffled-pileup.bed.gz.tbi"), - threads: 4 conda: DEFAULT_ENV + threads: 4 shell: """ - cat {input.beds} > {output.bed} + cat {input.beds} >{output.bed} tabix -p bed {output.bed} """ @@ -74,11 +74,11 @@ rule fdr_table: tbl="results/{sm}/additional-outputs-{v}/fire-peaks/{sm}-{v}-fire-score-to-fdr.tbl", conda: "../envs/python.yaml" - params: - script=workflow.source_path("../scripts/fdr-table.py"), threads: 16 resources: mem_mb=get_mem_mb_xl, + params: + script=workflow.source_path("../scripts/fdr-table.py"), shell: """ export POLARS_MAX_THREADS={threads} @@ -98,16 +98,16 @@ rule pileup_chromosome: bam=rules.fire.output.cram, output: bed=temp("temp/{sm}/{v}-{chrom}.pileup.bed.gz"), - threads: 4 conda: DEFAULT_ENV + threads: 4 shell: """ {FT_EXE} pileup -t {threads} \ --haps --fiber-coverage \ {input.bam} --rgn {wildcards.chrom} \ | bgzip -@ {threads} \ - > {output.bed} + >{output.bed} """ @@ -117,13 +117,13 @@ rule fdr_track_chromosome: fdr_tbl=rules.fdr_table.output.tbl, output: bed=temp("temp/{sm}/fire-peaks/{v}-{chrom}-FDR.track.bed"), - threads: 4 conda: "../envs/python.yaml" - params: - script=workflow.source_path("../scripts/fdr-table.py"), + threads: 4 resources: mem_mb=get_mem_mb_xl, + params: + script=workflow.source_path("../scripts/fdr-table.py"), shell: """ export POLARS_MAX_THREADS={threads} @@ -144,24 +144,24 @@ rule pileup: fofn=temp("temp/{sm}/fire/fire-{v}-pileup.fofn"), bed="results/{sm}/{sm}-fire-{v}-pileup.bed.gz", tbi="results/{sm}/{sm}-fire-{v}-pileup.bed.gz.tbi", - threads: 8 conda: DEFAULT_ENV + threads: 8 shell: """ printf '\nMaking FOFN\n' - echo {input.beds} > {output.fofn} - + echo {input.beds} >{output.fofn} + printf '\nMake header\n' - ((cat $(cat {output.fofn}) | grep "^#" | head -n 1) || true) \ + ( (cat $(cat {output.fofn}) | grep "^#" | head -n 1) || true) \ | bgzip -@ {threads} \ - > {output.bed} + >{output.bed} printf '\nConcatenating\n' cat $(cat {output.fofn}) \ | grep -v "^#" \ | bgzip -@ {threads} \ - >> {output.bed} + >>{output.bed} printf '\nIndexing\n' tabix -f -p bed {output.bed} @@ -176,9 +176,9 @@ rule helper_fdr_peaks_by_fire_elements: fire_tbi=rules.fire_sites_index.output.tbi, output: bed=temp("temp/{sm}/fire-peaks/{v}-{chrom}-fire-peaks.bed.gz"), - threads: 2 conda: DEFAULT_ENV + threads: 2 params: max_peak_fdr=MAX_PEAK_FDR, min_per_acc_peak=MIN_PER_ACC_PEAK, @@ -186,36 +186,37 @@ rule helper_fdr_peaks_by_fire_elements: """ HEADER=$(bgzip -cd {input.bed} | head -n 1 || true) NC=$(echo $HEADER | awk '{{print NF}}' || true) - FIRE_CT=$((NC+1)) - FIRE_ST=$((NC+2)) - FIRE_EN=$((NC+3)) - FIRE_SIZE=$((NC+4)) - FIRE_ID=$((NC+5)) + FIRE_CT=$((NC + 1)) + FIRE_ST=$((NC + 2)) + FIRE_EN=$((NC + 3)) + FIRE_SIZE=$((NC + 4)) + FIRE_ID=$((NC + 5)) OUT_HEADER=$(printf "$HEADER\\tpeak_chrom\\tpeak_start\\tpeak_end\\tFIRE_IDs\\tFIRE_size_mean\\tFIRE_size_ssd\\tFIRE_start_ssd\\tFIRE_end_ssd") echo $OUT_HEADER - - ( \ - printf "$OUT_HEADER\\n"; \ + + ( + printf "$OUT_HEADER\\n" tabix -h {input.bed} {wildcards.chrom} \ | bioawk -tc hdr '(NR==1)||($is_local_max=="true")' \ | csvtk filter -tT -C '$' -f "FDR<={params.max_peak_fdr}" \ | csvtk filter -tT -C '$' -f "fire_coverage>1" \ | bioawk -tc hdr '(NR==1)||(NF>0 && $fire_coverage/$coverage>={params.min_per_acc_peak})' \ | bedtools intersect -wa -wb -sorted -a - \ - -b <(tabix {input.fire} {wildcards.chrom} \ + -b <( + tabix {input.fire} {wildcards.chrom} \ | cut -f 1-3 \ - | awk -v OFMT="%f" '{{print $0"\t"$3-$2"\t"NR}}' \ - ) \ + | awk -v OFMT="%f" '{{print $0"\t"$3-$2"\t"NR}}' + ) \ | bedtools groupby -g 1-$NC \ -o first,median,median,collapse,mean,sstdev,sstdev,sstdev \ - -c $FIRE_CT,$FIRE_ST,$FIRE_EN,$FIRE_ID,$FIRE_SIZE,$FIRE_SIZE,$FIRE_ST,$FIRE_EN \ + -c $FIRE_CT,$FIRE_ST,$FIRE_EN,$FIRE_ID,$FIRE_SIZE,$FIRE_SIZE,$FIRE_ST,$FIRE_EN ) \ | hck -f 1,$FIRE_ST,$FIRE_EN,2-$NC,$FIRE_SIZE- \ | csvtk round -tT -C '$' -n 0 -f 2,3 \ | bedtools sort -header -i - \ | bgzip -@ {threads} \ - > {output.bed} + >{output.bed} """ @@ -226,9 +227,9 @@ rule fdr_peaks_by_fire_elements_chromosome: maximum=rules.coverage.output.maximum, output: bed=temp("temp/{sm}/fire-peaks/{v}-grouped-{chrom}-fire-peaks.bed.gz"), - threads: 4 conda: "../envs/python.yaml" + threads: 4 params: script=workflow.source_path("../scripts/merge_fire_peaks.py"), min_frac_accessible=MIN_FRAC_ACCESSIBLE, @@ -241,7 +242,7 @@ rule fdr_peaks_by_fire_elements_chromosome: --min-cov $(cat {input.minimum}) \ --min-frac-accessible {params.min_frac_accessible} \ | bgzip -@ {threads} \ - > {output.bed} + >{output.bed} """ @@ -256,22 +257,22 @@ rule fire_peaks: fofn=temp("temp/{sm}/fire-peaks/{sm}-fire-{v}-peaks.fofn"), bed="results/{sm}/{sm}-fire-{v}-peaks.bed.gz", tbi="results/{sm}/{sm}-fire-{v}-peaks.bed.gz.tbi", - threads: 4 conda: DEFAULT_ENV + threads: 4 shell: """ printf "\nMaking FOFN\n" - echo {input.beds} > {output.fofn} + echo {input.beds} >{output.fofn} - printf "\nMaking header\n" - ((cat $(cat {output.fofn}) | bgzip -cd | grep "^#" | head -n 1) || true) \ - | bgzip -@ {threads} > {output.bed} + printf "\nMaking header\n" + ( (cat $(cat {output.fofn}) | bgzip -cd | grep "^#" | head -n 1) || true) \ + | bgzip -@ {threads} >{output.bed} printf "\nConcatenating\n" cat $(cat {output.fofn}) | bgzip -cd -@ {threads} | grep -v "^#" \ - | bgzip -@ {threads} >> {output.bed} - + | bgzip -@ {threads} >>{output.bed} + printf "\nIndexing\n" tabix -f -p bed {output.bed} """ @@ -296,22 +297,22 @@ rule wide_fire_peaks: bed3_as=workflow.source_path("../templates/bed3.as"), shell: """ - ( \ - bgzip -cd {input.bed}; \ + ( + bgzip -cd {input.bed} bioawk -tc hdr 'NR==1 || $FDR<={params.max_peak_fdr}' {input.track} \ - | bioawk -tc hdr 'NR==1 || (NF>0 && $coverage>0 && $fire_coverage/$coverage>={params.min_frac_acc})' \ + | bioawk -tc hdr 'NR==1 || (NF>0 && $coverage>0 && $fire_coverage/$coverage>={params.min_frac_acc})' ) \ | cut -f 1-3 \ | bedtools sort \ | bedtools merge -d {params.nuc_size} \ | bgzip -@ {threads} \ - > {output.bed} - + >{output.bed} + bgzip -cd -@ 16 {output.bed} \ | bigtools bedtobigbed \ -s start -a {params.bed3_as} \ - - {input.fai} {output.bb} - + - {input.fai} {output.bb} + tabix -p bed {output.bed} """ @@ -325,9 +326,9 @@ rule one_percent_fire_peaks: tbi="results/{sm}/additional-outputs-{v}/fire-peaks/one-percent-FDR/{sm}-fire-{v}-01-fire-peaks.bed.gz.tbi", wide="results/{sm}/additional-outputs-{v}/fire-peaks/one-percent-FDR/{sm}-fire-{v}-01-fire-wide-peaks.bed.gz", wtbi="results/{sm}/additional-outputs-{v}/fire-peaks/one-percent-FDR/{sm}-fire-{v}-01-fire-wide-peaks.bed.gz.tbi", - threads: 4 conda: DEFAULT_ENV + threads: 4 params: nuc_size=config.get("nucleosome_size", 147), shell: @@ -335,18 +336,18 @@ rule one_percent_fire_peaks: bgzip -cd {input.bed} \ | csvtk filter -tT -C '$' -f "FDR<=0.01" \ | bgzip -@ {threads} \ - > {output.bed} + >{output.bed} tabix -f -p bed {output.bed} - ( \ - bgzip -cd {output.bed}; \ - bioawk -tc hdr '$FDR<=0.01' {input.track} \ + ( + bgzip -cd {output.bed} + bioawk -tc hdr '$FDR<=0.01' {input.track} ) \ | cut -f 1-3 \ | bedtools sort \ | bedtools merge -d {params.nuc_size} \ | bgzip -@ {threads} \ - > {output.wide} + >{output.wide} tabix -f -p bed {output.wide} """ @@ -359,8 +360,8 @@ rule peaks_vs_percent: "results/{sm}/additional-outputs-{v}/figures/{sm}-fire-{v}-peaks-vs-percent.pdf", category="Peak calls", ), - threads: 4 conda: "../envs/R.yaml" + threads: 4 script: "../scripts/peaks-vs-percent.R" diff --git a/workflow/rules/levio.smk b/workflow/rules/levio.smk index e6d4b1b16..409303abc 100644 --- a/workflow/rules/levio.smk +++ b/workflow/rules/levio.smk @@ -41,12 +41,12 @@ rule leviosam2: lifted=temp("temp/{sm}/leviosam2/{sm}-{chrom}-committed.bam"), deferred=temp("temp/{sm}/leviosam2/{sm}-{chrom}-deferred.bam"), unliftable=temp("temp/{sm}/leviosam2/{sm}-{chrom}-unliftable.bam"), + conda: + DEFAULT_ENV threads: MAX_THREADS resources: mem_mb=MAX_THREADS * 4 * 1024, runtime=16 * 60, - conda: - DEFAULT_ENV params: # maximum number of CIGAR opts to change, also the max gap size that can be spanned G=config.get("levio_G", 100_000), @@ -83,12 +83,12 @@ rule leviosam2_sorted: ref=REF, output: bam=temp("temp/{sm}/leviosam2/{sm}-{chrom}-sorted.bam"), + conda: + DEFAULT_ENV threads: SORT_THREADS resources: mem_mb=SORT_THREADS * 4 * 1024, runtime=16 * 60, - conda: - DEFAULT_ENV shell: """ samtools sort {input.lifted} \ diff --git a/workflow/rules/stats.smk b/workflow/rules/stats.smk index 366eca87e..f462c8592 100644 --- a/workflow/rules/stats.smk +++ b/workflow/rules/stats.smk @@ -9,19 +9,21 @@ rule clustering_vs_null: tmp=temp("temp/{sm}/tmp.pre.calls.bed"), null=temp("temp/{sm}/null.calls.bed"), bed="results/{sm}/clustering-vs-null.bed.gz", - threads: 4 conda: DEFAULT_ENV + threads: 4 shell: """ - bgzip -cd -@{threads} {input.bed} | cut -f 1-3 > {output.tmp} - bedtools shuffle -chrom -i {output.tmp} -g {input.fai} > {output.null} + bgzip -cd -@{threads} {input.bed} | cut -f 1-3 >{output.tmp} + bedtools shuffle -chrom -i {output.tmp} -g {input.fai} >{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.fai} | sed 's/$/\\tReal/g' + bedtools genomecov -bg -i {output.null} -g {input.fai} | sed 's/$/\\tNull/g' + ) \ | bedtools sort \ | bgzip -@ {threads} \ - > {output.bed} + >{output.bed} """ @@ -33,20 +35,20 @@ rule fires_in_peaks: 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", - threads: 4 conda: DEFAULT_ENV + threads: 4 params: 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 -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 "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} + echo "# of FIREs within peaks" >>{output.txt} + bedtools intersect -sorted -u -a {input.fire} -b {input.peaks} | wc -l >>{output.txt} """ @@ -78,8 +80,8 @@ rule hap_differences: ), bed="results/{sm}/{sm}-fire-{v}-hap-differences.bed.gz", bed9=temp("temp/{sm}/hap1-vs-hap2/FIRE-{v}.hap.differences.bed9"), - threads: 4 conda: "../envs/R.yaml" + threads: 4 script: "../scripts/hap-diffs.R" diff --git a/workflow/rules/track-hub.smk b/workflow/rules/track-hub.smk index 422fde256..1be0fcb55 100644 --- a/workflow/rules/track-hub.smk +++ b/workflow/rules/track-hub.smk @@ -5,9 +5,9 @@ rule percent_accessible: output: tmp=temp("temp/{sm}/{hp}/{v}-percent.accessible.bed"), bw="results/{sm}/trackHub-{v}/bw/{hp}.percent.accessible.bw", - threads: 4 conda: DEFAULT_ENV + threads: 4 resources: mem_mb=get_mem_mb, params: @@ -19,18 +19,17 @@ rule percent_accessible: bgzip -cd {input.bed} \ | bioawk -tc hdr '$coverage{params.suffix}>0' \ | bioawk -tc hdr \ - 'NR>1{{print $1,$2,$3,100*$fire_coverage{params.suffix}/$coverage{params.suffix}}}' \ - > {output.tmp} + 'NR>1{{print $1,$2,$3,100*$fire_coverage{params.suffix}/$coverage{params.suffix}}}' \ + >{output.tmp} # add fake if file is empty if [[ -s {output.tmp} ]]; then echo "File is not empty" else echo "File is empty" - printf "{params.chrom}\t0\t1\t0\\n" > {output.tmp} + printf "{params.chrom}\t0\t1\t0\\n" >{output.tmp} fi - bigtools bedgraphtobigwig \ --nzooms {params.nzooms} -s start \ {output.tmp} {input.fai} {output.bw} @@ -65,9 +64,9 @@ rule fdr_track_to_bw: fai=ancient(FAI), output: bw="results/{sm}/trackHub-{v}/bw/{col}.bw", - threads: 4 conda: DEFAULT_ENV + threads: 4 params: nzooms=NZOOMS, shell: @@ -86,9 +85,9 @@ rule fire_peaks_bb: fai=ancient(FAI), output: bb="results/{sm}/trackHub-{v}/bb/fire-peaks.bb", - threads: 4 conda: DEFAULT_ENV + threads: 4 params: bedfmt=workflow.source_path("../templates/fire_peak.as"), shell: @@ -109,19 +108,19 @@ rule hap_differences_track: fai=ancient(FAI), output: bb="results/{sm}/trackHub-{v}/bb/hap_differences.bb", + conda: + DEFAULT_ENV threads: 4 resources: mem_mb=get_mem_mb, - conda: - DEFAULT_ENV params: chrom=get_chroms()[0], bed9_as=workflow.source_path("../templates/bed9.as"), shell: """ - ( \ - printf "{params.chrom}\t0\t1\tfake\t100\t+\t0\t1\t230,230,230\\n"; \ - bedtools sort -i {input.bed9} \ + ( + printf "{params.chrom}\t0\t1\tfake\t100\t+\t0\t1\t230,230,230\\n" + bedtools sort -i {input.bed9} ) \ | bigtools bedtobigbed \ -s start -a {params.bed9_as} \ @@ -135,11 +134,11 @@ rule trackhub: output: hub="results/{sm}/trackHub-{v}/hub.txt", description="results/{sm}/trackHub-{v}/fire-description.html", - resources: - load=get_load, - threads: 4 conda: "../envs/python.yaml" + threads: 4 + resources: + load=get_load, params: ref=REF_NAME, script=workflow.source_path("../scripts/trackhub.py"), @@ -147,9 +146,9 @@ rule trackhub: shell: """ python {params.script} -v 2 \ - --trackhub-dir results/{wildcards.sm}/trackHub-{wildcards.v} \ - --reference {params.ref} \ - --sample {wildcards.sm} \ - --average-coverage $(cat {input.cov}) + --trackhub-dir results/{wildcards.sm}/trackHub-{wildcards.v} \ + --reference {params.ref} \ + --sample {wildcards.sm} \ + --average-coverage $(cat {input.cov}) cp {params.description} {output.description} """