Question: Use of if-else statement in snakemake rule
2
gravatar for Jokhe
22 months ago by
Jokhe90
Sweden
Jokhe90 wrote:

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!

snakemake python if-else • 2.8k views
ADD COMMENTlink modified 22 months ago • written 22 months ago by Jokhe90
2
if {wildcards.variant_status} == 'Somatic':

should be

if wildcards.variant_status == 'Somatic':

The run: statement interprets everything in the subsequent block as python code, so the {} used for templating variables in strings are no longer necessary.

ADD REPLYlink modified 22 months ago • written 22 months ago by Matt Shirley9.0k

Thank you Matt! Works like a charm!

ADD REPLYlink written 22 months ago by Jokhe90
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1622 users visited in the last hour