0
$\begingroup$

I am creating a bioinformatics pipeline using Snakemake doing everying from adapter trimming of raw sequencing data in .fq to mapping, SNP calling and SNP filtering.

The issue i am having is when i try to filter specific chromosomes of the bcf files, using the --regions. The rule i have created below works and does create a filtered bcf file.

rule bcftools_filter:
    input:
        bcf = "/dir/{sample}_aligned_call.bcf",
    output:
        "/dir/{sample}_aligned_call_filter.bcf"
    threads: 20
    shell:
        "bcftools view -Ob --threads {threads} -m2 -M2 -v snps {input.bcf} -o {output}"

However i would like to add the --regions parameter to the bcftools command. I have tried the shell part of the rule to several things

bcftools view -Ob --threads {threads} -m2 -M2 -v snps {input.bcf} --regions <chr1,chr2> -o {output}
bcftools view -Ob --threads {threads} -m2 -M2 -v snps {input.bcf} --regions {chr1,chr2} -o {output}
bcftools view -Ob --threads {threads} -m2 -M2 -v snps {input.bcf} --regions {{chr1,chr2}} -o {output}
bcftools view -Ob --threads {threads} -m2 -M2 -v snps {input.bcf} --regions {{chr1} {chr2}} 
 {{chr2}} -o {output}
bcftools view -Ob --threads {threads} -m2 -M2 -v snps {input.bcf} --regions ['chr1','chr2'] -o {output}

But it doesn't work, and i am having trouble finding any suggestions to this specific issue.

I know it works if i just run bcftools filter in my command line, so using --regions should work.

Could anyone have any suggestions how to add --regions to the snakemake rule?

$\endgroup$
2
  • $\begingroup$ Have you tried expand() function? snakemake.readthedocs.io/en/stable/snakefiles/… $\endgroup$
    – zorbax
    Aug 18 at 14:58
  • $\begingroup$ No i havent tried that, i didn't realize you could potentially use the expand() function for the shell part. But i will look into it $\endgroup$ Aug 18 at 16:29
1
$\begingroup$

Something like the following would work:

rule bcftools_filter:
    input:
        bcf = "/dir/{sample}_aligned_call.bcf",
    output:
        "/dir/{sample}_aligned_call_filter.bcf"
    params:
        regions="chr1,chr2"
    threads: 20
    shell:
        "bcftools view -Ob --threads {threads} -m2 -M2 --regions {params.regions} -v snps {input.bcf} -o {output}"

You could also include the regions as part of the output file name, rather than specifying it as a param.

$\endgroup$
2
  • $\begingroup$ Would a simple --regions \"chr1,chr2\" in the shell command not work? $\endgroup$
    – Ram RS
    Aug 18 at 22:12
  • 1
    $\begingroup$ You don't need to even quote it. However having it as a param allows you to set it elsewhere (e.g., the command line). $\endgroup$
    – Devon Ryan
    Aug 18 at 22:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.