feat: add savana for SV calling

.github/workflows/main.yml

@@ -25,42 +25,50 @@ jobs:
               --show-failed-logs -j 10 --conda-cleanup-pkgs cache --dryrun
             report: false
             cleanup: false
+            testcase_gen: false
           - name: local input
             args: >-
               --configfile .test/config-simple/config.yaml
               --show-failed-logs -j 10 --conda-cleanup-pkgs cache
             report: true
             cleanup: false
+            testcase_gen: false
           - name: testcase generation
             args: >-
               results/testcases/one/freebayes/IX:314200 --configfile .test/config-simple/config.yaml
             report: false
             cleanup: false
+            testcase_gen: true
           - name: target regions
             args: >-
               --configfile .test/config-target-regions/config.yaml
             report: true
             cleanup: false
+            testcase_gen: false
           - name: target regions, multiple BEDs
             args: >-
               --configfile .test/config-target-regions/config_multiple_beds.yaml
             report: true
             cleanup: false
+            testcase_gen: false
           - name: no candidate filtering
             args: >-
               --configfile .test/config-no-candidate-filtering/config.yaml
             report: true
             cleanup: false
+            testcase_gen: false
           - name: primers
             args: >-
               --configfile .test/config_primers/config.yaml
             report: true
             cleanup: false
+            testcase_gen: false
           - name: sra download
             args: >-
               --configfile .test/config-sra/config.yaml
             report: true
             cleanup: false
+            testcase_gen: false
     name: test ${{ matrix.case.name }}
 
     steps:
@@ -92,7 +100,17 @@ jobs:
           directory: .test
           snakefile: workflow/Snakefile
           args: >
-            --conda-cleanup-pkgs cache --sdm conda --cores 10 ${{ matrix.case.args }} 
+            --conda-cleanup-pkgs cache --sdm conda --cores 10 ${{ matrix.case.args }}
+          show-disk-usage-on-error: true
+
+      - name: test savana
+        uses: snakemake/snakemake-github-action@v2
+        if: ${{ !matrix.case.testcase_gen }}
+        with:
+          directory: .test
+          snakefile: workflow/Snakefile
+          args: >
+            only_savana --conda-cleanup-pkgs cache --sdm conda --cores 10 ${{ matrix.case.args }}
           show-disk-usage-on-error: true
 
       - name: generate report

.test/config-chm-eval/config.yaml

@@ -41,6 +41,10 @@ calling:
     activate: true
   freebayes:
     activate: true
+  savana:
+    activate: true
+    germline_events:
+      - present
   # See https://varlociraptor.github.io/docs/calling/#generic-variant-calling
   scenario: config-chm-eval/scenario.yaml
   # See http://snpeff.sourceforge.net/SnpSift.html#filter

.test/config-giab/config.yaml

@@ -85,6 +85,10 @@ calling:
     activate: false
   freebayes:
     activate: true
+  savana:
+    activate: true
+    germline_events:
+      - present
   # See https://varlociraptor.github.io/docs/calling/#generic-variant-calling
   scenario: config-giab/scenario.yaml
   filter:

.test/config-no-candidate-filtering/config.yaml

@@ -42,6 +42,8 @@ calling:
     activate: true
   freebayes:
     activate: true
+  savana:
+    activate: true
   # See https://varlociraptor.github.io/docs/calling/#generic-variant-calling
   scenario: config-no-candidate-filtering/scenario.yaml
   # See http://snpeff.sourceforge.net/SnpSift.html#filter

.test/config-simple/config.yaml

@@ -41,6 +41,8 @@ calling:
     activate: true
   freebayes:
     activate: true
+  savana:
+    activate: true
   # See https://varlociraptor.github.io/docs/calling/#generic-variant-calling
   scenario: config-simple/scenario.yaml
   # See http://snpeff.sourceforge.net/SnpSift.html#filter

.test/config-sra/config.yaml

@@ -39,6 +39,8 @@ calling:
     activate: true
   freebayes:
     activate: true
+  savana:
+    activate: true
   # See https://varlociraptor.github.io/docs/calling/#generic-variant-calling
   scenario: config-sra/scenario.yaml
   # See http://snpeff.sourceforge.net/SnpSift.html#filter

.test/config-target-regions/config.yaml

@@ -42,6 +42,8 @@ calling:
     activate: true
   freebayes:
     activate: true
+  savana:
+    activate: true
   # See https://varlociraptor.github.io/docs/calling/#generic-variant-calling
   scenario: config-target-regions/scenario.yaml
   # See http://snpeff.sourceforge.net/SnpSift.html#filter

.test/config_primers/config.yaml

@@ -43,6 +43,8 @@ calling:
     activate: true
   freebayes:
     activate: true
