1
$\begingroup$

I am learning to use snakemake. I am following the tutorial (https://snakemake.readthedocs.io/en/stable/tutorial/tutorial.html). In the section called "Basics", steps 1 to 4 worked well, processing one sample at a time, from mapping reads until indexing alignments.

snakemake --cores 2 sorted_reads/{A,B,C}.bam.bai

However, when I add the rule for joint variant calling (step 5):

snakemake --cores 2

I get the following error:

"WorkflowError:
Target rules may not contain wildcards. Please specify concrete files or a rule without wildcards at the command line, or have a rule without wildcards at the very top of your workflow (e.g. the typical "rule all" which just collects all results you want to generate in the end)."

It looks like I need to do this "rule all" step, any idea how to set this up?

This is how the snakemake file looks like:

SAMPLES = ["A", "B", "C"]

rule bwa_map:
    input:
        "data/genome.fa",
        "data/samples/{sample}.fastq"
    output:
        "mapped_reads/{sample}.bam"
    shell:
        "bwa mem {input} | samtools view -Sb - > {output}"

rule samtools_sort:
    input:
        "mapped_reads/{sample}.bam"
    output:
        "sorted_reads/{sample}.bam"
    shell:
        "samtools sort -T sorted_reads/{wildcards.sample} "
        "-O bam {input} > {output}"

rule samtools_index:
    input:
        "sorted_reads/{sample}.bam"
    output:
        "sorted_reads/{sample}.bam.bai"
    shell:
        "samtools index {input}"

rule bcftools_call:
    input:
        fa="data/genome.fa",
        bam=expand("sorted_reads/{sample}.bam", sample=SAMPLES),
        bai=expand("sorted_reads/{sample}.bam.bai", sample=SAMPLES)
    output:
        "calls/all.vcf"
    shell:
        "bcftools mpileup -f {input.fa} {input.bam} | "
        "bcftools call -mv - > {output}"
$\endgroup$

1 Answer 1

2
$\begingroup$

Read the snakemake documentation. all is a pseudo-rule that accepts a set of one or more final output files as input and does nothing with them.

$\endgroup$
2
  • $\begingroup$ One clarification: Snakemake defaults to trying the first rule in the file (whether called the idiomatic all or anything else) and since in your case that has a sample wildcard, it doesn't know what to fill in for that wildcard. (A typical all rule might use something like input: expand("mapped_reads/{sample}.bam", sample=SAMPLES) for example.) $\endgroup$
    – Jesse
    Jun 3, 2022 at 16:06
  • $\begingroup$ Even more to the point: move bcftools_call to the top of the file and it will behave as your all rule. $\endgroup$
    – Jesse
    Jun 3, 2022 at 19:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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