WIP: implementing second case: run CASE MODE using precomputed model...

snappy_pipeline/workflows/targeted_seq_cnv_calling/Snakefile

@@ -75,3 +75,8 @@ analsys_type = wf.substep_getattr("gcnv", "get_analysis_type")()
 if analsys_type == "cohort_mode":
 
     include: "gcnv_cohort_mode.rules"
+
+
+elif analsys_type == "case_mode_with_model":
+
+    include: "gcnv_case_mode_with_model.rules"

snappy_pipeline/workflows/targeted_seq_cnv_calling/gcnv_case_mode_with_model.rules

@@ -0,0 +1,86 @@
+# gCNV with CASE MODE using procomputed model.
+
+
+rule targeted_seq_cnv_calling_gcnv_preprocess_intervals:
+    input:
+        unpack(wf.get_input_files("gcnv", "preprocess_intervals")),
+    output:
+        **wf. get_output_files("gcnv","preprocess_intervals"),
+    log:
+        wf.get_log_file("gcnv", "preprocess_intervals"),
+    wrapper:
+        wf.wrapper_path("gcnv/preprocess_intervals")
+
+
+rule targeted_seq_cnv_calling_gcnv_coverage:
+    input:
+        unpack(wf.get_input_files("gcnv", "coverage")),
+    output:
+        **wf. get_output_files("gcnv","coverage"),
+    params:
+        args=wf.get_params("gcnv", "coverage"),
+    log:
+        wf.get_log_file("gcnv", "coverage"),
+    wrapper:
+        wf.wrapper_path("gcnv/coverage")
+
+
+rule targeted_seq_cnv_calling_gcnv_contig_ploidy:
+    input:
+        unpack(wf.get_input_files("gcnv", "contig_ploidy_case_mode")),
+    output:
+        **wf. get_output_files("gcnv","contig_ploidy_case_mode"),
+    log:
+        wf.get_log_file("gcnv", "contig_ploidy_case_mode"),
+    params:
+        args=wf.get_params("gcnv", "model"),
+    wrapper:
+        wf.wrapper_path("gcnv/contig_ploidy_case_mode")
+
+
+rule targeted_seq_cnv_calling_gcnv_call_cnvs:
+    input:
+        unpack(wf.get_input_files("gcnv", "call_cnvs_case_mode")),
+    output:
+        **wf. get_output_files("gcnv","call_cnvs_case_mode"),
+    log:
+        wf.get_log_file("gcnv", "call_cnvs_case_mode"),
+    params:
+        args=wf.get_params("gcnv", "model"),
+    wrapper:
+        wf.wrapper_path("gcnv/call_cnvs_case_mode")
+
+
+rule targeted_seq_cnv_calling_gcnv_post_germline_calls:
+    input:
+        unpack(wf.get_input_files("gcnv", "post_germline_calls_case_mode")),
+    output:
+        **wf. get_output_files("gcnv","post_germline_calls"),
+    log:
+        wf.get_log_file("gcnv", "post_germline_calls"),
+    wrapper:
+        wf.wrapper_path("gcnv/post_germline_calls")
+
+
+rule targeted_seq_cnv_calling_gcnv_merge_cohort_vcfs:
+    input:
+        wf.get_input_files("gcnv", "merge_cohort_vcfs"),
+    output:
+        **wf. get_output_files("gcnv","merge_cohort_vcfs"),
+    log:
+        wf.get_log_file("gcnv", "merge_cohort_vcfs"),
+    wrapper:
+        wf.wrapper_path("gcnv/merge_cohort_vcfs")
+
+
+rule targeted_seq_cnv_calling_gcnv_extract_ped:
+    input:
+        unpack(wf.get_input_files("gcnv", "extract_ped")),
+    output:
+        **( wf. get_output_files("gcnv","extract_ped")),
+    log:
+        wf.get_log_file("gcnv", "extract_ped"),
+    params:
+        ped_members=wf.substep_getattr("gcnv", "get_ped_members"),
+    wrapper:
+        wf.wrapper_path("gcnv/extract_ped")

snappy_pipeline/workflows/targeted_seq_cnv_calling/gcnv_cohort_mode.rules

@@ -1,5 +1,4 @@
-# COHORT mode
-# gCNV needs many rules, thus they are in their own file.
+# gCNV with COHORT MODE using procomputed model.
 
 
 rule targeted_seq_cnv_calling_gcnv_call_cnvs:
@@ -37,11 +36,11 @@ rule targeted_seq_cnv_calling_gcnv_post_germline_calls:
     input:
         unpack(targeted_seq_cnv_calling_gcnv_post_germline_calls_input),
     output:
-        **wf. get_output_files("gcnv","post_germline_calls_cohort_mode"),
+        **wf. get_output_files("gcnv","post_germline_calls"),
     log:
-        wf.get_log_file("gcnv", "post_germline_calls_cohort_mode"),
+        wf.get_log_file("gcnv", "post_germline_calls"),
     wrapper:
-        wf.wrapper_path("gcnv/post_germline_calls_cohort_mode")
+        wf.wrapper_path("gcnv/post_germline_calls")
 
 
 rule targeted_seq_cnv_calling_gcnv_merge_cohort_vcfs:

snappy_wrappers/wrappers/gcnv/call_cnvs_case_mode/wrapper.py

@@ -0,0 +1,24 @@
+# -*- coding: utf-8 -*-
+from snakemake.shell import shell
+
+paths_tsv = " ".join(snakemake.input.tsv)
+shell(
+    r"""
+set -x
+set -euo pipefail
+export TMPDIR=$(mktemp -d)
+trap "rm -rf $TMPDIR" ERR EXIT
+
+export MKL_NUM_THREADS=16
+export OMP_NUM_THREADS=16
+export THEANO_FLAGS="base_compiledir=$TMPDIR/theano_compile_dir"
+
+gatk GermlineCNVCaller \
+    --run-mode CASE \
+    $(for tsv in {paths_tsv}; do echo -I $tsv; done) \
+    --contig-ploidy-calls $(dirname {snakemake.input.ploidy})/ploidy-calls \
+    --model {snakemake.params.args[model]} \
+    --output $(dirname {snakemake.output.done}) \
+    --output-prefix cnv_calls
+"""
+)

snappy_wrappers/wrappers/gcnv/contig_ploidy_case_mode/environment.yaml

@@ -0,0 +1,8 @@
+channels:
+  - conda-forge
+  - bioconda
+dependencies:
+  - gatk4 ==4.1.7.0
+  - gcnvkernel ==0.7
+  - bcftools
+  - htslib

snappy_wrappers/wrappers/gcnv/contig_ploidy_case_mode/wrapper.py

@@ -0,0 +1,19 @@
+# -*- coding: utf-8 -*-
+
+from snakemake.shell import shell
+
+paths_tsv = " ".join(snakemake.input.tsv)
+
+shell(
+    r"""
+export THEANO_FLAGS="base_compiledir=$TMPDIR/theano_compile_dir"
+
+set -x
+
+gatk DetermineGermlineContigPloidy \
+    --model {snakemake.params.args[model]} \
+    $(for tsv in {paths_tsv}; do echo -I $tsv; done) \
+    --output $(dirname {snakemake.output}) \
+    --output-prefix ploidy
+"""
+)

snappy_wrappers/wrappers/gcnv/post_germline_calls/environment.yaml

@@ -0,0 +1 @@
+../environment.yaml