Skip to content

bgzip strelka region bed outside container #39

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 1 commit into
base: codex/add-snakemake-rule-for-strelka2-variant-caller
Choose a base branch
from
Open

bgzip strelka region bed outside container #39

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
90 changes: 70 additions & 20 deletions workflow/rules/strelka2.smk
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,11 +22,44 @@ rule strelka2_germline_chunkdirs:
"""


rule strelka2_germline_region_bed:
input:
ref_fa=lambda wc: config["supporting_files"]["files"]["huref"]["fasta"]["name"],
output:
bed=MDIR + "{sample}/align/{alnr}/snv/slk2g/work/{sample}.germline.{strelkachrm}/region.bed.gz",
tbi=MDIR + "{sample}/align/{alnr}/snv/slk2g/work/{sample}.germline.{strelkachrm}/region.bed.gz.tbi",
threads: 1
conda: "../envs/vanilla_v0.1.yaml"
log:
MDIR + "{sample}/align/{alnr}/snv/slk2g/vcfs/{strelkachrm}/log/{sample}.{alnr}.strelka2.{strelkachrm}.regionbed.log",
params:
schrm=get_strelka_chrm_day,
cpre="" if "b37" == config['genome_build'] else "chr",
mito_code="MT" if "b37" == config['genome_build'] else "M",
shell:
r"""
set -euo pipefail
mkdir -p $(dirname {output.bed}) >> {log} 2>&1

vchr=$(echo {params.cpre}{params.schrm} | sed 's/~/\:/g' | sed 's/23\:/X\:/' | sed 's/24\:/Y\:/' | sed 's/25\:/{params.mito_code}\:/' )
vchr=${vchr%:}
Comment on lines +44 to +45
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[P0] Escape braces in region bed shell commands

The new region‑bed rules use Bash parameter expansion (vchr=${vchr%:}) without escaping the braces for Snakemake’s .format call. When Snakemake renders the rule, it treats {vchr%:} as a placeholder and raises KeyError: 'vchr%', so neither the germline nor somatic Strelka jobs can run. The line should use doubled braces, e.g. vchr=${{vchr%:}}, and the same fix is needed in the somatic rule around the corresponding lines.

Useful? React with 👍 / 👎.

IFS=':' read -r vcontig vstart vend <<< "$vchr"
if [ -z "${vend:-}" ]; then
vstart=0
vend=$(awk -v c="$vcontig" '$1==c{print $2; exit}' {input.ref_fa}.fai)
fi
echo -e "$vcontig\t$vstart\t$vend" | bgzip -c 2>> {log} > {output.bed}
tabix -f -p bed {output.bed} >> {log} 2>&1
"""


rule strelka2_germline:
input:
cram=MDIR + "{sample}/align/{alnr}/{sample}.{alnr}.cram",
crai=MDIR + "{sample}/align/{alnr}/{sample}.{alnr}.cram.crai",
ref_fa=lambda wc: config["supporting_files"]["files"]["huref"]["fasta"]["name"],
region=MDIR + "{sample}/align/{alnr}/snv/slk2g/work/{sample}.germline.{strelkachrm}/region.bed.gz",
region_tbi=MDIR + "{sample}/align/{alnr}/snv/slk2g/work/{sample}.germline.{strelkachrm}/region.bed.gz.tbi",
d=MDIR + "{sample}/align/{alnr}/snv/slk2g/vcfs/{strelkachrm}/{sample}.ready",
output:
vcfgz=MDIR + "{sample}/align/{alnr}/snv/slk2g/vcfs/{strelkachrm}/{sample}.{alnr}.strelka2.{strelkachrm}.germline.vcf.gz",
Expand All @@ -52,19 +85,10 @@ rule strelka2_germline:
set -euo pipefail
mkdir -p {params.run_dir}

vchr=$(echo {params.cpre}{params.schrm} | sed 's/~/\:/g' | sed 's/23\:/X\:/' | sed 's/24\:/Y\:/' | sed 's/25\:/{params.mito_code}\:/' )
vchr=${{vchr%:}}
IFS=':' read -r vcontig vstart vend <<< "$vchr"
if [ -z "${{vend:-}}" ]; then
vstart=0
vend=$(awk -v c="$vcontig" '$1==c{{print $2; exit}}' {input.ref_fa}.fai)
fi
echo -e "$vcontig\t$vstart\t$vend" > {params.run_dir}/region.bed

/opt/strelka/bin/configureStrelkaGermlineWorkflow.py \
--bam {input.cram} \
--referenceFasta {input.ref_fa} \
--callRegions {params.run_dir}/region.bed \
--callRegions {input.region} \
--runDir {params.run_dir} >> {log} 2>&1
{params.run_dir}/runWorkflow.py -m local -j {threads} >> {log} 2>&1
cp {params.run_dir}/results/variants/variants.vcf.gz {output.vcfgz}
Expand Down Expand Up @@ -125,6 +149,39 @@ rule strelka2_somatic_chunkdirs:
"""