+  savana:
+    activate: true
   # See https://varlociraptor.github.io/docs/calling/#generic-variant-calling
   scenario: config_primers/scenario.yaml
   # See http://snpeff.sourceforge.net/SnpSift.html#filter

config/config.yaml

@@ -115,6 +115,12 @@ calling:
     # exclude_regions:
   freebayes:
     activate: true
+  savana:
+    activate: true
+    # Events from the varlociraptor scenario to use as germline event for CNV calling
+    # with savana. If not set, no CNV calling will be performed.
+    germline_events:
+      - germline
   # See https://varlociraptor.github.io/docs/calling/#generic-variant-calling
   scenario: config/scenario.yaml
   filter:

workflow/Snakefile

@@ -27,6 +27,7 @@ include: "rules/benchmarking.smk"
 include: "rules/ref.smk"
 include: "rules/trimming.smk"
 include: "rules/mapping.smk"
+include: "rules/germline_snvs.smk"
 include: "rules/candidate_calling.smk"
 include: "rules/calling.smk"
 include: "rules/annotation.smk"
@@ -61,6 +62,14 @@ rule only_alignment:
         expand("results/mapped/{sample}.cram", sample=samples["sample_name"]),
 
 
+rule only_savana:
+    input:
+        expand(
+            "results/candidate-calls/{sample}.savana.bcf",
+            sample=samples["sample_name"],
+        ),
+
+
 rule benchmark:
     input:
         expand(

workflow/envs/savana.yaml

@@ -0,0 +1,13 @@
+channels:
+  - conda-forge
+  - bioconda
+  - nodefaults
+dependencies:
+  - pybedtools==0.9.0
+  - pysam==0.20.0
+  - cyvcf2==0.30.16
+  - scikit-learn==1.2.2
+  - pandas==2.0.0
+  - matplotlib==3.7.1
+  - pip:
+    - git+https://github.com/johanneskoester/savana.git@fix/empty-calls

workflow/rules/candidate_calling.smk

@@ -23,6 +23,32 @@ rule freebayes:
         "v2.7.0/bio/freebayes"
 
 
+rule savana:
+    input:
+        ref=access.random(genome),
+        ref_idx=genome_fai,
+        aln=access.random("results/recal/{sample}.bam"),
+        index="results/recal/{sample}.bai",
+        germline_snvs="results/germline-snvs/{group}.bcf" if germline_events else [],
+    output:
+        bcf="results/candidate-calls/{sample}.savana.bcf",
+        outdir=directory("results/candidate-calls/savana/{sample}"),
+    log:
+        "logs/savana/{sample}.log",
+    params:
+        snvs=lambda w, input: (
+            f"--snp_vcf {input.germline_snvs}" if germline_events else ""
+        ),
+    conda:
+        "../envs/savana.yaml"
+    threads: 8
+    shell:
+        "(savana to --tumour {input.aln} --ref {input.ref} --outdir {output.outdir}"
+        " {params.snvs} &&"
+        " bcftools view -Ob -o {output.bcf} {output.outdir}/{wildcards.sample}.sv_breakpoints.vcf) "
+        "2>&1 > {log}"
+
+
 rule delly:
     input:
         ref=access.random(genome),

workflow/rules/common.smk

@@ -44,6 +44,8 @@ delly_excluded_regions = {
     ("homo_sapiens", "GRCh37"): "human.hg19",
 }
 
+germline_events = lookup("calling/savana/germline_events", within=config, default=[])
+
 
 def _group_or_sample(row):
     group = row.get("group", None)

workflow/rules/germline_snvs.smk

@@ -0,0 +1,27 @@
+rule gather_germline_calls:
+    input:
+        calls="results/calls/{group}.freebayes.{scatteritem}.bcf",
+        idx="results/calls/{group}.freebayes.{scatteritem}.bcf.csi",
+    output:
+        pipe("results/germline-snvs/{group}.germline_snv_candidates.bcf"),
+    log:
+        "logs/germline-snvs/gather-calls/{group}.log",
+    params:
+        extra="-a",
+    wrapper:
+        "v2.3.2/bio/bcftools/concat"
+
+
+rule control_fdr_germline_snvs:
+    input:
+        "results/germline-snvs/{group}.germline_snv_candidates.bcf",
+    output:
+        "results/germline-snvs/{group}.bcf",
+    params:
+        events=germline_events,
+    log:
+        "logs/germline-snvs/fdr-control/{group}.log",
+    conda:
+        "../envs/varlociraptor.yaml"
+    shell:
+        "varlociraptor filter-calls control-fdr {input} --events {params.events} --mode local-smart --var snv --fdr 0.05 > {output} 2> {log}"