Hi!
I'm building snakemake pipeline for VarScan2 variant calling. As a result of VarScan2, I have six VCF files that must be processed individually and bam-readcount tool is a part of this processing. The use of bam-readcount depend on the type of variant; if I have somatic variant calls I should use tumor BAM file as bam-readcount's input and reference BAM for germline and LOH variants correspondingly.
It is possible to do three separate rules for each status (germline, somatic, LOH) of variants when pipeline is working properly. However, I also tried to build if-else statement for my rule so that rule could select the right input file on the basis of wildcards. In other words: if wildcard variant_status is 'Somatic', use tumor BAM as input. If wildcard is not 'Somatic', use reference BAM instead. Seems to be easy but for some reason my if-else statement is not working; it gives FALSE result for all values of wildcard variant_status and therefore uses reference BAM for all input position files. I'm interested to know how this kind of statement could be used in snakemake as it would make pipeline more simply to follow.
# PART OF PIPELINE, SOME RULES ARE REMOVED TO SIMPLIFY POST
PATIENTS=['patientA', 'patientB', 'patientC'...]
TYPES=['snp', 'indel']
STATUSES=['Germline', 'Somatic', 'LOH']
rule all:
expand('{patient}/{patient}.{variant_type}.{variant_status}.readcounts', patient=PATIENTS,
variant_type=TYPES,
variant_status=STATUSES)
rule bam_readcount:
input:
positions='{patient}/{patient}.{variant_type}.{variant_status}.positions',
tumorbam='patient}/{patient}_tumor.bam',
referencebam='patient}/{patient}_reference.bam'
output:
'{patient}/{patient}.{variant_type}.{variant_status}.readcounts'
run:
if {wildcards.variant_status} == 'Somatic':
shell("bam-readcount -w1 -f {REF} {input.tumorbam} -l {input.positions} > {output}")
print("Tumor BAM used as input!")
else:
shell("bam-readcount -w1 -f {REF} {input.referencebam} -l {input.positions} > {output}")
print("Reference BAM used as input!")
I have read few tutorials about the use of if-else statements in snakemake (and Python in overall) and tried different kind of approaches to do if-else statements but none of them has worked so far and I believe that I'm missing something very obvious here. I'm learnign Python alongside my pipeline constructing, but could someone more experienced friendly point what could be wrong with my if-else statement and how could I fix it? Thank you in advance!
should be
The
run:statement interprets everything in the subsequent block as python code, so the{}used for templating variables in strings are no longer necessary.Thank you Matt! Works like a charm!