rule strelka2_somatic_region_bed:
wildcard_constraints:
sample=TUMORS_REGEX
input:
ref_fa=lambda wc: config["supporting_files"]["files"]["huref"]["fasta"]["name"],
output:
bed=MDIR + "{sample}/align/{alnr}/snv/slk2s/work/{sample}.somatic.{strelkachrm}/region.bed.gz",
tbi=MDIR + "{sample}/align/{alnr}/snv/slk2s/work/{sample}.somatic.{strelkachrm}/region.bed.gz.tbi",
threads: 1
conda: "../envs/vanilla_v0.1.yaml"
log:
MDIR + "{sample}/align/{alnr}/snv/slk2s/vcfs/{strelkachrm}/log/{sample}.{alnr}.strelka2.{strelkachrm}.regionbed.log",
params:
schrm=get_strelka_chrm_day,
cpre="" if "b37" == config['genome_build'] else "chr",
mito_code="MT" if "b37" == config['genome_build'] else "M",
shell:
r"""
set -euo pipefail
mkdir -p $(dirname {output.bed}) >> {log} 2>&1

vchr=$(echo {params.cpre}{params.schrm} | sed 's/~/\:/g' | sed 's/23\:/X\:/' | sed 's/24\:/Y\:/' | sed 's/25\:/{params.mito_code}\:/' )
vchr=${vchr%:}
IFS=':' read -r vcontig vstart vend <<< "$vchr"
if [ -z "${vend:-}" ]; then
vstart=0
vend=$(awk -v c="$vcontig" '$1==c{print $2; exit}' {input.ref_fa}.fai)
fi
echo -e "$vcontig\t$vstart\t$vend" | bgzip -c 2>> {log} > {output.bed}
tabix -f -p bed {output.bed} >> {log} 2>&1
"""


rule strelka2_somatic:
wildcard_constraints:
sample=TUMORS_REGEX
Expand All @@ -134,6 +191,8 @@ rule strelka2_somatic:
normal_cram=get_somcall_normal_cram,
normal_crai=get_somcall_normal_crai,
ref_fa=lambda wc: config["supporting_files"]["files"]["huref"]["fasta"]["name"],
region=MDIR + "{sample}/align/{alnr}/snv/slk2s/work/{sample}.somatic.{strelkachrm}/region.bed.gz",
region_tbi=MDIR + "{sample}/align/{alnr}/snv/slk2s/work/{sample}.somatic.{strelkachrm}/region.bed.gz.tbi",
d=MDIR + "{sample}/align/{alnr}/snv/slk2s/vcfs/{strelkachrm}/{sample}.ready",
output:
snv=MDIR + "{sample}/align/{alnr}/snv/slk2s/vcfs/{strelkachrm}/{sample}.{alnr}.strelka2.{strelkachrm}.somatic.snvs.vcf.gz",
Expand Down Expand Up @@ -161,20 +220,11 @@ rule strelka2_somatic:
set -euo pipefail
mkdir -p {params.run_dir}

vchr=$(echo {params.cpre}{params.schrm} | sed 's/~/\:/g' | sed 's/23\:/X\:/' | sed 's/24\:/Y\:/' | sed 's/25\:/{params.mito_code}\:/' )
vchr=${{vchr%:}}
IFS=':' read -r vcontig vstart vend <<< "$vchr"
if [ -z "${{vend:-}}" ]; then
vstart=0
vend=$(awk -v c="$vcontig" '$1==c{{print $2; exit}}' {input.ref_fa}.fai)
fi
echo -e "$vcontig\t$vstart\t$vend" > {params.run_dir}/region.bed

/opt/strelka/bin/configureStrelkaSomaticWorkflow.py \
--tumorBam {input.tumor_cram} \
--normalBam {input.normal_cram} \
--referenceFasta {input.ref_fa} \
--callRegions {params.run_dir}/region.bed \
--callRegions {input.region} \
--runDir {params.run_dir} >> {log} 2>&1
{params.run_dir}/runWorkflow.py -m local -j {threads} >> {log} 2>&1
cp {params.run_dir}/results/variants/somatic.snvs.vcf.gz {output.snv}
Expand Down