-
Notifications
You must be signed in to change notification settings - Fork 0
/
Snakefile
141 lines (111 loc) · 4.36 KB
/
Snakefile
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
# Copyright (C) 2018 Leiden University Medical Center, Sander Bollen
# This file is part of gnomad-merge
#
# pytest-workflow is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as
# published by the Free Software Foundation, either version 3 of the
# License, or (at your option) any later version.
#
# pytest-workflow is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with gnomad-merge. If not, see <https://www.gnu.org/licenses/
"""
Snakemake pipeline to merge gnomad exomes with genomes into one
big VCF containing only global allele frequencies.
:copyright: (c) 2018 Leiden University Medical Center
:copyright: (c) 2018 Sander Bollen
"""
from pathlib import Path
gnomad_data_dir = config.get("GNOMAD_DATA_DIR")
reference_fasta = config.get("REFERENCE_FASTA")
chunksize = config.get("CHUNKSIZE", 1000)
# genome set comes with additional vcf files for coding regions.
# these need to be pruned from from our inputs, lest we count acs and ans
# twice for these variants.
exclude_patterns = config.get("EXCLUDE_PATTERNS", "coding").split(",")
if gnomad_data_dir is None:
raise ValueError("--config GNOMAD_DATA_DIR must be set")
if reference_fasta is None:
raise ValueError("--config REFERENCE_FASTA must be set")
# can't use __file__ for location of Snakefile, have to use snakemake
# built-in scrdir
script_dir = Path(srcdir("src"))
create_py = str(script_dir / Path("create_db.py"))
export_py = str(script_dir / Path("db_to_vcf.py"))
# gnomad_data folder structure looks like:
# - vcf
# --- exomes
# ------ gnomad.exomes.r2.0.1.sites.vcf.gz
# --- genomes
# ------ gnomad.genomes.r2.0.1.sites.1.vcf.gz
# ------ gnomad.genomes.r2.0.1.sites.2.vcf.gz
vcf_dir = Path(gnomad_data_dir) / Path("vcf")
exome_dir = vcf_dir / Path("exomes")
genome_dir = vcf_dir / Path("genomes")
# get some sensible name for a file
# e.g. exomes.r2.0.1.sites
# so it can be used as a wildcard
def get_name_for_file(path: Path):
return path.name.split(".vcf.gz")[0].split("gnomad.")[-1]
def name_has_patterns(patterns, name):
return any(x in name for x in patterns)
# generate dictionary of names -> absolute paths
all_input_vcf_files = {}
for x in exome_dir.iterdir():
if x.name.endswith(".vcf.gz") and not name_has_patterns(
exclude_patterns, x.name):
all_input_vcf_files[get_name_for_file(x)] = x.absolute()
for g in genome_dir.iterdir():
if g.name.endswith(".vcf.gz") and not name_has_patterns(
exclude_patterns, g.name):
all_input_vcf_files[get_name_for_file(g)] = g.absolute()
fnames = list(all_input_vcf_files.keys())
def get_file_for_name(wildcards):
"""Input function to get absolute path pertaining to a file"""
return str(all_input_vcf_files[wildcards.fname])
rule all:
input: "exports/gnomad.all.vcf.gz"
rule decompose:
"""Decompose VCF files"""
input: get_file_for_name
output: temp("decomposed/{fname}.vcf.gz")
conda: "envs/vt.yml"
shell: "vt decompose -s {input} -o {output}"
rule normalize:
"""Normalize VCF files"""
input:
vcf="decomposed/{fname}.vcf.gz",
ref=reference_fasta
output: temp("normalized/{fname}.vcf.gz")
conda: "envs/vt.yml"
shell: "vt normalize {input.vcf} -r {input.ref} -o {output}"
rule fill_db:
"""File database with variants"""
input:
vcfs=expand("normalized/{fname}.vcf.gz", fname=fnames),
create_py=create_py
params:
chunksize=chunksize
output: "db/gnomad.db"
conda: "envs/create_db.yml"
shell: "python {input.create_py} --chunksize {params.chunksize} "
"-o {output} {input.vcfs}"
rule export_db:
"""Export db back to VCF format"""
input:
db="db/gnomad.db",
export_py=export_py
output: temp("exports/unsorted.vcf")
conda: "envs/create_db.yml"
shell: "python {input.export_py} {input.db} > {output}"
rule sort_bgzip:
"""Sort and bgzip VCF file"""
input: "exports/unsorted.vcf"
output: "exports/gnomad.all.vcf.gz"
conda: "envs/sort.yml"
shell: "bedtools sort -header -i {input} | bgzip -c > {output} "
"&& tabix -pvcf {output}"