Build a Reproducible Bioinformatics Pipeline with Snakemake

Most bioinformatics pipelines start as a folder of shell scripts run in whatever order the author remembers. That works until you add a third sample, hand the project to a labmate, or come back to it in six months. Snakemake fixes this by making the pipeline itself the documentation: a Snakefile full of declarative rules that says exactly which inputs produce which outputs, so the tool can figure out execution order, parallelize independent jobs, and rerun only what actually changed.

This tutorial builds a small but genuinely runnable pipeline — FastQC and fastp for quality control, then BWA-MEM alignment and a sorted, indexed BAM — wired together with wildcards, a config.yaml sample sheet, and per-rule conda environments. By the end you will have a Snakefile you can dry-run, visualize as a DAG, and execute for real.

Why Snakemake for reproducible pipelines

Snakemake extends Python with a Makefile-like rule language. Each rule declares an input:, an output:, and a shell: block that turns one into the other. Snakemake reads all the rules, works backward from whatever files you ask for, and builds a dependency graph — a directed acyclic graph, or DAG — of exactly the jobs needed to produce them. If an output already exists and is newer than its inputs, that job is skipped, so fixing one sample's trimming parameters reprocesses only that sample's downstream jobs, not the whole cohort.

If you have not installed it yet, see Install Snakemake with Conda for the isolated-environment setup used throughout this guide.

Project layout and the config.yaml sample sheet

Keep the pipeline's logic separate from the data it processes. A typical layout looks like this:

project/
├── Snakefile
├── config.yaml
├── envs/
│   ├── qc.yaml
│   └── align.yaml
├── data/
│   ├── reads/
│   │   ├── sample1_R1.fastq.gz
│   │   ├── sample1_R2.fastq.gz
│   │   ├── sample2_R1.fastq.gz
│   │   └── sample2_R2.fastq.gz
│   └── reference/
│       └── genome.fasta
└── results/   # created by Snakemake as rules run

config.yaml is the sample sheet: it lists every sample and its raw FASTQ pair, plus paths shared across the whole run, such as the reference genome. Nothing sample-specific is hard-coded into the Snakefile itself.

samples:
  sample1:
    - data/reads/sample1_R1.fastq.gz
    - data/reads/sample1_R2.fastq.gz
  sample2:
    - data/reads/sample2_R1.fastq.gz
    - data/reads/sample2_R2.fastq.gz

reference: "data/reference/genome.fasta"

To add a third sample later, you edit config.yaml, not the pipeline logic. That separation is most of what "reproducible" means in practice.

The Snakefile: rule all, wildcards, and per-rule conda envs

Load the config at the top of the Snakefile with configfile:, then declare rule all. rule all has no shell: block of its own — its only job is to list, as input:, every final file the pipeline should produce. Snakemake treats that list as the default target and works backward from it.

configfile: "config.yaml"

SAMPLES = list(config["samples"].keys())

rule all:
    input:
        expand("results/qc/{sample}_R1_fastqc.html", sample=SAMPLES),
        expand("results/qc/{sample}_R2_fastqc.html", sample=SAMPLES),
        expand("results/bam/{sample}.sorted.bam.bai", sample=SAMPLES)

rule bwa_index:
    input:
        config["reference"]
    output:
        multiext(config["reference"], ".amb", ".ann", ".bwt", ".pac", ".sa")
    conda:
        "envs/align.yaml"
    shell:
        "bwa index {input}"

rule fastqc_raw:
    input:
        r1=lambda wc: config["samples"][wc.sample][0],
        r2=lambda wc: config["samples"][wc.sample][1]
    output:
        html1="results/qc/{sample}_R1_fastqc.html",
        html2="results/qc/{sample}_R2_fastqc.html"
    conda:
        "envs/qc.yaml"
    threads: 2
    shell:
        "fastqc -o results/qc -t {threads} {input.r1} {input.r2}"

rule trim_reads:
    input:
        r1=lambda wc: config["samples"][wc.sample][0],
        r2=lambda wc: config["samples"][wc.sample][1]
    output:
        r1="results/trimmed/{sample}_R1.trim.fastq.gz",
        r2="results/trimmed/{sample}_R2.trim.fastq.gz",
        report="results/trimmed/{sample}.fastp.html"
    conda:
        "envs/qc.yaml"
    threads: 2
    shell:
        "fastp -i {input.r1} -I {input.r2} "
        "-o {output.r1} -O {output.r2} "
        "-h {output.report} -w {threads}"

