feat: fusion calling

.test/config-chm-eval/samples.tsv

@@ -1,2 +1,2 @@
-sample_name	group	alias	platform
-chm		chm	ILLUMINA
+sample_name	group	alias	platform	datatype
+chm		chm	ILLUMINA	dna

.test/config-giab/samples.tsv

@@ -1,2 +1,2 @@
-sample_name	alias	group	platform	purity
-NA12878	NA12878	NA12878	ILLUMINA	
+sample_name	alias	group	platform	purity	datatype
+NA12878	NA12878	NA12878	ILLUMINA		dna

.test/config-no-candidate-filtering/samples.tsv

@@ -1,3 +1,3 @@
-sample_name	group	alias	platform
-a		a	ILLUMINA
-b		b	ILLUMINA
+sample_name	group	alias	platform	datatype
+a		a	ILLUMINA	dna
+b		b	ILLUMINA	dna

.test/config-simple/samples.tsv

@@ -1,5 +1,5 @@
-sample_name	group	alias	platform
-a	one	x	ILLUMINA
-b	one	y	ILLUMINA
-b	two	x	ILLUMINA
-a	two	y	ILLUMINA
+sample_name	group	alias	platform	datatype
+a	one	x	ILLUMINA	dna
+b	one	y	ILLUMINA	dna
+b	two	x	ILLUMINA	dna
+a	two	y	ILLUMINA	dna

.test/config-sra/samples.tsv

@@ -1,3 +1,3 @@
-sample_name	group	alias	platform
-a		a	ILLUMINA
-b		b	ILLUMINA
+sample_name	group	alias	platform	datatype
+a		a	ILLUMINA	dna
+b		b	ILLUMINA	dna

.test/config-target-regions/samples.tsv

@@ -1,3 +1,3 @@
-sample_name	group	alias	platform
-a		a	ILLUMINA
-b		b	ILLUMINA
+sample_name	group	alias	platform	datatype
+a		a	ILLUMINA	dna
+b		b	ILLUMINA	dna

.test/config_primers/samples.tsv

@@ -1,3 +1,3 @@
-sample_name	group	alias	platform
-a		a	ILLUMINA
-b		b	ILLUMINA
+sample_name	group	alias	platform	datatype
+a		a	ILLUMINA	dna
+b		b	ILLUMINA	dna

config/config.yaml

@@ -117,6 +117,7 @@ calling:
       # only alphanumerics and underscores
       # ("somatic" below is just an example and can be modified as needed).
       some_id:
+        types: ["variants", "fusions"]
         # labels for the callset, displayed in the report. Will fall back to id if no labels specified
         labels:
           some-label: label text

config/samples.tsv

@@ -1 +1 @@
-sample_name	alias	group	platform	purity	panel	umi_read	umi_read_structure
+sample_name	alias	group	platform	purity	panel	datatype	umi_read	umi_read_structure

workflow/Snakefile

@@ -35,6 +35,7 @@ include: "rules/table.smk"
 include: "rules/regions.smk"
 include: "rules/plugins.smk"
 include: "rules/datavzrd.smk"
+include: "rules/fusion_calling.smk"
 include: "rules/testcase.smk"
 
 

workflow/envs/arriba.yaml

@@ -0,0 +1,5 @@
+channels:
+  - conda-forge
+  - bioconda
+dependencies:
+  - arriba =2.4

workflow/envs/bcftools.yaml

@@ -2,4 +2,4 @@ channels:
   - conda-forge
   - bioconda
 dependencies:
-  - bcftools =1.14
+  - bcftools =1.16

workflow/resources/datavzrd/fusion-calls-template.datavzrd.yaml

@@ -0,0 +1,16 @@
+name: ?f"Fusion calls {wildcards.event}"
+
+default-view: ?f"{params.groups[0]}-fusions"
+max-in-memory-rows: 1500
+
+datasets:
+  ?for group, path in zip(params.groups, params.fusion_calls):
+    ?f"{group}-fusions":
+      path: ?path
+      separator: "\t"
+
+views:
+  ?for group in params.groups:
+    ?f"{group}-fusions":
+          desc: ?f"Fusion calls.\n{config['calling']['fdr-control']['events'][wildcards.event]['desc']}"
+          dataset: ?f"{group}-fusions"

workflow/rules/annotation.smk

@@ -24,16 +24,16 @@ rule annotate_candidate_variants:
 
 rule annotate_variants:
     input:
-        calls="results/calls/{group}.{scatteritem}.bcf",
+        calls="results/calls/{group}.{analysis}.{scatteritem}.bcf",
         cache="resources/vep/cache",
         plugins="resources/vep/plugins",
         revel=lambda wc: get_plugin_aux("REVEL"),
         revel_tbi=lambda wc: get_plugin_aux("REVEL", True),
         fasta=genome,
         fai=genome_fai,
     output:
-        calls="results/calls/{group}.{scatteritem}.annotated.bcf",
-        stats="results/calls/{group}.{scatteritem}.stats.html",
+        calls="results/calls/{group}.{analysis}.{scatteritem}.annotated.bcf",
+        stats="results/calls/{group}.{analysis}.{scatteritem}.stats.html",
     params:
         # Pass a list of plugins to use, see https://www.ensembl.org/info/docs/tools/vep/script/vep_plugins.html
         # Plugin args can be added as well, e.g. via an entry "MyPlugin,1,FOO", see docs.
@@ -42,7 +42,7 @@ rule annotate_variants:
             config["annotations"]["vep"]["final_calls"]["params"]
         ),
     log:
-        "logs/vep/{group}.{scatteritem}.annotate.log",
+        "logs/vep/{group}.{analysis}.{scatteritem}.annotate.log",
     threads: get_vep_threads()
     wrapper:
         "v2.5.0/bio/vep/annotate"
@@ -89,9 +89,9 @@ rule gather_annotated_calls:
         calls=get_gather_annotated_calls_input(),
         idx=get_gather_annotated_calls_input(ext="bcf.csi"),
     output:
-        "results/final-calls/{group}.annotated.bcf",
+        "results/final-calls/{group}.{analysis}.annotated.bcf",
     log:
-        "logs/gather-annotated-calls/{group}.log",
+        "logs/gather-annotated-calls/{group}.{analysis}.log",
     params:
         extra="-a",
     wrapper:

workflow/rules/calling.smk

@@ -26,7 +26,8 @@ rule varlociraptor_alignment_properties:
     input:
         ref=genome,
         ref_idx=genome_fai,
-        bam="results/recal/{sample}.bam",
+        bam=get_sample_bam,
+        bai=lambda wc: get_sample_bam(wc, bai=True),
     output:
         "results/alignment-properties/{group}/{sample}.json",
     log:
@@ -42,8 +43,8 @@ rule varlociraptor_preprocess:
         ref=genome,
         ref_idx=genome_fai,
         candidates=get_candidate_calls(),
-        bam="results/recal/{sample}.bam",
-        bai="results/recal/{sample}.bai",
+        bam=get_sample_bam,
+        bai=lambda wc: get_sample_bam(wc, bai=True),
         alignment_props="results/alignment-properties/{group}/{sample}.json",
     output:
         "results/observations/{group}/{sample}.{caller}.{scatteritem}.bcf",
@@ -66,7 +67,7 @@ rule varlociraptor_call:
         obs=get_group_observations,
         scenario="results/scenarios/{group}.yaml",
     output:
-        temp("results/calls/{group}.{caller}.{scatteritem}.bcf"),
+        temp("results/calls/{group}.{caller}.{scatteritem}.unsorted.bcf"),
     log:
         "logs/varlociraptor/call/{group}.{caller}.{scatteritem}.log",
     params:
@@ -86,28 +87,30 @@ rule varlociraptor_call:
 
 rule sort_calls:
     input:
-        "results/calls/{group}.{caller}.{scatteritem}.bcf",
+        "results/calls/{group}.{caller}.{scatteritem}.unsorted.bcf",
     output:
         temp("results/calls/{group}.{caller}.{scatteritem}.bcf"),
+    params:
+        # Set to True, in case you want uncompressed BCF output
+        uncompressed_bcf=False,
+        # Extra arguments
+        extras="",
     log:
         "logs/bcf-sort/{group}.{caller}.{scatteritem}.log",
-    conda:
-        "../envs/bcftools.yaml"
     resources:
         mem_mb=8000,
-    shell:
-        "bcftools sort --max-mem {resources.mem_mb}M --temp-dir `mktemp -d` "
-        "-Ob {input} > {output} 2> {log}"
+    wrapper:
+        "v1.21.0/bio/bcftools/sort"
 
 
 rule bcftools_concat:
     input:
         calls=get_scattered_calls(),
         indexes=get_scattered_calls(ext="bcf.csi"),
     output:
-        "results/calls/{group}.{scatteritem}.bcf",
+        "results/calls/{group}.{analysis}.{scatteritem}.bcf",
     log:
-        "logs/concat-calls/{group}.{scatteritem}.log",
+        "logs/concat-calls/{group}.{analysis}.{scatteritem}.log",
     params:
         extra="-a",  # TODO Check this
     wrapper: