MissenseHMM/Snakefile at main · ernstlab/MissenseHMM · 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
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
import os

configfile: "config.yaml"

Z = config["Z"]
CHROMS = config["CHROMS"]
IN_LOC = config["IN_LOC"]
SCRIPT_LOC = config["SCRIPT_LOC"]
OUT_LOC = config["OUT_LOC"]
VERSION = config["VERSION"]
ENV = config["ENV"]
NUM_STATES = config["NUM_STATES"]
CHROMHMM_LOC = config["CHROMHMM_LOC"]


os.makedirs(OUT_LOC, exist_ok=True)
# os.makedirs(OUT_LOC + "/score_logs", exist_ok=True)


rule all:
    input:
        expand("{output}/dbnsfp_chr{chrm}.txt", chrm=CHROMS, output=OUT_LOC),
        expand("{output}/score_logs/dbnsfp_{z}.log", z=range(0, Z), output=OUT_LOC),
        OUT_LOC + "/thresholds.txt",
        expand("{output}/samples/log_chr{chrm}.log", chrm=CHROMS, output=OUT_LOC),
        expand("{output}/model.log", output=OUT_LOC),
        OUT_LOC + "/max_states/all.bed"


rule create_datafile:
    input:
        "%s/dbNSFP%s/dbNSFP%s_variant.chr{chrm}.gz" % (IN_LOC, VERSION, VERSION)
    output:
        "%s/dbnsfp_chr{chrm}.txt" % OUT_LOC
    params:
        script_loc = SCRIPT_LOC,
        version = VERSION,
        output_loc = OUT_LOC
    threads: 1
    shell:
        """
        {params.script_loc}/to_datafile.py -c {wildcards.chrm} -v {params.version} -g -l --output {params.output_loc}
        """


rule threshold:
    input:
        expand("{output}/dbnsfp_chr{chrm}.txt", chrm=CHROMS, output=OUT_LOC)
    output:
        "%s/score_logs/dbnsfp_{z}.log" % OUT_LOC
    threads: 1
    params:
        script_loc = SCRIPT_LOC,
        output_loc = OUT_LOC
    shell:
        """
        mkdir -p {params.output_loc}/score_logs
        {params.script_loc}/threshold.py {wildcards.z} {params.output_loc} {params.output_loc}/score_logs
        """


rule generate_threshold_file:
    input:
        expand("{output}/score_logs/dbnsfp_{z}.log", z=range(0, Z), output=OUT_LOC)
    output:
        "%s/thresholds.txt" % OUT_LOC
    threads: 1
    params:
        script_loc = SCRIPT_LOC,
        output_loc = OUT_LOC
    shell:
        """
        {params.script_loc}/generate_thresholds.py {params.output_loc}/score_logs {params.output_loc}
        """


rule generate_binary:
    input:
        "%s/thresholds.txt" % OUT_LOC,
        expand("{output}/dbnsfp_chr{chrm}.txt", chrm=CHROMS, output=OUT_LOC)
    output:
        "%s/samples/log_chr{chrm}.log" % OUT_LOC
    threads: 1
    params:
        script_loc = SCRIPT_LOC,
        output_loc = OUT_LOC
    shell:
        """
        mkdir -p {params.output_loc}/samples
        {params.script_loc}/to_binary.py -c {wildcards.chrm} -i {params.output_loc} -o {params.output_loc}/samples -t {params.output_loc}/thresholds.txt
        echo DONE >  {params.output_loc}/samples/log_chr{wildcards.chrm}.log
        """


rule train_hmm:
    input:
        expand("{output}/samples/log_chr{chrm}.log", chrm=CHROMS, output=OUT_LOC)
    output:
        "%s/model.log" % OUT_LOC
    threads: 4
    params:
        chromhmm_loc = CHROMHMM_LOC,
        output_loc = OUT_LOC,
        env = ENV,
        num_states = NUM_STATES
    shell:
        """
        if [ "{params.env}" = "cluster" ]; then
            . /u/local/Modules/default/init/modules.sh; module load java
        fi

        java -jar {params.chromhmm_loc}/ChromHMM.jar LearnModel -nobed -nobrowser -noenrich -noimage -pseudo -b 1 -n 2000 -d "-1" -lowmem -p 4 {params.output_loc}/samples {params.output_loc}/model_{params.num_states} {params.num_states} hg38

        echo DONE > {params.output_loc}/model.log
        """


rule apply_model:
    input:
        "%s/model.log" % OUT_LOC,
        "%s/thresholds.txt" % OUT_LOC,
        "%s/dbnsfp_chr{chrm}.txt" % OUT_LOC
    output:
        "%s/max_states/max_states_chr{chrm}.bed" % OUT_LOC
    threads: 1
    params:
        script_loc = SCRIPT_LOC,
        output_loc = OUT_LOC,
        num_states = NUM_STATES
    shell:
        """
        mkdir -p {params.output_loc}/max_states
        {params.script_loc}/max_pos_state.py -c {wildcards.chrm} -m {params.output_loc}/model_{params.num_states}/model_{params.num_states}.txt -o {params.output_loc}/max_states -t {params.output_loc}/thresholds.txt
        """


rule combine_states:
    input:
        expand("{output}/max_states/max_states_chr{chrm}.bed", chrm=CHROMS, output=OUT_LOC)
    output:
        "%s/max_states/all.bed" % OUT_LOC
    threads: 1
    params:
        output_loc = OUT_LOC
    shell:
        """
        cat {params.output_loc}/max_states/max_states_chr*.bed > {params.output_loc}/max_states/all.bed
        """