diff --git a/test/benchmark/caller_comparison_short_read/workflow/rules/callers.smk b/test/benchmark/caller_comparison_short_read/workflow/rules/callers.smk index 6cb6e2b2..268ff5fa 100644 --- a/test/benchmark/caller_comparison_short_read/workflow/rules/callers.smk +++ b/test/benchmark/caller_comparison_short_read/workflow/rules/callers.smk @@ -199,3 +199,29 @@ rule copy_Vaquita_LR_results: sed -i 's/S2L2/VaquitaLR_SL2/g' {output.res_SL2} sed -i 's/S2L3/VaquitaLR_SL3/g' {output.res_SL3} """) + +# Delly2 +rule run_Delly2: + output: + bcf = "results/caller_comparison_short_read/{dataset}/Delly2/variants.bcf" + log: + "logs/caller_comparison_short_read/Delly2_output.{dataset}.log" + threads: 8 + run: + if wildcards.dataset == 'Illumina_Paired_End': + short_bam = config["short_read_bam"]["s1"], + genome = config["reference_fa"]["Illumina_Paired_End"] + else: # wildcards.dataset == 'Illumina_Mate_Pair' + short_bam = config["short_read_bam"]["s2"], + genome = config["reference_fa"]["Illumina_Mate_Pair"] + shell(""" + export OMP_NUM_THREADS={threads} && \ + /usr/bin/time -v \ + delly call --outfile {output.bcf} --genome {genome} --map-qual 1 {short_bam} &>> {log} + """) + # -t [ --type ] arg (=DEL) SV type (DEL, DUP, INV, TRA, INS) + # -s [ --mad-cutoff ] arg (=9) insert size cutoff, median+s*MAD (deletions only) + # -n [ --noindels ] no small InDel calling + # Genotyping options: + # -v [ --vcffile ] arg input VCF/BCF file for re-genotyping + # -u [ --geno-qual ] arg (=5) min. mapping quality for genotyping diff --git a/test/benchmark/caller_comparison_short_read/workflow/rules/eval.smk b/test/benchmark/caller_comparison_short_read/workflow/rules/eval.smk index 46fd135a..7b2ded04 100644 --- a/test/benchmark/caller_comparison_short_read/workflow/rules/eval.smk +++ b/test/benchmark/caller_comparison_short_read/workflow/rules/eval.smk @@ -15,6 +15,16 @@ rule filter_vaquita_vcf: shell: "bcftools view -i 'INFO/SC>={wildcards.min_qual}' {input.vcf} > {output.vcf}" +rule filter_bcf: + input: + bcf = "results/caller_comparison_short_read/{dataset}/Delly2/variants.bcf" + output: + vcf ="results/caller_comparison_short_read/{dataset}/Delly2/variants.min_qual_{min_qual}.vcf" + conda: + "../../../envs/bcftools.yaml" + shell: + "bcftools convert {input.bcf} | bcftools view -i 'QUAL>={wildcards.min_qual}' > {output.vcf}" + rule bgzip: input: "{name}.vcf" @@ -37,6 +47,8 @@ rule truvari: index = "results/caller_comparison_short_read/{dataset}/{caller}/variants.min_qual_{min_qual}.vcf.gz.tbi" params: output_dir = "results/caller_comparison_short_read/{dataset}/eval/{caller}/min_qual_{min_qual}" + output: + summary = "results/caller_comparison_short_read/{dataset}/eval/{caller}/min_qual_{min_qual}/summary.txt" log: "logs/caller_comparison_short_read/truvari/truvari_output.{dataset}.{caller}.{min_qual}.log" run: @@ -89,7 +101,11 @@ rule cat_truvari_results_all: VaquitaLR_SL2 = expand("results/caller_comparison_short_read/{{dataset}}/eval/VaquitaLR_SL2/DUP_as_INS.min_qual_{min_qual}/pr_rec.txt", min_qual = min_qual_VaquitaLR), VaquitaLR_SL3 = expand("results/caller_comparison_short_read/{{dataset}}/eval/VaquitaLR_SL3/DUP_as_INS.min_qual_{min_qual}/pr_rec.txt", - min_qual = min_qual_VaquitaLR) + min_qual = min_qual_VaquitaLR), + Delly2 = expand("results/caller_comparison_short_read/{{dataset}}/eval/Delly2/min_qual_{min_qual}/pr_rec.txt", + min_qual=list(range(config["quality_ranges"]["Delly2"]["from"], + config["quality_ranges"]["Delly2"]["to"], + config["quality_ranges"]["Delly2"]["step"]))) output: iGenVar_S = temp("results/caller_comparison_short_read/{{dataset}}/eval/igenvar_s.all_results.txt"), iGenVar_SL1 = temp("results/caller_comparison_short_read/{{dataset}}/eval/igenvar_sl1.all_results.txt"), @@ -100,6 +116,7 @@ rule cat_truvari_results_all: VaquitaLR_SL1 = temp("results/caller_comparison_short_read/{{dataset}}/eval/vaquitaLR_SL1.all_results.txt"), VaquitaLR_SL2 = temp("results/caller_comparison_short_read/{{dataset}}/eval/vaquitaLR_SL2.all_results.txt"), VaquitaLR_SL3 = temp("results/caller_comparison_short_read/{{dataset}}/eval/vaquitaLR_SL3.all_results.txt"), + Delly2 = temp("results/caller_comparison_short_read/{{dataset}}/eval/Delly2.all_results.txt"), all = "results/caller_comparison_short_read/{dataset}/eval/all_results.txt" threads: 1 run: @@ -112,8 +129,9 @@ rule cat_truvari_results_all: shell("cat {input.VaquitaLR_SL1} > {output.VaquitaLR_SL1}") shell("cat {input.VaquitaLR_SL2} > {output.VaquitaLR_SL2}") shell("cat {input.VaquitaLR_SL3} > {output.VaquitaLR_SL3}") + shell("cat {input.Delly2} > {output.Delly2}") shell(""" cat {output.iGenVar_S} {output.iGenVar_SL1} {output.iGenVar_SL2} {output.iGenVar_SL3} \ {output.Vaquita} {output.VaquitaLR_S} {output.VaquitaLR_SL1} {output.VaquitaLR_SL2} {output.VaquitaLR_SL3} \ - > {output.all} + {output.Delly2} > {output.all} """) diff --git a/test/benchmark/caller_comparison_short_read/workflow/scripts/plot_all_results.R b/test/benchmark/caller_comparison_short_read/workflow/scripts/plot_all_results.R index 248bef7b..1d62c184 100644 --- a/test/benchmark/caller_comparison_short_read/workflow/scripts/plot_all_results.R +++ b/test/benchmark/caller_comparison_short_read/workflow/scripts/plot_all_results.R @@ -24,20 +24,25 @@ total <- rbind(res, f1) total$caller = factor(total$caller, levels=c('iGenVar_S', 'iGenVar_SL1', 'iGenVar_SL2', 'iGenVar_SL3', 'Vaquita', 'VaquitaLR_S', 'VaquitaLR_SL1', 'VaquitaLR_SL2', 'VaquitaLR_SL3', + 'Delly2', '0.1', '0.2', '0.3', '0.4', '0.5', '0.6', '0.7', '0.8', '0.9'), labels=c('iGenVar 0.0.3: S', 'iGenVar: SL1', 'iGenVar: SL2', 'iGenVar: SL3', 'Vaquita 0.4.0', 'Vaquita LR: S', 'Vaquita LR: SL1', 'Vaquita LR: SL2', 'Vaquita LR: SL3', + 'Delly2 1.0.3', '0.1', '0.2', '0.3', '0.4', '0.5', '0.6', '0.7', '0.8', '0.9')) scale_custom <- list( # https://stackoverflow.com/questions/46803260/assigning-40-shapes-or-more-in-scale-shape-manual scale_shape_manual(values = c(15, 16, 16, 16, 17, 18, 18, 18, 18, + 19, 46, 46, 46, 46, 46, 46, 46, 46, 46)), # https://www.farb-tabelle.de/de/farbtabelle.htm scale_color_manual(breaks = c('iGenVar 0.0.3: S', 'iGenVar: SL1', 'iGenVar: SL2', 'iGenVar: SL3', - 'Vaquita 0.4.0', 'Vaquita LR: S', 'Vaquita LR: SL1', 'Vaquita LR: SL2', 'Vaquita LR: SL3'), + 'Vaquita 0.4.0', 'Vaquita LR: S', 'Vaquita LR: SL1', 'Vaquita LR: SL2', 'Vaquita LR: SL3', + 'Delly2 1.0.3'), values = c("chartreuse3", "chartreuse1", "chartreuse1", "chartreuse1", "darkorchid1", "maroon1", "maroon3", "maroon3", "maroon3", + "dodgerblue", "tan", "tan", "tan", "tan", "tan", "tan", "tan", "tan", "tan"), na.value = "gray80")) diff --git a/test/benchmark/config/caller_comparison_config.yaml b/test/benchmark/config/caller_comparison_config.yaml index 1adac294..8a2c40be 100644 --- a/test/benchmark/config/caller_comparison_config.yaml +++ b/test/benchmark/config/caller_comparison_config.yaml @@ -96,6 +96,10 @@ quality_ranges: from: 1 to: 60 step: 1 + Delly2: + from: 1 + to: 1500 + step: 10 iGenVar: from: 1 to: 50