Skip to content

Commit

Permalink
[BENCHMARK] Add Delly2
Browse files Browse the repository at this point in the history 
Signed-off-by: Lydia Buntrock <[email protected]>
  • Loading branch information
@Irallia
Irallia committed Oct 25, 2022
1 parent 8647ce3 commit 1490d87
Show file tree
Hide file tree
Showing 4 changed files with 56 additions and 3 deletions.
26 changes: 26 additions & 0 deletions test/benchmark/caller_comparison_short_read/workflow/rules/callers.smk
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
22 changes: 20 additions & 2 deletions test/benchmark/caller_comparison_short_read/workflow/rules/eval.smk
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand All @@ -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:
Expand Down Expand Up @@ -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"),
Expand All @@ -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:
Expand All @@ -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}
""")
7 changes: 6 additions & 1 deletion test/benchmark/caller_comparison_short_read/workflow/scripts/plot_all_results.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"))

Expand Down
4 changes: 4 additions & 0 deletions test/benchmark/config/caller_comparison_config.yaml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -96,6 +96,10 @@ quality_ranges:
from: 1
to: 60
step: 1
Delly2:
from: 1
to: 1500
step: 10
iGenVar:
from: 1
to: 50
Expand Down

0 comments on commit 1490d87

Please sign in to comment.