From dff79fe05e1bc9defbe1547ff1e48fabb7bee392 Mon Sep 17 00:00:00 2001 From: Lydia Buntrock Date: Wed, 16 Jun 2021 17:25:27 +0200 Subject: [PATCH 1/5] [BENCHMARK] Copy files from https://github.com/eldariont/real_sv_calling_evaluation for a snakemake workflow. Signed-off-by: Lydia Buntrock --- test/benchmark/caller_comparison/Snakefile | 28 +++ .../caller_comparison/config/config.yaml | 24 +++ .../workflow/rules/callers.smk | 125 ++++++++++++ .../caller_comparison/workflow/rules/eval.smk | 147 ++++++++++++++ .../workflow/rules/plots.smk | 192 ++++++++++++++++++ .../workflow/scripts/plot_all_results.R | 27 +++ test/benchmark/envs/cyvcf2.yaml | 8 + .../{shell_scripts => envs}/environment.yml | 0 test/benchmark/envs/pbsv.yaml | 4 + test/benchmark/envs/samtools.yaml | 5 + test/benchmark/envs/sniffles.yaml | 4 + test/benchmark/envs/svim.yaml | 4 + test/benchmark/envs/truvari.yaml | 7 + test/benchmark/envs/truvari_environment.yaml | 8 - .../{shell_scripts => }/iGenVar_init.sh | 3 + .../{ => parameter_benchmarks}/Snakefile | 5 +- .../iGenVar_run_benchmark.sh | 2 +- .../plot_all.R | 0 18 files changed, 582 insertions(+), 11 deletions(-) create mode 100644 test/benchmark/caller_comparison/Snakefile create mode 100644 test/benchmark/caller_comparison/config/config.yaml create mode 100644 test/benchmark/caller_comparison/workflow/rules/callers.smk create mode 100644 test/benchmark/caller_comparison/workflow/rules/eval.smk create mode 100644 test/benchmark/caller_comparison/workflow/rules/plots.smk create mode 100644 test/benchmark/caller_comparison/workflow/scripts/plot_all_results.R create mode 100644 test/benchmark/envs/cyvcf2.yaml rename test/benchmark/{shell_scripts => envs}/environment.yml (100%) create mode 100644 test/benchmark/envs/pbsv.yaml create mode 100644 test/benchmark/envs/samtools.yaml create mode 100644 test/benchmark/envs/sniffles.yaml create mode 100644 test/benchmark/envs/svim.yaml create mode 100644 test/benchmark/envs/truvari.yaml delete mode 100644 test/benchmark/envs/truvari_environment.yaml rename test/benchmark/{shell_scripts => }/iGenVar_init.sh (90%) rename test/benchmark/{ => parameter_benchmarks}/Snakefile (96%) rename test/benchmark/{shell_scripts => parameter_benchmarks}/iGenVar_run_benchmark.sh (95%) rename test/benchmark/{scripts => parameter_benchmarks}/plot_all.R (100%) diff --git a/test/benchmark/caller_comparison/Snakefile b/test/benchmark/caller_comparison/Snakefile new file mode 100644 index 00000000..258ab9bb --- /dev/null +++ b/test/benchmark/caller_comparison/Snakefile @@ -0,0 +1,28 @@ +import os +import gzip +import sys +import shutil + +configfile: "config/config.yaml" + +VCFS=["giab", "giab.gt", "giab.seq", "giab.gt.seq"] +SVIM_THRESHOLDS = [1, 3, 4, 5, 6, 7, 8, 9, 10, 14, 18, 22, 26, 30, 40, 50, 60] + +include: "workflow/rules/callers.smk" +include: "workflow/rules/eval.smk" +include: "workflow/rules/plots.smk" + +##### Target rules ##### + +rule all: + input: + #SV lengths + expand("pipeline/SV-plots/SV-length_pbsv_{minscore}.png", minscore=[3, 5, 7]), + expand("pipeline/SV-plots/SV-length_Sniffles_{minscore}.png", minscore=[3, 5, 7]), + expand("pipeline/SV-plots/SV-length_SVIM_1000_900_1.0_0.5_{minscore}.png", minscore=[3, 5, 7]), + expand("pipeline/SV-plots/SV-counts.merged.png"), + #Evaluation + expand("pipeline/eval/results.all.png"), + expand("pipeline/eval/results.tools.{vcf}.png", vcf=VCFS), + expand("pipeline/eval/results.coverages.{vcf}.png", vcf=VCFS), + expand("pipeline/eval/results.svim.parameters.png"), diff --git a/test/benchmark/caller_comparison/config/config.yaml b/test/benchmark/caller_comparison/config/config.yaml new file mode 100644 index 00000000..26d3f0d3 --- /dev/null +++ b/test/benchmark/caller_comparison/config/config.yaml @@ -0,0 +1,24 @@ +long_bam: /data/long_reads/HG002.Sequel.10kb.pbmm2.hs37d5.whatshap.haplotag.RTG.10x.trio_sorted.bam +long_bai: /data/long_reads/HG002.Sequel.10kb.pbmm2.hs37d5.whatshap.haplotag.RTG.10x.trio_sorted.bam.bai + +reference: /path/to/database/genome.fa.gz + +truth: + giab: data/truth_set/HG002_SVs_Tier1_v0.6.vcf.gz + bed: data/truth_set/HG002_SVs_Tier1_v0.6.bed + +parameters: + pbmm_preset: CCS #SUBREAD (CLR), CCS (CCS) + survivor_distance: 500 + min_sv_size: 40 + +minimums: + sniffles_from: 1 + sniffles_to: 60 + sniffles_step: 2 + svim_from: 1 + svim_to: 100 + svim_step: 2 + pbsv_from: 10 + pbsv_to: 91 + pbsv_step: 10 diff --git a/test/benchmark/caller_comparison/workflow/rules/callers.smk b/test/benchmark/caller_comparison/workflow/rules/callers.smk new file mode 100644 index 00000000..33da51e4 --- /dev/null +++ b/test/benchmark/caller_comparison/workflow/rules/callers.smk @@ -0,0 +1,125 @@ +localrules: filter_svim, fix_sniffles, filter_insertions_and_deletions, filter_insertions_and_deletions_svim + +# SVIM +rule run_svim: + input: + bam = config["long_bam"], + bai = config["long_bai"], + genome = config["reference"] + output: + "pipeline/SVIM/{pmd,[0-9]+}_{pdn,[0-9]+}_{edn,[0-9\.]+}_{cmd,[0-9\.]+}/variants.vcf" + resources: + mem_mb = 20000, + time_min = 600, + io_gb = 100 + params: + working_dir = "pipeline/SVIM/{pmd}_{pdn}_{edn}_{cmd}/", + min_sv_size = config["parameters"]["min_sv_size"] + threads: 1 + conda: + "../../../envs/svim.yaml" + shell: + "svim alignment --sample {wildcards.data} \ + --partition_max_distance {wildcards.pmd} \ + --position_distance_normalizer {wildcards.pdn} \ + --edit_distance_normalizer {wildcards.edn} \ + --cluster_max_distance {wildcards.cmd} \ + --min_sv_size {params.min_sv_size} \ + --segment_gap_tolerance 20 \ + --segment_overlap_tolerance 20 \ + --interspersed_duplications_as_insertions \ + --tandem_duplications_as_insertions \ + --read_names \ + --max_sv_size 1000000 \ + --verbose \ + {params.working_dir} {input.bam} {input.genome}" + +rule filter_svim: + input: + "pipeline/SVIM/{parameters}/variants.vcf" + output: + temp("pipeline/SVIM/{parameters}/min_{minscore,[0-9]+}.vcf") + threads: 1 + shell: + "bcftools view -e \"GT=='0/0'\" {input} | \ + awk 'OFS=\"\\t\" {{ if($1 ~ /^#/) {{ print $0 }} \ + else {{ if($6>={wildcards.minscore}) {{ print $1, $2, $3, $4, $5, $6, \"PASS\", $8, $9, $10 }} }} }}' > {output}" + +# SNIFFLES +rule run_sniffles: + input: + bam = config["long_bam"], + bai = config["long_bai"] + output: + expand("pipeline/Sniffles/raw_{minsupport}.vcf", + minsupport=list(range(config["minimums"]["sniffles_from"], config["minimums"]["sniffles_to"]+1, config["minimums"]["sniffles_step"]))) + resources: + mem_mb = 400000, + time_min = 1200, + io_gb = 100 + params: + min_sv_size = config["parameters"]["min_sv_size"], + tmpdir = "300", + sniffles_from = config["minimums"]["sniffles_from"], + sniffles_to = config["minimums"]["sniffles_to"], + sniffles_step = config["minimums"]["sniffles_step"], + outdir = "pipeline/Sniffles/" + threads: 30 + conda: + "../../../envs/sniffles.yaml" + shell: + "bash workflow/scripts/run_sniffles.sh {input.bam} {input.bai} {params.sniffles_from} {params.sniffles_to} {params.sniffles_step} {params.min_sv_size} {threads} {params.outdir}" + +#see https://github.com/spiralgenetics/truvari/issues/43 +rule fix_sniffles: + input: + "pipeline/Sniffles/raw_{support}.vcf" + output: + "pipeline/Sniffles/min_{support,[0-9]+}.vcf" + shell: + "sed 's/##INFO= {output}" + +#PBSV +rule run_pbsv: + input: + bam = config["long_bam"], + bai = config["long_bai"] + genome = config["reference"], + output: + expand("pipeline/pbsv/min_{minsupport}.vcf", + minsupport=list(range(config["minimums"]["pbsv_from"], config["minimums"]["pbsv_to"]+1, config["minimums"]["pbsv_step"]))) + resources: + mem_mb = 400000, + time_min = 2000, + io_gb = 100 + params: + min_sv_size = config["parameters"]["min_sv_size"], + tmpdir = "1000", + pbsv_from = config["minimums"]["pbsv_from"], + pbsv_to = config["minimums"]["pbsv_to"], + pbsv_step = config["minimums"]["pbsv_step"], + outdir = "pipeline/pbsv/" + threads: 15 + conda: + "../../../envs/pbsv.yaml" + shell: + "bash workflow/scripts/run_pbsv.sh {input.bam} {input.bai} {input.genome} {params.pbsv_from} {params.pbsv_to} {params.pbsv_step} {params.min_sv_size} {threads} {params.outdir}" + +#Split to SV classes +rule filter_insertions_and_deletions: + input: + "pipeline/{caller}/min_{minscore}.vcf" + output: + "pipeline/{caller,Sniffles|pbsv}/min_{minscore,[0-9]+}.indel.vcf" + threads: 1 + shell: + "bcftools view -i 'SVTYPE=\"DEL\" | SVTYPE=\"INS\"' {input} | bcftools sort > {output}" + +rule filter_insertions_and_deletions_svim: + input: + "pipeline/SVIM/{parameters}/min_{minscore}.vcf" + output: + "pipeline/SVIM/{parameters}/min_{minscore,[0-9]+}.indel.vcf" + threads: 1 + shell: + "bcftools view -i 'SVTYPE=\"DEL\" | SVTYPE=\"INS\"' {input} | bcftools sort > {output}" diff --git a/test/benchmark/caller_comparison/workflow/rules/eval.smk b/test/benchmark/caller_comparison/workflow/rules/eval.smk new file mode 100644 index 00000000..c4af5b2a --- /dev/null +++ b/test/benchmark/caller_comparison/workflow/rules/eval.smk @@ -0,0 +1,147 @@ +localrules: bgzip, tabix, callset_eval_svim, callset_eval, reformat_truvari_results, reformat_truvari_results_svim, cat_truvari_results_all, cat_truvari_results_full, cat_truvari_results_svim_parameters + +def get_vcf(wildcards): + return config["truth"][wildcards.vcf.split(".")[0]] + +rule bgzip: + input: + "{name}.vcf" + output: + "{name}.vcf.gz" + shell: + "bgzip -c {input} > {output}" + + +rule tabix: + input: + "{name}.vcf.gz" + output: + "{name}.vcf.gz.tbi" + shell: + "tabix -p vcf {input}" + + +rule callset_eval_svim: + input: + genome = config["reference"], + truth_vcf = get_vcf, + truth_bed = config["truth"]["bed"], + calls = "pipeline/SVIM/{parameters}/min_{minscore}.indel.vcf.gz", + index = "pipeline/SVIM/{parameters}/min_{minscore}.indel.vcf.gz.tbi" + output: + summary="pipeline/SVIM_results/{parameters}/{minscore}/{vcf}/summary.txt" + params: + out_dir="pipeline/SVIM_results/{parameters}/{minscore}/{vcf}", + vcf="{vcf}" + threads: 1 + log: + log="logs/truvari/truvari.svim.{parameters}.{minscore}.{vcf}.log" + conda: + "../../../envs/truvari.yaml" + script: + "../scripts/run_truvari.py" + + +rule callset_eval: + input: + genome = config["reference"], + truth_vcf = get_vcf, + truth_bed = config["truth"]["bed"], + calls = "pipeline/{caller}/min_{minscore}.indel.vcf.gz", + index = "pipeline/{caller}/min_{minscore}.indel.vcf.gz.tbi" + output: + summary="pipeline/{caller,Sniffles|pbsv}_results/{minscore}/{vcf}/summary.txt" + params: + out_dir="pipeline/{caller}_results/{minscore}/{vcf}", + vcf="{vcf}" + threads: 1 + log: + log="logs/truvari/truvari.{caller}.{minscore}.{vcf}.log" + conda: + "../../../envs/truvari.yaml" + script: + "../scripts/run_truvari.py" + + +rule reformat_truvari_results: + input: + "pipeline/{caller}_results/{minscore}/{vcf}/summary.txt" + output: + "pipeline/{caller,Sniffles|pbsv}_results/{minscore}/{vcf}/pr_rec.txt" + threads: 1 + shell: + "cat {input} | grep 'precision\|recall' | tr -d ',' |sed 's/^[ \t]*//' | tr -d '\"' | tr -d ' ' | tr ':' '\t' | awk 'OFS=\"\\t\" {{ print \"{wildcards.caller}\", \"{wildcards.data}\", \"{wildcards.vcf}\", {wildcards.minscore}, $1, $2 }}' > {output}" + +rule reformat_truvari_results_svim: + input: + "pipeline/SVIM_results/{parameters}/{minscore}/{vcf}/summary.txt" + output: + "pipeline/SVIM_results/{parameters}/{minscore}/{vcf}/pr_rec.txt" + threads: 1 + shell: + "cat {input} | grep 'precision\|recall' | tr -d ',' |sed 's/^[ \t]*//' | tr -d '\"' | tr -d ' ' | tr ':' '\t' | awk 'OFS=\"\\t\" {{ print \"SVIM\", \"{wildcards.data}\", \"{wildcards.vcf}\", {wildcards.minscore}, $1, $2 }}' > {output}" + + +rule cat_truvari_results_all: + input: + svim = expand("pipeline/SVIM_results/1000_900_1.0_0.5/{minscore}/{vcf}/pr_rec.txt", + minscore=[0] + SVIM_THRESHOLDS, vcf=VCFS), + sniffles = expand("pipeline/Sniffles_results/{minscore}/{vcf}/pr_rec.txt", + minscore=list(range(config["minimums"]["sniffles_from"], config["minimums"]["sniffles_to"]+1, config["minimums"]["sniffles_step"])), + vcf=VCFS), + pbsv = expand("pipeline/pbsv_results/{minscore}/{vcf}/pr_rec.txt", + minscore=list(range(config["minimums"]["pbsv_from"], config["minimums"]["pbsv_to"]+1, config["minimums"]["pbsv_step"])), + vcf=VCFS) + output: + svim = temp("pipeline/eval/svim.all_results.txt"), + sniffles = temp("pipeline/eval/sniffles.all_results.txt"), + pbsv = temp("pipeline/eval/pbsv.all_results.txt"), + all = "pipeline/eval/all_results.txt" + threads: 1 + run: + shell("cat {input.svim} > {output.svim}") + shell("cat {input.sniffles} > {output.sniffles}") + shell("cat {input.pbsv} > {output.pbsv}") + shell("cat {output.svim} {output.sniffles} {output.pbsv} > {output.all}") + +rule cat_truvari_results_full: + input: + svim = expand("pipeline/SVIM_results/pooled/1000_900_1.0_0.5/{minscore}/{vcf}/pr_rec.txt", + minscore=[0] + SVIM_THRESHOLDS, vcf=VCFS), + sniffles = expand("pipeline/Sniffles_results/pooled/{minscore}/{vcf}/pr_rec.txt", + minscore=list(range(config["minimums"]["sniffles_from"], config["minimums"]["sniffles_to"]+1, config["minimums"]["sniffles_step"])), + vcf=VCFS), + pbsv = expand("pipeline/pbsv_results/pooled/{minscore}/{vcf}/pr_rec.txt", + minscore=list(range(config["minimums"]["pbsv_from"], config["minimums"]["pbsv_to"]+1, config["minimums"]["pbsv_step"])), + vcf=VCFS) + output: + svim = temp("pipeline/eval/svim.full_results.txt"), + sniffles = temp("pipeline/eval/sniffles.full_results.txt"), + pbsv = temp("pipeline/eval/pbsv.full_results.txt"), + all = "pipeline/eval/full_results.txt" + threads: 1 + run: + shell("cat {input.svim} > {output.svim}") + shell("cat {input.sniffles} > {output.sniffles}") + shell("cat {input.pbsv} > {output.pbsv}") + shell("cat {output.svim} {output.sniffles} {output.pbsv} > {output.all}") + +rule cat_truvari_results_svim_parameters: + input: + svim = expand("pipeline/SVIM_results/{pmd}_{pdn}_{edn}_{cmd}/{minscore}/{vcf}/pr_rec.txt", + pmd = [1000], + pdn = [900], + edn = [1.0, 1.5, 2.0, 2.5, 3.0, 4.0], + cmd = [0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5], + minscore= SVIM_THRESHOLDS, + vcf=VCFS) + output: + all = "pipeline/eval/svim_parameter_results.txt" + threads: 1 + run: + with open(output.all, 'w') as output_file: + for f in input.svim: + parameters = f.split("/")[4] + with open(f, 'r') as input_file: + for line in input_file: + print("%s\t%s" % (parameters, line.strip()), file=output_file) diff --git a/test/benchmark/caller_comparison/workflow/rules/plots.smk b/test/benchmark/caller_comparison/workflow/rules/plots.smk new file mode 100644 index 00000000..dce6f8f7 --- /dev/null +++ b/test/benchmark/caller_comparison/workflow/rules/plots.smk @@ -0,0 +1,192 @@ +localrules: plot_pr_all_results, plot_pr_tools, plot_pr_coverages, plot_pr_svim_parameters, run_pbsv1_all_types, run_pbsv2_all_types, SV_length_plot_pbsv, SV_length_plot_sniffles, SV_length_plot_svim, merge_counts, plot_counts + +rule plot_pr_all_results: + input: + "pipeline/eval/all_results.txt" + output: + "pipeline/eval/results.all.png" + threads: 1 + log: + "pipeline/logs/rplot.all.log" + shell: + "Rscript --vanilla workflow/scripts/plot_all_results.R {input} {output} > {log}" + +rule plot_pr_tools: + input: + "pipeline/eval/full_results.txt" + output: + "pipeline/eval/results.tools.{vcf}.png" + threads: 1 + log: + "pipeline/logs/rplot.tools.{vcf}.log" + shell: + "Rscript --vanilla workflow/scripts/plot_tools.R {input} {wildcards.vcf} {output} > {log}" + +rule plot_pr_coverages: + input: + "pipeline/eval/all_results.txt" + output: + png = "pipeline/eval/results.coverages.{vcf}.png", + tsv = "pipeline/eval/results.coverages.{vcf}.tsv" + threads: 1 + log: + "pipeline/logs/rplot.coverages.{vcf}.log" + shell: + "Rscript --vanilla workflow/scripts/plot_coverages.R {input} {wildcards.vcf} {output.png} {output.tsv} > {log}" + +rule plot_pr_svim_parameters: + input: + "pipeline/eval/svim_parameter_results.txt" + output: + "pipeline/eval/results.svim.parameters.png" + threads: 1 + log: + "pipeline/logs/rplot.all.log" + shell: + "Rscript --vanilla workflow/scripts/plot_svim_parameters.R {input} {output} > {log}" + +# rule plot_pr_coverages_bar: +# input: +# "pipeline/eval/all_results.txt" +# output: +# "pipeline/eval/results.coverages.bar.png" +# threads: 1 +# log: +# "pipeline/logs/rplot.coveragesbar.log" +# shell: +# "Rscript --vanilla workflow/scripts/plot_coverages_bar.R {input} {output} > {log}" + +#run pbsv for all SV types +rule run_pbsv1_all_types: + input: + bam = config["long_bam"], + bai = config["long_bai"], + genome = config["reference"], + output: + svsig = temp("pipeline/pbsv/svsig_all_types.svsig.gz") + resources: + mem_mb = 10000, + time_min = 600, + io_gb = 100 + threads: 1 + conda: + "../../../envs/pbsv.yaml" + shell: + """ + pbsv discover {input.bam} {output.svsig} + """ + +rule run_pbsv2_all_types: + input: + svsig = "pipeline/pbsv/svsig_all_types.svsig.gz", + genome = config["reference"], + output: + vcf = "pipeline/pbsv/min_{minsupport,[0-9]+}_all_types.vcf" + params: + min_sv_size = config["parameters"]["min_sv_size"] + resources: + mem_mb = 10000, + time_min = 600, + io_gb = 100 + threads: 1 + conda: + "../../../envs/pbsv.yaml" + shell: + """ + pbsv call -j 1 \ + --min-sv-length {params.min_sv_size} \ + --max-ins-length 100K \ + --call-min-reads-one-sample {wildcards.minsupport} \ + --call-min-reads-all-samples {wildcards.minsupport} \ + --call-min-reads-per-strand-all-samples 0 \ + --call-min-bnd-reads-all-samples 0 \ + --call-min-read-perc-one-sample 0 {input.genome} {input.svsig} {output.vcf} + """ + +rule SV_length_plot_pbsv: + input: + "pipeline/pbsv/min_{minimum}_all_types.vcf" + output: + plot = "pipeline/SV-plots/SV-length_pbsv_{minimum}.png", + counts = "pipeline/SV-plots/SV-counts_pbsv_{minimum}.txt", + log: + "logs/svplot/svlength_pbsv_{minimum}.log" + conda: + "../../../envs/cyvcf2.yaml" + shell: + "python workflow/scripts/SV-length-plot.py {input} --output {output.plot} --counts {output.counts} --filter 'hs37d5' --tool pbsv 2> {log}" + +rule SV_length_plot_sniffles: + input: + "pipeline/Sniffles/min_{minimum}.vcf" + output: + plot = "pipeline/SV-plots/SV-length_Sniffles_{minimum}.png", + counts = "pipeline/SV-plots/SV-counts_Sniffles_{minimum}.txt", + log: + "logs/svplot/svlength_Sniffles_{minimum}.log" + conda: + "../../../envs/cyvcf2.yaml" + shell: + "python workflow/scripts/SV-length-plot.py {input} --output {output.plot} --counts {output.counts} --filter 'hs37d5' --tool Sniffles 2> {log}" + +#run SVIM without converting DUP to INS +rule run_svim_all_types: + input: + genome = config["reference"], + bam = config["long_bam"], + bai = config["long_bai"] + output: + "pipeline/SVIM/{pmd,[0-9]+}_{pdn,[0-9]+}_{edn,[0-9\.]+}_{cmd,[0-9\.]+}_all_types/variants.vcf" + resources: + mem_mb = 10000, + time_min = 600, + io_gb = 100 + params: + working_dir = "pipeline/SVIM/{pmd}_{pdn}_{edn}_{cmd}_all_types/", + min_sv_size = config["parameters"]["min_sv_size"] + threads: 1 + shell: + "/home/heller_d/bin/anaconda3/envs/svim_test4/bin/python /home/heller_d/bin/anaconda3/envs/svim_test4/bin/svim alignment --sample {wildcards.data} \ + --partition_max_distance {wildcards.pmd} \ + --position_distance_normalizer {wildcards.pdn} \ + --edit_distance_normalizer {wildcards.edn} \ + --cluster_max_distance {wildcards.cmd} \ + --min_sv_size {params.min_sv_size} \ + --segment_gap_tolerance 20 \ + --segment_overlap_tolerance 20 \ + --read_names \ + --max_sv_size 1000000 \ + --verbose \ + {params.working_dir} {input.bam} {input.genome}" + +rule SV_length_plot_svim: + input: + "pipeline/SVIM/{parameters}_all_types/variants.vcf" + output: + plot = "pipeline/SV-plots/SV-length_SVIM_{parameters}_{minimum}.png", + counts = "pipeline/SV-plots/SV-counts_SVIM_{parameters}_{minimum}.txt", + log: + "logs/svplot/svlength_SVIM_{parameters}_{minimum}.log" + conda: + "../../../envs/cyvcf2.yaml" + shell: + "python workflow/scripts/SV-length-plot.py {input} --min_score {wildcards.minimum} --output {output.plot} --counts {output.counts} --filter 'hs37d5' --tool SVIM 2> {log}" + +rule merge_counts: + input: + svim = "pipeline/SV-plots/SV-counts_SVIM_1000_900_1.0_0.5_7.txt", + sniffles = "pipeline/SV-plots/SV-counts_Sniffles_5.txt", + pbsv = "pipeline/SV-plots/SV-counts_pbsv_5.txt", + output: + "pipeline/SV-plots/SV-counts.merged.txt" + shell: + "cat {input} | grep -v '#' > {output}" + +rule plot_counts: + input: + "pipeline/SV-plots/SV-counts.merged.txt" + output: + png = "pipeline/SV-plots/SV-counts.merged.png", + tsv = "pipeline/SV-plots/SV-counts.merged.tsv" + shell: + "Rscript --vanilla workflow/scripts/plot_counts.R {input} {output.png} {output.tsv}" diff --git a/test/benchmark/caller_comparison/workflow/scripts/plot_all_results.R b/test/benchmark/caller_comparison/workflow/scripts/plot_all_results.R new file mode 100644 index 00000000..a279e561 --- /dev/null +++ b/test/benchmark/caller_comparison/workflow/scripts/plot_all_results.R @@ -0,0 +1,27 @@ +library(tidyverse) +library(scales) + +args = commandArgs(trailingOnly=TRUE) + +res <- read_tsv(args[1], col_names = c("caller", "min_qual", "metric", "value")) +res$caller = factor(res$caller, levels=c('iGenVar', 'SVIM'), labels=c('iGenVar', 'SVIM')) +# res$caller = factor(res$caller, levels=c('pbsv', 'Sniffles', 'SVIM'), labels=c('pbsv', 'Sniffles', 'SVIM')) + +res %>% + filter(metric %in% c("recall", "precision")) %>% + pivot_wider(names_from=metric, values_from=value) %>% + filter(recall!=0 | precision!=0) %>% + mutate(precision = 100*precision, recall = 100*recall) %>% + ggplot(aes(recall, precision, color=caller, pch=caller)) + + geom_point(size=0.5) + + # scale_shape_manual(values=c(15,16,17)) + + # scale_color_manual(values=c("deepskyblue3", "goldenrod2", "firebrick2")) + + geom_path() + + # facet_wrap(~subsample+vcf) + + labs(y = "Precision", x = "Recall", color = "Tool", pch = "Tool") + + lims(x=c(0,100), y=c(0,100)) + + theme_bw() + + theme(panel.spacing = unit(0.75, "lines")) + + theme(text = element_text(size=14), axis.text.x = element_text(size=9), axis.text.y = element_text(size=9)) + +ggsave(args[2], width=20, height=12) diff --git a/test/benchmark/envs/cyvcf2.yaml b/test/benchmark/envs/cyvcf2.yaml new file mode 100644 index 00000000..ed48896a --- /dev/null +++ b/test/benchmark/envs/cyvcf2.yaml @@ -0,0 +1,8 @@ +channels: + - bioconda +dependencies: + - python>=3.6 + - pip + - pip: + - cyvcf2 + - matplotlib diff --git a/test/benchmark/shell_scripts/environment.yml b/test/benchmark/envs/environment.yml similarity index 100% rename from test/benchmark/shell_scripts/environment.yml rename to test/benchmark/envs/environment.yml diff --git a/test/benchmark/envs/pbsv.yaml b/test/benchmark/envs/pbsv.yaml new file mode 100644 index 00000000..3b5491c4 --- /dev/null +++ b/test/benchmark/envs/pbsv.yaml @@ -0,0 +1,4 @@ +channels: + - bioconda +dependencies: + - pbsv=2.3.0 diff --git a/test/benchmark/envs/samtools.yaml b/test/benchmark/envs/samtools.yaml new file mode 100644 index 00000000..1e13e14f --- /dev/null +++ b/test/benchmark/envs/samtools.yaml @@ -0,0 +1,5 @@ +channels: + - bioconda +#Sorting BAM files from NGMLR fails in samtools version>=1.10, therefore fix to version 1.9 +dependencies: + - samtools=1.9 diff --git a/test/benchmark/envs/sniffles.yaml b/test/benchmark/envs/sniffles.yaml new file mode 100644 index 00000000..01e2fdce --- /dev/null +++ b/test/benchmark/envs/sniffles.yaml @@ -0,0 +1,4 @@ +channels: + - bioconda +dependencies: + - sniffles=1.0.11 diff --git a/test/benchmark/envs/svim.yaml b/test/benchmark/envs/svim.yaml new file mode 100644 index 00000000..0d692987 --- /dev/null +++ b/test/benchmark/envs/svim.yaml @@ -0,0 +1,4 @@ +channels: + - bioconda +dependencies: + - svim=1.4.2 diff --git a/test/benchmark/envs/truvari.yaml b/test/benchmark/envs/truvari.yaml new file mode 100644 index 00000000..68bcaed8 --- /dev/null +++ b/test/benchmark/envs/truvari.yaml @@ -0,0 +1,7 @@ +channels: + - bioconda +dependencies: + - python>=3.6 + - pip + - pip: + - Truvari diff --git a/test/benchmark/envs/truvari_environment.yaml b/test/benchmark/envs/truvari_environment.yaml deleted file mode 100644 index 807f9d14..00000000 --- a/test/benchmark/envs/truvari_environment.yaml +++ /dev/null @@ -1,8 +0,0 @@ -channels: - - defaults - - conda-forge - - bioconda -dependencies: - - python=3 - - pip: - - truvari diff --git a/test/benchmark/shell_scripts/iGenVar_init.sh b/test/benchmark/iGenVar_init.sh similarity index 90% rename from test/benchmark/shell_scripts/iGenVar_init.sh rename to test/benchmark/iGenVar_init.sh index 9caff0ad..75f5f902 100755 --- a/test/benchmark/shell_scripts/iGenVar_init.sh +++ b/test/benchmark/iGenVar_init.sh @@ -46,6 +46,9 @@ mkdir -p long_reads && cd long_reads wget --retry-connrefused --waitretry=30 --read-timeout=30 --timeout=30 --tries=20 --no-clobber --no-verbose \ https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/data/AshkenazimTrio/HG002_NA24385_son/PacBio_CCS_10kb/alignment/HG002.Sequel.10kb.pbmm2.hs37d5.whatshap.haplotag.RTG.10x.trio.bam +wget --retry-connrefused --waitretry=30 --read-timeout=30 --timeout=30 --tries=20 --no-clobber --no-verbose \ + https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/data/AshkenazimTrio/HG002_NA24385_son/PacBio_CCS_10kb/alignment/HG002.Sequel.10kb.pbmm2.hs37d5.whatshap.haplotag.RTG.10x.trio.bam.bai + echo "$(tput setaf 1)$(tput setab 7)------- long reads downloaded (4/5) --------$(tput sgr 0)" 1>&3 diff --git a/test/benchmark/Snakefile b/test/benchmark/parameter_benchmarks/Snakefile similarity index 96% rename from test/benchmark/Snakefile rename to test/benchmark/parameter_benchmarks/Snakefile index 9fdf45fa..23a693e0 100644 --- a/test/benchmark/Snakefile +++ b/test/benchmark/parameter_benchmarks/Snakefile @@ -55,7 +55,8 @@ rule truvari: output_dir = "results/truvari/{parameter_name}/{parameter_value}_min_qual_{min_qual}" output: summary = "results/truvari/{parameter_name}/{parameter_value}_min_qual_{min_qual}/summary.txt" - # conda: "envs/truvari_environment.yaml" + conda: + "../envs/truvari.yaml" shell: """ rm -rf {params.output_dir} && truvari bench -b data/truth_set/HG002_SVs_Tier1_v0.6.vcf.gz \ @@ -110,4 +111,4 @@ rule plot_pr_all_results: output: "results/plots/{parameter_name}.results.all.png" shell: - "Rscript --vanilla Repos/iGenVar/test/benchmark/scripts/plot_all.R {input} {wildcards.parameter_name} {output}" + "Rscript --vanilla Repos/iGenVar/test/benchmark/parameter_benchmarks/plot_all.R {input} {wildcards.parameter_name} {output}" diff --git a/test/benchmark/shell_scripts/iGenVar_run_benchmark.sh b/test/benchmark/parameter_benchmarks/iGenVar_run_benchmark.sh similarity index 95% rename from test/benchmark/shell_scripts/iGenVar_run_benchmark.sh rename to test/benchmark/parameter_benchmarks/iGenVar_run_benchmark.sh index ac0a094a..eead516c 100755 --- a/test/benchmark/shell_scripts/iGenVar_run_benchmark.sh +++ b/test/benchmark/parameter_benchmarks/iGenVar_run_benchmark.sh @@ -38,7 +38,7 @@ else echo "Use one core." fi -/usr/bin/time -v snakemake --snakefile Repos/iGenVar/test/benchmark/Snakefile --cores $cores \ # --quiet --forcerun run_igenvar \ +/usr/bin/time -v snakemake --snakefile Repos/iGenVar/test/benchmark/parameter_benchmarks/Snakefile --cores $cores \ # --quiet --forcerun run_igenvar \ --rerun-incomplete --stats logs/${current_date}_snakemake_stats.txt --timestamp --use-conda echo "$(tput setaf 1)$(tput setab 7)------- done: iGenVar with cigar string & split read method --------$(tput sgr 0)" 1>&3 diff --git a/test/benchmark/scripts/plot_all.R b/test/benchmark/parameter_benchmarks/plot_all.R similarity index 100% rename from test/benchmark/scripts/plot_all.R rename to test/benchmark/parameter_benchmarks/plot_all.R From 26f53505c21410fbc7753b407503d2a46c481148 Mon Sep 17 00:00:00 2001 From: Lydia Buntrock Date: Fri, 20 Aug 2021 13:19:07 +0200 Subject: [PATCH 2/5] [BENCHMARK] Update files and add iGenVar to workflow Signed-off-by: Lydia Buntrock --- test/benchmark/caller_comparison/Snakefile | 28 +-- .../caller_comparison/config/config.yaml | 31 ++- .../workflow/rules/callers.smk | 151 ++++---------- .../caller_comparison/workflow/rules/eval.smk | 158 +++++---------- .../workflow/rules/plots.smk | 190 +----------------- test/benchmark/iGenVar_init.sh | 27 ++- 6 files changed, 138 insertions(+), 447 deletions(-) diff --git a/test/benchmark/caller_comparison/Snakefile b/test/benchmark/caller_comparison/Snakefile index 258ab9bb..bc2ce45d 100644 --- a/test/benchmark/caller_comparison/Snakefile +++ b/test/benchmark/caller_comparison/Snakefile @@ -1,12 +1,4 @@ -import os -import gzip -import sys -import shutil - -configfile: "config/config.yaml" - -VCFS=["giab", "giab.gt", "giab.seq", "giab.gt.seq"] -SVIM_THRESHOLDS = [1, 3, 4, 5, 6, 7, 8, 9, 10, 14, 18, 22, 26, 30, 40, 50, 60] +configfile: "Repos/iGenVar/test/benchmark/caller_comparison/config/config.yaml" include: "workflow/rules/callers.smk" include: "workflow/rules/eval.smk" @@ -16,13 +8,11 @@ include: "workflow/rules/plots.smk" rule all: input: - #SV lengths - expand("pipeline/SV-plots/SV-length_pbsv_{minscore}.png", minscore=[3, 5, 7]), - expand("pipeline/SV-plots/SV-length_Sniffles_{minscore}.png", minscore=[3, 5, 7]), - expand("pipeline/SV-plots/SV-length_SVIM_1000_900_1.0_0.5_{minscore}.png", minscore=[3, 5, 7]), - expand("pipeline/SV-plots/SV-counts.merged.png"), - #Evaluation - expand("pipeline/eval/results.all.png"), - expand("pipeline/eval/results.tools.{vcf}.png", vcf=VCFS), - expand("pipeline/eval/results.coverages.{vcf}.png", vcf=VCFS), - expand("pipeline/eval/results.svim.parameters.png"), + # SV calling + expand("results/caller_comparison/SVIM/variants.vcf"), + # SV lengths + # expand("results/caller_comparison/SV-plots/SV-length_SVIM_{minscore}.png", minscore=[3, 5, 7]), + # Evaluation + expand("results/caller_comparison/eval/results.all.png"), + # expand("pipeline/eval/results.tools.{vcf}.png", vcf=VCFS), + # expand("pipeline/eval/results.coverages.{vcf}.png", vcf=VCFS), diff --git a/test/benchmark/caller_comparison/config/config.yaml b/test/benchmark/caller_comparison/config/config.yaml index 26d3f0d3..850085a5 100644 --- a/test/benchmark/caller_comparison/config/config.yaml +++ b/test/benchmark/caller_comparison/config/config.yaml @@ -1,24 +1,23 @@ -long_bam: /data/long_reads/HG002.Sequel.10kb.pbmm2.hs37d5.whatshap.haplotag.RTG.10x.trio_sorted.bam -long_bai: /data/long_reads/HG002.Sequel.10kb.pbmm2.hs37d5.whatshap.haplotag.RTG.10x.trio_sorted.bam.bai +long_bam: data/long_reads/HG002.Sequel.10kb.pbmm2.hs37d5.whatshap.haplotag.RTG.10x.trio_sorted.bam +long_bai: data/long_reads/HG002.Sequel.10kb.pbmm2.hs37d5.whatshap.haplotag.RTG.10x.trio_sorted.bam.bai -reference: /path/to/database/genome.fa.gz - -truth: - giab: data/truth_set/HG002_SVs_Tier1_v0.6.vcf.gz - bed: data/truth_set/HG002_SVs_Tier1_v0.6.bed +reference: data/reference/hs37d5.fa.gz parameters: - pbmm_preset: CCS #SUBREAD (CLR), CCS (CCS) - survivor_distance: 500 - min_sv_size: 40 + sample: HG002 + min_var_length: 40 + max_var_length: 1000000 minimums: - sniffles_from: 1 - sniffles_to: 60 - sniffles_step: 2 + igenvar_from: 1 + igenvar_to: 100 + igenvar_step: 2 svim_from: 1 svim_to: 100 svim_step: 2 - pbsv_from: 10 - pbsv_to: 91 - pbsv_step: 10 + # sniffles_from: 1 + # sniffles_to: 60 + # sniffles_step: 2 + # pbsv_from: 10 + # pbsv_to: 91 + # pbsv_step: 10 diff --git a/test/benchmark/caller_comparison/workflow/rules/callers.smk b/test/benchmark/caller_comparison/workflow/rules/callers.smk index 33da51e4..cc13cbcf 100644 --- a/test/benchmark/caller_comparison/workflow/rules/callers.smk +++ b/test/benchmark/caller_comparison/workflow/rules/callers.smk @@ -1,4 +1,25 @@ -localrules: filter_svim, fix_sniffles, filter_insertions_and_deletions, filter_insertions_and_deletions_svim +rule run_igenvar: + input: + bam = config["long_bam"] + output: + vcf = "results/caller_comparison/iGenVar/variants.vcf" + params: + sample = config["parameters"]["sample"], + min_var_length = config["parameters"]["min_var_length"], + max_var_length = config["parameters"]["max_var_length"] + shell: + """ + ./build/iGenVar/bin/iGenVar -t 1 -j {input.bam} -o {output.vcf} \ + --vcf_sample_name {params.sample} \ + --method cigar_string \ + --method split_read \ + --min_var_length {params.min_var_length} \ + --max_var_length {params.max_var_length} \ + --min_qual 2 + """ + # Defaults: + # --clustering_methods hierarchical_clustering --refinement_methods no_refinement + # --max_tol_inserted_length 5 --max_overlap 10 --hierarchical_clustering_cutoff 10 # SVIM rule run_svim: @@ -7,119 +28,33 @@ rule run_svim: bai = config["long_bai"], genome = config["reference"] output: - "pipeline/SVIM/{pmd,[0-9]+}_{pdn,[0-9]+}_{edn,[0-9\.]+}_{cmd,[0-9\.]+}/variants.vcf" + "results/caller_comparison/SVIM/variants.vcf" resources: mem_mb = 20000, time_min = 600, io_gb = 100 params: - working_dir = "pipeline/SVIM/{pmd}_{pdn}_{edn}_{cmd}/", - min_sv_size = config["parameters"]["min_sv_size"] + working_dir = "results/caller_comparison/SVIM/", + sample = config["parameters"]["sample"], + min_var_length = config["parameters"]["min_var_length"], + max_var_length = config["parameters"]["max_var_length"] threads: 1 conda: "../../../envs/svim.yaml" shell: - "svim alignment --sample {wildcards.data} \ - --partition_max_distance {wildcards.pmd} \ - --position_distance_normalizer {wildcards.pdn} \ - --edit_distance_normalizer {wildcards.edn} \ - --cluster_max_distance {wildcards.cmd} \ - --min_sv_size {params.min_sv_size} \ - --segment_gap_tolerance 20 \ - --segment_overlap_tolerance 20 \ - --interspersed_duplications_as_insertions \ - --tandem_duplications_as_insertions \ - --read_names \ - --max_sv_size 1000000 \ - --verbose \ - {params.working_dir} {input.bam} {input.genome}" - -rule filter_svim: - input: - "pipeline/SVIM/{parameters}/variants.vcf" - output: - temp("pipeline/SVIM/{parameters}/min_{minscore,[0-9]+}.vcf") - threads: 1 - shell: - "bcftools view -e \"GT=='0/0'\" {input} | \ - awk 'OFS=\"\\t\" {{ if($1 ~ /^#/) {{ print $0 }} \ - else {{ if($6>={wildcards.minscore}) {{ print $1, $2, $3, $4, $5, $6, \"PASS\", $8, $9, $10 }} }} }}' > {output}" - -# SNIFFLES -rule run_sniffles: - input: - bam = config["long_bam"], - bai = config["long_bai"] - output: - expand("pipeline/Sniffles/raw_{minsupport}.vcf", - minsupport=list(range(config["minimums"]["sniffles_from"], config["minimums"]["sniffles_to"]+1, config["minimums"]["sniffles_step"]))) - resources: - mem_mb = 400000, - time_min = 1200, - io_gb = 100 - params: - min_sv_size = config["parameters"]["min_sv_size"], - tmpdir = "300", - sniffles_from = config["minimums"]["sniffles_from"], - sniffles_to = config["minimums"]["sniffles_to"], - sniffles_step = config["minimums"]["sniffles_step"], - outdir = "pipeline/Sniffles/" - threads: 30 - conda: - "../../../envs/sniffles.yaml" - shell: - "bash workflow/scripts/run_sniffles.sh {input.bam} {input.bai} {params.sniffles_from} {params.sniffles_to} {params.sniffles_step} {params.min_sv_size} {threads} {params.outdir}" - -#see https://github.com/spiralgenetics/truvari/issues/43 -rule fix_sniffles: - input: - "pipeline/Sniffles/raw_{support}.vcf" - output: - "pipeline/Sniffles/min_{support,[0-9]+}.vcf" - shell: - "sed 's/##INFO= {output}" - -#PBSV -rule run_pbsv: - input: - bam = config["long_bam"], - bai = config["long_bai"] - genome = config["reference"], - output: - expand("pipeline/pbsv/min_{minsupport}.vcf", - minsupport=list(range(config["minimums"]["pbsv_from"], config["minimums"]["pbsv_to"]+1, config["minimums"]["pbsv_step"]))) - resources: - mem_mb = 400000, - time_min = 2000, - io_gb = 100 - params: - min_sv_size = config["parameters"]["min_sv_size"], - tmpdir = "1000", - pbsv_from = config["minimums"]["pbsv_from"], - pbsv_to = config["minimums"]["pbsv_to"], - pbsv_step = config["minimums"]["pbsv_step"], - outdir = "pipeline/pbsv/" - threads: 15 - conda: - "../../../envs/pbsv.yaml" - shell: - "bash workflow/scripts/run_pbsv.sh {input.bam} {input.bai} {input.genome} {params.pbsv_from} {params.pbsv_to} {params.pbsv_step} {params.min_sv_size} {threads} {params.outdir}" - -#Split to SV classes -rule filter_insertions_and_deletions: - input: - "pipeline/{caller}/min_{minscore}.vcf" - output: - "pipeline/{caller,Sniffles|pbsv}/min_{minscore,[0-9]+}.indel.vcf" - threads: 1 - shell: - "bcftools view -i 'SVTYPE=\"DEL\" | SVTYPE=\"INS\"' {input} | bcftools sort > {output}" - -rule filter_insertions_and_deletions_svim: - input: - "pipeline/SVIM/{parameters}/min_{minscore}.vcf" - output: - "pipeline/SVIM/{parameters}/min_{minscore,[0-9]+}.indel.vcf" - threads: 1 - shell: - "bcftools view -i 'SVTYPE=\"DEL\" | SVTYPE=\"INS\"' {input} | bcftools sort > {output}" + """ + svim alignment --sample {params.sample} \ + --partition_max_distance 1000 \ + --cluster_max_distance 0.5 \ + --min_sv_size {params.min_var_length} \ + --segment_gap_tolerance 20 \ + --segment_overlap_tolerance 20 \ + --interspersed_duplications_as_insertions \ + --tandem_duplications_as_insertions \ + --read_names \ + --max_sv_size {params.max_var_length} \ + --verbose \ + {params.working_dir} {input.bam} {input.genome} + """ + # Defaults: + # --position_distance_normalizer 900 --edit_distance_normalizer 1.0 diff --git a/test/benchmark/caller_comparison/workflow/rules/eval.smk b/test/benchmark/caller_comparison/workflow/rules/eval.smk index c4af5b2a..d45a77d7 100644 --- a/test/benchmark/caller_comparison/workflow/rules/eval.smk +++ b/test/benchmark/caller_comparison/workflow/rules/eval.smk @@ -1,7 +1,10 @@ -localrules: bgzip, tabix, callset_eval_svim, callset_eval, reformat_truvari_results, reformat_truvari_results_svim, cat_truvari_results_all, cat_truvari_results_full, cat_truvari_results_svim_parameters - -def get_vcf(wildcards): - return config["truth"][wildcards.vcf.split(".")[0]] +rule filter_vcf: + input: + vcf = "results/caller_comparison/{caller}/variants.vcf" + output: + "results/caller_comparison/{caller}/variants.min_qual_{min_qual}.vcf" + shell: + "bcftools view -i 'QUAL>={wildcards.min_qual}' {input.vcf} > {output}" rule bgzip: input: @@ -11,7 +14,6 @@ rule bgzip: shell: "bgzip -c {input} > {output}" - rule tabix: input: "{name}.vcf.gz" @@ -20,128 +22,60 @@ rule tabix: shell: "tabix -p vcf {input}" - -rule callset_eval_svim: +rule truvari: input: - genome = config["reference"], - truth_vcf = get_vcf, - truth_bed = config["truth"]["bed"], - calls = "pipeline/SVIM/{parameters}/min_{minscore}.indel.vcf.gz", - index = "pipeline/SVIM/{parameters}/min_{minscore}.indel.vcf.gz.tbi" - output: - summary="pipeline/SVIM_results/{parameters}/{minscore}/{vcf}/summary.txt" + vcf = "results/caller_comparison/{caller}/variants.min_qual_{min_qual}.vcf.gz", + index = "results/caller_comparison/{caller}/variants.min_qual_{min_qual}.vcf.gz.tbi" params: - out_dir="pipeline/SVIM_results/{parameters}/{minscore}/{vcf}", - vcf="{vcf}" - threads: 1 - log: - log="logs/truvari/truvari.svim.{parameters}.{minscore}.{vcf}.log" - conda: - "../../../envs/truvari.yaml" - script: - "../scripts/run_truvari.py" - - -rule callset_eval: - input: - genome = config["reference"], - truth_vcf = get_vcf, - truth_bed = config["truth"]["bed"], - calls = "pipeline/{caller}/min_{minscore}.indel.vcf.gz", - index = "pipeline/{caller}/min_{minscore}.indel.vcf.gz.tbi" + output_dir = "results/caller_comparison/eval/{caller}/min_qual_{min_qual}" output: - summary="pipeline/{caller,Sniffles|pbsv}_results/{minscore}/{vcf}/summary.txt" - params: - out_dir="pipeline/{caller}_results/{minscore}/{vcf}", - vcf="{vcf}" - threads: 1 - log: - log="logs/truvari/truvari.{caller}.{minscore}.{vcf}.log" + summary = "results/caller_comparison/eval/{caller}/min_qual_{min_qual}/summary.txt" conda: "../../../envs/truvari.yaml" - script: - "../scripts/run_truvari.py" - - -rule reformat_truvari_results: - input: - "pipeline/{caller}_results/{minscore}/{vcf}/summary.txt" - output: - "pipeline/{caller,Sniffles|pbsv}_results/{minscore}/{vcf}/pr_rec.txt" - threads: 1 shell: - "cat {input} | grep 'precision\|recall' | tr -d ',' |sed 's/^[ \t]*//' | tr -d '\"' | tr -d ' ' | tr ':' '\t' | awk 'OFS=\"\\t\" {{ print \"{wildcards.caller}\", \"{wildcards.data}\", \"{wildcards.vcf}\", {wildcards.minscore}, $1, $2 }}' > {output}" + """ + rm -rf {params.output_dir} && truvari bench -b data/truth_set/HG002_SVs_Tier1_v0.6.vcf.gz \ + -c {input.vcf} -o {params.output_dir} --passonly --includebed data/truth_set/HG002_SVs_Tier1_v0.6.bed -p 0 + """ -rule reformat_truvari_results_svim: +rule reformat_truvari_results: input: - "pipeline/SVIM_results/{parameters}/{minscore}/{vcf}/summary.txt" + "results/caller_comparison/eval/{caller}/min_qual_{min_qual}/summary.txt" output: - "pipeline/SVIM_results/{parameters}/{minscore}/{vcf}/pr_rec.txt" + "results/caller_comparison/eval/{caller}/min_qual_{min_qual}/pr_rec.txt" threads: 1 shell: - "cat {input} | grep 'precision\|recall' | tr -d ',' |sed 's/^[ \t]*//' | tr -d '\"' | tr -d ' ' | tr ':' '\t' | awk 'OFS=\"\\t\" {{ print \"SVIM\", \"{wildcards.data}\", \"{wildcards.vcf}\", {wildcards.minscore}, $1, $2 }}' > {output}" - + """ + cat {input} | grep '\\|\' | tr -d ',' |sed 's/^[ \t]*//' | tr -d '\"' | tr -d ' ' \ + | tr ':' '\t' | awk 'OFS=\"\\t\" {{ print \"{wildcards.caller}\", \"{wildcards.min_qual}\", $1, $2 }}' > {output} + """ rule cat_truvari_results_all: input: - svim = expand("pipeline/SVIM_results/1000_900_1.0_0.5/{minscore}/{vcf}/pr_rec.txt", - minscore=[0] + SVIM_THRESHOLDS, vcf=VCFS), - sniffles = expand("pipeline/Sniffles_results/{minscore}/{vcf}/pr_rec.txt", - minscore=list(range(config["minimums"]["sniffles_from"], config["minimums"]["sniffles_to"]+1, config["minimums"]["sniffles_step"])), - vcf=VCFS), - pbsv = expand("pipeline/pbsv_results/{minscore}/{vcf}/pr_rec.txt", - minscore=list(range(config["minimums"]["pbsv_from"], config["minimums"]["pbsv_to"]+1, config["minimums"]["pbsv_step"])), - vcf=VCFS) + igenvar = expand("results/caller_comparison/eval/iGenVar/min_qual_{min_qual}/pr_rec.txt", + min_qual=list(range(config["minimums"]["igenvar_from"], + config["minimums"]["igenvar_to"]+1, + config["minimums"]["igenvar_step"]))), + svim = expand("results/caller_comparison/eval/SVIM/min_qual_{min_qual}/pr_rec.txt", + min_qual=list(range(config["minimums"]["svim_from"], + config["minimums"]["svim_to"]+1, + config["minimums"]["svim_step"]))), + # sniffles = expand("results/caller_comparison/eval/Sniffles/min_qual_{min_qual}/pr_rec.txt", + # min_qual=list(range(config["minimums"]["sniffles_from"], config["minimums"]["sniffles_to"]+1, config["minimums"]["sniffles_step"]))), + # pbsv = expand("results/caller_comparison/eval/pbsv/min_qual_{min_qual}/pr_rec.txt", + # min_qual=list(range(config["minimums"]["pbsv_from"], config["minimums"]["pbsv_to"]+1, config["minimums"]["pbsv_step"]))) output: - svim = temp("pipeline/eval/svim.all_results.txt"), - sniffles = temp("pipeline/eval/sniffles.all_results.txt"), - pbsv = temp("pipeline/eval/pbsv.all_results.txt"), - all = "pipeline/eval/all_results.txt" + igenvar = temp("results/caller_comparison/eval/igenvar.all_results.txt"), + # svim = temp("results/caller_comparison/eval/svim.all_results.txt"), + svim = temp("results/caller_comparison/eval/svim.all_results.txt"), + # sniffles = temp("results/caller_comparison/eval/sniffles.all_results.txt"), + # pbsv = temp("results/caller_comparison/eval/pbsv.all_results.txt"), + all = "results/caller_comparison/eval/all_results.txt" threads: 1 run: + shell("cat {input.igenvar} > {output.igenvar}") shell("cat {input.svim} > {output.svim}") - shell("cat {input.sniffles} > {output.sniffles}") - shell("cat {input.pbsv} > {output.pbsv}") - shell("cat {output.svim} {output.sniffles} {output.pbsv} > {output.all}") - -rule cat_truvari_results_full: - input: - svim = expand("pipeline/SVIM_results/pooled/1000_900_1.0_0.5/{minscore}/{vcf}/pr_rec.txt", - minscore=[0] + SVIM_THRESHOLDS, vcf=VCFS), - sniffles = expand("pipeline/Sniffles_results/pooled/{minscore}/{vcf}/pr_rec.txt", - minscore=list(range(config["minimums"]["sniffles_from"], config["minimums"]["sniffles_to"]+1, config["minimums"]["sniffles_step"])), - vcf=VCFS), - pbsv = expand("pipeline/pbsv_results/pooled/{minscore}/{vcf}/pr_rec.txt", - minscore=list(range(config["minimums"]["pbsv_from"], config["minimums"]["pbsv_to"]+1, config["minimums"]["pbsv_step"])), - vcf=VCFS) - output: - svim = temp("pipeline/eval/svim.full_results.txt"), - sniffles = temp("pipeline/eval/sniffles.full_results.txt"), - pbsv = temp("pipeline/eval/pbsv.full_results.txt"), - all = "pipeline/eval/full_results.txt" - threads: 1 - run: - shell("cat {input.svim} > {output.svim}") - shell("cat {input.sniffles} > {output.sniffles}") - shell("cat {input.pbsv} > {output.pbsv}") - shell("cat {output.svim} {output.sniffles} {output.pbsv} > {output.all}") - -rule cat_truvari_results_svim_parameters: - input: - svim = expand("pipeline/SVIM_results/{pmd}_{pdn}_{edn}_{cmd}/{minscore}/{vcf}/pr_rec.txt", - pmd = [1000], - pdn = [900], - edn = [1.0, 1.5, 2.0, 2.5, 3.0, 4.0], - cmd = [0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5], - minscore= SVIM_THRESHOLDS, - vcf=VCFS) - output: - all = "pipeline/eval/svim_parameter_results.txt" - threads: 1 - run: - with open(output.all, 'w') as output_file: - for f in input.svim: - parameters = f.split("/")[4] - with open(f, 'r') as input_file: - for line in input_file: - print("%s\t%s" % (parameters, line.strip()), file=output_file) + # shell("cat {input.sniffles} > {output.sniffles}") + # shell("cat {input.pbsv} > {output.pbsv}") + shell("cat {output.igenvar} {output.svim} > {output.all}") + # shell("cat {output.igenvar} {output.svim} {output.sniffles} {output.pbsv} > {output.all}") diff --git a/test/benchmark/caller_comparison/workflow/rules/plots.smk b/test/benchmark/caller_comparison/workflow/rules/plots.smk index dce6f8f7..2f951bb6 100644 --- a/test/benchmark/caller_comparison/workflow/rules/plots.smk +++ b/test/benchmark/caller_comparison/workflow/rules/plots.smk @@ -1,192 +1,12 @@ -localrules: plot_pr_all_results, plot_pr_tools, plot_pr_coverages, plot_pr_svim_parameters, run_pbsv1_all_types, run_pbsv2_all_types, SV_length_plot_pbsv, SV_length_plot_sniffles, SV_length_plot_svim, merge_counts, plot_counts - rule plot_pr_all_results: input: - "pipeline/eval/all_results.txt" + "results/caller_comparison/eval/all_results.txt" output: - "pipeline/eval/results.all.png" - threads: 1 + "results/caller_comparison/eval/results.all.png" log: - "pipeline/logs/rplot.all.log" + "logs/rplot.all.log" shell: - "Rscript --vanilla workflow/scripts/plot_all_results.R {input} {output} > {log}" - -rule plot_pr_tools: - input: - "pipeline/eval/full_results.txt" - output: - "pipeline/eval/results.tools.{vcf}.png" - threads: 1 - log: - "pipeline/logs/rplot.tools.{vcf}.log" - shell: - "Rscript --vanilla workflow/scripts/plot_tools.R {input} {wildcards.vcf} {output} > {log}" - -rule plot_pr_coverages: - input: - "pipeline/eval/all_results.txt" - output: - png = "pipeline/eval/results.coverages.{vcf}.png", - tsv = "pipeline/eval/results.coverages.{vcf}.tsv" - threads: 1 - log: - "pipeline/logs/rplot.coverages.{vcf}.log" - shell: - "Rscript --vanilla workflow/scripts/plot_coverages.R {input} {wildcards.vcf} {output.png} {output.tsv} > {log}" - -rule plot_pr_svim_parameters: - input: - "pipeline/eval/svim_parameter_results.txt" - output: - "pipeline/eval/results.svim.parameters.png" - threads: 1 - log: - "pipeline/logs/rplot.all.log" - shell: - "Rscript --vanilla workflow/scripts/plot_svim_parameters.R {input} {output} > {log}" - -# rule plot_pr_coverages_bar: -# input: -# "pipeline/eval/all_results.txt" -# output: -# "pipeline/eval/results.coverages.bar.png" -# threads: 1 -# log: -# "pipeline/logs/rplot.coveragesbar.log" -# shell: -# "Rscript --vanilla workflow/scripts/plot_coverages_bar.R {input} {output} > {log}" - -#run pbsv for all SV types -rule run_pbsv1_all_types: - input: - bam = config["long_bam"], - bai = config["long_bai"], - genome = config["reference"], - output: - svsig = temp("pipeline/pbsv/svsig_all_types.svsig.gz") - resources: - mem_mb = 10000, - time_min = 600, - io_gb = 100 - threads: 1 - conda: - "../../../envs/pbsv.yaml" - shell: - """ - pbsv discover {input.bam} {output.svsig} """ - -rule run_pbsv2_all_types: - input: - svsig = "pipeline/pbsv/svsig_all_types.svsig.gz", - genome = config["reference"], - output: - vcf = "pipeline/pbsv/min_{minsupport,[0-9]+}_all_types.vcf" - params: - min_sv_size = config["parameters"]["min_sv_size"] - resources: - mem_mb = 10000, - time_min = 600, - io_gb = 100 - threads: 1 - conda: - "../../../envs/pbsv.yaml" - shell: + Rscript --vanilla Repos/iGenVar/test/benchmark/caller_comparison/workflow/scripts/plot_all_results.R \ + {input} {output} > {log} """ - pbsv call -j 1 \ - --min-sv-length {params.min_sv_size} \ - --max-ins-length 100K \ - --call-min-reads-one-sample {wildcards.minsupport} \ - --call-min-reads-all-samples {wildcards.minsupport} \ - --call-min-reads-per-strand-all-samples 0 \ - --call-min-bnd-reads-all-samples 0 \ - --call-min-read-perc-one-sample 0 {input.genome} {input.svsig} {output.vcf} - """ - -rule SV_length_plot_pbsv: - input: - "pipeline/pbsv/min_{minimum}_all_types.vcf" - output: - plot = "pipeline/SV-plots/SV-length_pbsv_{minimum}.png", - counts = "pipeline/SV-plots/SV-counts_pbsv_{minimum}.txt", - log: - "logs/svplot/svlength_pbsv_{minimum}.log" - conda: - "../../../envs/cyvcf2.yaml" - shell: - "python workflow/scripts/SV-length-plot.py {input} --output {output.plot} --counts {output.counts} --filter 'hs37d5' --tool pbsv 2> {log}" - -rule SV_length_plot_sniffles: - input: - "pipeline/Sniffles/min_{minimum}.vcf" - output: - plot = "pipeline/SV-plots/SV-length_Sniffles_{minimum}.png", - counts = "pipeline/SV-plots/SV-counts_Sniffles_{minimum}.txt", - log: - "logs/svplot/svlength_Sniffles_{minimum}.log" - conda: - "../../../envs/cyvcf2.yaml" - shell: - "python workflow/scripts/SV-length-plot.py {input} --output {output.plot} --counts {output.counts} --filter 'hs37d5' --tool Sniffles 2> {log}" - -#run SVIM without converting DUP to INS -rule run_svim_all_types: - input: - genome = config["reference"], - bam = config["long_bam"], - bai = config["long_bai"] - output: - "pipeline/SVIM/{pmd,[0-9]+}_{pdn,[0-9]+}_{edn,[0-9\.]+}_{cmd,[0-9\.]+}_all_types/variants.vcf" - resources: - mem_mb = 10000, - time_min = 600, - io_gb = 100 - params: - working_dir = "pipeline/SVIM/{pmd}_{pdn}_{edn}_{cmd}_all_types/", - min_sv_size = config["parameters"]["min_sv_size"] - threads: 1 - shell: - "/home/heller_d/bin/anaconda3/envs/svim_test4/bin/python /home/heller_d/bin/anaconda3/envs/svim_test4/bin/svim alignment --sample {wildcards.data} \ - --partition_max_distance {wildcards.pmd} \ - --position_distance_normalizer {wildcards.pdn} \ - --edit_distance_normalizer {wildcards.edn} \ - --cluster_max_distance {wildcards.cmd} \ - --min_sv_size {params.min_sv_size} \ - --segment_gap_tolerance 20 \ - --segment_overlap_tolerance 20 \ - --read_names \ - --max_sv_size 1000000 \ - --verbose \ - {params.working_dir} {input.bam} {input.genome}" - -rule SV_length_plot_svim: - input: - "pipeline/SVIM/{parameters}_all_types/variants.vcf" - output: - plot = "pipeline/SV-plots/SV-length_SVIM_{parameters}_{minimum}.png", - counts = "pipeline/SV-plots/SV-counts_SVIM_{parameters}_{minimum}.txt", - log: - "logs/svplot/svlength_SVIM_{parameters}_{minimum}.log" - conda: - "../../../envs/cyvcf2.yaml" - shell: - "python workflow/scripts/SV-length-plot.py {input} --min_score {wildcards.minimum} --output {output.plot} --counts {output.counts} --filter 'hs37d5' --tool SVIM 2> {log}" - -rule merge_counts: - input: - svim = "pipeline/SV-plots/SV-counts_SVIM_1000_900_1.0_0.5_7.txt", - sniffles = "pipeline/SV-plots/SV-counts_Sniffles_5.txt", - pbsv = "pipeline/SV-plots/SV-counts_pbsv_5.txt", - output: - "pipeline/SV-plots/SV-counts.merged.txt" - shell: - "cat {input} | grep -v '#' > {output}" - -rule plot_counts: - input: - "pipeline/SV-plots/SV-counts.merged.txt" - output: - png = "pipeline/SV-plots/SV-counts.merged.png", - tsv = "pipeline/SV-plots/SV-counts.merged.tsv" - shell: - "Rscript --vanilla workflow/scripts/plot_counts.R {input} {output.png} {output.tsv}" diff --git a/test/benchmark/iGenVar_init.sh b/test/benchmark/iGenVar_init.sh index 75f5f902..1065b1f7 100755 --- a/test/benchmark/iGenVar_init.sh +++ b/test/benchmark/iGenVar_init.sh @@ -20,7 +20,7 @@ git clone https://github.com/seqan/iGenVar.git cd iGenVar git submodule update --recursive --init -echo "$(tput setaf 1)$(tput setab 7)------- iGenVar downloaded (1/5) --------$(tput sgr 0)" 1>&3 +echo "$(tput setaf 1)$(tput setab 7)------- iGenVar downloaded (1/6) --------$(tput sgr 0)" 1>&3 cd ../.. mkdir -p build && cd build && mkdir -p iGenVar && cd iGenVar @@ -30,7 +30,7 @@ make -j 16 make test && make doc -echo "$(tput setaf 1)$(tput setab 7)------- iGenVar built (2/5) --------$(tput sgr 0)" 1>&3 +echo "$(tput setaf 1)$(tput setab 7)------- iGenVar built (2/6) --------$(tput sgr 0)" 1>&3 # -------- -------- get data and unzip -------- -------- # cd ../.. @@ -40,7 +40,7 @@ mkdir -p data && cd data # wget ... # cd .. -echo "$(tput setaf 1)$(tput setab 7)------- short reads downloaded (3/5) --------$(tput sgr 0)" 1>&3 +echo "$(tput setaf 1)$(tput setab 7)------- short reads downloaded (3/6) --------$(tput sgr 0)" 1>&3 mkdir -p long_reads && cd long_reads wget --retry-connrefused --waitretry=30 --read-timeout=30 --timeout=30 --tries=20 --no-clobber --no-verbose \ @@ -48,9 +48,16 @@ wget --retry-connrefused --waitretry=30 --read-timeout=30 --timeout=30 --tries=2 wget --retry-connrefused --waitretry=30 --read-timeout=30 --timeout=30 --tries=20 --no-clobber --no-verbose \ https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/data/AshkenazimTrio/HG002_NA24385_son/PacBio_CCS_10kb/alignment/HG002.Sequel.10kb.pbmm2.hs37d5.whatshap.haplotag.RTG.10x.trio.bam.bai +cd .. +echo "$(tput setaf 1)$(tput setab 7)------- long reads downloaded (4/6) --------$(tput sgr 0)" 1>&3 -echo "$(tput setaf 1)$(tput setab 7)------- long reads downloaded (4/5) --------$(tput sgr 0)" 1>&3 +# -------- -------- get reference ../../data -------- -------- # +mkdir -p reference && cd reference +wget --retry-connrefused --waitretry=30 --read-timeout=30 --timeout=30 --tries=20 --no-clobber --no-verbose \ + ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/technical/reference/phase2_reference_assembly_sequence/hs37d5.fa.gz +cd .. +echo "$(tput setaf 1)$(tput setab 7)------- reference downloaded (5/6) --------$(tput sgr 0)" 1>&3 # -------- -------- get truth set ../../data -------- -------- # mkdir -p truth_set && cd truth_set @@ -61,12 +68,18 @@ wget --retry-connrefused --waitretry=30 --read-timeout=30 --timeout=30 --tries=2 https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/../../data/AshkenazimTrio/analysis/NIST_SVs_Integration_v0.6/HG002_SVs_Tier1_v0.6.vcf.gz.tbi wget --retry-connrefused --waitretry=30 --read-timeout=30 --timeout=30 --tries=20 --no-clobber --no-verbose \ https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/../../data/AshkenazimTrio/analysis/NIST_SVs_Integration_v0.6/HG002_SVs_Tier1_v0.6.bed +cd ../.. +echo "$(tput setaf 1)$(tput setab 7)------- truth set downloaded (7/6) --------$(tput sgr 0)" 1>&3 -echo "$(tput setaf 1)$(tput setab 7)------- truth set downloaded (5/5) --------$(tput sgr 0)" 1>&3 +mkdir -p results +# -------- -------- pre installation steps -------- -------- # +# We do our benchmarks with snakemake, and need several tools like +# python3, samtools, snakemake, tabix, pip and truvari +# we recommend to use miniconda: https://docs.conda.io/en/latest/miniconda.html +# and run -cd ../.. -mkdir -p results +conda env create -f Repos/iGenVar/test/benchmark/envs/environment.yml # ---------------------------------------- echo "$(tput setaf 1)$(tput setab 7)------- Initial steps - done --------$(tput sgr 0)" 1>&3 From fa1d374790909c98c59716781f5e24e98e629a86 Mon Sep 17 00:00:00 2001 From: Lydia Buntrock Date: Mon, 23 Aug 2021 15:38:11 +0200 Subject: [PATCH 3/5] [BENCHMARK] Add sniffles to plot (loop over quals with sniffles) Signed-off-by: Lydia Buntrock --- doc/CallerComparisonPlot.pdf | Bin 0 -> 38242 bytes .../caller_comparison/config/config.yaml | 8 ++- .../workflow/rules/callers.smk | 53 +++++++++++++++++- .../caller_comparison/workflow/rules/eval.smk | 39 ++++++------- .../workflow/scripts/plot_all_results.R | 3 +- test/benchmark/iGenVar_init.sh | 19 +++++-- 6 files changed, 91 insertions(+), 31 deletions(-) create mode 100644 doc/CallerComparisonPlot.pdf diff --git a/doc/CallerComparisonPlot.pdf b/doc/CallerComparisonPlot.pdf new file mode 100644 index 0000000000000000000000000000000000000000..e54eebb521c44143a9765a52b1c1564b78f09795 GIT binary patch literal 38242 zcmeIb2Y6J~+V@RIqI3{YkfDP>lIatA??sAqkjZ2MNJt=s-jt@aLlF^?CJHDZARq#Q z6lv0nAW8@6MLGiF`}?m+Lb054zUO-(-NMl!SaUVE*3-Q|CmwRU!~*J{)%YYlb; z6?^N{Tc_XJ{nq|=u|ZZ#s3j)+&7cYuf@;X{=P=IVF}H^_^~c2|#QNyZ@l1beki?D#{^(Beoh>$x&1MOy)g>yP z8HCh|f;lz(zL*HV=#GnrXWfE&f3&mB_Oba_}j*)*uAbp z-VS+3rmKIV&O6JtUOV?uzioxP9(d{Yt&&w<+&_Ev+7bWS{-EfV;L5>0`>c5Q;MVS# zTQ`r+oiV8Ql5g(Ee>vY8Ii+&tHGPk-9oQ%C?4_xb2c7%l>@St)9v}F|jvI?6wjcQ2 z?$r~W*;)-Ow=LwicZxUow$j=T9d4yBv;K?dw*6ZzZQK8=t|jZ&8CbUPol7?!buHb0 z+#iuq>lRFn+x~D{<`10rm!yu1Nn7qKr+-Y(m1*O+Vg9+UwgJ) zQ!q{DtZ$uI{HX3sZ%`{=i{tLzhrB;7&MnpQqgGcLGI8*~1~%-I zr*rq|r@HT0btdm2TdA5Se;nWNp0jkr>z^!p@1VncUT$GGpC28{R{YaabE=e?U-rT9 zTst3Rif-6*UW1bhY?H$dePJ86^!qYp%4VKa>ZA0ok`eja4w;xLDqD%DDia$|-t^%c zN2afRaBRbtnB!FrttgmhOOKN^?`NIiozkjzZ2miUqWiR%5whZ1ja#=4{k&*jvkql* zrFS>D@Y1T*WqV%mzI|v<^+F-(mY(gM{%G#KJ+@xzx#@EB_O6=?ZCO0LfBFKix;<4s zKk?G4QS+aFY0CGfFZ7JNRI6{BD?QU!ob-b=_*AMdT75Id(|ul+w@$w_e`%$Cx6?j! zx9@TF;QH5Hv&s)je>5WO>oV3;dv7(YwsOJI3e($%mbx`P?dcz8AD-E#$mQ5!>o*?C z*W=5Qe_W4>%9~;GqNRtves^Wz0;@|m%74AhviA@69(J(sooO>RWVu=@s@(UtqwcO- zJLW)z=cX?%wzS;b+gY2v?R@TwTBV;mSY`1`*Pi$Fuh_@)dt}B(IctuuoxfDc)YZbP ztj(Ek$jzYexo2vux;S=H&Z2b-Kde`(dcIt(3$0rl|H`)+$}f&xaDK>Jm4mh{i9Pkr zk0rxet$w+cXK(D=gWuR(&v~lV>>i`XY&%*c@8L|(m*~B1cebNJVeLGdhn1xsME#QINbAifrWbgaD^r&(8>${2{I;dafU2H%W;Phmy>v+aYpus$Y2|s**5I|V zKd0I9@$`K|f0y)24go3{CNC=WqY%H5INMBuiUGPn>5&yL<(4GM(*(D3Z_cfI|$Z0OP9?ow;R!9(`F{QSAR!zX>1 z(8%6v;|A-?V+u!nHey}X-#X1c)3Qk4>LbSOUw0srW*U1)lXcj}*H+blAWV)k0(2wJY|o3ZKpGH2cbP6^af`UH8PG z)>GE+YrM7P#L)OrKQ35MwNS1s1!8Bd$+W)nd%N1DedubqID3uyV+LocS%1ih-#evy zadXZ>o728k?(4?gYrkDS$B|CyE_{%oQ2x7?V*8q9x?aa~b^Uj{Gu^4Pxb}uA6(7cA z`gPof7pBf!+hgO2*Yg}}{nGJi5oITq`yy?>SK1c1k-2M&o2MGoX;Aaa5>3+;S$L}4 z;&k)Moq73SnHBH;a&qz)X)iSR^4`INdqZt=znHzP*8*SOe3R?8uzy;s|E7eF*H_c> zeQUc3FYMS+GT%!_uSb^uB5lULwdi$w-dgrFo5xHTFlO?-xnF0=?`(7BMvodHk&Pl> z89lZB_66&Ew#(RP)gSFEWgg(U*1UETdj|Xbh;u)F8MmeN!idScH@=zw;MvzV+m@AI zaCG?COY4TsTA1ax;m-MQHt)D#+>S}^;d6&y8r!V=_TU+FmM(8ndf=7j)(l&}+0vw4 z#ye%kT^bkDYUvkOY*R0z_-mDyzDS#ARR@T#zd3W8=0y*->@~OSl0i|6uP*K!c3{}{d1J>lTH5Gj z(P<;gyIYmt(e&UC!AE9Av^+HV{Kt`pTIBT9bJR(bd2Oa%C8v~qvq+JHaa)em&-z-c zUU?5(Sk*o6q|Ebchomd@$&{#uBN9gJY`v=ckg`oSR~VW9;Eh=uS5=?y@15tX*{^n* z@aB~U2fpdzqfO12YOA{s*|59miG)g>e$4Rg#$7QRYOj5-`jF%OmTta#XI-9?MUKx+ zdwjV5|Mt^g(;hGQ+u5z|+H3##zH{9diD8uHtY7wfci;5j9b1b}Kas6$o1-(jwC=fo$JV~>y06}K zHssylI{ro#OVxYtW^AWgt7Ga>3XW)dk(((%11+EcD2iJYvW1( z;NyKa9)0QenHQpBPWigV_8T&N-o#bz;~smLV|OcYfMwK-=18`}AY)LZ^4=!&M2b+Qsk5`^s|fgm#Bk&u%f{*0C)&-kf!HTEqGo*TmF)bJCj9 zOWurMbFOo~#l3Q!?(RRcEyviNg_l2?ocZ+jp8HCl7045kDQg}m4DUpZ0nXi zXhF}L6IYLjyH>LD?alA5C|m1d)h3xQRh_-O{fk9xjc)Cj+WY6x#ot_5qW{5Gw_eNp z`rf|Vw#`}G&RTesWz)~o4xC+jWo>AfGq%#X8fzNvDVOo~up?=fWh{O$ORb+9SwAjY zuXx(`=6vZ{v*@s6xBs)T>AqiIJmJ*Ly>(;Sb}kaTvc!E5^5h_l6b9`q}FZS1{+^| zdjtEeYNwVbuHSi{;gr}f$5^(D0pn>IroV>zn`qzdtjRA;{`hmJ(uUi zqu(D|D$GCeZg9?D>SgY;`@xo#FT|CYaJFNm4gQbv-=l_^cqDoB3O82sjS+FHX!s(l54nAshc53P1nY(I- zdWUr`>UnhM@}01Pag#e0taRYqw;?MkA3Od-%U&~!6z{Wl*^O3nW{fD%VEfurA!Xj4 zygk#FJ133|xi)28{1^56dP1*@HSF%t zyqy!?xzprWj%*9}R+>37p+TN9^X5&=mGRNdrL*1`_*tQq4YEcRyzxoKO)bx-KJA?~ z`TWTvHAC*Mw7O=6+CLt#q1G2CSDo(oLEVeN$FJwP75l}UwGBc%89T?Fui2wwv$+*F zx4%-cSMD7b+dbOh_tnYPYrunzm+pRB;f+1LD?T`Vy7h?P&n}&PqyGIo#ogaNyfSx0 z;p6V9`|95gX|{IIgR>0>{C;8I;C^?`-FoBmAL9lP#3#nXiQ&AW^*}V86wU{t_>lTt zBH}DJU=1cYN(;@k(-NMi9^NQ8BU2 zx_f=*Oh{bQ+SM%1q~?#44V}{1ZA0o-uU^d?=Z~7a#-pnOD?6$SK~mB?zvumzi!0B ze8<1(b8Y9;;|1>y>^fu3?xUTq4LGX0^(7^4z9(#*FIJF=Of% zeod7&D%FYZOP5xA@m{SoQ_}SRadgy_ocHpqS(U}|+k;p4S024S)gLq7TiLp>HPve$ zj6Qec-r;un@@(@|Z23{K6^*xj@Rg(CuF!hB3b%gS)?#AZjWTWf-s)=h<@>U5`5!%Dp3rV(yR9vGU{L#qQzzW* z`F@9?gRbVQQSjczJ1qxJ$((BO9RKpe_eXT=k~>xVOQ~O5{riZW9Ya!${BmtVmH|Cp zS=GNo)9+4|=~R0Dkt6Nur=HQZRLJ>j_h)1oRQ-+qyPK{(dhHi)@0$}n6XLDaJGi4G z`_8HH-PT>Fy8rgV*%>4570mc3;guJ&KPnO7m%+bi_tm=W zU#H4?B5j(gCErQcqoei2+quT2u3^cOYw+@!dEa|`VaF2lG9=8*?ik|go;Ut2-wB5& zWAIzsJEl95x@l;Z>4QuC8nmNIwxPqeX3O1lM7*W?*8E=%Yc=lWGilRT&;DBe_ns?O zy;_|#xvN~7*>qn1&xc=}6?wkx`I7N3txA=-`ihxNI_J2byKlxKp_%5aZ?iaiRQjbI z8+)D`J2cPwZNJ()1%uySv2IfLwEeSP9ol%~=^DTXO|Jbl$`8JRC&RNoGLAitZFOAH+DRyyql=Z8ie%bov z=$F3y*a1U|oQm)MmFJCi6&k%$`n%UpwCP@=d&}-Qzs~nfj$_%bjn2C1=C9kY`|s82 zpLO7KmB-E~R5Iu2OkWJGQ|(o6*;y0z)biB-x_MHeC(L&t{T`4iG@X^9+EEBiYipUyUp>=~})uM~fD0IBUw@qfg+UnJG zC4bG;qV$21@vpXP(tFN3Egj8oHE*}zgGH@o^{&yQ&h^i7bl6ts-G$#ZY4d6G+40RU z4BxlR_WI#kVa>y;_gNY{>PGrsYUe7m>77&6Pc|7@@wBf|-CLRNm-)Tu?*;C!x!){9xeSFeq|dN0!-~w;hnLxyYh(N2zYkwD zX7rfw(XPA+`6}94*t3N0az{JYcCF^>;cgwduVd`aT)yvZ4edv}^*#K8YpJ`wb-wjw z+lZY7cBDVMdgq9;dtPkxVzp5jM!no1$Nq^^Hb(i!w%J#9NZIGht}Zj9b-q@c+ICxd zamD@4-D0LiAMHF~%kc1=ott&u@4GpF;gqt6zS?}}r*=2m9bY;BT+iJDqB})zS@E;4 z*Z%2m77YJ({n|s1h76fKWb3S*-hF%5S@$kEcs;|Y%|X6yk<)rrE8izew>h4Y*4t6NKD9fTlQJwD8yHhIWyPtQj{;uDi5vOwg((Fpy z@s_*Yhgba2>-*iuVul?rHE7wOaeuh}csX^O)JF!N9sF69E$b$|HK6JzGt*c7p=#E* zcfMV7eZ@lE>vY_DXWZ(Wna@3cuGg@)_YB^1W!Bvdqc^@kY(T-WnX*2=DAS<4-keRc z?frJ}s*c;T%>8FxxmEWV~<+R5WT=+V4R_Bojz*p}OG+D0th zv2^6pPqu$F`AR!iv7j7IfAzzozFfV@ad`KE!Qa*RvV-lz)9*K%c(b^>+}%3UW2cA4 zbXwO>ACOkzrN@5Pv`e-d1+B_n~G~U$LxG- z`OxKEx(+J2EVxWa|CPa!?HiA8vG>D|m)stBsqLuTBWrGQyS|Q^_wzeHmHg53b(gDi zK5tmCe(z#=>MyKcy58}av(ansP1@OV(i;PQny}2jc*m@$-7R*sIC=2W&qr1iDA_W2 zcSymXU-@!PTt;8*ZDS{wpM3e7*Ssy8tQfVp{+IRJ0n3fIu zEsVS}=chxf`>bAc`TCFhTg-jEWt+Rvi#lH0+4tM~yE4okK61bZ1HvA>;;QAoGic}i zRV(Yx{%qNWjRQ{(yztgLX$R~ayK`sR56k)&cAx)A!i|m_)^9)Fd(i%)o1^B(%v~}0 zl{6!Yj(BCz)O%NMZ>&0^-@g`*T{fW8z4qNQE}#8L&Wr7ih5ayTU$HA=FTGOfrQXX{ zZ25HV>3-X<-7UDkVvD$Tzw|zwX~Vl8=bX{#?2L~pO-i_U;li2yzBTpxUi;|Xig61= zzKQuYYR}=D=lbuScH!>j%6*RAez@fMU3tcKh-iO#)tT;}#jP7S^K`}s1^^JMo{`kkftNv^8Ki>MiWxhp`9VUO%WKi6RvX{F3vHF8vIqfU2ExA^@@&^OY z|CaaG^{aD7H=Z#t?cJd_AN(oo4oeP0gA%LgHKQq5cZsHquCVp-QOhpTRiop&B%KHf zsnI1e(ogAJwehwVhcl#Gm*|8zOPDL9M?y@zKgu5&A7phnLLy?Q_Q%C_v4nJr_4Xvj z_4*PN>iW9G`VzWDM)`YNLgKqbMR2N{*B2WTZ3zjF^{XLwdVM~s@j)RGT^K5^OPnPn zI65Xi!XIe~i8Vy}$O59Gyic6(oDki~8=KHA%9{{x35n?x6YcNHXkOEEoFe!3Pb?#t zusB*&x-m(73piqm|78IKNpT*^L+MSk2OY%7zX^j4} zL?@rKy2!JhK4-Uj$VHwxXLB*+Gw1AKf4`A}4%ZX+#>aZQMEPS`gR+D)>(bW`4=f>dqa$O;2u=FdG$tnAqO~Ezuef+`Y`ozX z+zAURR;*^@T0xH&f|7LJECjLm^yEsAQoW~^f`b0E78DdxH5zSC%Fxu&H`0b&N2nYW;oerzh>Ima23!@yl-FCYh8p0f2 zhsWu4`4ElY;dh6+Y<8O+9p${+9wwo?BAhOl*X3}~pVO=J4wpO38D_V+JhXS{I>#=D z-D$Jyemf)Zk012N73m1Kdl=j6!04gC@AkUfj1USHekeda_~e5huhSzvff|Q14F0%K z96RIEzl$#E4+c8I+#rn02h}K)JJb>GWE8iRGfr2y%jWWPpWS6+_8QmCxZ#d)huzJ8 z4!?uW%z3BN<_L3g=zs=|>h`$Y%+JlN>>8EXKnfJw6GvraXf>m1R?y&fTOHxdT~zrv zj&#!=Lc`o@8|I3%TOEFvFOd%Thf%>F(Z@*4*zIM`P{1+sRa+Mn+r!Lt&0l=w6H1^o zk!~oXol_J-jolvV^dbjjikzTIpKi`cK3wy-J;+1%doV;ddO_##2>mefvb!Sa4FlIP z-|W!rg(G|qjHhw_VLTV(_~1Lzae3gE6D;N{`lNpSiQ{>Y2W(;yr9(F zFQ7;o*JJpA{vln5^wJd}(K#K^AI3brPV5Pt3zMEp_I4yrKgQtTf#GMOKGN3ct>5sO z_VC@qeO|_sR=H>=zDTnY`4i{|TImP16U}N*Pv|#xCHiS^ll(*<@)VyU6Zy@t=5O>6 zJ=eSf^MOw61bt2%huKBYui?M+BXAuFnRz6SZ`HBr6y49tgZ=~dZgY!H>^dL^^a`hC z2Wv#1k%w3O(s97fk%>8me*rzQA^!k<+{fKMCtxjoGak%9CYyT@=)1?j6teB>B9pwoq3B;rY;en`&3gGkW{PM|;HgY4XH zXMFTO8J{>$Kj=M{?BlKo_(1=%1CFsbqtEh0emhuX=80VXTl+J}AXV19G+wCDAK@c> z5&m)!ej9qFKcK$!D3SD#Lnagh=*+9PQpH+Mxs7Q zIPejU1pJQ2@Kbhb@WSYh6yV%yf%2CXb21^7V=0hWw+QIR12@{PGS`#37+6r0`ZaLjfeBY7v_Z|gz_Hf zl{~?DgHuWUNK*c!{)D56?eP1|TX+cW!+C7M^oK%qM;b z(?q^}qOWSNKA=)~Dy|Yi$u2>dfZUOSaFYHEzg!+iBmhY4lsKND-|d0l+>Z{2;cpD1 z67e;hd4k;726W)>!lbuw2CR5|UUX_+gj2{Kn{X!a)!;8Gs`G|ke`3G#|KN}Kjh(P@ zkA4H~A~ilb1U(WMl#0(L?xDTpA$%9E%f6HSYZAW&C=^WPs^hO|eTfm}-<%z>ZzkD@qykM`< z6Z9Yg`oj#pXr^=>$^!TZzq$U8b{Mb|V|SXL@DG6p>`wDx{@}H6QsE{38$5H_h>87* z`_L2S0$QPmfQ=g~bTdzAgkC|7jW+U%(o?V)Tv1dO=-6 zdPx|?vCx-KpRvQo`l?)ycnYcED-*9Sc)UWV^PckMv3OGkw z6SKf=`6Vu4%S8NW9r%vPkXzyb1AI5}ARfi|Q8b46Di0wCk=*5}kRE!?B)N~7AXXu} znGIL*XZS(Q(!|L$wsNHN(ir%{-Pji*K%*7ifToHQmE)pK%t|da4)K=m;U7KnTi<1G zz#nlBnWu@>0#|eujfGxjg{?{MiSjmA0@wfkH#UIp)L4O-EYMPRD>Kj>+=iyV@8SRc z73I2q#oy?4(rC;!d43Vf7vLoLA}lgLke_rNJR?`bj+7xM-37iSU1yf^S0LjP*NO3- zxGuS2y~G~@DIq5ly99dUer)ua>wrYkSf>AHu2Wk`{^2_IiXDT`NoSRFBz_lO>$>nd z`MU2heds3j5aWj)Uq^z;_j}-*;(V}*zD%#CmXSP`@E9!++6CfsE3pqYl=uUETqSM< zK@xj~j^zG{SB0yJE2$qP_Afh3x<8=z7#h7N-)~Q{4OzYT_m9`%|36(f<71hR<;&QK zc_0A}bes7r-jMC#p^1M~4*@UaPlYeUMzUY@2%9i|8@jQnMEnl;fkeM4J50u1NtWm> zATMHR@fm;eq=b_=q4@c+eQLab-71eUbeOn`d_@%_kR|E9zy72q0vbHV6F^%aW5yOe z&_KV!O0pb-Wq;kz|K?fhcL6&}(v)N^vD(qWfL1>-Z)_hFC+mi*~qWoV*(D5e7&^<%W2Eb#B>!viK^27V)PBz!a$D0@~t?Jx6X zxRfN%FziEFrtHW{xJq%G6$QTVk8oA~24wOoZ$suuqb1JinVGBY|CN599q0emW7UVA z#`h;@r&?Z;|1jqcD?${|>i+UNUzhSK+ScQa%%q7xno>zn55# zC5FI98aXG{Z;+_sHQ}x747!aU_@ChNvpoHqW97F2em&tU)t0h6KqODBxCVQdo~h~u z_f226aAkSYHL5mfvAN%5Z%=q_)k%`$%O~|xQ99mQ{y7oD%t&O3pdz7=(+wXQ4^_ApSssq4p$szfidX*mo(1^5EXF!g`a?kSZ=`)6QgLVq+RRe@y2A>tN zs}4ydZAPIYjGq_YD=VksMEgKrI85$TEAW^rCaOfc)sy;DPGL1UAMqE*9(2PBk7#e= z7ZW8CGbj`Eaoz9?-fCn~?^P=@U5?FfMHyC*hxSO2T9_yYVHD|$9}$WgiBLF(8_M`~ ztdW(Yh;kh(aw0~7@=%WTDGqWL%?{<*l>hNA%ALtHiJiFaRy)NwHj`;MrC*Avn7R7b zVik2G_@@;%j)^T~7i37tgJX@#l?aaIuT-U>@&W&q$-*hkk4Rg&8!Hk}!+1!U{>WB1 z7CrEk>#7?&Jdy`FbR?WptQe>wNWTM7xN?1}&Wvg|^s3@ZjSqsPCvX<@W}Odt;sd$M z3!GDh(EZ9XiS6j0ydKQf%vsUnH+NHnbc=qf?Ho(SN~$>q(*rSUz@8#xAJ7PMXp3Yd z1-MI&sijEph4$nfM6X7d{o*GH6S*li?WU^9_$sv^SYczR&g$b!a;7zS*jh(Pb8>@Mteo+Ojj#TrbbCSFq1p@ z)cwj)lQaW)gTv5A{H--S(d$-K-tC48GzA;OreG$wsCJ=at+A+zgytYMa#`9_U*Wu& z4@CtuNAwZX)4!1%6+o^_e{Ixku>(m0qJ3~v`VY@(uNs=>7D?s}zmX1_qpXfdo7OJs z0UY~jP34UKsIO{;0V-KPhMsVRTX0R8pfuTspHijHWERpyT_-XSa={;$U-vVrlNAr} zk9AdL9nvKHAOtg8VFGl9lX=1u_{%NGhqY&lAHowf6?%N+C(1%-PyP}f@CFD_(TFgE zXb??BYCwYYDMDJqf5H@U2y_C?Q+~%Wc>oxLei(FDra*g*2cM~$fWfMQO0$J}symr# zFu@YCGw6wDmB9$UI95IwDT?5Rmt*N)B4j-thi=S|V>ujr_ z-mOJtB!@CEp9m}ie%XX{SSRge>p~NwD=45)gSi^4A;O|B<#6=HHCaI5pQ$5(Us~M| zc8DL8w`h;naIDylHE?1(q|7~3Ryl^#W=-4!*BK9*Rby5E)UIS%Iu7J&x^7m$$gJpr z+P%Cf{lN*Hfgd;uRr$b6q9&&25{eidQq>7_&EdOx#K zwoB~~y9?!*^*gR>HjD*5aEAUObU!MpRY5FMyB8EK!87?o4>lw-RV-$ZOW!zFbsZhx zy6gZgko74GN6%&H>^~5kaO~llhc+NOx+%+-&(MAbG+)}w9n6d;!2?r$U^_pVkD$B zAL@`C3-7$xqjmzRm_QfU3DhV;q*ia%deD5WJrh?WV~|GN(^{RW;YUjEXs<}hfRTJ6 z{KsD5v#O%uSP5&7a2%?kDO|P=U$l0j$OIciM4{@zRN7R#WloCCI7@{epYut&oA7N~(YjaR3GQLRw-|j0Ri9Z*8aoH}J?(1T9*i=U6k>WdLD< zQjiTUpmbi4PP`Z9iY3GjY}69cq3%dX*q}%QY@$EewQP@LgYM`a$BH+ogP^zAJQ~Jo zys2?&Eg7JJJ75o}hIEBgpb)tqfMoWnFswpytx{15g9cfPYA?v&RNk2#b0>e{SSu{V zxb)&@WJT&iMj5MLo`o?O6Ro4J&1uC9x-P2|`bhLzhxV~f$<7<`5o{!`fM1c4ijgVR zO8CNDSj)h-ie}^lPoPg8Nt)mTX@!m+c!(cGPeb9c6HaO+mTs^MJO{lK!+_TuLo+dn z_88I1Lk*CNi|Jl1JIQOyvzq~6JX3XX*(=nVH6 zCBnjp0-#ozs-g;Xl-c9;pjOC_))OIxQuRT)%qM{KkgVkxtbsLHqlb#f;3ZljDTay^ z$ih$RNR=0)qPz=#6_r7Yzji*1F#546?&TNl399n4CGT=&~?@g zawm=8GY5O}+>agen&{ zBack1D?~&aROP1%WiSTqU?zAL&{TH|{T1Ph2G|RK1?u>9{YEbp6M{pM1ksRCgS`*j zua#E(4)SA-7d@m~j?pF)Gb(shtg4Iy4<@@{eu|2<#wG0(La|Oy|40(b6j_1|0E6^M zHs*sDq1wj;Kj1ywLvf*)by(5D8DbvjVO^YKc75hsA5j#G^$}}{uX1X5G$_Nf!BJvbOdH%oHh4`Pfh6cG zOi&+0tHdOV$h6)7|KOzbOvf_)9|UhAq!Y|VQ4rix8!SgY1Nt-p)QcJ#P>xTL zA7lU|COkosXbGZ3ii);5mK7sehgbN4fE7(~4Z38f0vKf>%E-VAGh?t6D&(6~A(aMc zl>wu}-atZTqu5=H)Gy{N-ISnFMadFPktjZjTKR9TiH zd+nuEZUP@c5wmuQB|{B5r;J3#Lo9-(YA233WLDx0M)};K9AnAopwVPTCn9GSL=emo z>iJQm3L2tO(o<%OrGO*w5sYP)imTv`+Mo#n4N)VXuV^2FjRs?<^dVkzkEFtwk!YE4 zUlB49zx=jlMEs%bO*srUfKB0Rz-mQU=n%8g+CNA1AZyWTgHVE5Vi)M1A`3-P>O~rd zPAQTNr7d2W-zH9B)>>~S^yexOx5)tIfrXe@EVL4*V0mZ z1Wea;G>hwG`@~j6-rx?d8B1m72;2zt%zjtoAymY2;Sg(y*dH}abVM~c6PaK!#MIa| zGZ%?af=19cT%a2t92IsWG-K%zI2MDxPO?h%p*fkz1zw{k+=o{+5iOC6;$`ipkd3h$ z3yBK>sA39Z(Me@CWHpRWG)r5|lG$S$@JLY%w55NNQCNM=iJ&85JVR(iv0Kd=+ zE<-ofL=z=p+a~g%Pr4+^F;NM52ht>pf<0UZGl-tiVroGskhaxB2Q-4n)F9cH!LBt~ zlHw<>Lx;(1m@HKj^pcAVt{n{8TY&dAepzTHuuwgo32q?27(DR%EvwA1n4jgxvE?mqG^U9Dwn__RqKME7@x2aZNtixVZd#! zvvz}xYnKn+UD#;qpnQie?S<10X+`CzlTM)f%pFep;1NS2OJodvXae*Jll%seq#@MS z5UL`xNMt9!EWF3&aC=AyGa-h9C(0|y3m8HAsl9rLU4B>67wpTzu_z!D(eMeDdYM(Y ztP$$LBttdafGV(v)>H$C+Q4gYgSrdYqIgm@60Xpte&dn!$OF?5N+<5K!B&l5hl$$n z22at2z&vRV-r< zA9BD%kqx6MW|A&I9}t~9E-@wkRI;HzG)I10t8MaAcrJn@vKl{`jAtA|3q9m zFX_W1=nm}PMAq<9*{141rmle=s&*%ruRRtDvX$#9O)J>UMAp ziv>N{%bYl#sury4K)5W0+7wAqG)9&_ z8dj-a(QoAaq&^z@uuIX0)goW*d!RaK^qoqXxes=V9w0DzeJe0NEH&dRl7I~KL~=qy zw6jk7B7MTniLq35p#Q(0XOcc4M*IqBFD`Ju=JW69ko}7fy46H1$PO7JyMG!Fc^SE& zNB_hh)fde^E*PNwax^ymVsD^|says0399AyrCW%==!^J#Rp%XC^PESF>J+MR6SjB_|O*_)vj{ zFj83>IgjWej)I4L(|Qbh1C`;xQ`K2Xl`j@3ru^O;tTahzpnr%!kt4Ix4mj0n znJLHAHSv|0JH07iQQZZ5;h2b$<4E{M#f;-{9Wz%Q6XDT^>gCF4<$E{=CBS(iDB>uN zS?@u|(NPdrwLt2H+zq8#<=~hZDkS9o0O#Z*w8|k)igQ$6#7TOf>dCR_MgGudY(@dH z1jzkb)qs--kNkjR?Hi^kmVUd`V|`X4;C zX&!dn#Wg@8l^n7nex=8>U>n#h?+1c<4uj1SJh2kpKB>kZE8IkPl%^{ z;^{47R?zV;qA?Ra5_OqfX+RxO7;u#G)TW7wlcPdq&ZhZ;!2E& ziLMl7D0Wb6$$5$rDNk)GB7Vy5rMHRMF4Y?JOSMar9TAD>X*XqWDNk*xPD1SsJyLBR z)({N_ek+dWJocON)Mm<4o9x6PuS|TJQ~0iG1(}HUE8!i;5du3T^mL^%KGk!zijeZu zCh3yWDAg0m4U?Yc3+%T{d1~{i>b8k1P$yHZB+doDRh?B8(^P4x&8VtGEkJc)@f95m zRLND%nZKgc6nex=8_P8leOL=OO7?y2%DNk*tJhhqf)Mm<4o8~!;l&3Z+52rk}nex=8 z{*Fb;Q=2JIZKgc6nexMf*tU8mkcTJ<$vxS9Ap@inJ5y;4)x0! z^+Hlx(hEr|R7iRiY|=YU|6@KRz0wsv*ppu9n#c!x(wkZhAM8mlQ~lrJ1CK&`cpjaH zqP@JNHPY?pPm@GMhI%}HtB1F|JMfBQ6T(%V*P^yDMJM%sbrdH>TFx7w0k zr}{s6ajPxq(c^z2+K41k)1M^Th~~c|+Q=_(%m4GY&6;-+yDiqlmo)2d9Xx*FF)zin zg(h9$om7vn1m1_Mzr^$BbDlrFmYKgo^yhYt&_B0xIMCRqUbn0l1V2HMxoERlL!W&u zGv}V9Jh7ckFG7B*KX$u{)lZ*uhlZhW&zxgL`_JdXo_)Qv%=uZGIp_ZK7|uV!X3%i~v0^U7tta@uYO%$|7XG&uM8rPI7cnwOXBW!zzjPs+z7zk^zT LOUat}66*g08k@L$ literal 0 HcmV?d00001 diff --git a/test/benchmark/caller_comparison/config/config.yaml b/test/benchmark/caller_comparison/config/config.yaml index 850085a5..3234e18f 100644 --- a/test/benchmark/caller_comparison/config/config.yaml +++ b/test/benchmark/caller_comparison/config/config.yaml @@ -1,4 +1,5 @@ long_bam: data/long_reads/HG002.Sequel.10kb.pbmm2.hs37d5.whatshap.haplotag.RTG.10x.trio_sorted.bam +long_bam_md: data/long_reads/HG002.Sequel.10kb.pbmm2.hs37d5.whatshap.haplotag.RTG.10x.trio_sorted.md.bam long_bai: data/long_reads/HG002.Sequel.10kb.pbmm2.hs37d5.whatshap.haplotag.RTG.10x.trio_sorted.bam.bai reference: data/reference/hs37d5.fa.gz @@ -7,6 +8,7 @@ parameters: sample: HG002 min_var_length: 40 max_var_length: 1000000 + min_qual: 2 minimums: igenvar_from: 1 @@ -15,9 +17,9 @@ minimums: svim_from: 1 svim_to: 100 svim_step: 2 - # sniffles_from: 1 - # sniffles_to: 60 - # sniffles_step: 2 + sniffles_from: 1 + sniffles_to: 60 + sniffles_step: 2 # pbsv_from: 10 # pbsv_to: 91 # pbsv_step: 10 diff --git a/test/benchmark/caller_comparison/workflow/rules/callers.smk b/test/benchmark/caller_comparison/workflow/rules/callers.smk index cc13cbcf..da65169d 100644 --- a/test/benchmark/caller_comparison/workflow/rules/callers.smk +++ b/test/benchmark/caller_comparison/workflow/rules/callers.smk @@ -5,6 +5,7 @@ rule run_igenvar: vcf = "results/caller_comparison/iGenVar/variants.vcf" params: sample = config["parameters"]["sample"], + min_qual = config["parameters"]["min_qual"], min_var_length = config["parameters"]["min_var_length"], max_var_length = config["parameters"]["max_var_length"] shell: @@ -15,7 +16,7 @@ rule run_igenvar: --method split_read \ --min_var_length {params.min_var_length} \ --max_var_length {params.max_var_length} \ - --min_qual 2 + --min_qual {params.min_qual} """ # Defaults: # --clustering_methods hierarchical_clustering --refinement_methods no_refinement @@ -58,3 +59,53 @@ rule run_svim: """ # Defaults: # --position_distance_normalizer 900 --edit_distance_normalizer 1.0 + +# SNIFFLES (we have to loop over min_support, because sniffles does not write a quality score into the vcf) +rule run_sniffles: + input: + bam = config["long_bam_md"], + output: + expand("results/caller_comparison/Sniffles/raw_variants.{minsupport}.vcf", + minsupport=list(range(config["minimums"]["sniffles_from"], + config["minimums"]["sniffles_to"]+1, + config["minimums"]["sniffles_step"]))) + resources: + mem_mb = 400000, + time_min = 1200, + io_gb = 100 + params: + min_support = config["parameters"]["min_qual"], + min_length = config["parameters"]["min_var_length"], + qual_from = config["minimums"]["sniffles_from"], + qual_to = config["minimums"]["sniffles_to"]+1, + qual_step = config["minimums"]["sniffles_step"] + threads: 10 + conda: + "../../../envs/sniffles.yaml" + shell: + """ + for i in $(seq {params.qual_from} {params.qual_step} {params.qual_to}) + do + sniffles --mapped_reads {input.bam} --vcf results/caller_comparison/Sniffles/raw_variants.$i.vcf \ + --min_support $i --min_length {params.min_length} --threads {threads} --genotype + done + """ + +#see https://github.com/spiralgenetics/truvari/issues/43 +rule fix_sniffles: + input: + "results/caller_comparison/Sniffles/raw_variants.{support}.vcf" + output: + "results/caller_comparison/Sniffles/variants.unsorted.min_qual_{support,[0-9]+}.vcf" + shell: + "sed 's/##INFO= {output}" + +# Split to SV classes +# Since iGenVar can only find INS and DEL so far, we filter these out for better comparability. +rule fix_sniffles_2_and_filter_insertions_and_deletions: + input: + "results/caller_comparison/Sniffles/variants.unsorted.min_qual_{support,[0-9]+}.vcf" + output: + "results/caller_comparison/Sniffles/variants.min_qual_{support,[0-9]+}.vcf" + shell: + "bcftools view -i 'SVTYPE=\"DEL\" | SVTYPE=\"INS\"' {input} | bcftools sort > {output}" diff --git a/test/benchmark/caller_comparison/workflow/rules/eval.smk b/test/benchmark/caller_comparison/workflow/rules/eval.smk index d45a77d7..f0f0bba6 100644 --- a/test/benchmark/caller_comparison/workflow/rules/eval.smk +++ b/test/benchmark/caller_comparison/workflow/rules/eval.smk @@ -1,10 +1,10 @@ rule filter_vcf: input: - vcf = "results/caller_comparison/{caller}/variants.vcf" + "results/caller_comparison/{caller,iGenVar|SVIM}/variants.vcf" output: - "results/caller_comparison/{caller}/variants.min_qual_{min_qual}.vcf" + "results/caller_comparison/{caller,iGenVar|SVIM}/variants.min_qual_{min_qual}.vcf" shell: - "bcftools view -i 'QUAL>={wildcards.min_qual}' {input.vcf} > {output}" + "bcftools view -i 'QUAL>={wildcards.min_qual}' {input} > {output}" rule bgzip: input: @@ -52,30 +52,31 @@ rule reformat_truvari_results: rule cat_truvari_results_all: input: - igenvar = expand("results/caller_comparison/eval/iGenVar/min_qual_{min_qual}/pr_rec.txt", - min_qual=list(range(config["minimums"]["igenvar_from"], - config["minimums"]["igenvar_to"]+1, - config["minimums"]["igenvar_step"]))), - svim = expand("results/caller_comparison/eval/SVIM/min_qual_{min_qual}/pr_rec.txt", - min_qual=list(range(config["minimums"]["svim_from"], - config["minimums"]["svim_to"]+1, - config["minimums"]["svim_step"]))), - # sniffles = expand("results/caller_comparison/eval/Sniffles/min_qual_{min_qual}/pr_rec.txt", - # min_qual=list(range(config["minimums"]["sniffles_from"], config["minimums"]["sniffles_to"]+1, config["minimums"]["sniffles_step"]))), + igenvar = expand("results/caller_comparison/eval/iGenVar/min_qual_{min_qual}/pr_rec.txt", + min_qual=list(range(config["minimums"]["igenvar_from"], + config["minimums"]["igenvar_to"]+1, + config["minimums"]["igenvar_step"]))), + svim = expand("results/caller_comparison/eval/SVIM/min_qual_{min_qual}/pr_rec.txt", + min_qual=list(range(config["minimums"]["svim_from"], + config["minimums"]["svim_to"]+1, + config["minimums"]["svim_step"]))), + sniffles = expand("results/caller_comparison/eval/Sniffles/min_qual_{min_qual}/pr_rec.txt", + min_qual=list(range(config["minimums"]["sniffles_from"], + config["minimums"]["sniffles_to"]+1, + config["minimums"]["sniffles_step"]))), # pbsv = expand("results/caller_comparison/eval/pbsv/min_qual_{min_qual}/pr_rec.txt", # min_qual=list(range(config["minimums"]["pbsv_from"], config["minimums"]["pbsv_to"]+1, config["minimums"]["pbsv_step"]))) output: - igenvar = temp("results/caller_comparison/eval/igenvar.all_results.txt"), - # svim = temp("results/caller_comparison/eval/svim.all_results.txt"), - svim = temp("results/caller_comparison/eval/svim.all_results.txt"), - # sniffles = temp("results/caller_comparison/eval/sniffles.all_results.txt"), + igenvar = temp("results/caller_comparison/eval/igenvar.all_results.txt"), + svim = temp("results/caller_comparison/eval/svim.all_results.txt"), + sniffles = temp("results/caller_comparison/eval/sniffles.all_results.txt"), # pbsv = temp("results/caller_comparison/eval/pbsv.all_results.txt"), all = "results/caller_comparison/eval/all_results.txt" threads: 1 run: shell("cat {input.igenvar} > {output.igenvar}") shell("cat {input.svim} > {output.svim}") - # shell("cat {input.sniffles} > {output.sniffles}") + shell("cat {input.sniffles} > {output.sniffles}") # shell("cat {input.pbsv} > {output.pbsv}") - shell("cat {output.igenvar} {output.svim} > {output.all}") + shell("cat {output.igenvar} {output.svim} {output.sniffles} > {output.all}") # shell("cat {output.igenvar} {output.svim} {output.sniffles} {output.pbsv} > {output.all}") diff --git a/test/benchmark/caller_comparison/workflow/scripts/plot_all_results.R b/test/benchmark/caller_comparison/workflow/scripts/plot_all_results.R index a279e561..fdf397fa 100644 --- a/test/benchmark/caller_comparison/workflow/scripts/plot_all_results.R +++ b/test/benchmark/caller_comparison/workflow/scripts/plot_all_results.R @@ -4,7 +4,7 @@ library(scales) args = commandArgs(trailingOnly=TRUE) res <- read_tsv(args[1], col_names = c("caller", "min_qual", "metric", "value")) -res$caller = factor(res$caller, levels=c('iGenVar', 'SVIM'), labels=c('iGenVar', 'SVIM')) +res$caller = factor(res$caller, levels=c('iGenVar', 'SVIM', 'Sniffles'), labels=c('iGenVar', 'SVIM', 'Sniffles')) # res$caller = factor(res$caller, levels=c('pbsv', 'Sniffles', 'SVIM'), labels=c('pbsv', 'Sniffles', 'SVIM')) res %>% @@ -17,7 +17,6 @@ res %>% # scale_shape_manual(values=c(15,16,17)) + # scale_color_manual(values=c("deepskyblue3", "goldenrod2", "firebrick2")) + geom_path() + - # facet_wrap(~subsample+vcf) + labs(y = "Precision", x = "Recall", color = "Tool", pch = "Tool") + lims(x=c(0,100), y=c(0,100)) + theme_bw() + diff --git a/test/benchmark/iGenVar_init.sh b/test/benchmark/iGenVar_init.sh index 1065b1f7..f3d28928 100755 --- a/test/benchmark/iGenVar_init.sh +++ b/test/benchmark/iGenVar_init.sh @@ -20,7 +20,7 @@ git clone https://github.com/seqan/iGenVar.git cd iGenVar git submodule update --recursive --init -echo "$(tput setaf 1)$(tput setab 7)------- iGenVar downloaded (1/6) --------$(tput sgr 0)" 1>&3 +echo "$(tput setaf 1)$(tput setab 7)------- iGenVar downloaded (1/7) --------$(tput sgr 0)" 1>&3 cd ../.. mkdir -p build && cd build && mkdir -p iGenVar && cd iGenVar @@ -30,7 +30,7 @@ make -j 16 make test && make doc -echo "$(tput setaf 1)$(tput setab 7)------- iGenVar built (2/6) --------$(tput sgr 0)" 1>&3 +echo "$(tput setaf 1)$(tput setab 7)------- iGenVar built (2/7) --------$(tput sgr 0)" 1>&3 # -------- -------- get data and unzip -------- -------- # cd ../.. @@ -40,7 +40,7 @@ mkdir -p data && cd data # wget ... # cd .. -echo "$(tput setaf 1)$(tput setab 7)------- short reads downloaded (3/6) --------$(tput sgr 0)" 1>&3 +echo "$(tput setaf 1)$(tput setab 7)------- short reads downloaded (3/7) --------$(tput sgr 0)" 1>&3 mkdir -p long_reads && cd long_reads wget --retry-connrefused --waitretry=30 --read-timeout=30 --timeout=30 --tries=20 --no-clobber --no-verbose \ @@ -49,15 +49,22 @@ wget --retry-connrefused --waitretry=30 --read-timeout=30 --timeout=30 --tries=2 wget --retry-connrefused --waitretry=30 --read-timeout=30 --timeout=30 --tries=20 --no-clobber --no-verbose \ https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/data/AshkenazimTrio/HG002_NA24385_son/PacBio_CCS_10kb/alignment/HG002.Sequel.10kb.pbmm2.hs37d5.whatshap.haplotag.RTG.10x.trio.bam.bai cd .. -echo "$(tput setaf 1)$(tput setab 7)------- long reads downloaded (4/6) --------$(tput sgr 0)" 1>&3 +echo "$(tput setaf 1)$(tput setab 7)------- long reads downloaded (4/7) --------$(tput sgr 0)" 1>&3 # -------- -------- get reference ../../data -------- -------- # mkdir -p reference && cd reference wget --retry-connrefused --waitretry=30 --read-timeout=30 --timeout=30 --tries=20 --no-clobber --no-verbose \ ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/technical/reference/phase2_reference_assembly_sequence/hs37d5.fa.gz +gzip --decompress --keep hs37d5.fa.gz cd .. -echo "$(tput setaf 1)$(tput setab 7)------- reference downloaded (5/6) --------$(tput sgr 0)" 1>&3 +echo "$(tput setaf 1)$(tput setab 7)------- reference downloaded (5/7) --------$(tput sgr 0)" 1>&3 + +samtools calmd -b data/long_reads/HG002.Sequel.10kb.pbmm2.hs37d5.whatshap.haplotag.RTG.10x.trio_sorted.bam \ +> data/long_reads/HG002.Sequel.10kb.pbmm2.hs37d5.whatshap.haplotag.RTG.10x.trio_sorted.md.bam \ +data/reference/GCA_000001405.28_GRCh38.p13_genomic.fna.gz + +echo "$(tput setaf 1)$(tput setab 7)------- missing MD tags added (6/7) --------$(tput sgr 0)" 1>&3 # -------- -------- get truth set ../../data -------- -------- # mkdir -p truth_set && cd truth_set @@ -69,7 +76,7 @@ wget --retry-connrefused --waitretry=30 --read-timeout=30 --timeout=30 --tries=2 wget --retry-connrefused --waitretry=30 --read-timeout=30 --timeout=30 --tries=20 --no-clobber --no-verbose \ https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/../../data/AshkenazimTrio/analysis/NIST_SVs_Integration_v0.6/HG002_SVs_Tier1_v0.6.bed cd ../.. -echo "$(tput setaf 1)$(tput setab 7)------- truth set downloaded (7/6) --------$(tput sgr 0)" 1>&3 +echo "$(tput setaf 1)$(tput setab 7)------- truth set downloaded (7/7) --------$(tput sgr 0)" 1>&3 mkdir -p results From 766b72de78573b50c89833441ae4e23d94255666 Mon Sep 17 00:00:00 2001 From: Lydia Buntrock Date: Tue, 31 Aug 2021 13:42:48 +0200 Subject: [PATCH 4/5] [BENCHMARK] Add pbsv to plot Signed-off-by: Lydia Buntrock --- doc/CallerComparisonPlot.pdf | Bin 38242 -> 7667 bytes .../caller_comparison/config/config.yaml | 15 ++--- .../workflow/rules/callers.smk | 56 +++++++++++++++++- .../caller_comparison/workflow/rules/eval.smk | 15 ++--- .../workflow/scripts/plot_all_results.R | 6 +- test/benchmark/envs/pbsv.yaml | 2 +- 6 files changed, 75 insertions(+), 19 deletions(-) diff --git a/doc/CallerComparisonPlot.pdf b/doc/CallerComparisonPlot.pdf index e54eebb521c44143a9765a52b1c1564b78f09795..a30922a2a4dc7275fc57e393b57296a1db52a04c 100644 GIT binary patch delta 3825 zcmZXVX*d)N_r|mDvS%4<$(BL}GeX8tV#pq{mG!YEjD#r$(b#vPvhQm{wy{ozF_wg} zjvvz*x$pbaS*oBuR8LZo1StX)Re&nWDl*KsV7|^!rvMi* zX((7pK}k{Jzgn{G9L;SN#Z)lWb%c*_Kj%EDHmmp(u85Ra9YuUY#Sp8n#?s4jFWFH* z2bBD3<=AZ9T+0<_uz6t0)WA{RTb{F7Q&R)Bm=8J$s>j#jwTuxbJs05@`@elLHui`+ z(-ZuS+lcVcGmJ)f;PbV8LZueD_QlEZX8(J>#y@AynZB!}S}9Mc1mjCgjs z=Xz3c*^wzIFmT_A;Iltn^YMoG3Hsr7E&lgN{9fv*?>g;}O`h#9Ek&wEJuNB69oSux zH_W=wZiluNxSr5F4&(Nl#Iuy%?KSI~Z4J?)6U?XX9-zRf*>t)D8&5U6+t_u|6*r{3 zm^`1E-6b1!6}Psk(F)EXm?#@lX)$*re;WfbtYZMnHbHtOYHrS{e`bVW`@`-s)gHde zSoK1e=A4c(o?TSNp`79K7uybPAA?I&eA#{}{%Eb%$GG!iX3)MfOs)xYR;Si(6@*AW?|0ib6bSYY{_$h} z;K|FDx$_5O`QHnYHT8S6RAlV(Ds4F4KjGX94@=zS(_@*RlW3>i`GV zTJwRTfh%&~uq%3e9M?bV&0{C%;tf|fx0+hf5RK|4`;Ym9TO3zhsA~Az;~N&AfSCYO z!FfM9ZMAUcFe^k=TVR8Q>_lECWZwDaV^>sgLx_v?I*~Arej94Oolmngk&4zO2>p{O z2eAJImM|*We_tcd%@>28I`OOuu+s0}iQ~--hE)~=?VZrlhr=Ud>^qQgwknUu{d zu}+kv(6dGN;Cs6M$6Tmwb~CBqH8;*`IC5gTy zVU%8+mTIk+mk2bC|NtgKPzM75!o!K4zRsarI(uU@u`dEn1Y=C z173-19?@XFk59S1%uAiN?n|@Y*X{ADE3t$iXeH>V=_g0GN<_232h^0*&w1Uak6tid zq{8|AAn#b$rxXMG#!$sZsB0ck)=+2aZD;IQ*K_rfgps3xw28}AY8s|7pZ{5fheG)O zUKsu!D*F|=!T2-d^|CO%Apu;P-R88#AM2}hZm4+;47u=f4hMdEXU+(MvF6!m`+h_e zyyp~gsES)eLs&WCV5HuFAbtGTMbpKyT%pj)H+F+a^P~w@a~wZ!it(ZxH@%6o{4(-k zQgLIBqijaQx=;+EA#b^Un(2PaZt z;|-9zlO}ODiJ7nFjTYUuLh_|$9t}jIb=><`8}R#dy{kRa^(*@D5EYvQ96OT}M0(LU z!=a8SHENZqK;)}`VI65M9UO?GJJ+KtGyy+lZSWmuSD~{hm%%Z36qhfa(a>=$=Rr^m zbi_qWgu3kK!ibk2y5c-!k;{Pl$@!HWiOl_=EXu=f5#DVR@(|#>jXB2xi|FrjtB!d#06+z|_Ghvjk_Hp`IeT!AA5o&rm(y(z6k^^*4z@)a6(5)H&#gycIba#PrQ zr7MEp5i1%mzfjK;Ri@fa@03iq>fUG74%;4cA(n*phQ=5SbtR9{WZ#Ah;7{OyI^QN4 zk^;#|`RW^^wZ|y8KW6O+7yxZAa=_;4N*KDM5qY&SI(e))HWl+6r&^UgSsrhSA5O8( z&Wml1ZKUrOfA97%KuM)miCr8wK=j@Se)qWv9`rX-R{)X4+1vlc)C&h|Y+hLE>K1ooco5K(CAh{kdwMKD8*dPtzZo% zq&`RQ#WEQM|ITn6U?%ENX4ctTLqaB=7U`IEPPxKJvyqxniKm*vN8V~1vwow#GJiho zCZ4X9`6is6Ul=$jIclpwCqTC6^zV*F4~`#{c4W#Uq!SuG<#mG^O)`T?J}e=DD^5r&oRE&zbls;udMp~{l(%OiH{wm6u`3~CQds<@&IB{v6pmHh7deW0D; zibBZziQ~w4q;AR7-<6dmDmLZ^88BknPqhEGC+tpG`CG*^o(uON<283i7FWU#kw0{K z>2wxV2EvudUWujZ`=a=+7Tft8+}Z0+M{g5!PbzOE^oeD+E$C|k%Sbb{V5{Yx1VXs= zX|pmMw|=M0ohEagEhTH|#n=;Oms?F$!P3>%U!Hr@JFV`Wv&$LAt-A+#JkNAZTCs%1 zgzhLC*Z?6hnAuBqtb}?2Zo_Jtd${%3c*%ASmNBP|g+w!XSR4JJj-jmZLCg5yPKBkf(3=TkiT{t%)ElefeI}3Ku5Zl8{%v|qN z?_B0WysSR0zSP>JqV8#?d|#ns0Qhy7`YKf9b!oJ{704S(J6W*7(;$pTxD#5bon#gB zlP^GIT=gv025VW(Fcu}@6Z4r;JhxD>E~v*Cub7@cB=`-C&tv*ErEX3U34#=M>SA@w z{SCBGgTQCeNleyzgL7}WuBN8xOQ2R%zBz}?zL^kf6&)`o8pbkge5Y!C=_eKXyo4@c zbxYn)98e9MpP4m@2iYK|>UpO+2*Vj-~R?rpj^5lzVwdbky`sm$!lMsz={6_Io#hrU`SdER9#pt2MIYQcfmA z40T1DF%tB`pZsHIE{-^MG_-4@f_5p^qp!S=4sAHUwFT;TG#n$As=#);pL9F$qrcss zKP%J<=@~O7l0)VcMUGQ3_*)UGZmpUD75AfGuSEqt7{w^8gw3on3d%hzxhvHp5b6J! z*e0f@F1Qu&Wo%fnEVbKX>3q9OeP^@`b-V&ykhwChoB3fijpS$ItuASBCxRFpD_%0) z0<|D4z`wqUjGIQdH>Mxfn0V+6b2wePaBeihIfK1iR2(={8`WKOR7|&P)Sq5wTv5I9^6T2dqpdV6 z{r#fG2YymOiw3Yq3D=0eKswauv=r%Q3X+7UX2_TbT|-p)0S8x z3Q#V4@jiQBu87zT;^E?a_-F44@7HUS!VdJN2#vn&WR4~WncjUqrGY;K2EGpJ8im!% z)Itc35S3jC9o1oXwqf^v@{SabWK!^4@!~X(vR!D7tXjlWB$?7z8&Sj4(V2AXqGI)C z_bOe-U4$cnvb}%0Sc&I5KK}_We^~_0xiYYZR$BP5G(NA;nFjh7#n0Ot0jt)~nzNZDnR<>~&Yb*}`R(-Z zEDznml;7eaVXFc$ZAZ*xXp8Yg-v+*~a{*jIgjCpLT-uRVkEaDIGJORch+lFm+?s|; zmyZKSySRE*gy~Lp$@Gvtg!>ypXz*1lsV85z&SA z@Z^`NAo33 z6S747Wu&74q(E>XD3$_Ah5mYVQ*ofh(}>Eal&rb*>IP#AcmDU$9G$udm}xF0lc6^J z=cAM0Y~hsrId*kvOb@$b@lOf$Ty+%UdTjbgOjS(R27-mj`|m>s3Py+aCK6X4(*-(P71qew!b!8D;cna(1_d331g{GUYUzf~^8GjKV zmTB$HT4e3bz5>!7USRR>E4!OKTojIy$94Q(Bm5S%ziGPMPC4nFRyf9XvM<~&+1y-< z(F#8lIatA??sAqkjZ2MNJt=s-jt@aLlF^?CJHDZARq#Q z6lv0nAW8@6MLGiF`}?m+Lb054zUO-(-NMl!SaUVE*3-Q|CmwRU!~*J{)%YYlb; z6?^N{Tc_XJ{nq|=u|ZZ#s3j)+&7cYuf@;X{=P=IVF}H^_^~c2|#QNyZ@l1beki?D#{^(Beoh>$x&1MOy)g>yP z8HCh|f;lz(zL*HV=#GnrXWfE&f3&mB_Oba_}j*)*uAbp z-VS+3rmKIV&O6JtUOV?uzioxP9(d{Yt&&w<+&_Ev+7bWS{-EfV;L5>0`>c5Q;MVS# zTQ`r+oiV8Ql5g(Ee>vY8Ii+&tHGPk-9oQ%C?4_xb2c7%l>@St)9v}F|jvI?6wjcQ2 z?$r~W*;)-Ow=LwicZxUow$j=T9d4yBv;K?dw*6ZzZQK8=t|jZ&8CbUPol7?!buHb0 z+#iuq>lRFn+x~D{<`10rm!yu1Nn7qKr+-Y(m1*O+Vg9+UwgJ) zQ!q{DtZ$uI{HX3sZ%`{=i{tLzhrB;7&MnpQqgGcLGI8*~1~%-I zr*rq|r@HT0btdm2TdA5Se;nWNp0jkr>z^!p@1VncUT$GGpC28{R{YaabE=e?U-rT9 zTst3Rif-6*UW1bhY?H$dePJ86^!qYp%4VKa>ZA0ok`eja4w;xLDqD%DDia$|-t^%c zN2afRaBRbtnB!FrttgmhOOKN^?`NIiozkjzZ2miUqWiR%5whZ1ja#=4{k&*jvkql* zrFS>D@Y1T*WqV%mzI|v<^+F-(mY(gM{%G#KJ+@xzx#@EB_O6=?ZCO0LfBFKix;<4s zKk?G4QS+aFY0CGfFZ7JNRI6{BD?QU!ob-b=_*AMdT75Id(|ul+w@$w_e`%$Cx6?j! zx9@TF;QH5Hv&s)je>5WO>oV3;dv7(YwsOJI3e($%mbx`P?dcz8AD-E#$mQ5!>o*?C z*W=5Qe_W4>%9~;GqNRtves^Wz0;@|m%74AhviA@69(J(sooO>RWVu=@s@(UtqwcO- zJLW)z=cX?%wzS;b+gY2v?R@TwTBV;mSY`1`*Pi$Fuh_@)dt}B(IctuuoxfDc)YZbP ztj(Ek$jzYexo2vux;S=H&Z2b-Kde`(dcIt(3$0rl|H`)+$}f&xaDK>Jm4mh{i9Pkr zk0rxet$w+cXK(D=gWuR(&v~lV>>i`XY&%*c@8L|(m*~B1cebNJVeLGdhn1xsME#QINbAifrWbgaD^r&(8>${2{I;dafU2H%W;Phmy>v+aYpus$Y2|s**5I|V zKd0I9@$`K|f0y)24go3{CNC=WqY%H5INMBuiUGPn>5&yL<(4GM(*(D3Z_cfI|$Z0OP9?ow;R!9(`F{QSAR!zX>1 z(8%6v;|A-?V+u!nHey}X-#X1c)3Qk4>LbSOUw0srW*U1)lXcj}*H+blAWV)k0(2wJY|o3ZKpGH2cbP6^af`UH8PG z)>GE+YrM7P#L)OrKQ35MwNS1s1!8Bd$+W)nd%N1DedubqID3uyV+LocS%1ih-#evy zadXZ>o728k?(4?gYrkDS$B|CyE_{%oQ2x7?V*8q9x?aa~b^Uj{Gu^4Pxb}uA6(7cA z`gPof7pBf!+hgO2*Yg}}{nGJi5oITq`yy?>SK1c1k-2M&o2MGoX;Aaa5>3+;S$L}4 z;&k)Moq73SnHBH;a&qz)X)iSR^4`INdqZt=znHzP*8*SOe3R?8uzy;s|E7eF*H_c> zeQUc3FYMS+GT%!_uSb^uB5lULwdi$w-dgrFo5xHTFlO?-xnF0=?`(7BMvodHk&Pl> z89lZB_66&Ew#(RP)gSFEWgg(U*1UETdj|Xbh;u)F8MmeN!idScH@=zw;MvzV+m@AI zaCG?COY4TsTA1ax;m-MQHt)D#+>S}^;d6&y8r!V=_TU+FmM(8ndf=7j)(l&}+0vw4 z#ye%kT^bkDYUvkOY*R0z_-mDyzDS#ARR@T#zd3W8=0y*->@~OSl0i|6uP*K!c3{}{d1J>lTH5Gj z(P<;gyIYmt(e&UC!AE9Av^+HV{Kt`pTIBT9bJR(bd2Oa%C8v~qvq+JHaa)em&-z-c zUU?5(Sk*o6q|Ebchomd@$&{#uBN9gJY`v=ckg`oSR~VW9;Eh=uS5=?y@15tX*{^n* z@aB~U2fpdzqfO12YOA{s*|59miG)g>e$4Rg#$7QRYOj5-`jF%OmTta#XI-9?MUKx+ zdwjV5|Mt^g(;hGQ+u5z|+H3##zH{9diD8uHtY7wfci;5j9b1b}Kas6$o1-(jwC=fo$JV~>y06}K zHssylI{ro#OVxYtW^AWgt7Ga>3XW)dk(((%11+EcD2iJYvW1( z;NyKa9)0QenHQpBPWigV_8T&N-o#bz;~smLV|OcYfMwK-=18`}AY)LZ^4=!&M2b+Qsk5`^s|fgm#Bk&u%f{*0C)&-kf!HTEqGo*TmF)bJCj9 zOWurMbFOo~#l3Q!?(RRcEyviNg_l2?ocZ+jp8HCl7045kDQg}m4DUpZ0nXi zXhF}L6IYLjyH>LD?alA5C|m1d)h3xQRh_-O{fk9xjc)Cj+WY6x#ot_5qW{5Gw_eNp z`rf|Vw#`}G&RTesWz)~o4xC+jWo>AfGq%#X8fzNvDVOo~up?=fWh{O$ORb+9SwAjY zuXx(`=6vZ{v*@s6xBs)T>AqiIJmJ*Ly>(;Sb}kaTvc!E5^5h_l6b9`q}FZS1{+^| zdjtEeYNwVbuHSi{;gr}f$5^(D0pn>IroV>zn`qzdtjRA;{`hmJ(uUi zqu(D|D$GCeZg9?D>SgY;`@xo#FT|CYaJFNm4gQbv-=l_^cqDoB3O82sjS+FHX!s(l54nAshc53P1nY(I- zdWUr`>UnhM@}01Pag#e0taRYqw;?MkA3Od-%U&~!6z{Wl*^O3nW{fD%VEfurA!Xj4 zygk#FJ133|xi)28{1^56dP1*@HSF%t zyqy!?xzprWj%*9}R+>37p+TN9^X5&=mGRNdrL*1`_*tQq4YEcRyzxoKO)bx-KJA?~ z`TWTvHAC*Mw7O=6+CLt#q1G2CSDo(oLEVeN$FJwP75l}UwGBc%89T?Fui2wwv$+*F zx4%-cSMD7b+dbOh_tnYPYrunzm+pRB;f+1LD?T`Vy7h?P&n}&PqyGIo#ogaNyfSx0 z;p6V9`|95gX|{IIgR>0>{C;8I;C^?`-FoBmAL9lP#3#nXiQ&AW^*}V86wU{t_>lTt zBH}DJU=1cYN(;@k(-NMi9^NQ8BU2 zx_f=*Oh{bQ+SM%1q~?#44V}{1ZA0o-uU^d?=Z~7a#-pnOD?6$SK~mB?zvumzi!0B ze8<1(b8Y9;;|1>y>^fu3?xUTq4LGX0^(7^4z9(#*FIJF=Of% zeod7&D%FYZOP5xA@m{SoQ_}SRadgy_ocHpqS(U}|+k;p4S024S)gLq7TiLp>HPve$ zj6Qec-r;un@@(@|Z23{K6^*xj@Rg(CuF!hB3b%gS)?#AZjWTWf-s)=h<@>U5`5!%Dp3rV(yR9vGU{L#qQzzW* z`F@9?gRbVQQSjczJ1qxJ$((BO9RKpe_eXT=k~>xVOQ~O5{riZW9Ya!${BmtVmH|Cp zS=GNo)9+4|=~R0Dkt6Nur=HQZRLJ>j_h)1oRQ-+qyPK{(dhHi)@0$}n6XLDaJGi4G z`_8HH-PT>Fy8rgV*%>4570mc3;guJ&KPnO7m%+bi_tm=W zU#H4?B5j(gCErQcqoei2+quT2u3^cOYw+@!dEa|`VaF2lG9=8*?ik|go;Ut2-wB5& zWAIzsJEl95x@l;Z>4QuC8nmNIwxPqeX3O1lM7*W?*8E=%Yc=lWGilRT&;DBe_ns?O zy;_|#xvN~7*>qn1&xc=}6?wkx`I7N3txA=-`ihxNI_J2byKlxKp_%5aZ?iaiRQjbI z8+)D`J2cPwZNJ()1%uySv2IfLwEeSP9ol%~=^DTXO|Jbl$`8JRC&RNoGLAitZFOAH+DRyyql=Z8ie%bov z=$F3y*a1U|oQm)MmFJCi6&k%$`n%UpwCP@=d&}-Qzs~nfj$_%bjn2C1=C9kY`|s82 zpLO7KmB-E~R5Iu2OkWJGQ|(o6*;y0z)biB-x_MHeC(L&t{T`4iG@X^9+EEBiYipUyUp>=~})uM~fD0IBUw@qfg+UnJG zC4bG;qV$21@vpXP(tFN3Egj8oHE*}zgGH@o^{&yQ&h^i7bl6ts-G$#ZY4d6G+40RU z4BxlR_WI#kVa>y;_gNY{>PGrsYUe7m>77&6Pc|7@@wBf|-CLRNm-)Tu?*;C!x!){9xeSFeq|dN0!-~w;hnLxyYh(N2zYkwD zX7rfw(XPA+`6}94*t3N0az{JYcCF^>;cgwduVd`aT)yvZ4edv}^*#K8YpJ`wb-wjw z+lZY7cBDVMdgq9;dtPkxVzp5jM!no1$Nq^^Hb(i!w%J#9NZIGht}Zj9b-q@c+ICxd zamD@4-D0LiAMHF~%kc1=ott&u@4GpF;gqt6zS?}}r*=2m9bY;BT+iJDqB})zS@E;4 z*Z%2m77YJ({n|s1h76fKWb3S*-hF%5S@$kEcs;|Y%|X6yk<)rrE8izew>h4Y*4t6NKD9fTlQJwD8yHhIWyPtQj{;uDi5vOwg((Fpy z@s_*Yhgba2>-*iuVul?rHE7wOaeuh}csX^O)JF!N9sF69E$b$|HK6JzGt*c7p=#E* zcfMV7eZ@lE>vY_DXWZ(Wna@3cuGg@)_YB^1W!Bvdqc^@kY(T-WnX*2=DAS<4-keRc z?frJ}s*c;T%>8FxxmEWV~<+R5WT=+V4R_Bojz*p}OG+D0th zv2^6pPqu$F`AR!iv7j7IfAzzozFfV@ad`KE!Qa*RvV-lz)9*K%c(b^>+}%3UW2cA4 zbXwO>ACOkzrN@5Pv`e-d1+B_n~G~U$LxG- z`OxKEx(+J2EVxWa|CPa!?HiA8vG>D|m)stBsqLuTBWrGQyS|Q^_wzeHmHg53b(gDi zK5tmCe(z#=>MyKcy58}av(ansP1@OV(i;PQny}2jc*m@$-7R*sIC=2W&qr1iDA_W2 zcSymXU-@!PTt;8*ZDS{wpM3e7*Ssy8tQfVp{+IRJ0n3fIu zEsVS}=chxf`>bAc`TCFhTg-jEWt+Rvi#lH0+4tM~yE4okK61bZ1HvA>;;QAoGic}i zRV(Yx{%qNWjRQ{(yztgLX$R~ayK`sR56k)&cAx)A!i|m_)^9)Fd(i%)o1^B(%v~}0 zl{6!Yj(BCz)O%NMZ>&0^-@g`*T{fW8z4qNQE}#8L&Wr7ih5ayTU$HA=FTGOfrQXX{ zZ25HV>3-X<-7UDkVvD$Tzw|zwX~Vl8=bX{#?2L~pO-i_U;li2yzBTpxUi;|Xig61= zzKQuYYR}=D=lbuScH!>j%6*RAez@fMU3tcKh-iO#)tT;}#jP7S^K`}s1^^JMo{`kkftNv^8Ki>MiWxhp`9VUO%WKi6RvX{F3vHF8vIqfU2ExA^@@&^OY z|CaaG^{aD7H=Z#t?cJd_AN(oo4oeP0gA%LgHKQq5cZsHquCVp-QOhpTRiop&B%KHf zsnI1e(ogAJwehwVhcl#Gm*|8zOPDL9M?y@zKgu5&A7phnLLy?Q_Q%C_v4nJr_4Xvj z_4*PN>iW9G`VzWDM)`YNLgKqbMR2N{*B2WTZ3zjF^{XLwdVM~s@j)RGT^K5^OPnPn zI65Xi!XIe~i8Vy}$O59Gyic6(oDki~8=KHA%9{{x35n?x6YcNHXkOEEoFe!3Pb?#t zusB*&x-m(73piqm|78IKNpT*^L+MSk2OY%7zX^j4} zL?@rKy2!JhK4-Uj$VHwxXLB*+Gw1AKf4`A}4%ZX+#>aZQMEPS`gR+D)>(bW`4=f>dqa$O;2u=FdG$tnAqO~Ezuef+`Y`ozX z+zAURR;*^@T0xH&f|7LJECjLm^yEsAQoW~^f`b0E78DdxH5zSC%Fxu&H`0b&N2nYW;oerzh>Ima23!@yl-FCYh8p0f2 zhsWu4`4ElY;dh6+Y<8O+9p${+9wwo?BAhOl*X3}~pVO=J4wpO38D_V+JhXS{I>#=D z-D$Jyemf)Zk012N73m1Kdl=j6!04gC@AkUfj1USHekeda_~e5huhSzvff|Q14F0%K z96RIEzl$#E4+c8I+#rn02h}K)JJb>GWE8iRGfr2y%jWWPpWS6+_8QmCxZ#d)huzJ8 z4!?uW%z3BN<_L3g=zs=|>h`$Y%+JlN>>8EXKnfJw6GvraXf>m1R?y&fTOHxdT~zrv zj&#!=Lc`o@8|I3%TOEFvFOd%Thf%>F(Z@*4*zIM`P{1+sRa+Mn+r!Lt&0l=w6H1^o zk!~oXol_J-jolvV^dbjjikzTIpKi`cK3wy-J;+1%doV;ddO_##2>mefvb!Sa4FlIP z-|W!rg(G|qjHhw_VLTV(_~1Lzae3gE6D;N{`lNpSiQ{>Y2W(;yr9(F zFQ7;o*JJpA{vln5^wJd}(K#K^AI3brPV5Pt3zMEp_I4yrKgQtTf#GMOKGN3ct>5sO z_VC@qeO|_sR=H>=zDTnY`4i{|TImP16U}N*Pv|#xCHiS^ll(*<@)VyU6Zy@t=5O>6 zJ=eSf^MOw61bt2%huKBYui?M+BXAuFnRz6SZ`HBr6y49tgZ=~dZgY!H>^dL^^a`hC z2Wv#1k%w3O(s97fk%>8me*rzQA^!k<+{fKMCtxjoGak%9CYyT@=)1?j6teB>B9pwoq3B;rY;en`&3gGkW{PM|;HgY4XH zXMFTO8J{>$Kj=M{?BlKo_(1=%1CFsbqtEh0emhuX=80VXTl+J}AXV19G+wCDAK@c> z5&m)!ej9qFKcK$!D3SD#Lnagh=*+9PQpH+Mxs7Q zIPejU1pJQ2@Kbhb@WSYh6yV%yf%2CXb21^7V=0hWw+QIR12@{PGS`#37+6r0`ZaLjfeBY7v_Z|gz_Hf zl{~?DgHuWUNK*c!{)D56?eP1|TX+cW!+C7M^oK%qM;b z(?q^}qOWSNKA=)~Dy|Yi$u2>dfZUOSaFYHEzg!+iBmhY4lsKND-|d0l+>Z{2;cpD1 z67e;hd4k;726W)>!lbuw2CR5|UUX_+gj2{Kn{X!a)!;8Gs`G|ke`3G#|KN}Kjh(P@ zkA4H~A~ilb1U(WMl#0(L?xDTpA$%9E%f6HSYZAW&C=^WPs^hO|eTfm}-<%z>ZzkD@qykM`< z6Z9Yg`oj#pXr^=>$^!TZzq$U8b{Mb|V|SXL@DG6p>`wDx{@}H6QsE{38$5H_h>87* z`_L2S0$QPmfQ=g~bTdzAgkC|7jW+U%(o?V)Tv1dO=-6 zdPx|?vCx-KpRvQo`l?)ycnYcED-*9Sc)UWV^PckMv3OGkw z6SKf=`6Vu4%S8NW9r%vPkXzyb1AI5}ARfi|Q8b46Di0wCk=*5}kRE!?B)N~7AXXu} znGIL*XZS(Q(!|L$wsNHN(ir%{-Pji*K%*7ifToHQmE)pK%t|da4)K=m;U7KnTi<1G zz#nlBnWu@>0#|eujfGxjg{?{MiSjmA0@wfkH#UIp)L4O-EYMPRD>Kj>+=iyV@8SRc z73I2q#oy?4(rC;!d43Vf7vLoLA}lgLke_rNJR?`bj+7xM-37iSU1yf^S0LjP*NO3- zxGuS2y~G~@DIq5ly99dUer)ua>wrYkSf>AHu2Wk`{^2_IiXDT`NoSRFBz_lO>$>nd z`MU2heds3j5aWj)Uq^z;_j}-*;(V}*zD%#CmXSP`@E9!++6CfsE3pqYl=uUETqSM< zK@xj~j^zG{SB0yJE2$qP_Afh3x<8=z7#h7N-)~Q{4OzYT_m9`%|36(f<71hR<;&QK zc_0A}bes7r-jMC#p^1M~4*@UaPlYeUMzUY@2%9i|8@jQnMEnl;fkeM4J50u1NtWm> zATMHR@fm;eq=b_=q4@c+eQLab-71eUbeOn`d_@%_kR|E9zy72q0vbHV6F^%aW5yOe z&_KV!O0pb-Wq;kz|K?fhcL6&}(v)N^vD(qWfL1>-Z)_hFC+mi*~qWoV*(D5e7&^<%W2Eb#B>!viK^27V)PBz!a$D0@~t?Jx6X zxRfN%FziEFrtHW{xJq%G6$QTVk8oA~24wOoZ$suuqb1JinVGBY|CN599q0emW7UVA z#`h;@r&?Z;|1jqcD?${|>i+UNUzhSK+ScQa%%q7xno>zn55# zC5FI98aXG{Z;+_sHQ}x747!aU_@ChNvpoHqW97F2em&tU)t0h6KqODBxCVQdo~h~u z_f226aAkSYHL5mfvAN%5Z%=q_)k%`$%O~|xQ99mQ{y7oD%t&O3pdz7=(+wXQ4^_ApSssq4p$szfidX*mo(1^5EXF!g`a?kSZ=`)6QgLVq+RRe@y2A>tN zs}4ydZAPIYjGq_YD=VksMEgKrI85$TEAW^rCaOfc)sy;DPGL1UAMqE*9(2PBk7#e= z7ZW8CGbj`Eaoz9?-fCn~?^P=@U5?FfMHyC*hxSO2T9_yYVHD|$9}$WgiBLF(8_M`~ ztdW(Yh;kh(aw0~7@=%WTDGqWL%?{<*l>hNA%ALtHiJiFaRy)NwHj`;MrC*Avn7R7b zVik2G_@@;%j)^T~7i37tgJX@#l?aaIuT-U>@&W&q$-*hkk4Rg&8!Hk}!+1!U{>WB1 z7CrEk>#7?&Jdy`FbR?WptQe>wNWTM7xN?1}&Wvg|^s3@ZjSqsPCvX<@W}Odt;sd$M z3!GDh(EZ9XiS6j0ydKQf%vsUnH+NHnbc=qf?Ho(SN~$>q(*rSUz@8#xAJ7PMXp3Yd z1-MI&sijEph4$nfM6X7d{o*GH6S*li?WU^9_$sv^SYczR&g$b!a;7zS*jh(Pb8>@Mteo+Ojj#TrbbCSFq1p@ z)cwj)lQaW)gTv5A{H--S(d$-K-tC48GzA;OreG$wsCJ=at+A+zgytYMa#`9_U*Wu& z4@CtuNAwZX)4!1%6+o^_e{Ixku>(m0qJ3~v`VY@(uNs=>7D?s}zmX1_qpXfdo7OJs z0UY~jP34UKsIO{;0V-KPhMsVRTX0R8pfuTspHijHWERpyT_-XSa={;$U-vVrlNAr} zk9AdL9nvKHAOtg8VFGl9lX=1u_{%NGhqY&lAHowf6?%N+C(1%-PyP}f@CFD_(TFgE zXb??BYCwYYDMDJqf5H@U2y_C?Q+~%Wc>oxLei(FDra*g*2cM~$fWfMQO0$J}symr# zFu@YCGw6wDmB9$UI95IwDT?5Rmt*N)B4j-thi=S|V>ujr_ z-mOJtB!@CEp9m}ie%XX{SSRge>p~NwD=45)gSi^4A;O|B<#6=HHCaI5pQ$5(Us~M| zc8DL8w`h;naIDylHE?1(q|7~3Ryl^#W=-4!*BK9*Rby5E)UIS%Iu7J&x^7m$$gJpr z+P%Cf{lN*Hfgd;uRr$b6q9&&25{eidQq>7_&EdOx#K zwoB~~y9?!*^*gR>HjD*5aEAUObU!MpRY5FMyB8EK!87?o4>lw-RV-$ZOW!zFbsZhx zy6gZgko74GN6%&H>^~5kaO~llhc+NOx+%+-&(MAbG+)}w9n6d;!2?r$U^_pVkD$B zAL@`C3-7$xqjmzRm_QfU3DhV;q*ia%deD5WJrh?WV~|GN(^{RW;YUjEXs<}hfRTJ6 z{KsD5v#O%uSP5&7a2%?kDO|P=U$l0j$OIciM4{@zRN7R#WloCCI7@{epYut&oA7N~(YjaR3GQLRw-|j0Ri9Z*8aoH}J?(1T9*i=U6k>WdLD< zQjiTUpmbi4PP`Z9iY3GjY}69cq3%dX*q}%QY@$EewQP@LgYM`a$BH+ogP^zAJQ~Jo zys2?&Eg7JJJ75o}hIEBgpb)tqfMoWnFswpytx{15g9cfPYA?v&RNk2#b0>e{SSu{V zxb)&@WJT&iMj5MLo`o?O6Ro4J&1uC9x-P2|`bhLzhxV~f$<7<`5o{!`fM1c4ijgVR zO8CNDSj)h-ie}^lPoPg8Nt)mTX@!m+c!(cGPeb9c6HaO+mTs^MJO{lK!+_TuLo+dn z_88I1Lk*CNi|Jl1JIQOyvzq~6JX3XX*(=nVH6 zCBnjp0-#ozs-g;Xl-c9;pjOC_))OIxQuRT)%qM{KkgVkxtbsLHqlb#f;3ZljDTay^ z$ih$RNR=0)qPz=#6_r7Yzji*1F#546?&TNl399n4CGT=&~?@g zawm=8GY5O}+>agen&{ zBack1D?~&aROP1%WiSTqU?zAL&{TH|{T1Ph2G|RK1?u>9{YEbp6M{pM1ksRCgS`*j zua#E(4)SA-7d@m~j?pF)Gb(shtg4Iy4<@@{eu|2<#wG0(La|Oy|40(b6j_1|0E6^M zHs*sDq1wj;Kj1ywLvf*)by(5D8DbvjVO^YKc75hsA5j#G^$}}{uX1X5G$_Nf!BJvbOdH%oHh4`Pfh6cG zOi&+0tHdOV$h6)7|KOzbOvf_)9|UhAq!Y|VQ4rix8!SgY1Nt-p)QcJ#P>xTL zA7lU|COkosXbGZ3ii);5mK7sehgbN4fE7(~4Z38f0vKf>%E-VAGh?t6D&(6~A(aMc zl>wu}-atZTqu5=H)Gy{N-ISnFMadFPktjZjTKR9TiH zd+nuEZUP@c5wmuQB|{B5r;J3#Lo9-(YA233WLDx0M)};K9AnAopwVPTCn9GSL=emo z>iJQm3L2tO(o<%OrGO*w5sYP)imTv`+Mo#n4N)VXuV^2FjRs?<^dVkzkEFtwk!YE4 zUlB49zx=jlMEs%bO*srUfKB0Rz-mQU=n%8g+CNA1AZyWTgHVE5Vi)M1A`3-P>O~rd zPAQTNr7d2W-zH9B)>>~S^yexOx5)tIfrXe@EVL4*V0mZ z1Wea;G>hwG`@~j6-rx?d8B1m72;2zt%zjtoAymY2;Sg(y*dH}abVM~c6PaK!#MIa| zGZ%?af=19cT%a2t92IsWG-K%zI2MDxPO?h%p*fkz1zw{k+=o{+5iOC6;$`ipkd3h$ z3yBK>sA39Z(Me@CWHpRWG)r5|lG$S$@JLY%w55NNQCNM=iJ&85JVR(iv0Kd=+ zE<-ofL=z=p+a~g%Pr4+^F;NM52ht>pf<0UZGl-tiVroGskhaxB2Q-4n)F9cH!LBt~ zlHw<>Lx;(1m@HKj^pcAVt{n{8TY&dAepzTHuuwgo32q?27(DR%EvwA1n4jgxvE?mqG^U9Dwn__RqKME7@x2aZNtixVZd#! zvvz}xYnKn+UD#;qpnQie?S<10X+`CzlTM)f%pFep;1NS2OJodvXae*Jll%seq#@MS z5UL`xNMt9!EWF3&aC=AyGa-h9C(0|y3m8HAsl9rLU4B>67wpTzu_z!D(eMeDdYM(Y ztP$$LBttdafGV(v)>H$C+Q4gYgSrdYqIgm@60Xpte&dn!$OF?5N+<5K!B&l5hl$$n z22at2z&vRV-r< zA9BD%kqx6MW|A&I9}t~9E-@wkRI;HzG)I10t8MaAcrJn@vKl{`jAtA|3q9m zFX_W1=nm}PMAq<9*{141rmle=s&*%ruRRtDvX$#9O)J>UMAp ziv>N{%bYl#sury4K)5W0+7wAqG)9&_ z8dj-a(QoAaq&^z@uuIX0)goW*d!RaK^qoqXxes=V9w0DzeJe0NEH&dRl7I~KL~=qy zw6jk7B7MTniLq35p#Q(0XOcc4M*IqBFD`Ju=JW69ko}7fy46H1$PO7JyMG!Fc^SE& zNB_hh)fde^E*PNwax^ymVsD^|says0399AyrCW%==!^J#Rp%XC^PESF>J+MR6SjB_|O*_)vj{ zFj83>IgjWej)I4L(|Qbh1C`;xQ`K2Xl`j@3ru^O;tTahzpnr%!kt4Ix4mj0n znJLHAHSv|0JH07iQQZZ5;h2b$<4E{M#f;-{9Wz%Q6XDT^>gCF4<$E{=CBS(iDB>uN zS?@u|(NPdrwLt2H+zq8#<=~hZDkS9o0O#Z*w8|k)igQ$6#7TOf>dCR_MgGudY(@dH z1jzkb)qs--kNkjR?Hi^kmVUd`V|`X4;C zX&!dn#Wg@8l^n7nex=8>U>n#h?+1c<4uj1SJh2kpKB>kZE8IkPl%^{ z;^{47R?zV;qA?Ra5_OqfX+RxO7;u#G)TW7wlcPdq&ZhZ;!2E& ziLMl7D0Wb6$$5$rDNk)GB7Vy5rMHRMF4Y?JOSMar9TAD>X*XqWDNk*xPD1SsJyLBR z)({N_ek+dWJocON)Mm<4o9x6PuS|TJQ~0iG1(}HUE8!i;5du3T^mL^%KGk!zijeZu zCh3yWDAg0m4U?Yc3+%T{d1~{i>b8k1P$yHZB+doDRh?B8(^P4x&8VtGEkJc)@f95m zRLND%nZKgc6nex=8_P8leOL=OO7?y2%DNk*tJhhqf)Mm<4o8~!;l&3Z+52rk}nex=8 z{*Fb;Q=2JIZKgc6nexMf*tU8mkcTJ<$vxS9Ap@inJ5y;4)x0! z^+Hlx(hEr|R7iRiY|=YU|6@KRz0wsv*ppu9n#c!x(wkZhAM8mlQ~lrJ1CK&`cpjaH zqP@JNHPY?pPm@GMhI%}HtB1F|JMfBQ6T(%V*P^yDMJM%sbrdH>TFx7w0k zr}{s6ajPxq(c^z2+K41k)1M^Th~~c|+Q=_(%m4GY&6;-+yDiqlmo)2d9Xx*FF)zin zg(h9$om7vn1m1_Mzr^$BbDlrFmYKgo^yhYt&_B0xIMCRqUbn0l1V2HMxoERlL!W&u zGv}V9Jh7ckFG7B*KX$u{)lZ*uhlZhW&zxgL`_JdXo_)Qv%=uZGIp_ZK7|uV!X3%i~v0^U7tta@uYO%$|7XG&uM8rPI7cnwOXBW!zzjPs+z7zk^zT LOUat}66*g08k@L$ diff --git a/test/benchmark/caller_comparison/config/config.yaml b/test/benchmark/caller_comparison/config/config.yaml index 3234e18f..01e3a2a6 100644 --- a/test/benchmark/caller_comparison/config/config.yaml +++ b/test/benchmark/caller_comparison/config/config.yaml @@ -2,7 +2,8 @@ long_bam: data/long_reads/HG002.Sequel.10kb.pbmm2.hs37d5.whatshap.haplotag.RTG.1 long_bam_md: data/long_reads/HG002.Sequel.10kb.pbmm2.hs37d5.whatshap.haplotag.RTG.10x.trio_sorted.md.bam long_bai: data/long_reads/HG002.Sequel.10kb.pbmm2.hs37d5.whatshap.haplotag.RTG.10x.trio_sorted.bam.bai -reference: data/reference/hs37d5.fa.gz +reference_fa_gz: data/reference/hs37d5.fa.gz +reference_fa: data/reference/hs37d5.fa parameters: sample: HG002 @@ -12,14 +13,14 @@ parameters: minimums: igenvar_from: 1 - igenvar_to: 100 + igenvar_to: 80 igenvar_step: 2 svim_from: 1 - svim_to: 100 + svim_to: 80 svim_step: 2 sniffles_from: 1 - sniffles_to: 60 + sniffles_to: 80 sniffles_step: 2 - # pbsv_from: 10 - # pbsv_to: 91 - # pbsv_step: 10 + pbsv_from: 1 + pbsv_to: 80 + pbsv_step: 2 diff --git a/test/benchmark/caller_comparison/workflow/rules/callers.smk b/test/benchmark/caller_comparison/workflow/rules/callers.smk index da65169d..fe619bdf 100644 --- a/test/benchmark/caller_comparison/workflow/rules/callers.smk +++ b/test/benchmark/caller_comparison/workflow/rules/callers.smk @@ -27,7 +27,7 @@ rule run_svim: input: bam = config["long_bam"], bai = config["long_bai"], - genome = config["reference"] + genome = config["reference_fa_gz"] output: "results/caller_comparison/SVIM/variants.vcf" resources: @@ -109,3 +109,57 @@ rule fix_sniffles_2_and_filter_insertions_and_deletions: "results/caller_comparison/Sniffles/variants.min_qual_{support,[0-9]+}.vcf" shell: "bcftools view -i 'SVTYPE=\"DEL\" | SVTYPE=\"INS\"' {input} | bcftools sort > {output}" + +#PBSV +rule run_pbsv_dicsover: + input: + bam = config["long_bam"] + output: + svsig_gz = dynamic("results/caller_comparison/pbsv/signatures.{region}.svsig.gz") + resources: + mem_mb = 400000, + time_min = 2000, + io_gb = 100 + threads: 1 + conda: + "../../../envs/pbsv.yaml" + shell: + # "pbsv discover {input.bam} {output.svsig_gz}" + """ + for i in $(samtools view -H {input.bam} | grep '^@SQ' | cut -f2 | cut -d':' -f2); do + pbsv discover --region $i {input.bam} results/caller_comparison/pbsv/signatures.$i.svsig.gz + done + """ + +rule run_pbsv_call: + input: + genome = config["reference_fa"], + svsig_gz = dynamic("results/caller_comparison/pbsv/signatures.{region}.svsig.gz") + output: + vcf = expand("results/caller_comparison/pbsv/variants.min_qual_{minsupport}.vcf", + minsupport=list(range(config["minimums"]["pbsv_from"], + config["minimums"]["pbsv_to"]+1, + config["minimums"]["pbsv_step"]))) + resources: + mem_mb = 400000, + time_min = 2000, + io_gb = 100 + params: + min_sv_length = config["parameters"]["min_var_length"], + qual_from = config["minimums"]["pbsv_from"], + qual_to = config["minimums"]["pbsv_to"]+1, + qual_step = config["minimums"]["pbsv_step"] + threads: 1 + conda: + "../../../envs/pbsv.yaml" + shell: + # pbsv call --types DEL,INS,DUP --min-sv-length {params.min_sv_length} --max-ins-length 100K \ + """ + for i in $(seq {params.qual_from} {params.qual_step} {params.qual_to}) + do + pbsv call --types DEL,INS --min-sv-length {params.min_sv_length} --max-ins-length 100K \ + --call-min-reads-all-samples $i --call-min-reads-one-sample $i \ + --call-min-reads-per-strand-all-samples 0 --call-min-bnd-reads-all-samples 0 --call-min-read-perc-one-sample 0 \ + --num-threads {threads} {input.genome} {input.svsig_gz} results/caller_comparison/pbsv/variants.min_qual_$i.vcf + done + """ diff --git a/test/benchmark/caller_comparison/workflow/rules/eval.smk b/test/benchmark/caller_comparison/workflow/rules/eval.smk index f0f0bba6..139e32cc 100644 --- a/test/benchmark/caller_comparison/workflow/rules/eval.smk +++ b/test/benchmark/caller_comparison/workflow/rules/eval.smk @@ -64,19 +64,20 @@ rule cat_truvari_results_all: min_qual=list(range(config["minimums"]["sniffles_from"], config["minimums"]["sniffles_to"]+1, config["minimums"]["sniffles_step"]))), - # pbsv = expand("results/caller_comparison/eval/pbsv/min_qual_{min_qual}/pr_rec.txt", - # min_qual=list(range(config["minimums"]["pbsv_from"], config["minimums"]["pbsv_to"]+1, config["minimums"]["pbsv_step"]))) + pbsv = expand("results/caller_comparison/eval/pbsv/min_qual_{min_qual}/pr_rec.txt", + min_qual=list(range(config["minimums"]["pbsv_from"], + config["minimums"]["pbsv_to"]+1, + config["minimums"]["pbsv_step"]))) output: igenvar = temp("results/caller_comparison/eval/igenvar.all_results.txt"), svim = temp("results/caller_comparison/eval/svim.all_results.txt"), sniffles = temp("results/caller_comparison/eval/sniffles.all_results.txt"), - # pbsv = temp("results/caller_comparison/eval/pbsv.all_results.txt"), - all = "results/caller_comparison/eval/all_results.txt" + pbsv = temp("results/caller_comparison/eval/pbsv.all_results.txt"), + all = "results/caller_comparison/eval/all_results.txt" threads: 1 run: shell("cat {input.igenvar} > {output.igenvar}") shell("cat {input.svim} > {output.svim}") shell("cat {input.sniffles} > {output.sniffles}") - # shell("cat {input.pbsv} > {output.pbsv}") - shell("cat {output.igenvar} {output.svim} {output.sniffles} > {output.all}") - # shell("cat {output.igenvar} {output.svim} {output.sniffles} {output.pbsv} > {output.all}") + shell("cat {input.pbsv} > {output.pbsv}") + shell("cat {output.igenvar} {output.svim} {output.sniffles} {output.pbsv} > {output.all}") diff --git a/test/benchmark/caller_comparison/workflow/scripts/plot_all_results.R b/test/benchmark/caller_comparison/workflow/scripts/plot_all_results.R index fdf397fa..d88fb701 100644 --- a/test/benchmark/caller_comparison/workflow/scripts/plot_all_results.R +++ b/test/benchmark/caller_comparison/workflow/scripts/plot_all_results.R @@ -4,9 +4,9 @@ library(scales) args = commandArgs(trailingOnly=TRUE) res <- read_tsv(args[1], col_names = c("caller", "min_qual", "metric", "value")) -res$caller = factor(res$caller, levels=c('iGenVar', 'SVIM', 'Sniffles'), labels=c('iGenVar', 'SVIM', 'Sniffles')) -# res$caller = factor(res$caller, levels=c('pbsv', 'Sniffles', 'SVIM'), labels=c('pbsv', 'Sniffles', 'SVIM')) - +res$caller = factor(res$caller, + levels=c('iGenVar', 'SVIM', 'Sniffles', 'pbsv'), + labels=c('iGenVar', 'SVIM', 'Sniffles', 'pbsv')) res %>% filter(metric %in% c("recall", "precision")) %>% pivot_wider(names_from=metric, values_from=value) %>% diff --git a/test/benchmark/envs/pbsv.yaml b/test/benchmark/envs/pbsv.yaml index 3b5491c4..028d22e5 100644 --- a/test/benchmark/envs/pbsv.yaml +++ b/test/benchmark/envs/pbsv.yaml @@ -1,4 +1,4 @@ channels: - bioconda dependencies: - - pbsv=2.3.0 + - pbsv=2.6.2 From 7c72242024e4a92d665fe6af40a5c0d68ad20a4c Mon Sep 17 00:00:00 2001 From: Lydia Buntrock Date: Thu, 16 Sep 2021 17:18:48 +0200 Subject: [PATCH 5/5] [MISC] Clean up parameters in Snakemake workflow Signed-off-by: Lydia Buntrock --- .../caller_comparison/config/config.yaml | 21 ++----- .../workflow/rules/callers.smk | 60 +++++++++---------- .../caller_comparison/workflow/rules/eval.smk | 24 ++++---- 3 files changed, 46 insertions(+), 59 deletions(-) diff --git a/test/benchmark/caller_comparison/config/config.yaml b/test/benchmark/caller_comparison/config/config.yaml index 01e3a2a6..cd726088 100644 --- a/test/benchmark/caller_comparison/config/config.yaml +++ b/test/benchmark/caller_comparison/config/config.yaml @@ -1,5 +1,5 @@ long_bam: data/long_reads/HG002.Sequel.10kb.pbmm2.hs37d5.whatshap.haplotag.RTG.10x.trio_sorted.bam -long_bam_md: data/long_reads/HG002.Sequel.10kb.pbmm2.hs37d5.whatshap.haplotag.RTG.10x.trio_sorted.md.bam +long_md_bam: data/long_reads/HG002.Sequel.10kb.pbmm2.hs37d5.whatshap.haplotag.RTG.10x.trio_sorted.md.bam long_bai: data/long_reads/HG002.Sequel.10kb.pbmm2.hs37d5.whatshap.haplotag.RTG.10x.trio_sorted.bam.bai reference_fa_gz: data/reference/hs37d5.fa.gz @@ -9,18 +9,9 @@ parameters: sample: HG002 min_var_length: 40 max_var_length: 1000000 - min_qual: 2 -minimums: - igenvar_from: 1 - igenvar_to: 80 - igenvar_step: 2 - svim_from: 1 - svim_to: 80 - svim_step: 2 - sniffles_from: 1 - sniffles_to: 80 - sniffles_step: 2 - pbsv_from: 1 - pbsv_to: 80 - pbsv_step: 2 +quality_ranges: + igenvar: {from: 1, to: 80, step: 2} + svim: {from: 1, to: 80, step: 2} + sniffles: {from: 1, to: 80, step: 2} + pbsv: {from: 1, to: 80, step: 2} diff --git a/test/benchmark/caller_comparison/workflow/rules/callers.smk b/test/benchmark/caller_comparison/workflow/rules/callers.smk index fe619bdf..14866156 100644 --- a/test/benchmark/caller_comparison/workflow/rules/callers.smk +++ b/test/benchmark/caller_comparison/workflow/rules/callers.smk @@ -1,22 +1,24 @@ +wildcard_constraints: + sample = config["parameters"]["sample"], + min_var_length = config["parameters"]["min_var_length"], + max_var_length = config["parameters"]["max_var_length"] + rule run_igenvar: input: bam = config["long_bam"] output: vcf = "results/caller_comparison/iGenVar/variants.vcf" params: - sample = config["parameters"]["sample"], - min_qual = config["parameters"]["min_qual"], - min_var_length = config["parameters"]["min_var_length"], - max_var_length = config["parameters"]["max_var_length"] + min_qual = config["parameters"]["min_qual"] shell: """ ./build/iGenVar/bin/iGenVar -t 1 -j {input.bam} -o {output.vcf} \ - --vcf_sample_name {params.sample} \ + --vcf_sample_name {sample} \ --method cigar_string \ --method split_read \ - --min_var_length {params.min_var_length} \ - --max_var_length {params.max_var_length} \ - --min_qual {params.min_qual} + --min_var_length {min_var_length} \ + --max_var_length {max_var_length} \ + --min_qual 2 """ # Defaults: # --clustering_methods hierarchical_clustering --refinement_methods no_refinement @@ -36,24 +38,21 @@ rule run_svim: io_gb = 100 params: working_dir = "results/caller_comparison/SVIM/", - sample = config["parameters"]["sample"], - min_var_length = config["parameters"]["min_var_length"], - max_var_length = config["parameters"]["max_var_length"] threads: 1 conda: "../../../envs/svim.yaml" shell: """ - svim alignment --sample {params.sample} \ + svim alignment --sample {sample} \ --partition_max_distance 1000 \ --cluster_max_distance 0.5 \ - --min_sv_size {params.min_var_length} \ + --min_sv_size {min_var_length} \ --segment_gap_tolerance 20 \ --segment_overlap_tolerance 20 \ --interspersed_duplications_as_insertions \ --tandem_duplications_as_insertions \ --read_names \ - --max_sv_size {params.max_var_length} \ + --max_sv_size {max_var_length} \ --verbose \ {params.working_dir} {input.bam} {input.genome} """ @@ -63,22 +62,20 @@ rule run_svim: # SNIFFLES (we have to loop over min_support, because sniffles does not write a quality score into the vcf) rule run_sniffles: input: - bam = config["long_bam_md"], + bam = config["long_md_bam"], output: expand("results/caller_comparison/Sniffles/raw_variants.{minsupport}.vcf", - minsupport=list(range(config["minimums"]["sniffles_from"], - config["minimums"]["sniffles_to"]+1, - config["minimums"]["sniffles_step"]))) + minsupport=list(range(config["quality_ranges"]["sniffles"]["from"], + config["quality_ranges"]["sniffles"]["to"]+1, + config["quality_ranges"]["sniffles"]["step"]))) resources: mem_mb = 400000, time_min = 1200, io_gb = 100 params: - min_support = config["parameters"]["min_qual"], - min_length = config["parameters"]["min_var_length"], - qual_from = config["minimums"]["sniffles_from"], - qual_to = config["minimums"]["sniffles_to"]+1, - qual_step = config["minimums"]["sniffles_step"] + qual_from = config["quality_ranges"]["sniffles"]["from"], + qual_to = config["quality_ranges"]["sniffles"]["to"]+1, + qual_step = config["quality_ranges"]["sniffles"]["step"] threads: 10 conda: "../../../envs/sniffles.yaml" @@ -87,7 +84,7 @@ rule run_sniffles: for i in $(seq {params.qual_from} {params.qual_step} {params.qual_to}) do sniffles --mapped_reads {input.bam} --vcf results/caller_comparison/Sniffles/raw_variants.$i.vcf \ - --min_support $i --min_length {params.min_length} --threads {threads} --genotype + --min_support $i --min_length {min_var_length} --threads {threads} --genotype done """ @@ -137,18 +134,17 @@ rule run_pbsv_call: svsig_gz = dynamic("results/caller_comparison/pbsv/signatures.{region}.svsig.gz") output: vcf = expand("results/caller_comparison/pbsv/variants.min_qual_{minsupport}.vcf", - minsupport=list(range(config["minimums"]["pbsv_from"], - config["minimums"]["pbsv_to"]+1, - config["minimums"]["pbsv_step"]))) + minsupport=list(range(config["quality_ranges"]["pbsv"]["from"], + config["quality_ranges"]["pbsv"]["to"]+1, + config["quality_ranges"]["pbsv"]["step"]))) resources: mem_mb = 400000, time_min = 2000, io_gb = 100 params: - min_sv_length = config["parameters"]["min_var_length"], - qual_from = config["minimums"]["pbsv_from"], - qual_to = config["minimums"]["pbsv_to"]+1, - qual_step = config["minimums"]["pbsv_step"] + qual_from = config["quality_ranges"]["pbsv"]["from"], + qual_to = config["quality_ranges"]["pbsv"]["to"]+1, + qual_step = config["quality_ranges"]["pbsv"]["step"] threads: 1 conda: "../../../envs/pbsv.yaml" @@ -157,7 +153,7 @@ rule run_pbsv_call: """ for i in $(seq {params.qual_from} {params.qual_step} {params.qual_to}) do - pbsv call --types DEL,INS --min-sv-length {params.min_sv_length} --max-ins-length 100K \ + pbsv call --types DEL,INS --min-sv-length {min_var_length} --max-ins-length 100K \ --call-min-reads-all-samples $i --call-min-reads-one-sample $i \ --call-min-reads-per-strand-all-samples 0 --call-min-bnd-reads-all-samples 0 --call-min-read-perc-one-sample 0 \ --num-threads {threads} {input.genome} {input.svsig_gz} results/caller_comparison/pbsv/variants.min_qual_$i.vcf diff --git a/test/benchmark/caller_comparison/workflow/rules/eval.smk b/test/benchmark/caller_comparison/workflow/rules/eval.smk index 139e32cc..6f921fe3 100644 --- a/test/benchmark/caller_comparison/workflow/rules/eval.smk +++ b/test/benchmark/caller_comparison/workflow/rules/eval.smk @@ -53,21 +53,21 @@ rule reformat_truvari_results: rule cat_truvari_results_all: input: igenvar = expand("results/caller_comparison/eval/iGenVar/min_qual_{min_qual}/pr_rec.txt", - min_qual=list(range(config["minimums"]["igenvar_from"], - config["minimums"]["igenvar_to"]+1, - config["minimums"]["igenvar_step"]))), + min_qual=list(range(config["quality_ranges"]["igenvar"]["from"], + config["quality_ranges"]["igenvar"]["to"]+1, + config["quality_ranges"]["igenvar"]["step"]))), svim = expand("results/caller_comparison/eval/SVIM/min_qual_{min_qual}/pr_rec.txt", - min_qual=list(range(config["minimums"]["svim_from"], - config["minimums"]["svim_to"]+1, - config["minimums"]["svim_step"]))), + min_qual=list(range(config["quality_ranges"]["svim"]["from"], + config["quality_ranges"]["svim"]["to"]+1, + config["quality_ranges"]["svim"]["step"]))), sniffles = expand("results/caller_comparison/eval/Sniffles/min_qual_{min_qual}/pr_rec.txt", - min_qual=list(range(config["minimums"]["sniffles_from"], - config["minimums"]["sniffles_to"]+1, - config["minimums"]["sniffles_step"]))), + min_qual=list(range(config["quality_ranges"]["sniffles"]["from"], + config["quality_ranges"]["sniffles"]["to"]+1, + config["quality_ranges"]["sniffles"]["step"]))), pbsv = expand("results/caller_comparison/eval/pbsv/min_qual_{min_qual}/pr_rec.txt", - min_qual=list(range(config["minimums"]["pbsv_from"], - config["minimums"]["pbsv_to"]+1, - config["minimums"]["pbsv_step"]))) + min_qual=list(range(config["quality_ranges"]["pbsv"]["from"], + config["quality_ranges"]["pbsv"]["to"]+1, + config["quality_ranges"]["pbsv"]["step"]))) output: igenvar = temp("results/caller_comparison/eval/igenvar.all_results.txt"), svim = temp("results/caller_comparison/eval/svim.all_results.txt"),