Skip to content

Commit 765897c

Browse files
Smedsfgvieira
authored and
tgroth97
committed
feat: gatk modelsegments wrapper (snakemake#1321)
### Description gatk ModelSegments wrapper ### QC <!-- Make sure that you can tick the boxes below. --> * [ ] I confirm that: For all wrappers added by this PR, * there is a test case which covers any introduced changes, * `input:` and `output:` file paths in the resulting rule can be changed arbitrarily, * either the wrapper can only use a single core, or the example rule contains a `threads: x` statement with `x` being a reasonable default, * rule names in the test case are in [snake_case](https://en.wikipedia.org/wiki/Snake_case) and somehow tell what the rule is about or match the tools purpose or name (e.g., `map_reads` for a step that maps reads), * all `environment.yaml` specifications follow [the respective best practices](https://stackoverflow.com/a/64594513/2352071), * wherever possible, command line arguments are inferred and set automatically (e.g. based on file extensions in `input:` or `output:`), * all fields of the example rules in the `Snakefile`s and their entries are explained via comments (`input:`/`output:`/`params:` etc.), * `stderr` and/or `stdout` are logged correctly (`log:`), depending on the wrapped tool, * temporary files are either written to a unique hidden folder in the working directory, or (better) stored where the Python function `tempfile.gettempdir()` points to (see [here](https://docs.python.org/3/library/tempfile.html#tempfile.gettempdir); this also means that using any Python `tempfile` default behavior works), * the `meta.yaml` contains a link to the documentation of the respective tool or command, * `Snakefile`s pass the linting (`snakemake --lint`), * `Snakefile`s are formatted with [snakefmt](https://github.com/snakemake/snakefmt), * Python wrapper scripts are formatted with [black](https://black.readthedocs.io). * Conda environments use a minimal amount of channels, in recommended ordering. E.g. for bioconda, use (conda-forge, bioconda, nodefaults, as conda-forge should have highest priority and defaults channels are usually not needed because most packages are in conda-forge nowadays). --------- Co-authored-by: Filipe G. Vieira <1151762+fgvieira@users.noreply.github.com>
1 parent 71373f2 commit 765897c

File tree

6 files changed

+140
-0
lines changed

6 files changed

+140
-0
lines changed

‎bio/gatk/modelsegments/environment.yaml

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,7 @@
1+
channels:
2+
- conda-forge
3+
- bioconda
4+
- nodefaults
5+
dependencies:
6+
- gatk4 =4.4.0.0
7+
- snakemake-wrapper-utils =0.5.3

‎bio/gatk/modelsegments/meta.yaml

Lines changed: 19 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,19 @@
1+
name: gatk ModelSegments
2+
url: https://gatk.broadinstitute.org/hc/en-us/articles/13832747657883-ModelSegments
3+
description: |
4+
Models segmented copy ratios from denoised copy ratios and segmented minor-allele fractions from allelic counts
5+
authors:
6+
- Patrik Smeds
7+
input:
8+
- denoised_copy_ratios: denoised_copy_ratios file (optional)
9+
- allelic_counts: allelic_counts file (optional)
10+
- normal_allelic_counts: matched_normal allelic-counts (optional)
11+
- segments: segments Picard interval-list file containing a multisample segmentation output by a previous run (optional)
12+
output:
13+
- list of files ending with either '.modelFinal.seq', '.cr.seg', '.af.igv.seg', '.cr.igv.seg', '.hets.tsv', '.modelBegin.cr.param', '.modelBegin.af.param', '.modelBegin.seg', '.modelFinal.af.param' or '.modelFinal.cr.param'
14+
params:
15+
- java_opts: additional arguments to be passed to the java compiler, e.g. "-XX:ParallelGCThreads=10" (not for `-XmX` or `-Djava.io.tmpdir`, since they are handled automatically).
16+
- extra: additional program arguments.
17+
note: |
18+
* Expected input is either `denoised_copy_ratios` or `allelic_counts`
19+
* `normal_allelic_counts` must be provided together with `allelic_counts`

‎bio/gatk/modelsegments/test/Snakefile

Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,17 @@
1+
rule modelsegments_denoise_input:
2+
input:
3+
denoised_copy_ratios="a.denoisedCR.tsv",
4+
output:
5+
"a.den.modelFinal.seg",
6+
"a.n.cr.seg",
7+
log:
8+
"logs/gatk/modelsegments_denoise.log",
9+
params:
10+
#prefix="a.den.test",
11+
extra="", # optional
12+
java_opts="", # optional
13+
resources:
14+
mem_mb=1024,
15+
wrapper:
16+
"master/bio/gatk/modelsegments"
17+

‎bio/gatk/modelsegments/test/a.denoisedCR.tsv

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,6 @@
1+
@HD VN:1.6
2+
@SQ SN:ref LN:45
3+
@SQ SN:ref2 LN:40
4+
@RG ID:GATKCopyNumber SM:mysample
5+
CONTIG START END LOG2_COPY_RATIO
6+
ref2 1 30 0.000000

‎bio/gatk/modelsegments/wrapper.py

Lines changed: 83 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,83 @@
1+
__author__ = "Patrik Smeds"
2+
__copyright__ = "Copyright 2023, Patrik Smed"
3+
__email__ = "patrik.smeds@gmail.com"
4+
__license__ = "MIT"
5+
6+
import os
7+
import tempfile
8+
from snakemake.shell import shell
9+
from snakemake_wrapper_utils.java import get_java_opts
10+
11+
12+
denoised_copy_ratios = ""
13+
if snakemake.input.get("denoised_copy_ratios", None):
14+
denoised_copy_ratios = (
15+
f"--denoised-copy-ratios {snakemake.input.denoised_copy_ratios}"
16+
)
17+
18+
allelic_counts = ""
19+
if snakemake.input.get("allelic_counts", None):
20+
allelic_counts = f"--allelic-counts {snakemake.input.allelic_counts}"
21+
22+
normal_allelic_counts = ""
23+
if snakemake.input.get("normal_allelic_counts", None):
24+
matched_normal_allelic_counts = (
25+
f"--normal-allelic-counts {snakemake.inut.normal_allelic_counts}"
26+
)
27+
28+
segments = ""
29+
if snakemake.input.get("segments", None):
30+
interval_list = f"--segments {snakemake.input.segments}"
31+
32+
if not allelic_counts and not denoised_copy_ratios:
33+
raise Exception(
34+
"wrapper input requires either 'allelic_counts' or 'denoise_copy_ratios' to be set"
35+
)
36+
37+
if normal_allelic_counts and not allelic_counts:
38+
raise Exception(
39+
"'allelica_counts' is required when 'normal-allelic-counts' is an input to the rule!"
40+
)
41+
42+
extra = snakemake.params.get("extra", "")
43+
java_opts = get_java_opts(snakemake)
44+
45+
log = snakemake.log_fmt_shell(stdout=True, stderr=True)
46+
47+
48+
with tempfile.TemporaryDirectory() as tmpdir:
49+
output_folder = os.path.join(tmpdir, "output_folder")
50+
shell(
51+
"gatk --java-options '{java_opts}' ModelSegments"
52+
" {segments}"
53+
" {denoised_copy_ratios}"
54+
" {allelic_counts}"
55+
" {normal_allelic_counts}"
56+
" --output-prefix temp_name__"
57+
" -O {output_folder}"
58+
" --tmp-dir {tmpdir}"
59+
" {extra}"
60+
" {log}"
61+
)
62+
63+
created_files = {}
64+
# Find all created files
65+
for new_file in os.listdir(output_folder):
66+
file_path = os.path.join(output_folder, new_file)
67+
if os.path.isfile(file_path):
68+
file_end = os.path.basename(file_path).split("__")[1]
69+
created_files[file_end] = file_path
70+
71+
# Match expected output with found files
72+
for output in snakemake.output:
73+
file_found = False
74+
for file_ending in created_files:
75+
if output.endswith(file_ending):
76+
shell(f"cp {created_files[file_ending]} {output}")
77+
file_found = True
78+
break
79+
if not file_found:
80+
created_files_list = [f"{e}" for e in created_files]
81+
raise IOError(
82+
f"Could not create file {output}, possible files ends with {created_files_list}"
83+
)

‎test.py

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -4131,6 +4131,14 @@ def test_gatk_haplotypecaller_gvcf():
41314131
)
41324132

41334133

4134+
@skip_if_not_modified
4135+
def test_gatk_modelsegments():
4136+
run(
4137+
"bio/gatk/modelsegments",
4138+
["snakemake", "--cores", "1", "a.den.modelFinal.seg", "--use-conda", "-F"],
4139+
)
4140+
4141+
41344142
@skip_if_not_modified
41354143
def test_gatk_variantrecalibrator():
41364144
def check_log(log):

0 commit comments

Comments
 (0)