Bio applied workflow engines

Workflow Engines: Snakemake and Nextflow/nf-core

When to Use

• Turning a tangle of shell scripts into a reproducible, resumable pipeline with dependency tracking (FASTQ → QC → align → call → annotate).
• Writing or debugging a Snakefile (rules, wildcards, configfile:, conda:, container:, checkpoints).
• Writing or debugging a Nextflow DSL2 script (process, channel, workflow blocks).
• Running/customizing a curated nf-core pipeline (nf-core/rnaseq, nf-core/sarek, nf-core/atacseq, ...) and building its CSV samplesheet.
• Scaling a pipeline to SLURM, AWS Batch, or Google Batch/Life Sciences.
• Deciding between Snakemake and Nextflow for a new project.

Version Compatibility

• Snakemake ≥ 8.0 (executor-plugin architecture: --executor slurm, not the old --cluster)
• Nextflow ≥ 24.04, DSL2 is default (no nextflow.enable.dsl=2 needed, but harmless to keep)
• nf-core tools ≥ 2.14, Python ≥ 3.10
• Java ≥ 17 (required by Nextflow)

Prerequisites

• Comfortable with bash scripting and YAML.
• A container or conda toolchain available (Docker, Singularity/Apptainer, or conda/mamba) — both engines assume per-step environment isolation.
• Know the tools being wrapped (bwa, fastp, GATK, etc.) — the engine only orchestrates them, it does not replace understanding the commands.

Snakemake: rules, wildcards, config

Goal: replace a linear shell script with a DAG that Snakemake resumes on failure and parallelizes across samples. Approach: define rule all listing final targets; each other rule declares input/output with {wildcards}; Snakemake infers execution order backward from rule all by matching output patterns.

# Snakefile
configfile: "config/config.yaml"
SAMPLES = config["samples"]

rule all:
    input:
        expand("results/vcf/{sample}.vcf.gz", sample=SAMPLES),
        "results/qc/multiqc_report.html"

rule fastqc:
    input:  "data/{sample}_R1.fastq.gz"
    output: "results/qc/{sample}_fastqc.html"
    conda:  "envs/qc.yaml"
    log:    "logs/fastqc/{sample}.log"
    shell:  "fastqc {input} --outdir results/qc/ 2> {log}"

rule bwa_mem:
    input:
        r1  = "trimmed/{sample}_R1.fq.gz",
        r2  = "trimmed/{sample}_R2.fq.gz",
        ref = config["reference"]
    output: "aligned/{sample}.bam"
    params:
        rg = r"@RG\tID:{sample}\tSM:{sample}\tPL:ILLUMINA"
    log:       "logs/bwa_mem/{sample}.log"
    benchmark: "benchmarks/bwa_mem/{sample}.tsv"
    threads: 8
    resources: mem_mb = 16000
    shell:
        "(bwa mem -t {threads} -R '{params.rg}' {input.ref} "
        "{input.r1} {input.r2} | samtools sort -o {output}) 2> {log}"

# config/config.yaml
samples: [sample1, sample2, sample3]
reference: ref/hg38.fa

snakemake -n --cores 8            # dry-run: preview the job DAG first
snakemake --use-conda --cores 8   # real run, per-rule conda envs
snakemake --dag | dot -Tpdf > dag.pdf

Checkpoint for dynamic output counts

checkpoint split_by_chromosome:
    input:  "assembly.fa"
    output: directory("chromosomes/")
    shell:  "csplit assembly.fa /^>/ '{{*}}' -f chromosomes/ -b '%d.fa'"

def aggregate_chromosomes(wildcards):
    out = checkpoints.split_by_chromosome.get(**wildcards).output[0]
    chroms = glob_wildcards(f"{out}/{{chrom}}.fa").chrom
    return expand("annotated/{chrom}.gff3", chrom=chroms)

rule all:
    input: aggregate_chromosomes

Validating rule/wildcard matching before you run

Goal: catch a bad output pattern or unmatched target before burning cluster time on a dry-run. Approach: Snakemake matches a target filename against each rule's output pattern by turning {wildcard} into a regex group — reproduce that logic to sanity-check rule design.

import re

def match_rule_output(output_pattern: str, target_file: str) -> dict | None:
    """
    Check whether target_file matches a Snakemake-style output pattern
    (e.g. "results/vcf/{sample}.vcf.gz") and return the wildcard values.
    Mirrors Snakemake's own {name} -> (?P<name>[^/]+) substitution.
    """
    regex = re.sub(r"\{(\w+)\}", r"(?P<\1>[^/]+)", re.escape(output_pattern))
    regex = regex.replace(r"\/", "/")
    m = re.fullmatch(regex, target_file)
    return m.groupdict() if m else None

rules = [
    ("fastqc",  "results/qc/{sample}_fastqc.html"),
    ("bwa_mem", "aligned/{sample}.bam"),
    ("call",    "results/vcf/{sample}.vcf.gz"),
]
for target in ["aligned/sample2.bam", "some/unknown/file.txt"]:
    hit = next((n for n, p in rules if match_rule_output(p, target)), None)
    print(f"{target:<35} -> {hit or 'NO MATCH (fix pattern or target path)'}")

