Bio Workflow Management Snakemake Workflows

Version Compatibility

Reference examples tested with: Snakemake 8.0+, Python 3.11+

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show <package> then help(module.function) to check signatures
• CLI: <tool> --version then <tool> --help to confirm flags

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Note: Snakemake 8 (Jan 2024) removed --cluster, --drmaa, and the *RemoteProvider classes from core and moved them to pip-installable EXECUTOR plugins (--executor slurm, package snakemake-executor-plugin-slurm) and STORAGE plugins (storage.s3(...), snakemake-storage-plugin-s3). --use-conda/--use-singularity became --software-deployment-method / --sdm conda apptainer. A Snakemake 7 command line does not run unchanged on 8/9. Run snakemake --version FIRST and branch all execution guidance on 7 vs 8/9.

Snakemake Workflows

"Build a reproducible bioinformatics pipeline with Snakemake" -> Declare each step as a rule that says "a file matching THIS output pattern is produced FROM those inputs", let the engine resolve the DAG backward from requested targets, fan out over samples with wildcards, and pin the software environment so the result reproduces next year.

• Python: Snakefile rule/checkpoint blocks with expand(), wildcards, config, resources, and conda:/container: (Snakemake)

The governing principle: Snakemake is pull/goal-oriented - it builds a STATIC DAG backward from requested target files

Snakemake is a pull, make-like engine. An author does NOT describe a forward flow of data. Each rule is a pattern-matched recipe ("a file that looks like THIS can be produced FROM that"), and the engine takes the requested target files and works BACKWARD, unifying wildcards by string-matching output filename patterns, until it reaches files already on disk. The whole plan - a static DAG - is computed at parse time, before a single job runs (Köster & Rahmann 2012 Bioinformatics 28:2520-2522). Almost every Snakemake bug a biologist hits is a downstream consequence of this one model:

• Rules are wired by OUTPUT-FILE PATTERN, not call order. A missing or typo'd output path silently drops a rule from the DAG - there is no error, the job just never runs. Debugging "why didn't it run" means tracing the backward dependency from the target, not reading top-to-bottom.
• Because the plan is fully known up front, snakemake -n (dry run), --dag, and --report are first-class. This is the payoff of the static model. Nextflow's reactive-dataflow model (processes connected by asynchronous channels, DAG emerges at runtime) has no true dry-run - hold both models in mind and most "why did/didn't it run" questions answer themselves.
• Data-dependent branching is impossible in the base model. If the NUMBER or identity of outputs is unknown until a step runs (split into one file per detected cluster, scatter over however many contigs an assembler emits), the static DAG cannot represent it -> that is exactly what CHECKPOINTS exist for. A biologist who thinks "the pipeline decides at runtime how many chunks" is fighting the paradigm and needs a checkpoint, not a clever run: block.
• A workflow manager buys reproducible LOGIC and nothing else automatically. The DAG being deterministic says nothing about tool versions. A rule with no conda:/container: runs against whatever is on $PATH; "reproducible" is unearned until the software environment is pinned (Grüning et al. 2018 Cell Syst 6:631-635). Pin containers by DIGEST and conda by LOCKFILE - see Software Deployment below.

Decision: Snakemake vs Nextflow (pick by team and infrastructure, not benchmarks)

Reuse before authoring: for a mainstream analysis (RNA-seq, variant calling, ATAC-seq), a curated community pipeline already encodes years of QC and edge cases. Adopting one means RUNNING it (e.g. nf-core/rnaseq via workflow-management/nf-core-pipelines), not authoring Groovy - so a Python-shop preference for Snakemake only decides the authoring case, not whether to build at all. Author from scratch only for a novel method or an unsupported combination of steps.

Decision: rerun triggers - why a job reran, or did not

Since Snakemake 7.8 the default is NOT pure mtime. A rerun fires on a SET of triggers: {mtime, params, input, code, software-env} (Mölder et al. 2021 F1000Research 10:33). This surprises everyone upgrading from old Snakemake.

The code trigger catches shell/script/run body changes - reformatting whitespace or editing a comment counts as a code change and reruns the job. On very large DAGs or multi-TB inputs the provenance triggers add a hashing/stat storm; --rerun-triggers mtime skips it.

Decision: run vs script vs shell vs notebook vs wrapper

run: executes IN the main Snakemake process - it shares the interpreter and GIL, cannot be conda/container-isolated (the conda: directive is disallowed with run:), blocks the scheduler, and an OOM in it takes down the whole workflow. Move anything beyond trivial glue to script:.

Decision: execution backend (Snakemake 8/9)

Porting a v7 --cluster "sbatch --account=X --partition=Y --mem=Z --time=T" command: for the native slurm executor, map those sbatch flags to resource keys (slurm_account, slurm_partition, mem_mb, runtime in minutes) set per-rule in resources: or globally via --default-resources; for a drop-in port keep the old string under cluster-generic (its own plugin, above). --cores = local cores; --jobs/-j = number of concurrent cluster/cloud jobs (in v8 these are separate). Profiles are versioned: the file is config/config.v8+.yaml, every long option becomes a YAML key.

Rules, wildcards, and expand

expand() returns a LIST of strings by combinatorial substitution - it does NOT touch the filesystem. Use it to enumerate targets in rule all.

configfile: 'config/config.yaml'
SAMPLES = config['samples']

rule all:                                              # the requested targets; the DAG is built backward from here
    input:
        expand('results/{sample}.bam', sample=SAMPLES)

rule align:
    input:
        r1 = 'data/{sample}_R1.fq.gz',
        r2 = 'data/{sample}_R2.fq.gz',
        index = 'ref/genome.fa'
    output:
        bam = 'aligned/{sample}.bam'                   # this OUTPUT PATTERN, matched against the target, wires the rule in
    threads: 8
    log:
        'logs/align/{sample}.log'
    shell:
        'bwa mem -t {threads} {input.index} {input.r1} {input.r2} | '
        'samtools sort -@ {threads} -o {output.bam} 2> {log}'

Wildcard constraints - stop silent mis-routing

Wildcards are greedy regex string-unification ({sample} compiles to .+), not typed parameters. An unconstrained wildcard swallows path separators and adjacent tokens: data/{sample}.txt matches data/a/b.txt as sample=a/b, and {a}.{b}.txt on 101.B.normal.txt has no unique parse. The failure is silent mis-routing, not an error. Constrain whenever a value can contain /, ., or _, or a filename has multiple variable tokens.

wildcard_constraints:
    sample = '[^/]+',                                  # no path separators
    chrom = r'\d+|X|Y|MT'                              # only real chromosome tokens

# two rules whose output patterns can both produce a requested file raise AmbiguousRuleException;
# prefer non-overlapping constraints to disambiguate, and fall back to `ruleorder: a > b` only if needed.