Skip to content

Commit

Permalink
fix: add rule to restrict callsets to contigs in reference (needed fo…
Browse files Browse the repository at this point in the history
…r wgs callsets where no intersection with target region is possible).
  • Loading branch information
BiancaStoecker committed Sep 12, 2024
1 parent 0f73007 commit eb3632d
Show file tree
Hide file tree
Showing 2 changed files with 24 additions and 2 deletions.
6 changes: 6 additions & 0 deletions workflow/envs/bcftools.yaml
Original file line number Diff line number Diff line change
@@ -0,0 +1,6 @@
channels:
- conda-forge
- bioconda
- nodefaults
dependencies:
- bcftools =1.20
20 changes: 18 additions & 2 deletions workflow/rules/eval.smk
Original file line number Diff line number Diff line change
Expand Up @@ -106,12 +106,28 @@ rule intersect_calls_with_target_regions:
"<(bcftools view {input.bcf}) -wa -f 1.0 -header > {output}) 2> {log}"


rule restrict_to_reference_contigs:
input:
calls="results/filtered-variants/{callset}.bcf",
calls_index="results/filtered-variants/{callset}.bcf.csi",
ref_index="resources/reference/genome.fasta.fai",
output:
"results/filtered-variants/{callset}_restricted.bcf",
log:
"logs/restrict-to-reference-contigs/{callset}.log",
conda:
"../envs/bcftools.yaml"
shell:
"bcftools view --regions $(cut -f1 {input.ref_index} | tr '\\n' ',') {input.calls} | "
"bcftools reheader -f {input.ref_index} > {output} 2> {log}"


rule normalize_calls:
input:
calls=branch(
intersect_calls,
then="results/normalized-variants/{callset}_intersected.vcf",
otherwise="results/filtered-variants/{callset}.bcf",
otherwise="results/filtered-variants/{callset}_restricted.bcf",
),
ref="resources/reference/genome.fasta",
ref_index="resources/reference/genome.fasta.fai",
Expand All @@ -122,7 +138,7 @@ rule normalize_calls:
log:
"logs/normalize-calls/{callset}.log",
conda:
"../envs/tools.yaml"
"../envs/bcftools.yaml"
shell:
"(bcftools norm {params.extra} --fasta-ref {input.ref} {input.calls} | "
"bcftools view -Oz > {output}) 2> {log}"
Expand Down

0 comments on commit eb3632d

Please sign in to comment.