FIRE/workflow/rules/apply-model.smk at 30018c0a61117e4cdefb4f8fe5bde277579a04ad · fiberseq/FIRE · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
#
# Applying the model
#
rule fire:
    input:
        bam=ancient(get_input_bam),
        ref=ancient(get_ref)
        #ref=ancient(REF),
    output:
        cram="results/{sm}/{sm}-fire-{v}-filtered.cram",
        crai="results/{sm}/{sm}-fire-{v}-filtered.cram.crai",
    threads: 32
    resources:
        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 -@ {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 \
                - - \
            | 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}

        # 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"
            exit 1
        fi
        """


rule fire_sites_chrom:
    input:
        cram=rules.fire.output.cram,
    output:
        bed=temp("temp/{sm}/chrom/{v}-{chrom}.sorted.bed.gz"),
    threads: 4
    conda:
        DEFAULT_ENV
    resources:
        mem_mb=16 * 1024,
    params:
        min_fdr=MIN_FIRE_FDR,
    shell:
        """
        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}
        """


rule fire_sites:
    input:
        beds=expand(
            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",
    threads: 1
    conda:
        DEFAULT_ENV
    shell:
        """
        cat {input.beds} > {output.bed}
        """


rule fire_sites_index:
    input:
        bed=rules.fire_sites.output.bed,
    output:
        tbi=rules.fire_sites.output.bed + ".tbi",
    threads: 1
    conda:
        DEFAULT_ENV
    shell:
        """
        tabix -p bed {input.bed}
        """