Generate-ERs/snakefile at main · sid-sethi/Generate-ERs · 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
import os
from os import path

from snakemake.utils import min_version
min_version("5.3")


configfile: "config.yml"
workdir: config["workdir"]

WORKDIR = config["workdir"]
BIGWIGDIR = config["bigwig"]
SNAKEDIR = path.dirname(workflow.snakefile)

sample = config["sample_name"]
gtf = config["gtf"]


rule all:
    input:
        "ERs/" + sample + "_ers_raw.txt",
        "Intergenic_ERs/" + sample + "_3prime_intergenic_ers.txt",
        "Intergenic_ERs/" + sample + "_5prime_intergenic_ers.txt"

#########################################################################


rule generate_ERs:
    """
    Call ERs using ODER
    """
    input:
        gtf_file = gtf

    output:
        opt_ers = "ERs/" + sample + "_optimal_ers.rda"

    params:
        bw_dir = BIGWIGDIR,
        resPath = "ERs",
        prefix = sample,
        script = SNAKEDIR + "/scripts/oder.R"

    conda: "envs/r_v4.yml"

    shell:
        """
        Rscript {params.script} {params.bw_dir} {input.gtf_file} {params.prefix} {params.resPath}
        """


rule generate_txdb:
    """
    Generate TxDb object from provided GTF
    """
    input:
        gtf_file = gtf,
        chrom_len = SNAKEDIR + "/data/chromInfo_hg38.txt"

    output:
        txdb = "TxDb/gtf_TxDb.Rsqlite"

    params:
        resPath = "TxDb",
        script = SNAKEDIR + "/scripts/txdb_from_gtf.R"

    conda: "envs/r_v3.yml"

    shell:
        """
        Rscript {params.script} {input.gtf_file} {input.chrom_len} {params.resPath}
        """


rule annotate_ERs:
    """
    Compare ERs with known annotation and annotate ER type
    """
    input:
        opt_ers = rules.generate_ERs.output.opt_ers,
        gtf_txdb = rules.generate_txdb.output.txdb

    output:
        ers = "ERs/" + sample + "_ers_raw.txt"

    params:
        resPath = "ERs",
        tissue = sample,
        script = SNAKEDIR + "/scripts/annotate_ERs.R",
        scriptPath = SNAKEDIR + "/scripts/"

    conda: "envs/r_v3.yml"

    shell:
        """
        Rscript {params.script} {input.opt_ers} {input.gtf_txdb} {params.tissue} {params.resPath} {params.scriptPath}
        """


rule generate_intergenic_data:
    """
    Select optimal intergenic ERs (within 10 kb of a protein-coding gene)
    """
    input:
        ers = rules.annotate_ERs.output.ers,
        gtf_file = gtf

    output:
        data_three = "Intergenic_ERs/" + sample + "_3prime_intergenic_ers.txt",
        data_five = "Intergenic_ERs/" + sample + "_5prime_intergenic_ers.txt"

    params:
        resPath = "Intergenic_ERs",
        prefix = sample,
        script = SNAKEDIR + "/scripts/select_intergenic_ers.R"

    conda: "envs/r_v3.yml"

    shell:
        """
        Rscript {params.script} {input.ers} {input.gtf_file} {params.resPath} {params.prefix}
        """