rule align:
    input:
        r1="results/trimmed/{sample}_R1.trim.fastq.gz",
        r2="results/trimmed/{sample}_R2.trim.fastq.gz",
        ref=config["reference"],
        idx=multiext(config["reference"], ".amb", ".ann", ".bwt", ".pac", ".sa")
    output:
        "results/bam/{sample}.sorted.bam"
    conda:
        "envs/align.yaml"
    threads: 4
    shell:
        "bwa mem -t {threads} {input.ref} {input.r1} {input.r2} "
        "| samtools sort -@ {threads} -o {output} -"

rule index_bam:
    input:
        "results/bam/{sample}.sorted.bam"
    output:
        "results/bam/{sample}.sorted.bam.bai"
    conda:
        "envs/align.yaml"
    shell:
        "samtools index {input}"

How the wildcard {sample} does the work

Notice that no rule mentions sample1 or sample2 by name. Every output: path uses {sample} as a wildcard, and expand() in rule all fills it in with each entry from SAMPLES. When Snakemake needs results/bam/sample1.sorted.bam.bai, it matches {sample} to sample1, then walks backward through index_bam, align, and trim_reads to the raw FASTQ pair looked up in config["samples"]. Add a third sample to config.yaml and the same rules produce a third set of outputs — nothing in the Snakefile changes.

Wiring QC into alignment

The dependency that actually makes this a pipeline, rather than a pile of independent scripts, is that align's input: is trim_reads's output:. Snakemake sees results/trimmed/{sample}_R1.trim.fastq.gz is both an output and an input, so it inserts an edge between the two rules in the DAG. Change fastp's trimming parameters and every downstream alignment is marked stale on the next run. fastqc_raw runs on the untrimmed reads in parallel as a sanity check only — it has no downstream consumer, so it never blocks alignment.

Per-rule conda environments

Each rule above declares a conda: directive pointing at a small environment file, so QC tools and alignment tools never share a namespace:

# envs/qc.yaml
channels:
  - conda-forge
  - bioconda
dependencies:
  - fastqc=0.12.1
  - fastp=0.23.4

# envs/align.yaml
channels:
  - conda-forge
  - bioconda
dependencies:
  - bwa=0.7.17
  - samtools=1.20

Pinning exact versions here, and channels to conda-forge/bioconda only, is what makes the pipeline portable: anyone who clones the repository gets identical tool versions, not whatever happened to be on their PATH. See Conda environments and Bioconda channels for why channel order matters.

Dry runs before you run anything real

Before executing a single job, ask Snakemake what it would do. The -n (or --dry-run) flag builds the full DAG and prints the job plan without running any shell command:

snakemake -n -p --cores 4

-p additionally prints the resolved shell command for each job, which is invaluable for spotting a malformed path before it fails halfway through a real run. For two samples, you should see something like:

Building DAG of jobs...
Job stats:
job              count
-------------  -------
align                2
all                   1
bwa_index             1
fastqc_raw            2
index_bam             2
trim_reads            2
total                10

This was a dry-run (flag -n). The order of jobs does not reflect the order of execution.

Ten jobs, not six rules — fastqc_raw, trim_reads, align, and index_bam each expand once per sample, while bwa_index and all run once regardless of sample count. A dry run that reports MissingInputException at this stage is the single most useful debugging signal Snakemake gives you, and it costs nothing to run.

Visualizing the DAG

--dag prints the fully wildcard-expanded job graph in Graphviz's dot format — one node per sample per rule, exactly the jobs a real run would execute:

snakemake --dag | dot -Tpng > dag.png

For a simpler picture with one node per rule rather than per job, use --rulegraph instead:

snakemake --rulegraph | dot -Tpng > rulegraph.png

The rule graph for this pipeline shows a straight chain — raw reads feeding both fastqc_raw and trim_reads, trim_reads feeding align, align feeding index_bam — with bwa_index joining in as a separate branch that only align depends on.

Running the pipeline for real

Once the dry run and DAG look right, drop -n and add --use-conda so each rule actually gets its declared environment:

snakemake --cores 4 --use-conda

On the first run, Snakemake solves and builds each envs/*.yaml before executing any job that needs it, so that run takes longer than later ones — the environments are cached under .snakemake/conda/ and reused after that. Interrupt a run and restart it, and Snakemake reruns only what's missing or stale, which is the entire point of describing a pipeline as inputs and outputs instead of an imperative script.

Common pitfalls

Next steps

The pattern here — config.yaml as the sample sheet, wildcards for per-sample rules, conda: for isolated per-rule dependencies — scales to real pipelines with dozens of samples and many more stages without the Snakefile itself growing more complex per sample. From here, a natural extension is swapping BWA-MEM for a splice-aware aligner and adding a counting rule, following the same rule-chaining approach used for quantification and differential expression in a standard RNA-seq workflow. The same wildcard-and-conda pattern applies regardless of which aligner or QC tool sits inside each rule.