Skip to content

Commit

Permalink
Fix: Split intersection with target region and normalize calls in sep…
Browse files Browse the repository at this point in the history
…arate rules.
  • Loading branch information
BiancaStoecker committed Jun 27, 2024
1 parent dce88d8 commit b5ecb09
Show file tree
Hide file tree
Showing 2 changed files with 29 additions and 6 deletions.
9 changes: 8 additions & 1 deletion workflow/rules/common.smk
Original file line number Diff line number Diff line change
Expand Up @@ -93,7 +93,7 @@ def get_plot_cov_labels(): # TODO check if ever used anywhere
def label(name):
lower, upper = get_cov_interval(name)
if upper:
return f"{lower}-{upper - 1}"
return f"{lower}-{upper-1}"
return f"≥{lower}"

return {name: label(name) for name in low_coverages}
Expand Down Expand Up @@ -260,6 +260,13 @@ def get_target_regions(wildcards):
return []


def intersect_calls(wildcards):
if get_target_regions(wildcards) == []:
return False
else:
return True


def get_target_regions_intersect_statement(wildcards, input):
if input.target:
return f"bedtools intersect -a /dev/stdin -b {input.target} |"
Expand Down
26 changes: 21 additions & 5 deletions workflow/rules/eval.smk
Original file line number Diff line number Diff line change
Expand Up @@ -91,12 +91,30 @@ rule remove_non_pass:
"v3.3.6/bio/bcftools/view"


rule normalize_calls:
rule intersect_calls_with_target_regions:
input:
bcf="results/filtered-variants/{callset}.bcf",
regions=get_target_regions,
output:
"results/normalized-variants/{callset}_intersected.vcf.gz",
log:
"logs/intersect-calls/{callset}.log",
conda:
"../envs/tools.yaml"
shell:
"(bedtools intersect -b {input.regions} -a "
"<(bcftools view {input.bcf}) -wa -f 1.0 -header > {output}) 2> {log}"


rule normalize_calls:
input:
calls=branch(
intersect_calls,
then="results/normalized-variants/{callset}_intersected.vcf.gz",
otherwise="results/filtered-variants/{callset}.bcf",
),
ref="resources/reference/genome.fasta",
ref_index="resources/reference/genome.fasta.fai",
regions=get_target_regions,
output:
"results/normalized-variants/{callset}.vcf.gz",
params:
Expand All @@ -106,9 +124,7 @@ rule normalize_calls:
conda:
"../envs/tools.yaml"
shell:
"(bedtools intersect -b {input.regions} -a "
"<(bcftools view {input.bcf}) -wa -f 1.0 -header | "
"bcftools norm {params.extra} --fasta-ref {input.ref} | "
"(bcftools norm {params.extra} --fasta-ref {input.ref} {input.calls} | "
"bcftools view -Oz > {output}) 2> {log}"


Expand Down

0 comments on commit b5ecb09

Please sign in to comment.