Nextflow DSL2: processes, channels, workflow

Goal: wire tool invocations together with explicit data channels instead of Snakemake's filename-pattern matching. Approach: a process declares input/output and a script: block; a workflow {} block connects processes via Channel objects.

#!/usr/bin/env nextflow
nextflow.enable.dsl = 2

params.reads     = "data/*_{R1,R2}.fastq.gz"
params.reference = "ref/hg38.fa"
params.outdir    = "results"

process BWA_MEM {
    tag "$sample_id"
    cpus 8
    memory '16 GB'
    conda "bioconda::bwa=0.7.17 bioconda::samtools=1.19"
    publishDir "${params.outdir}/aligned", mode: 'copy'

    input:
        tuple val(sample_id), path(reads)
        path  ref
    output:
        tuple val(sample_id), path("${sample_id}.bam")
    script:
    """
    bwa mem -t ${task.cpus} ${ref} ${reads[0]} ${reads[1]} \
        | samtools sort -o ${sample_id}.bam
    """
}

workflow {
    reads_ch = Channel.fromFilePairs(params.reads)
    BWA_MEM(reads_ch, file(params.reference))
}

nextflow run main.nf -profile docker --outdir results/
nextflow run main.nf -resume            # hash-based incremental re-run

nf-core: running curated pipelines + samplesheet validation

Goal: run a production-grade, CI-tested pipeline (e.g. nf-core/rnaseq, nf-core/sarek) instead of writing one from scratch, and catch samplesheet errors before Nextflow does. Approach: install nf-core tools, build the CSV samplesheet the pipeline expects, validate it locally, then launch with a profile (docker/singularity/conda).

pip install nf-core
nf-core list                                   # browse available pipelines
nf-core launch nf-core/rnaseq                  # interactive samplesheet + params builder

nextflow run nf-core/rnaseq \
    -profile docker -r 3.14.0 \
    --input samplesheet.csv --genome GRCh38 --outdir results/

import re
from dataclasses import dataclass
from typing import Optional

@dataclass
class RNAseqSample:
    """One row of an nf-core/rnaseq samplesheet."""
    sample: str
    fastq_1: str
    fastq_2: Optional[str] = None
    strandedness: str = "auto"

    def validate(self) -> list[str]:
        """Return a list of validation errors (empty list = valid row)."""
        errors = []
        if not re.match(r"^[A-Za-z0-9_\-]+$", self.sample):
            errors.append(f"sample name '{self.sample}' has invalid characters")
        if not self.fastq_1.endswith((".fastq.gz", ".fq.gz")):
            errors.append("fastq_1 must end with .fastq.gz or .fq.gz")
        if self.strandedness not in ("auto", "forward", "reverse", "unstranded"):
            errors.append(f"strandedness '{self.strandedness}' is invalid")
        return errors

samples = [
    RNAseqSample("CTRL_1", "data/ctrl_1_R1.fastq.gz", "data/ctrl_1_R2.fastq.gz"),
    RNAseqSample("bad sample!", "data/bad.txt"),
]
for s in samples:
    errs = s.validate()
    print(f"{s.sample}: {'OK' if not errs else '; '.join(errs)}")

Cluster / cloud execution

# Snakemake — SLURM via executor plugin
pip install snakemake-executor-plugin-slurm
snakemake --executor slurm --jobs 100 --default-resources slurm_account=my_account --use-conda

# Snakemake — Google Batch / AWS Batch
pip install snakemake-executor-plugin-googlebatch
snakemake --executor googlebatch --default-storage-provider gcs \
    --default-storage-prefix gs://my-bucket/results --jobs 200

# Nextflow — set executor in nextflow.config, then just run normally
# process.executor = 'slurm'; process.queue = 'long'
nextflow run main.nf -profile cluster

Pitfalls

• Snakemake runs from the Snakefile's directory: use relative paths or workflow.basedir for portability.
• Wildcards are greedy by default: {sample} matches slashes; restrict with {sample,[^/]+}.
• rule all must list every final output: Snakemake builds the DAG backward from this target only.
• Missing log: directive: without it, a failed job leaves no captured stderr — always add one.
• benchmark: is overwritten per run: include {sample} in the path or successive runs clobber the TSV.
• Nextflow channels are consumed once: reusing a channel across processes needs .multiMap{} or duplicating the source, or the second consumer gets nothing.
• nf-core samplesheets are pipeline-specific: column names/order vary between rnaseq and sarek — validate with nf-core launch or the pipeline's schema, don't assume one CSV shape fits all.
• Coordinate systems inside custom rules: BED is 0-based half-open, VCF/GFF are 1-based inclusive — mixing them causes off-by-one errors in scripts a rule/process shells out to.

See Also

• bio-workflow-management-snakemake-workflows
• bio-workflow-management-nextflow-pipelines
• bio-workflow-management-cwl-workflows
• bio-workflows-fastq-to-variants