Pipelines and reproducibility

≈ 15 min

Turn a variant-calling analysis into a pipeline anyone can re-run with one command – with Snakemake and with Galaxy – and find out what an AI agent can and cannot prove about its results.

Claerbout’s principle

“An article about computational science in a scientific publication is not the scholarship itself, it is merely advertising of the scholarship. The actual scholarship is the complete software development environment and the complete set of instructions which generated the figures.”

J. B. Buckheit and D. L. Donoho (1995), WaveLab and reproducible research, summarising the principle of Jon Claerbout of the Stanford Exploration Project.

At the Stanford Exploration Project, Jon Claerbout and Martin Karrenbach set four basic rules for reproducible research: every result can be burned (deleted), built again, viewed, and its folder cleaned of intermediate files – each with one command. By the end of this practical your own analysis will pass the same test.

In the reads-to-variants practical you turned sequencing reads into variant calls, typing one command at a time. Could someone else – or you, in six months – get exactly the same calls again? Only if they know which data you used, every command and option, in which order, with which settings, and with which versions of the programs. A pipeline writes all of that down in a form a computer can run.

0.1 The route

Chapter	You will…
1 · A project folder	set up a project the data-science way: raw data read-only and checksummed, a README
2 · By hand, then a script	run the analysis by hand, then as a bash script – and see what a script cannot do
3 · Python for Snakemake	learn the Python you need to read and write a Snakefile, in the notebook and as a script
4 · Your first Snakefile	write the pipeline as Snakemake rules, with help (and mistakes) from an AI assistant
5 · Any sample, any setting	wildcards, a config file, parameters, logs and a Python script in the pipeline
6 · Software environments	find out why the same pipeline can give different numbers on another computer, and pin the software
7 · The same pipeline in Galaxy	build the workflow graphically in Galaxy – and get the identical result
8 · An AI agent	let an AI agent “do the analysis” and ask what it can prove
9 · Which is best?	burn and rebuild your results, compare the four approaches, and export your project

The data

The same sample and genes as in the reads-to-variants practical, made small enough for every step to run in seconds in your browser:

Reads: 3,519 pairs of real 76-base Illumina exome reads from NA12878 (1000 Genomes run SRR098401) that map around CYP2C19 and CYP2C9, two genes that metabolise many drugs.
Reference: two 26–30 kb slices of the human genome (hg19), each covering part of one gene – including the sites of CYP2C19*2 and CYP2C9*2.
Programs: minimap2 2.22, samtools 1.17 and bcftools 1.10 – the real programs, compiled to WebAssembly and running in your browser. (minimap2 replaces BWA, which is not available in a browser.)

NA12878 is a public reference sample whose donor consented to open data sharing. Nothing in this practical is a clinical result.

0.2 The workbench

Activity 0.2 · Look around

The Terminal is a Linux-like shell. List the course data: ls -l /data/course
Read what the data is: cat /data/course/README.md
The Files tab shows your folders and has an editor for Snakefiles, scripts and config files.
The Notebook runs Python in your browser (the first time, it takes 10–30 seconds to start).
Galaxy is a practice Galaxy server running in the page (chapter 7).
The AI assistant has a Chat tab (a coding assistant) and an Agent tab (chapter 8). Its prepared answers contain mistakes on purpose – spotting them is part of the practical.

The practical takes about 5½ hours – for example chapters 0–5 in one session and 6–9 in another; stop whenever you like and carry on later. Your own files (Snakefile, config, scripts …), conda environments, Galaxy workflows, the notebook, ticks and answers are saved in this browser – use the same browser next time. To continue on another computer (lab computers may clear the browser when you log out), use My answers → Save all my work to a file before you leave, and Open my work from a file on the next computer. Results are not saved: rebuilding them is one command, which is rather the point. My answers also downloads your answers; the last chapter shows how to download your whole project.

Think back to the reads-to-variants practical. If a colleague wanted to regenerate your filtered VCF file exactly, what would they need from you? List as many things as you can.

At least: the raw data – the exact read files and the exact reference genome (version and any changes) – ideally with checksums; every command, in order, with all its options; the parameters (e.g. the filter thresholds); the program versions (and the libraries they were built with); the computing environment (operating system, how the software was installed); anything done by hand (renamed files, a step repeated after an error, a file edited in a text editor). A shell history or lab-book notes rarely record all of this – that is what this practical fixes.

In Claerbout’s principle, a paper is only “advertising”. What is the “scholarship”?

the figures and tables in the paper the complete software environment and the complete set of instructions that produced the results the methods section of the paper the peer reviews

A project folder

≈ 20 min

Reproducibility starts before the first command: with a folder in which data, code, settings and results each have their place – and raw data that cannot change unnoticed.

Before the first command, decide where everything goes. A good layout keeps apart what you were given (data), what you write (code and settings) and what the computer makes (results) – so that anyone can see at a glance what can be regenerated and what cannot.

1.1 Make the folders

Activity 1.1 · A project skeleton

Make sure you are in your home folder: cd ~
Make the folders in one go: mkdir -p cyp2c19-pipeline/{config,data/raw,logs,results,workflow/{envs,scripts}}
The braces are brace expansion: the shell turns {a,b} into two words, so one command makes the project folder and all eight folders inside it.
Go into the project and look at it: cd cyp2c19-pipeline && tree

Folder	What goes in it	Can it be deleted and made again?
`data/raw/`	the input data, exactly as received	No – protect it
`config/`	settings: sample names, file paths, thresholds	no – you write it
`workflow/`	the pipeline (`Snakefile`), software environments (`envs/`), scripts (`scripts/`)	no – you write it
`results/`, `logs/`	everything the pipeline makes	Yes – so never edit anything in here by hand

This is close to the layout the Snakemake documentation recommends, and to W. S. Noble’s classic guide, A quick guide to organizing computational biology projects (PLoS Comput Biol 2009).

1.2 Bring in the data – and prove it is the same data

Activity 1.2 · Copy and checksum

Copy the reads, giving them names that say which sample they are: cp /data/course/SRR098401_1.fastq data/raw/NA12878_R1.fastq && cp /data/course/SRR098401_2.fastq data/raw/NA12878_R2.fastq
Copy the reference and its index: cp /data/course/hg19_CYP2C_slices.fa data/raw/reference.fa && cp /data/course/hg19_CYP2C_slices.fa.fai data/raw/reference.fa.fai
Record a checksum of each file: md5sum data/raw/NA12878_R1.fastq data/raw/NA12878_R2.fastq data/raw/reference.fa > data/raw/MD5SUMS
Compare your checksums with the ones published with the course data: cat data/raw/MD5SUMS /data/course/MD5SUMS
Anyone can now check that their copies are identical to yours: md5sum -c data/raw/MD5SUMS

The reads now have new names. How do the checksums show that data/raw/NA12878_R1.fastq is still exactly the file SRR098401_1.fastq?

Both files have the same MD5 checksum, d2b8ddb488975fa5ef18e6e59edb5681. A checksum is calculated from the contents only – not the name or the date – and changing even one byte gives a completely different checksum. Same checksum = same data, whatever the file is called.

1.3 Protect the raw data

Activity 1.3 · Read-only

Take away write permission from everyone (a-w): chmod a-w data/raw/*
Look at the permissions: ls -l data/raw – -r--r--r-- means readable, not writable. (Compare ls -l /data/course: those files belong to root, so you could never change them.)
Try an accident – writing over the reference: echo oops > data/raw/reference.fa

Why protect the raw data, but not the results?

Everything else is derived from the raw data: if it changes, every result and conclusion silently loses its foundation, and the data often cannot be obtained again. Results, on the other hand, are disposable: the pipeline can make them again at any time – so they are never edited by hand, and they need no protection. (Write protection prevents accidents; as the owner you could still change the permissions back.)

1.4 Say what it is: a README

The README is the first thing anyone (including you, later) reads. Start it now and complete it as the project grows.

README.md

# CYP2C19 / CYP2C9 variant calling for NA12878

## What this project does
Maps short exome reads of sample NA12878 to two slices of the human
reference genome (hg19) that cover parts of the genes CYP2C19 and CYP2C9,
calls the variants and summarises them.

## Data (data/raw – read-only, never edited)
| File | Where it came from |
| --- | --- |
| NA12878_R1.fastq | /data/course/SRR098401_1.fastq – run SRR098401, 1000 Genomes |
| NA12878_R2.fastq | /data/course/SRR098401_2.fastq |
| reference.fa (+ .fai) | /data/course/hg19_CYP2C_slices.fa – hg19 chr10:96,530,001-96,560,000 and chr10:96,699,001-96,725,000 |

Checksums: data/raw/MD5SUMS – check them with  md5sum -c data/raw/MD5SUMS

## How to run
(to be written when the pipeline exists)

## Author and date
(your name, today's date)

Activity 1.4 · README and .gitignore

Press Create this file above. It opens in the Files tab: fill in your name and the date, and save with Ctrl+S.
Real projects are kept under version control with git, so every change to the code is recorded. Results are not put in git – they can be regenerated. Create the .gitignore file below that says so (git is not available in this browser, but the file will be part of your project when you download it).

.gitignore

# Everything the pipeline can make again is not kept in version control
results/
logs/
.snakemake/

The README says where the data came from. What do the checksums add that the README text cannot?

The README claims which files were used; the checksums let anyone verify the claim – that their copy is byte-for-byte the same as yours, and that nothing has changed since. Provenance you can check is worth much more than provenance you have to trust.

By hand, then a script

≈ 30 min

First the familiar way: type each step. Then put the steps in a bash script – the simplest kind of pipeline – and find out what a script cannot do.

Make sure the terminal is in your project folder (the prompt ends in cyp2c19-pipeline; if not: cd ~/cyp2c19-pipeline ). All paths from now on are relative to the project folder.

2.1 By hand

These are the steps of the reads-to-variants practical, with minimap2 for mapping. Run them one after the other and read what each one prints.

Activity 2.1 · Seven commands

Folders for the results: mkdir -p results/mapped results/qc results/variants
Map the read pairs to the reference and sort the alignments by position (about 5 s): minimap2 -ax sr -R '@RG\tID:NA12878\tSM:NA12878' data/raw/reference.fa data/raw/NA12878_R1.fastq data/raw/NA12878_R2.fastq | samtools sort -o results/mapped/NA12878.sorted.bam -
-ax sr: short reads, SAM output; -R adds a read group naming the sample (\t means a tab); the final - tells samtools to read from the pipe.
Index the BAM file: samtools index results/mapped/NA12878.sorted.bam
Count the mapped reads: samtools flagstat results/mapped/NA12878.sorted.bam > results/qc/NA12878.flagstat.txt then head -5 results/qc/NA12878.flagstat.txt
Call the variants: bcftools mpileup -f data/raw/reference.fa results/mapped/NA12878.sorted.bam -Ou | bcftools call -mv -Oz -o results/variants/NA12878.raw.vcf.gz
Keep the confident calls: bcftools filter -i 'QUAL>=20 && INFO/DP>=5' -Oz -o results/variants/NA12878.filtered.vcf.gz results/variants/NA12878.raw.vcf.gz
How many variants passed? bcftools view -H results/variants/NA12878.filtered.vcf.gz | wc -l

How many variants passed the filter?

15 of 60 raw calls.

A number is easy to compare, but it does not prove two files are the same. A checksum of the variant records does. (Checksum the records, not the file: the header of a VCF file records the date and time it was made, so two identical analyses never give identical files.)

Activity 2.2 · A fingerprint of the result

bcftools view -H results/variants/NA12878.filtered.vcf.gz | md5sum – -H prints the variant records without the header.

Copy the checksum here. You will compare every later run – Snakemake, Galaxy, a colleague’s computer – with it.

09dd51431f385ed3863928eee0064125 – if yours differs, check that you used exactly the commands above.

2.2 A script – the simplest pipeline

Typing commands is slow and error-prone, and the only record is your shell history. Put them in a file instead: a bash script. Lines starting with # are comments; a \ at the end of a line continues the command on the next line.

workflow/run_all.sh

#!/bin/bash
# Variant calling for NA12878 – every step, in order (first version)
mkdir -p results/mapped results/qc results/variants

# 1. map the reads and sort the alignments
minimap2 -ax sr -R '@RG\tID:NA12878\tSM:NA12878' \
    data/raw/reference.fa data/raw/NA12878_R1.fastq data/raw/NA12878_R2.fastq \
    | samtools sort -o results/mapped/NA12878.sorted.bam -

# 2. index the BAM file, and count the mapped reads
samtools index results/mapped/NA12878.sorted.bam
samtools flagstat results/mapped/NA12878.sorted.bam > results/qc/NA12878.flagstat.txt

# 3. call the variants and keep the good ones
bcftools mpileup -f data/raw/reference.fa results/mapped/NA12878.sorted.bam -Ou \
    | bcftools call -mv -Oz -o results/variants/NA12878.raw.vcf.gz
bcftools filter -i 'QUAL>=20 && INFO/DP>=5' -Oz -o results/variants/NA12878.filtered.vcf.gz \
    results/variants/NA12878.raw.vcf.gz

echo "All steps finished."

Activity 2.3 · Run the script

Press Create this file. Read it in the editor: it is exactly the commands you typed.
Run it: bash workflow/run_all.sh – watch the time in the status bar under the terminal.
Run it again: bash workflow/run_all.sh . Nothing has changed – what does the script do?

The second run redid every step. Why is that a problem for a real project – say 100 samples, when one new sample arrives, or when only the filter threshold changes?

A script has no idea what is already done: it cannot tell that the BAM files are up to date, so it repeats hours of mapping to change one threshold or add one sample. It lacks two things a workflow engine provides: dependencies (which file is made from which) and re-entrancy (continuing from where it stopped, redoing only what is needed). You could add checks by hand (“if the file exists, skip”) – but then an outdated file would never be updated.

2.3 When a step fails

A typo is all it takes. Break the script on purpose and see what it reports.

Activity 2.4 · A silent failure

In the editor, change NA12878_R1.fastq in the minimap2 command to NA12878_R3.fastq (a file that does not exist) and save with Ctrl+S.
Run it: bash workflow/run_all.sh – read all the output. Then ask the shell for the script’s exit status: echo $? (0 means success).
Look at the results: ls -l results/mapped results/variants – compare the sizes with before. How many reads are in the new BAM file? samtools view -c results/mapped/NA12878.sorted.bam

minimap2 could not open its input – yet the script made new results and said “All steps finished.” Why did the later steps still “work”, and why is this dangerous?

minimap2 writes the SAM header before it opens the reads. So samtools sort received a valid header with no reads and wrote an empty BAM file (309 bytes, 0 reads, instead of about 500 kB). Indexing, flagstat, calling and filtering all ran without complaint on it and made empty results, and the script exited with status 0 – the status of its last command, echo. Nothing reports the failure: the results are simply wrong. Here they are conspicuously empty, but a failure halfway through a file – or one of two read files missing – gives results that look normal, and a report would say “no variants” where there are some.

Bash can be told to stop at the first error. Three settings, known together as strict mode:

`set -e`	stop as soon as a command fails
`set -u`	treat a variable that was never set as an error
`set -o pipefail`	a pipe (`a \| b`) fails if any program in it fails, not only the last one

Activity 2.5 · Strict mode

Add the line set -euo pipefail as the second line of the script (after #!/bin/bash), and save.
Run the broken script again: bash workflow/run_all.sh then echo $? . Where does it stop now?
Fix the typo (back to NA12878_R1.fastq), save, and run the script once more: bash workflow/run_all.sh

Without set -o pipefail, what is the exit status of minimap2 … | samtools sort … when minimap2 fails but samtools sort succeeds?

the exit status of minimap2 (the first program) the exit status of samtools sort – success always 1

What a workflow engine adds

Strict mode stops a script at the first error – but it still cannot tell which files are up to date, it still redoes everything, and it leaves half-made files behind. A workflow engine such as Snakemake adds exactly what is missing:

dependencies – it knows which file is made from which;
re-entrancy – it redoes only what is needed;
clean failures – it deletes the output of a failed step so that nothing half-made is mistaken for a result;
a record – of what was run, when, with which code and settings.

Python for Snakemake

≈ 35 min

A Snakefile is Python with a few extra keywords. Learn the handful of Python ideas you need – values, lists, dictionaries, file-name patterns, indentation – in the notebook, then run the same code as a script.

Snakemake is written in Python, and so is a Snakefile: a Python file with a few extra keywords. You do not need to be a programmer to write one – but you need to read a handful of Python ideas. Try each of them in the notebook, where you can change a line and run it again.

How the notebook works

Each grey box is a cell of Python. Run a cell with ▶ or Shift+Enter; its output appears below it. Python remembers the variables made by the cells you have run – in the order you ran them. The buttons in these instructions put code into a new cell and run it; you can then edit it and run it again.

The notebook works in your project folder, so config/config.yaml means ~/cyp2c19-pipeline/config/config.yaml.

3.1 Values, variables and lists

Python

sample = "NA12878"          # text (a "string") needs quotes
min_qual = 20               # a number does not
samples = ["NA12878", "NA12891", "NA12892"]   # a list: items in [ ], separated by commas

print(sample, min_qual)
print(samples[0], len(samples))

`name = value`	makes a variable: a name for a value
`"NA12878"`	a string (text) – always in quotes; `20` is a number
`[a, b, c]`	a list; `samples[0]` is the first item (Python counts from 0)
`len(x)`, `print(x)`	calling a function: its name, then the inputs (arguments) in brackets
`# …`	a comment, ignored by Python

Activity 3.1 · Your first Python

Run the cell above (the button puts it in the notebook and runs it).
In the notebook, add the line print(samples[1]) at the end of the cell, and run it again.

What does samples[2] give?

NA12892 – the third item, because counting starts at 0.

3.2 The comma trap

Python

bad = ["NA12878" "NA12891"]     # oops: no comma between the two names
print(bad, len(bad))

What is len(bad), and why?

1. With no comma between them, Python joins two strings that stand next to each other into one: the list holds a single item, 'NA12878NA12891'. There is no error message – it is a silent bug. Remember it: in a Snakefile, a missing comma between two input files glues their names together.

3.3 Dictionaries and a config file

Settings belong in a config file, not in the code. YAML is a simple text format for them: key: value pairs, and lists with -. Create the pipeline’s config file now – Snakemake will read it in chapter 5.

config/config.yaml

# Settings for the pipeline. Change them here, not in the Snakefile.
samples:
  - NA12878

reference: data/raw/reference.fa

# variant filters
min_qual: 20
min_depth: 5

Python

import yaml                               # a package that reads YAML files

with open("config/config.yaml") as f:     # open the file; the indented block uses it
    config = yaml.safe_load(f)            # YAML text -> Python dictionaries and lists

print(config)
print(config["min_qual"], config["samples"])

`{"min_qual": 20, ...}`	a dictionary: values looked up by a key – `config["min_qual"]` gives 20
`import yaml`	loads a package (a library of ready-made functions)
`with open(...) as f:`	opens a file for the indented block below it, and closes it afterwards. Indentation (4 spaces) is how Python groups lines – it is not decoration.

Activity 3.2 · Read the settings

Create config/config.yaml (button above), then run the notebook cell that reads it.

What kind of Python value is config["samples"]?

A list (['NA12878']) – because the YAML file lists the samples with -.

3.4 Making file names – the heart of Snakemake

Python

for s in config["samples"]:                 # repeat the indented line for each sample
    print(f"data/raw/{s}_R1.fastq")           # an f-string fills in {s}

pattern = "results/variants/{sample}.filtered.vcf.gz"
print(pattern.format(sample="NA12878"))       # .format fills in {sample} later – as Snakemake does

An f-string (f"…{s}…") fills in a value now; .format() fills in a pattern later. Snakemake works like .format(): your rules contain patterns such as results/{sample}.bam, and Snakemake fills in the sample when it knows which file you want. Snakemake’s expand() makes one name per value:

Python

from snakemake.io import expand

print(expand("results/variants/{sample}.filtered.vcf.gz", sample=config["samples"]))
print(expand("results/qc/{sample}.flagstat.txt", sample=["A", "B", "C"]))

What does expand("results/qc/{sample}.flagstat.txt", sample=["A", "B", "C"]) return?

A list of three file names: ['results/qc/A.flagstat.txt', 'results/qc/B.flagstat.txt', 'results/qc/C.flagstat.txt'].

3.5 Functions and keyword arguments

Python

def vcf_name(sample, kind="filtered"):      # def: a reusable, named piece of code
    return f"results/variants/{sample}.{kind}.vcf.gz"

print(vcf_name("NA12878"))                   # kind uses its default value
print(vcf_name("NA12878", kind="raw"))       # a keyword argument: name=value

A rule in a Snakefile looks much like this: a name followed by a colon, an indented body, and keyword=value pairs such as ref="data/raw/reference.fa".

3.6 Notebook or script?

A notebook remembers what you ran, in the order you ran it – not the order of the cells on the page. That makes it easy to fool yourself.

Activity 3.3 · Hidden state

Put this in a cell and run it: threshold = 20
And this in the next cell: print("keeping variants with QUAL >=", threshold)
Now edit the first cell to threshold = 30 – but do not run it. Run only the second cell. What does the notebook claim now?
Press Restart (notebook toolbar) and then Run all. Which value is printed now?

After the edit, the page showed threshold = 30 but printed 20. Why is that a reproducibility problem, and how do scripts avoid it?

The notebook’s memory (state) no longer matched the code on the page, so the output could not be reproduced from the code you can see. A script always runs from the top, in a fresh Python, so the code is the record. Good practice: explore in a notebook, check with Restart and Run all, and move the code that makes results into scripts that the pipeline runs.

Now the same idea as a script: save it, then run it from the terminal.

workflow/scripts/hello.py

# List the files the pipeline will make for each sample in config/config.yaml
import yaml

with open("config/config.yaml") as f:
    config = yaml.safe_load(f)

print("Reference:", config["reference"])
for sample in config["samples"]:
    print(f"results/variants/{sample}.filtered.vcf.gz")
    print(f"results/qc/{sample}.flagstat.txt")

Activity 3.4 · From notebook to script

Create workflow/scripts/hello.py (button above).
In the Terminal, in your project folder: python workflow/scripts/hello.py

Claerbout on interactive programs

“Interactive programs should always be able to save their state so they can restart. Otherwise, dependence on an interactive program can be a form of slavery (nonreproducible research).” – Jon Claerbout, Reproducible computational research

Your first Snakefile

≈ 45 min

Describe each step as a rule – what it needs, what it makes and the command in between – and let Snakemake work out the order, run only what is needed, and keep a record.

A rule

A Snakefile is a set of rules. Each rule says how to make its output files from its input files – with a shell command, a Python script or a run block – and refers to its files as {input} and {output} in the command. Snakemake starts from the files you ask for (by default those of the first rule: its outputs – or, for a target rule such as all, its inputs), finds the rules that make them, then the rules that make their inputs, and so on – building a DAG (directed acyclic graph) of jobs.

4.1 The first rule

workflow/Snakefile

rule map_reads:
    input:
        "data/raw/reference.fa",
        "data/raw/NA12878_R1.fastq",
        "data/raw/NA12878_R2.fastq",
    output:
        "results/mapped/NA12878.sorted.bam",
    shell:
        "minimap2 -ax sr -R '@RG\\tID:NA12878\\tSM:NA12878' {input} | samtools sort -o {output} -"

The \\t inside the Python string becomes \t in the command, which minimap2 turns into a tab.

`rule map_reads:`	a rule and its name (just a label), then an indented body
`input:`	the files the rule needs – a list of strings, so commas between them
`output:`	the files it makes
`shell:`	the command. `{input}` becomes all the input files, in order, separated by spaces; `{output}` the output file

Activity 4.1 · Run a Snakefile

Create workflow/Snakefile with the button above. (Snakemake finds workflow/Snakefile by itself.)
“Burn” the results you made by hand – from now on the pipeline makes them: rm -r results
Ask Snakemake what it would do – a dry run: snakemake -n . Find the job, and the reason it gives for running it.
Run it, printing the command it runs (-p), on one CPU core: snakemake --cores 1 -p
Run it again: snakemake --cores 1

The second run said “Nothing to be done”. Why?

Snakemake only runs a rule once per session the output file exists and is newer than its inputs, and nothing it was made with has changed the dry run already made the file

Snakemake runs a job only if its output is missing, older than an input, or was made with different code, parameters, inputs or software than now (it keeps records in .snakemake/metadata) – or if you force it. This is re-entrancy: re-run a pipeline at any time and only what is needed happens.

4.2 Two rules of your own

Now write the rules for the next two steps yourself. Here is the skeleton of flagstat: add it below map_reads in the editor, and replace each ____ with a file name. Its input is the sorted BAM file; its output should be results/qc/NA12878.flagstat.txt.

rule flagstat:
    input:
        "____",
    output:
        "____",
    shell:
        "samtools flagstat {input} > {output}"

Then write index_bam the same way. samtools index writes its output, the index results/mapped/NA12878.sorted.bam.bai, next to the BAM file by itself – so the command does not mention {output}, but the rule still needs an output: line to tell Snakemake which file it makes.

rule index_bam:
    input:
        "____",
    output:
        "____",
    shell:
        "samtools index {input}"

Activity 4.2 · Write, check, run

Add rule flagstat and save. Ask for its output file – Snakemake works out the rest: snakemake --cores 1 -p results/qc/NA12878.flagstat.txt
Add rule index_bam, save, and ask for the index: snakemake --cores 1 -p results/mapped/NA12878.sorted.bam.bai

Stuck? – and compare its answer with yours.

4.3 A rule from the AI assistant

Variant calling is a longer command. Ask the AI assistant to write the rule – and treat its answer the way you would treat code from a stranger.

Activity 4.3 · Check the AI’s code

Read the answer: does the explanation match the code?
Press Add rule call_variants to Snakefile under the code (it goes at the end of your Snakefile), then check it with a dry run: snakemake -n results/variants/NA12878.raw.vcf.gz
Read the error message carefully. What file is Snakemake looking for? Fix the rule yourself in the editor – or press ✦ Ask the AI assistant about this error under the error in the terminal.
When the dry run works, run it: snakemake --cores 1 results/variants/NA12878.raw.vcf.gz

What was wrong with the assistant’s rule, how did you find out, and why did Python not complain?

A comma was missing after "data/raw/reference.fa". Python joins two adjacent strings into one, so Snakemake looked for a single input file called data/raw/reference.faresults/mapped/NA12878.sorted.bam – the dry run’s MissingInputException showed the glued name. Python raises no error because joining adjacent strings is legal (it is how long commands are split over lines). The assistant also explained the code confidently – an explanation is not evidence that code works; running it is.

4.4 Filter – and one rule to rule them all

Write the last step yourself. The command is the one you ran by hand in chapter 2 – with the file names replaced by {input} and {output}:

rule filter_variants:
    input:
        "results/variants/NA12878.raw.vcf.gz",
    output:
        "results/variants/NA12878.filtered.vcf.gz",
    shell:
        "____"

Finally, a target rule. Snakemake builds the first rule of the Snakefile when you give no target, so pipelines start with a rule – by convention called all – whose inputs are the files you want at the end. It has no output and no command. Put it at the very top of the Snakefile:

rule all:
    input:
        "results/variants/NA12878.filtered.vcf.gz",
        "results/qc/NA12878.flagstat.txt",

Activity 4.4 · The whole pipeline

Add rule filter_variants (with your command) below the other rules, and save.
Add rule all at the top, and save.
Burn and build: rm -r results then snakemake --cores 1
Is the result the same as by hand? bcftools view -H results/variants/NA12878.filtered.vcf.gz | md5sum

Look at the job list of the last run. Which of your rules did not run, and why not?

index_bam did not run. No rule needs its output: rule all does not ask for the .bai file, and call_variants does not list it as an input (bcftools reads the whole BAM file without the index). Snakemake runs only the jobs needed for the files you ask for: a job runs if its output is a target or is needed by another job – and only if that output is missing or out of date.

If something does not work, compare your Snakefile with this reference version (it is what you should have now):

workflow/Snakefile

# Variant calling for NA12878 – first version (one sample, names written out)

rule all:
    input:
        "results/variants/NA12878.filtered.vcf.gz",
        "results/qc/NA12878.flagstat.txt",


rule map_reads:
    input:
        "data/raw/reference.fa",
        "data/raw/NA12878_R1.fastq",
        "data/raw/NA12878_R2.fastq",
    output:
        "results/mapped/NA12878.sorted.bam",
    shell:
        "minimap2 -ax sr -R '@RG\\tID:NA12878\\tSM:NA12878' {input} | samtools sort -o {output} -"


rule index_bam:
    input:
        "results/mapped/NA12878.sorted.bam",
    output:
        "results/mapped/NA12878.sorted.bam.bai",
    shell:
        "samtools index {input}"


rule flagstat:
    input:
        "results/mapped/NA12878.sorted.bam",
    output:
        "results/qc/NA12878.flagstat.txt",
    shell:
        "samtools flagstat {input} > {output}"


rule call_variants:
    input:
        "data/raw/reference.fa",
        "results/mapped/NA12878.sorted.bam",
    output:
        "results/variants/NA12878.raw.vcf.gz",
    shell:
        "bcftools mpileup -f {input} -Ou | bcftools call -mv -Oz -o {output}"


rule filter_variants:
    input:
        "results/variants/NA12878.raw.vcf.gz",
    output:
        "results/variants/NA12878.filtered.vcf.gz",
    shell:
        "bcftools filter -i 'QUAL>=20 && INFO/DP>=5' -Oz -o {output} {input}"

4.5 See the DAG

Activity 4.5 · Draw the graph

Snakemake describes its DAG in the dot language; the Graphviz program dot draws it: snakemake --dag | dot -Tsvg > dag.svg
Open the picture: open dag.svg

On your own computer, snakemake --dag | dot -Tpdf > dag.pdf makes a PDF; this browser’s Graphviz makes SVG and PNG pictures, but not PDF. snakemake --rulegraph draws one box per rule instead of per job – simpler when there are many samples.

4.6 Re-entrancy and provenance

Now test the promises of a workflow engine. Run each dry run and read the reason Snakemake gives for each job.

Activity 4.6 · What does Snakemake re-run?

Delete one result: rm results/qc/NA12878.flagstat.txt , then snakemake -n . Which jobs would run?
Change QUAL>=20 to QUAL>=30 in filter_variants, save, and dry-run again: snakemake -n
Change it back to 20, save, and dry-run once more: snakemake -n – then build what is missing: snakemake --cores 1
Now make a step fail: in rule flagstat change samtools flagstat to samtools flagstatt (a typo), save, and force that rule to run: snakemake --cores 1 -R flagstat . Read the whole message.
Is there a half-made flagstat file left? ls -l results/qc Then fix the typo, save, and run snakemake --cores 1 .

What reasons did Snakemake give (a) after you deleted the flagstat file and (b) after you changed the threshold? Why did changing it back need no re-run?

(a) Missing output files: only flagstat (and all) had to run – not the mapping or calling. (b) Code has changed since last execution for filter_variants: Snakemake stores the code each output was made with (in .snakemake/metadata) and notices that the rule is no longer the same – provenance information, beyond comparing the times of the files. After changing it back, the code matched the record of the last run again, so the existing file was still valid.

Compare what Snakemake did when the flagstat step failed with what your bash script did in chapter 2.

Snakemake stopped at the failing job with a clear Error in rule flagstat block (and the program’s own error message), exited with a non-zero status, and removed the output of the failed job – the shell had already created an empty flagstat file with > – so no half-made file could be mistaken for a result. The bash script without strict mode carried on, reported success and left wrong (empty) results behind.

One Snakefile for any sample

≈ 50 min

Make the pipeline general – wildcards, a config file, named inputs, parameters, logs and a Python script – and watch Snakemake keep track of what each result was made with.

Your pipeline works – for one sample, with its name and the thresholds written into every rule. Now make it general: wildcards for samples, the config file for settings, named inputs, parameters, logs, and a Python script for a report.

5.1 Wildcards

Replace the sample name with a wildcard, {sample}. When Snakemake needs results/mapped/NA12878.sorted.bam, it matches it against the output pattern results/mapped/{sample}.sorted.bam, sets sample = NA12878, and fills the same value into the inputs. Inside the command, the value is {wildcards.sample}. Here is map_reads with wildcards – replace yours with it:

rule map_reads:
    input:
        "data/raw/reference.fa",
        "data/raw/{sample}_R1.fastq",
        "data/raw/{sample}_R2.fastq",
    output:
        "results/mapped/{sample}.sorted.bam",
    shell:
        "minimap2 -ax sr -R '@RG\\tID:{wildcards.sample}\\tSM:{wildcards.sample}' {input} "
        "| samtools sort -o {output} -"

Activity 5.1 · Any sample

Replace your map_reads rule with the one above, and save. (Two strings next to each other in shell: are joined into one command – on purpose this time.)
Now do the same in index_bam, flagstat, call_variants and filter_variants yourself: replace NA12878 by {sample} in their input: and output:. Leave rule all as it is – it names the files you want. Save.
Dry run: snakemake -n . Which job has a different reason from the others? Then run it: snakemake --cores 1

Which job re-ran because of its own change, and why did the jobs after it re-run too?

map_reads: Code has changed since last execution – its command now uses {wildcards.sample}. The jobs downstream re-ran with the reason Input files updated by another job: their input, the BAM file, was about to be remade. Rules whose command did not change (e.g. samtools flagstat {input} > {output}) were only re-run because their input changed – replacing NA12878 by {sample} in a file pattern changes no command.

5.2 Settings from the config file

You created config/config.yaml in chapter 3. configfile: reads it into a Python dictionary called config – exactly what you did with yaml.safe_load in the notebook. Replace your rule all (at the top) with this block, which uses the sample list from the config file and expand():

configfile: "config/config.yaml"

SAMPLES = config["samples"]
REF = config["reference"]


rule all:
    input:
        expand("results/variants/{sample}.filtered.vcf.gz", sample=SAMPLES),
        expand("results/qc/{sample}.flagstat.txt", sample=SAMPLES),

Activity 5.2 · One list of samples

Replace rule all with the block above, save, and dry-run: snakemake -n . Nothing to do – the same files, now listed by expand().
Pretend a second sample arrives: in config/config.yaml add a line - NA12891 under - NA12878, save, and dry-run: snakemake -n . Read the error: which files would the new sample need? Then remove NA12891 again (we have no reads for it) and save.

What did Snakemake tell you about the new sample, and what would you need to do to analyse it?

A MissingInputException for rule map_reads: to make results/mapped/NA12891.sorted.bam it needs data/raw/NA12891_R1.fastq and data/raw/NA12891_R2.fastq, which do not exist. With the two read files in data/raw/, one command would analyse the new sample – and leave NA12878’s results untouched, because they are up to date.

5.3 Named inputs and logs

With several inputs, names are clearer than positions: ref=… in input:, {input.ref} in the command. A log: file keeps each program’s messages. Replace your map_reads with this version – it also takes the reference from the config file:

rule map_reads:
    input:
        ref=REF,
        r1="data/raw/{sample}_R1.fastq",
        r2="data/raw/{sample}_R2.fastq",
    output:
        "results/mapped/{sample}.sorted.bam",
    log:
        "logs/map_reads/{sample}.log",
    threads: 1
    shell:
        "minimap2 -ax sr -t {threads} -R '@RG\\tID:{wildcards.sample}\\tSM:{wildcards.sample}' "
        "{input.ref} {input.r1} {input.r2} 2> {log} "
        "| samtools sort -o {output} - 2>> {log}"

`ref=REF`	a named input, used as `{input.ref}` in the command. `REF` is the Python variable set from the config file at the top of the Snakefile.
`threads: 1`	how many CPU cores the job may use; `{threads}` puts the number into the command (`-t` is minimap2’s option for threads).
`log:`	the job’s log file, `{log}` in the command.
`2> {log}`	sends a program’s messages – its standard error, stream 2 – to the log file. `2>>` adds them to the end of the file instead of replacing it, so both programs’ messages are kept.

Activity 5.3 · Names, logs – and the AI again

Replace map_reads with the version above, save, and run snakemake --cores 1 . Then read the log: cat logs/map_reads/NA12878.log
Ask the assistant to do the same for variant calling: Press Replace rule call_variants in Snakefile under its code.
Run it: snakemake --cores 1 . If it fails, read the error and the log file it points to – then fix the rule, or ask the assistant about the error.

What was wrong with the assistant’s call_variants, and why did index_bam run this time?

The command still used {input}, which – with named inputs too – means all the inputs: the reference, the BAM and the index. bcftools tried to read the .bai index as a BAM file (“Exec format error”) and failed. The fix is -f {input.ref} {input.bam}. index_bam ran because the new rule lists the index as an input (bai=…): now a job needs it.

5.4 Parameters

The filter thresholds are settings, so they belong in the config file. params: passes values to the command as {params.name} – and Snakemake records them.

Activity 5.4 · Thresholds from the config file

Press Replace rule filter_variants in Snakefile under the rule, and run snakemake --cores 1 .
If it fails, find out why (the program’s help text is printed with the error) and fix the command – keeping {params.min_qual} and {params.min_depth}. Run again until it works.
Same result as ever? bcftools view -H results/variants/NA12878.filtered.vcf.gz | md5sum
Change min_qual: 20 to min_qual: 30 in config/config.yaml, save, and dry-run: snakemake -n . Then run it and count: snakemake --cores 1 bcftools view -H results/variants/NA12878.filtered.vcf.gz | wc -l Finally set it back to 20, save, and run snakemake --cores 1 once more.

What was wrong with the assistant’s filter_variants? What does this tell you about AI-written code?

It used options that do not exist: bcftools filter has no --min-qual or --min-depth (bcftools printed “unrecognized option” and its real options). The fix is an expression, -i 'QUAL>={params.min_qual} && INFO/DP>={params.min_depth}'. AI assistants can invent plausible options, functions or even tools; check against --help or the manual, and test.

With min_qual: 30, how many variants passed – and what reason did the dry run give?

13 (two calls have QUAL between 20 and 30). The reason was Params have changed since last execution, showing the old and new value – Snakemake compared the recorded parameters of the existing file with the config file. Only filter_variants and what comes after it re-ran.

5.5 A Python script in the pipeline

The last step makes a table and a plot with Python. With script:, Snakemake runs a Python script and gives it an object called snakemake, so the script reads snakemake.input and writes snakemake.output.table instead of file names written into it. Read the script with the table below it: some ideas are from chapter 3, the others are new – you need to read them, not write them.

workflow/scripts/summarise_variants.py

# Summarise the filtered variants: a table of counts and a plot of call quality.
# Snakemake runs this script and gives it an object called `snakemake`
# with the rule's input, output, params, wildcards, ...
import gzip

import pandas as pd
import matplotlib

matplotlib.use("Agg")  # draw plots into files, not on a screen
import matplotlib.pyplot as plt

rows = []
for path in snakemake.input:
    sample = path.split("/")[-1].split(".")[0]  # results/variants/NA12878.filtered.vcf.gz -> NA12878
    with gzip.open(path, "rt") as vcf:
        for line in vcf:
            if line.startswith("#"):  # skip the header lines
                continue
            fields = line.rstrip("\n").split("\t")
            chrom, pos, ref, alt, qual = fields[0], int(fields[1]), fields[3], fields[4], float(fields[5])
            genotype = fields[9].split(":")[0]
            kind = "SNP" if len(ref) == 1 and len(alt) == 1 else "indel"
            rows.append({"sample": sample, "chrom": chrom, "pos": pos, "ref": ref, "alt": alt,
                         "qual": qual, "genotype": genotype, "type": kind})

variants = pd.DataFrame(rows)
summary = variants.groupby(["sample", "chrom", "type", "genotype"]).size().reset_index(name="variants")
summary.to_csv(snakemake.output.table, sep="\t", index=False)

fig, ax = plt.subplots(figsize=(6, 4))
for sample, group in variants.groupby("sample"):
    ax.hist(group["qual"], bins=20, alpha=0.7, label=sample)
ax.set_xlabel("Variant quality (QUAL)")
ax.set_ylabel("Number of variants")
ax.set_title("Quality of the filtered variant calls")
ax.legend()
fig.tight_layout()
fig.savefig(snakemake.output.plot, dpi=100)

`for path in snakemake.input:`	repeat the indented lines for each input file of the rule
`path.split("/")[-1].split(".")[0]`	`.split("/")` cuts text into a list at every `/`; `[-1]` is the last item (the file name); the part before the first `.` is the sample name
`with gzip.open(path, "rt") as vcf:`	open a compressed file as text (`"rt"`) for the indented block
`if …: continue`	if the line is a header line, skip to the next line
`line.rstrip("\n").split("\t")`	remove the line break, then cut the line into its tab-separated columns
`chrom, pos, … = fields[0], int(fields[1]), …`	several variables at once; `int()` and `float()` turn text into numbers
`"SNP" if … else "indel"`	one of two values, chosen by a condition
`rows.append({…})`	add one row – a dictionary – to the list `rows`
`pd.DataFrame(rows)`, `.groupby([…]).size()`	pandas: a table made from the rows; count the rows in each group (sample, chromosome, type, genotype)
`plt.subplots()`, `ax.hist(…)`, `fig.savefig(…)`	matplotlib: make a figure, draw a histogram of the QUAL values, save it as a file

rule summarise:
    input:
        expand("results/variants/{sample}.filtered.vcf.gz", sample=SAMPLES),
    output:
        table="results/report/variant_summary.tsv",
        plot="results/report/variant_quality.png",
    script:
        "scripts/summarise_variants.py"

Activity 5.5 · A report

Create the script (button above).
Add rule summarise (above) at the end of the Snakefile, and add "results/report/variant_summary.tsv", to the inputs of rule all. Save.
Run it: snakemake --cores 1 – then cat results/report/variant_summary.tsv and open results/report/variant_quality.png

According to the summary table, how many homozygous (1/1) variants are there on human_CYP2C9?

3 – and 8 heterozygous (0/1) ones; CYP2C19 has 4, all heterozygous.

Your Snakefile should now be equivalent to this reference version – compare, or replace yours if it will not work:

workflow/Snakefile

# =====================================================================
#  Variant calling for CYP2C19 / CYP2C9 – a reproducible pipeline
#  Run from the project folder:   snakemake --cores 1
# =====================================================================

configfile: "config/config.yaml"

SAMPLES = config["samples"]
REF = config["reference"]


# the target rule: the files we want at the end
rule all:
    input:
        expand("results/variants/{sample}.filtered.vcf.gz", sample=SAMPLES),
        expand("results/qc/{sample}.flagstat.txt", sample=SAMPLES),
        "results/report/variant_summary.tsv",


# 1. align the reads to the reference and sort the alignments
rule map_reads:
    input:
        ref=REF,
        r1="data/raw/{sample}_R1.fastq",
        r2="data/raw/{sample}_R2.fastq",
    output:
        "results/mapped/{sample}.sorted.bam",
    log:
        "logs/map_reads/{sample}.log",
    threads: 1
    shell:
        "minimap2 -ax sr -t {threads} -R '@RG\\tID:{wildcards.sample}\\tSM:{wildcards.sample}' "
        "{input.ref} {input.r1} {input.r2} 2> {log} "
        "| samtools sort -o {output} - 2>> {log}"


# 2. index the sorted BAM file
rule index_bam:
    input:
        "results/mapped/{sample}.sorted.bam",
    output:
        "results/mapped/{sample}.sorted.bam.bai",
    shell:
        "samtools index {input}"


# 3. mapping statistics
rule flagstat:
    input:
        "results/mapped/{sample}.sorted.bam",
    output:
        "results/qc/{sample}.flagstat.txt",
    shell:
        "samtools flagstat {input} > {output}"


# 4. call variants
rule call_variants:
    input:
        ref=REF,
        bam="results/mapped/{sample}.sorted.bam",
        bai="results/mapped/{sample}.sorted.bam.bai",
    output:
        "results/variants/{sample}.raw.vcf.gz",
    log:
        "logs/call_variants/{sample}.log",
    shell:
        "bcftools mpileup -f {input.ref} {input.bam} -Ou 2> {log} "
        "| bcftools call -mv -Oz -o {output} 2>> {log}"


# 5. keep good-quality calls
rule filter_variants:
    input:
        "results/variants/{sample}.raw.vcf.gz",
    output:
        "results/variants/{sample}.filtered.vcf.gz",
    params:
        min_qual=config["min_qual"],
        min_depth=config["min_depth"],
    shell:
        "bcftools filter -i 'QUAL>={params.min_qual} && INFO/DP>={params.min_depth}' "
        "-Oz -o {output} {input}"


# 6. a table and a plot, made by a Python script
rule summarise:
    input:
        expand("results/variants/{sample}.filtered.vcf.gz", sample=SAMPLES),
    output:
        table="results/report/variant_summary.tsv",
        plot="results/report/variant_quality.png",
    script:
        "scripts/summarise_variants.py"

Extension: a one-liner from the assistant Extension

Add it, run snakemake --cores 1 results/qc/NA12878.per_chrom.txt, and work out from the error message how to fix it.

Software environments

≈ 40 min

Same data, same Snakefile, same settings – different results. Find out why, pin the software, and record the environment.

Same data, same Snakefile, same settings – different results? It happens when the software differs. In Claerbout’s principle, the scholarship includes “the complete software development environment”.

6.1 Which software are you using?

This practical uses the package manager conda – here Miniforge, with the conda-forge and bioconda channels. An environment is a folder of programs with particular versions; the active one decides which program a command runs.

Activity 6.1 · Look at your environment

conda env list – the * marks the active environment (also shown in the prompt).
conda list – every package in it, with version, build and channel.
which bcftools then bcftools --version | head -2
Python has its own package list: python --version pip list

Which version of the htslib library does bcftools --version say it uses – and which htslib version does conda list show?

bcftools says htslib 1.10, conda lists htslib 1.17 (the one samtools uses). In this browser each program was compiled with its own built-in copies of its libraries (htslib, zlib), so bcftools always uses htslib 1.10 – and one environment can hold every program. On a real computer the programs in an environment share its single copy of each library – whichever version conda chose. Keep this in mind for the next activities.

6.2 The same pipeline on another computer

A second run of this pipeline – the same rules, config and data, with the variant-calling steps in a conda environment – was made on a Linux computer with snakemake --use-conda. Its result is in /data/course/colleague/.

Activity 6.2 · Spot the difference

cat /data/course/colleague/README.md
Compare checksums: bcftools view -H /data/course/colleague/NA12878.filtered.vcf.gz | md5sum (yours was 09dd51431f385ed3863928eee0064125).
Save both sets of records as text and compare them line by line – in a scratch folder outside the project, because nothing in results/ is made by hand: mkdir -p ~/compare && bcftools view -H results/variants/NA12878.filtered.vcf.gz > ~/compare/mine.txt bcftools view -H /data/course/colleague/NA12878.filtered.vcf.gz > ~/compare/colleague.txt
diff ~/compare/mine.txt ~/compare/colleague.txt | head -20 – < lines are yours, > lines the colleague’s.
Just position, genotype and quality: cut -f 1,2,6,10 ~/compare/mine.txt > ~/compare/mine.short cut -f 1,2,6,10 ~/compare/colleague.txt > ~/compare/colleague.short diff ~/compare/mine.short ~/compare/colleague.short
The VCF header records the software that made it: bcftools view -h /data/course/colleague/NA12878.filtered.vcf.gz | grep Version bcftools view -h results/variants/NA12878.filtered.vcf.gz | grep Version
And the colleague’s environment: grep -E 'bcftools|htslib' /data/course/colleague/calling-environment.yml

What is the same in the two runs, and what differs? What caused the difference?

The same: all 15 variant sites, their alleles and their genotypes. Different: the QUAL of three calls (e.g. human_CYP2C9 23059: 111 in yours, 132 in theirs; 1402: 50 vs 56; 20827: 167 vs 165), the genotype likelihoods (PL) of several calls, and statistics in INFO such as DP4. Cause: the software. The colleague’s environment file asked for bcftools=1.10, which conda read as “any 1.10.x” and resolved to bcftools 1.10.2 – a build from 2020 – installed with htslib 1.21; yours is bcftools 1.10 with htslib 1.10. The ##bcftoolsVersion header lines record the htslib each run actually used. Same code, same data, different environment – different numbers.

Activity 6.3 · Does it matter?

Imagine a stricter filter, QUAL>=120. Count what each run keeps – yours, by changing the setting on the command line: snakemake --cores 1 --config min_qual=120 bcftools view -H results/variants/NA12878.filtered.vcf.gz | wc -l
…and the colleague’s: bcftools filter -i 'QUAL>=120' /data/course/colleague/NA12878.filtered.vcf.gz | bcftools view -H | wc -l
Go back to the settings in the config file: snakemake --cores 1

How many variants does each run keep with QUAL>=120?

6 in yours, 7 in the colleague’s: the call at human_CYP2C9 23059 has QUAL 111 in one run and 132 in the other. A version difference that looks harmless can change which variants are reported.

6.3 One environment for everything?

Activity 6.4 · Ask the assistant

Press Save as, then build it: conda env create -f workflow/envs/pipeline.yaml Read the whole message.
Press ✦ Ask the AI assistant about this error under the message.
Ask its follow-up question: without Python, would that environment give exactly the same bcftools and htslib as your pipeline? Compare with the colleague’s environment.
Try to pin bcftools and its library exactly, in one environment with samtools: conda create -n exact samtools=1.17 bcftools==1.10 htslib==1.10

The assistant called its single environment “fully reproducible”. What was wrong with it – and why does pinning bcftools and htslib exactly in one environment not help either?

It could not even be built. The bioconda builds of minimap2 2.22, samtools 1.17 and bcftools 1.10 were made against zlib 1.2 and need zlib older than 1.3, while Python 3.13 needs zlib 1.3.1 or newer – and an environment holds only one version of each library, so conda refuses (LibMambaUnsatisfiableError).

Without Python its pins are loose: bcftools=1.10 installs the newest 1.10.x (1.10.2), htslib is not mentioned so the solver picks the newest one that fits (1.21 today), and pandas and matplotlib have no version at all – exactly the colleague’s environment, with its different QUAL values.

Exact pins conflict too: htslib==1.10 for bcftools rules out samtools 1.17, which needs htslib 1.17 or newer. One environment for everything is fragile; pipelines give each step its own small, exactly pinned environment.

6.4 One environment per rule

Create three environment files. The tools and their htslib are pinned exactly (==) – the lesson of the colleague’s run. For the Python step, minor versions are enough here; conda env export records the exact ones (6.5).

workflow/envs/mapping.yaml

channels:
  - conda-forge
  - bioconda
dependencies:
  - minimap2==2.22
  - samtools==1.17
  - htslib==1.17

workflow/envs/calling.yaml

channels:
  - conda-forge
  - bioconda
dependencies:
  - bcftools==1.10
  - htslib==1.10

workflow/envs/python.yaml

channels:
  - conda-forge
dependencies:
  - python=3.13
  - pandas=2.3
  - matplotlib=3.9

In this browser, Python and its packages come from Pyodide (Python 3.13, pandas 2.3, matplotlib 3.8): Snakemake says so when it builds this environment. The file is written for a real computer, where matplotlib 3.9 is the first version available for Python 3.13.

Then give every rule a conda: line with its environment (the path is relative to the Snakefile): envs/mapping.yaml for map_reads, index_bam and flagstat; envs/calling.yaml for call_variants and filter_variants; envs/python.yaml for summarise. Put the conda: line above shell: or script: – Snakemake allows nothing after them. For example:

rule call_variants:
    input:
        ...
    output:
        ...
    log:
        ...
    conda:
        "envs/calling.yaml"
    shell:
        ...

Activity 6.5 · Environments in the Snakefile

The assistant’s environment file cannot be built – delete it, so that nobody mistakes it for part of the pipeline: rm workflow/envs/pipeline.yaml
Create the three environment files (buttons above).
Add a conda: line to each of the six rules, and save.
Run with the environments: snakemake --cores 1 --use-conda . Snakemake builds each environment once (in .snakemake/conda/). What reason does it give for re-running the jobs?
Still the same result? bcftools view -H results/variants/NA12878.filtered.vcf.gz | md5sum

What is the difference between bcftools=1.10 and bcftools==1.10, and why pin htslib as well?

bcftools=1.10 means “any version starting with 1.10” – 1.10, 1.10.1 or 1.10.2 – and conda picks the newest, together with whatever htslib that build needs; bcftools==1.10 means exactly 1.10. The numbers come from the library as much as from bcftools itself, and the dependency ranges of old packages are sometimes widened later – which is how a 2020 build of bcftools came to be installed with htslib 1.21 – so pinning htslib==1.10 too fixes the part that actually changed the colleague’s results. The most exact records are a lock file of the package files (conda list --explicit --md5; Snakemake uses one saved as envs/calling.linux-64.pin.txt) or a container image; conda env export lists every package with version and build.

Always the same way

Once your rules have conda: environments, run the pipeline with --use-conda every time. Without it, Snakemake notices that the software would not be the one recorded for your results (Software environment definition has changed) and wants to re-run everything – the provenance record at work.

6.5 Record the environment

Activity 6.6 · Records

Make the calling environment as a named environment: conda env create -n calling -f workflow/envs/calling.yaml
The packages you asked for: conda env export -n calling --from-history ; every package, version and build: conda env export -n calling (on a real computer this lists about thirty packages: every library bcftools needs).
Python packages in the format of a requirements.txt file: pip freeze
Snakemake’s own record of every result – which rule made it, when, from which inputs, with which command: snakemake --summary --use-conda and snakemake -D --use-conda | head -4

A reader wants the exact software you used, years from now. Which record helps most?

the environment files with loose pins (bcftools=1.10) conda env export (every package with version and build), a lock file, or a container image the version numbers in the README

Your finished Snakefile – for comparison, or to replace yours if it will not work:

workflow/Snakefile

# =====================================================================
#  Variant calling for CYP2C19 / CYP2C9 – a reproducible pipeline
#  Run from the project folder:   snakemake --cores 1 --use-conda
# =====================================================================

configfile: "config/config.yaml"

SAMPLES = config["samples"]
REF = config["reference"]


# the target rule: the files we want at the end
rule all:
    input:
        expand("results/variants/{sample}.filtered.vcf.gz", sample=SAMPLES),
        expand("results/qc/{sample}.flagstat.txt", sample=SAMPLES),
        "results/report/variant_summary.tsv",


# 1. align the reads to the reference and sort the alignments
rule map_reads:
    input:
        ref=REF,
        r1="data/raw/{sample}_R1.fastq",
        r2="data/raw/{sample}_R2.fastq",
    output:
        "results/mapped/{sample}.sorted.bam",
    log:
        "logs/map_reads/{sample}.log",
    threads: 1
    conda:
        "envs/mapping.yaml"
    shell:
        "minimap2 -ax sr -t {threads} -R '@RG\\tID:{wildcards.sample}\\tSM:{wildcards.sample}' "
        "{input.ref} {input.r1} {input.r2} 2> {log} "
        "| samtools sort -o {output} - 2>> {log}"


# 2. index the sorted BAM file
rule index_bam:
    input:
        "results/mapped/{sample}.sorted.bam",
    output:
        "results/mapped/{sample}.sorted.bam.bai",
    conda:
        "envs/mapping.yaml"
    shell:
        "samtools index {input}"


# 3. mapping statistics
rule flagstat:
    input:
        "results/mapped/{sample}.sorted.bam",
    output:
        "results/qc/{sample}.flagstat.txt",
    conda:
        "envs/mapping.yaml"
    shell:
        "samtools flagstat {input} > {output}"


# 4. call variants
rule call_variants:
    input:
        ref=REF,
        bam="results/mapped/{sample}.sorted.bam",
        bai="results/mapped/{sample}.sorted.bam.bai",
    output:
        "results/variants/{sample}.raw.vcf.gz",
    log:
        "logs/call_variants/{sample}.log",
    conda:
        "envs/calling.yaml"
    shell:
        "bcftools mpileup -f {input.ref} {input.bam} -Ou 2> {log} "
        "| bcftools call -mv -Oz -o {output} 2>> {log}"


# 5. keep good-quality calls
rule filter_variants:
    input:
        "results/variants/{sample}.raw.vcf.gz",
    output:
        "results/variants/{sample}.filtered.vcf.gz",
    params:
        min_qual=config["min_qual"],
        min_depth=config["min_depth"],
    conda:
        "envs/calling.yaml"
    shell:
        "bcftools filter -i 'QUAL>={params.min_qual} && INFO/DP>={params.min_depth}' "
        "-Oz -o {output} {input}"


# 6. a table and a plot, made by a Python script
rule summarise:
    input:
        expand("results/variants/{sample}.filtered.vcf.gz", sample=SAMPLES),
    output:
        table="results/report/variant_summary.tsv",
        plot="results/report/variant_quality.png",
    conda:
        "envs/python.yaml"
    script:
        "scripts/summarise_variants.py"

The same pipeline in Galaxy

≈ 40 min

Galaxy runs pipelines without programming: tools in forms, a history of every dataset, and workflows you build by drawing connections. Build yours – and compare the result with Snakemake’s.

A graphical workflow engine

Galaxy (usegalaxy.org, usegalaxy.eu and many institutional servers) runs bioinformatics tools through web forms, keeps every dataset in a history with a full record of how it was made, and lets you build workflows by connecting tools on a canvas. James Taylor and Anton Nekrutenko began Galaxy at Penn State; it was first described in 2005. In James Taylor’s words: “Good ideas don’t have owners – they belong to everyone.”

The Galaxy tab is a practice server that works like Galaxy and runs the same WebAssembly programs as your terminal and your Snakefile.

7.1 Data into a history

Activity 7.1 · Upload

Open the Galaxy tab. Tools are on the left, your history on the right.
Press Upload Data (top left), tick SRR098401_1.fastq, SRR098401_2.fastq and hg19_CYP2C_slices.fa (or your copies in data/raw/), and press Start. Galaxy detects the format of each file (fastqsanger, fasta).
Click the name of a dataset in the history to expand it; the eye shows its contents.

This practice server keeps your history only while the page is open (your workflows are saved). A real Galaxy server keeps every history until you delete it.

7.2 One tool, by hand

Activity 7.2 · Map with minimap2

Click Map with minimap2 in the tool panel.
Check the inputs – the reference, the forward reads (_1) and the reverse reads (_2) – and the preset (sr), then press Run tool. Watch the new dataset: grey = queued, yellow = running, green = done.
Expand the new BAM dataset and click its ⓘ (dataset details). Scroll through: tool ID and version, dependencies, command line, parameters.

Compare Galaxy’s command line with your map_reads rule. What is the same? What did Galaxy add, and what does it record for every dataset?

The same programs and options: minimap2 -x sr -a with a read group, piped into samtools sort. Galaxy first links the history datasets into a job working folder under simple names (ln -s … input_f.fastq), writes to output.bam and moves the result into its dataset store; it also indexes BAM datasets itself. For every dataset it records the tool ID and version (…/iuc/minimap2/minimap2/2.22+galaxy0), the dependencies with versions (minimap2 2.22, samtools 1.17), every parameter, the input datasets and the exact command line – provenance, automatically, without the user writing anything down.

7.3 Build the workflow

Now draw the whole pipeline. Each box is a step; each connection carries a dataset from an output (dot on the right of a box) to an input (dot on the left). Only matching data types can be connected.

Activity 7.3 · The workflow editor

Click Workflows (top bar of the Galaxy tab), then Create a new workflow.
Press + Input dataset three times. Click each box to set its label (forward reads, reverse reads, reference) and format (fastqsanger for the reads, fasta for the reference).
With the editor open, click these tools in the tool panel to add them: Map with minimap2, Samtools flagstat, bcftools mpileup, bcftools call, bcftools filter.
Try a wrong connection: drag from the forward reads output to the reference input of bcftools mpileup. What happens?
Now connect everything correctly: the three inputs into minimap2; minimap2’s BAM into flagstat and into mpileup (and the reference into mpileup); mpileup into call; call into filter. Click an empty part of the canvas: the Checks panel on the right says when every input is connected. (Tidy layout arranges the boxes.)
Click the bcftools filter box and check its expression: QUAL>=20 && INFO/DP>=5 – the same as min_qual and min_depth in your config file.

Why did Galaxy refuse the first connection? How does that help you?

Every dataset has a datatype and every tool input accepts only some: mpileup’s reference input needs fasta, and the reads are fastqsanger. The editor checks while you build, so mistakes surface immediately rather than as a failed (or silently wrong) job later. On the command line and in a Snakefile nothing checks file types for you – you find out from an error message, or from strange results.

7.4 Run it – and compare with Snakemake

Activity 7.4 · One click, five jobs

Press Run (in the editor). Check that each input has the right dataset, then press Run workflow.
Watch the jobs finish in the history and on the invocation page (about 15 seconds).
On the history dataset made by bcftools filter, click copy to my files. It is copied into ~/galaxy/.
In the terminal: ls ~/galaxy then compare the checksum of its records with your Snakemake result: bcftools view -H ~/galaxy/*filter*.vcf.gz | md5sum (if you ran bcftools filter more than once, use the full name of one file – press Tab to complete it)

Is Galaxy’s result the same as Snakemake’s? Why?

Yes – the variant records have the same checksum, 09dd51431f385ed3863928eee0064125. Two very different engines – one written as code, one drawn on a canvas – ran the same programs, in the same versions, with the same parameters, on the same data. The result depends on data, software and settings, not on the engine. (Galaxy’s VCF files differ from yours in their headers, which record dates and file names – the reason we compare the records.)

Activity 7.5 · A workflow from the history

In the history, press Extract workflow. Galaxy turns the jobs in your history into a workflow. Compare it with the one you drew – what is different?
On the Workflows page, save your workflow into your project with to project (it becomes workflow/galaxy/….ga), or download it with .ga. Open the .ga file in the Files tab: it is JSON text describing every step.

What would a colleague need to re-run your Galaxy workflow and get the same result?

The .ga file (it lists every tool by Tool Shed ID and version, every parameter and every connection), a Galaxy server with those tool versions installed, and the same input data, uploaded with the right datatypes. A shared history adds the full record of each dataset. The weak point is the server: re-running depends on a server that still offers those exact tool versions. (The extracted workflow also contains any extra jobs you ran – here the separate minimap2 run – so check it before sharing.)

An AI agent: results without provenance

≈ 30 min

An AI agent can do the whole analysis for you and hand you a file. Convenient – but can you check how it was made?

Before you start

The agent in this practical is a scripted simulation: it does not run any analysis. It was written to show, in miniature, problems people meet with AI tools – confident answers that are not based on your data, invented details, and answers that change when you ask again. Your job is to find out what its result is worth from the evidence alone, as if it were real. Real agents can run tools – and when they do, the same questions apply: which tool, which version, which data, and where is the record?

If you have switched the assistant to live mode (⚙), use the suggested tasks under the box for 8.1–8.4: a task you type goes to a real agent instead (see the extension at the end of this chapter).

8.1 Give the agent the job

Activity 8.1 · Delegate

Open the AI assistant and choose the Agent tab.
Give it the task “Call the variants for NA12878 from my reads and give me a VCF file” (the first suggestion under the box). Read what it says it did, and what it concludes.

Before you look at the file: what would you need to know to trust this result?

Which data it used (your reads? which reference?), which programs and versions, which parameters and filters, what happened at each step (logs), and ideally a way to re-run it and get the same file. In short: its provenance – the same things your Snakefile and Galaxy history record.

8.2 Check the file

Activity 8.2 · Evidence

The records: grep -v '^##' ~/ai-agent/NA12878_variants.vcf | head -5
Which software and reference does the file’s header name? grep -E '^##(source|reference|contig|GATKCommandLine)' ~/ai-agent/NA12878_variants.vcf
Is that software available to you at all? type gatk bwa conda list | grep -i gatk
CYP2C19*2 (rs4244285) is at hg19 chr10:96,541,616 – position 11,616 of your human_CYP2C19 slice. Your pipeline’s call: bcftools view -H results/variants/NA12878.filtered.vcf.gz | grep -w 11616 The agent’s: grep rs4244285 ~/ai-agent/NA12878_variants.vcf
The agent also reports CYP2C19*17 (rs12248560, GRCh38 chr10:94,761,900 = hg19 chr10:96,521,657). Look up the region your reads come from: grep -A2 human_CYP2C19 /data/course/provenance.json | head -8

List at least three pieces of evidence that the agent’s file was not made from your reads with your reference.

The positions are on GRCh38 chr10 (header: Homo_sapiens_assembly38.fasta, contig length 133,797,422) – you gave it hg19 slices called human_CYP2C19 and human_CYP2C9.
The header names GATK HaplotypeCaller 4.5.0.0 – not installed anywhere in your environment.
It calls rs4244285 homozygous (1/1) with all 104 reads showing A (AD 0,104), but your reads show both alleles: DP4 40,4,46,2 – about 44 high-quality reference bases and 48 alternative ones, a heterozygote (0/1).
It reports rs12248560 at GRCh38 chr10:94,761,900 – that is hg19 96,521,657, outside the region the reads were collected from (96,530,001–96,560,000): there are no reads there to call anything from.

The position of rs4244285 on GRCh38 (chr10:94,781,859) and hg19 (chr10:96,541,616) can be checked in ClinVar; rs12248560 in ClinVar too.

8.3 Ask again – and ask for proof

Activity 8.3 · Reproducible?

Ask the agent to run the analysis again. Read its new conclusion.
md5sum ~/ai-agent/*.vcf grep rs4244285 ~/ai-agent/*.vcf
Ask it which exact commands and software versions it used. Compare the answer with the two file headers.
Ask it to show the log files and the BAM file.
Finally ask: did you actually run these tools on my reads?

The two runs disagree (1/1 and 0/1 for rs4244285, “poor” and “intermediate” metaboliser). What would you have concluded if you had asked only once – and how could you have known which run to believe?

You would have reported whichever answer you happened to get – and with a single run nothing looks wrong. From the agent’s output alone there was no way to decide between them: no record of the inputs, programs, versions or parameters, no logs, no intermediate files, and no way to re-run the same process. Only your own reads (8.2) could settle it. Its excuse – that variant calling is “stochastic” – does not hold either: a genotype changing from 1/1 to 0/1 at about 100 reads is not a small difference, and your pinned pipeline gave the same checksum in the browser, in Galaxy and on another computer. Only an analysis you can re-run, and that gives the same answer when you do, can be checked.

When asked, the agent listed commands and versions. Why is that not provenance?

It is a story told afterwards, not a record made while running. Nothing links it to the file: it contradicts the second run (BWA-MEM2 2.2.1 and GATK 4.4.0.0), names files that do not exist (the BAM, the logs – “no longer available”), and in the end the agent admits it did not run anything. Provenance has to be produced by the process itself – like Snakemake’s metadata or Galaxy’s job records – and be checkable, for example by re-running.

8.4 An agent that works in the open

AI agents can be useful here – if they work through your pipeline, where every step leaves a record. (Agents that really run tools usually show the commands they ran; that transcript is part of the provenance, but only part: it does not record the software versions or prove which files were used. A pipeline does.)

Activity 8.4 · A glass box

Give the agent the task “Run my Snakemake pipeline for me and tell me what you did”. Then look at the Terminal: its commands are there, with their complete output.

In Claerbout’s terms, what is the difference between what the two agents gave you?

The first agent gave you advertising: a file and a confident story, with no code, data record or environment behind it – nothing that can be inspected or re-run. The second ran your pipeline, so the scholarship exists: the Snakefile, the config, the environment files, the data with checksums, and Snakemake’s record of every job – visible in your terminal and in .snakemake/metadata. It does not matter who (or what) typed the commands, as long as the record is made by the process and can be checked.

Extension: a real agent Extension

With a live model connected (⚙ in the assistant – for example a free Google Gemini key), a task you type in the Agent tab goes to a real agent: the model chooses a command, runs it in your terminal, reads the output and decides what to do next. It works in a folder of its own, ~/ai-agent/live-1 (then live-2, …), with copies of the course data; anything it changes elsewhere is undone. While it works, the terminal is locked for typing; Stop ends the task. Each step is a request to the AI service, so a task uses up to 16 requests – and on Google’s free tier each model answers only a limited number per day.

Extension · A real agent

In the Agent tab, type the task “Call the variants for NA12878 from the reads in input/ and give me a filtered VCF file”. Watch its commands in the Terminal, then read its report.
Give it exactly the same task again, and compare the two runs: the commands and options it chose, the files it made, and the number of variants. For each run’s VCF, a checksum of the records: bcftools view -H FILE | md5sum.

This agent really ran the tools on your reads. Is its result as trustworthy as your pipeline’s? What record is there of how it was made, and what is missing?

It is far better than the simulation: the commands, their output and the files are real, and you saw them happen. But it is not a pipeline. The commands and options – mapping preset, read groups, calling model, filters – were chosen during the run; a second run often chooses differently, and its results then differ. The transcript stays in the browser tab, not with the results; nothing pins the software versions; and giving the agent the task again does not repeat the same steps. To make the result reproducible, turn its commands into rules in your Snakefile, with pinned environments: use the agent to help write the pipeline, not to produce the result.

Using AI in an analysis

Use AI to write and explain code you can read, run, test and keep – a Snakefile rule, a script – not to produce results you cannot check.
Never accept a result without provenance: which data, which software and version, which parameters – recorded by the process, not described afterwards.
Check AI output: run it, read the errors, look at file headers, compare checksums with a trusted run.
Say where AI was used (tool, version, for what) in your methods or README.
Never paste personal or patient data into an AI tool.

Clean up Tidy

The agent’s files should not be mistaken for results later: rm -r ~/ai-agent

Burn, build – and which is best?

≈ 30 min

Put your pipeline through Claerbout’s test, compare four ways of running an analysis, and take your reproducible project with you.

9.1 Burn, build, view

Claerbout’s test: can every result be deleted and rebuilt, with one command, identically?

Activity 9.1 · The test

Burn – let Snakemake delete every output file it made (log files and empty folders stay): snakemake --delete-all-output then ls -R results
Build – one command: snakemake --cores 1 --use-conda
View – the same result? bcftools view -H results/variants/NA12878.filtered.vcf.gz | md5sum and open results/report/variant_quality.png

Extension: clean Extension

The fourth command, clean, removes intermediate files. In Snakemake you mark them with temp(): change the output of call_variants to temp("results/variants/{sample}.raw.vcf.gz") and save. Then burn and build again – snakemake --delete-all-output and snakemake --cores 1 --use-conda – and see Snakemake delete the raw VCF as soon as filter_variants has used it. (A file that already existed is only treated as temporary once it is made again. Keep intermediates you might want to inspect – such as the BAM file – until you are sure.)

9.2 Which is best?

	Bash script	Snakemake	Galaxy	AI agent (black box)
The pipeline is…	a text file of commands	a text file of rules (Python)	a graph on a server (.ga file)	a conversation
Record of how each result was made	none (unless you add it)	per output: code, params, inputs, environment, logs	per dataset: tool, version, parameters, command	none – a story told afterwards
Software versions	whatever is installed	declared per rule (conda / containers)	fixed by each tool version on the server	unknown
Re-runs only what is needed	no	yes	per job, by hand	no – and may answer differently
Stops cleanly on errors	only with strict mode	yes, and removes broken outputs	yes (red datasets)	you cannot tell
Many samples, clusters, cloud	hard	built in	possible on a large server	–
Getting started	easy	some Python and command line	easiest – forms and menus	easiest of all
Sharing	the file (if it runs elsewhere)	a folder or git repository; runs wherever conda has the pinned packages	a .ga file, histories – needs a server with the same tools	–

Which approach would you choose – and why – for (a) a diagnostic lab running the same analysis on hundreds of exomes a week, under audit; (b) a biologist with no programming experience analysing three samples and sharing the analysis with a collaborator; (c) trying out a new variant caller quickly?

(a) A code-based workflow engine (Snakemake, or Nextflow) with pinned environments or containers, under version control – it scales, runs unattended and records the provenance of every result; a validated Galaxy server is a possible alternative. (b) Galaxy: no programming, a complete history for every dataset, and workflows and histories that can be shared on the same server. (c) Ask an AI assistant to help write the new rule – then add it to your pipeline, so that it runs with a record; not an agent that returns a file. An agent-produced result without provenance is not acceptable in any of the three.

Which did you find easiest to use, and which gives you the most confidence in a result? Are they the same?

Personal answers. Often the easiest (the agent, or Galaxy) is not the one that gives the most confidence; confidence comes from a record you can check and a result you can regenerate.

9.3 Export your project

Finish the README: replace its How to run section – the placeholder – with these two sections, and fill in your name and the date. Then download the project. It runs with conda on Linux, on a Mac, or on Windows with WSL.

## How to run
Needs conda (e.g. Miniforge) and Snakemake 9, on Linux, a Mac or Windows with WSL:

    conda create -n snakemake -c conda-forge -c bioconda snakemake=9.27
    conda activate snakemake
    snakemake --cores 1 --use-conda

On a Mac with an Apple processor (M1 or later) the pinned tool versions exist
only for Intel Macs, so run  CONDA_SUBDIR=osx-64 snakemake --cores 1 --use-conda
(macOS runs them with Rosetta 2).

Settings are in config/config.yaml; the software of each step is in workflow/envs/.

## Checking a re-run
The variant records of a correct run have this checksum (on a Mac: md5 instead of md5sum):

    gunzip -c results/variants/NA12878.filtered.vcf.gz | grep -v '^#' | md5sum
    09dd51431f385ed3863928eee0064125  -

Activity 9.2 · Take it with you

Complete README.md and save it.
– everything needed to run the pipeline, without the results (they can be rebuilt).

Run it on your own computer

Install Miniforge (conda with the conda-forge channel).
Make an environment with Snakemake: conda create -n snakemake -c conda-forge -c bioconda snakemake=9.27, then conda activate snakemake.
Unzip the project, cd cyp2c19-pipeline, and run snakemake --cores 1 --use-conda. Snakemake downloads and builds the three environments the first time. On a Mac with an Apple processor (M1 or later), run CONDA_SUBDIR=osx-64 snakemake --cores 1 --use-conda instead: these older tool versions were only ever built for Intel Macs, and macOS runs them with Rosetta 2. (The computer is part of the environment, too.)
Check: bcftools is not on your computer outside the environments, so look at the result with gunzip -c results/variants/NA12878.filtered.vcf.gz | grep -v '^#' | md5sum (on a Mac: md5 instead of md5sum) – it should print 09dd51431f385ed3863928eee0064125.

We ran the exported project this way on Linux (Snakemake 9.27.0, conda 26.7): same checksum. With CONDA_SUBDIR=osx-64 the same three environments resolve for macOS, with the same versions (checked with conda; not run on a Mac).

Does your exported project meet Claerbout’s principle? What does it contain – and what is still missing?

It contains the data (with checksums and its source in the README), the code (Snakefile and script), the settings (config file), the software environment as pinned environment files, and instructions. Missing or weaker: an exact lock of every package and build (conda env export or a container image), the provenance records themselves (.snakemake/metadata, logs) if you do not include results, the operating system and processor type (these old tool versions exist only for Linux and Intel Macs), a version-control history of how the code changed, and tests. Good enough for a reader to regenerate and check the result – which is the point.

Reference

≈ 5 min

Commands, syntax and words from this practical, in one place.

Snakemake

`snakemake -n`	dry run: show the jobs and the reason for each, run nothing
`snakemake --cores 1`	run the pipeline (build the first rule – usually `all`)
`snakemake --cores 1 FILE`	build one file (and what it needs)
`-p`	print the shell commands
`--use-conda`	run each rule in its `conda:` environment
`-R RULE` / `-F`	force a rule (and everything after it) / everything to run
`--config key=value`	override a setting of the config file
`--dag \| dot -Tsvg > dag.svg`	draw the graph of jobs (`--rulegraph`: of rules)
`--summary`, `-D`	the record of every output: rule, date, status (and inputs, command)
`--delete-all-output`	delete every output file the pipeline made (not logs)
`--lint`	check the Snakefile for common problems

A rule

rule call_variants:                        # a name, a colon, an indented body
    input:
        ref=REF,                           # named inputs: {input.ref} in the command
        bam="results/mapped/{sample}.sorted.bam",   # {sample} is a wildcard
    output:
        "results/variants/{sample}.raw.vcf.gz",
    log:
        "logs/call_variants/{sample}.log",
    params:
        depth=250,                         # {params.depth}
    conda:
        "envs/calling.yaml"                # with --use-conda
    shell:
        "bcftools mpileup -d {params.depth} -f {input.ref} {input.bam} -Ou 2> {log} "
        "| bcftools call -mv -Oz -o {output} 2>> {log}"
# braces meant for the shell are doubled:  awk '{{print $1}}'
# script: "scripts/x.py"  runs Python with a `snakemake` object (snakemake.input, snakemake.output ...)

Python

`x = "text"`, `n = 20`	variables: a string (quotes) and a number
`[a, b]`, `xs[0]`, `len(xs)`	a list, its first item, its length – mind the commas
`{"k": v}`, `d["k"]`	a dictionary and a lookup
`f"{x}.bam"`, `"{s}.bam".format(s=x)`	fill values into text now / later
`for x in xs:` + indented lines	repeat for each item
`def f(a, b=1):` … `return`	a function with a default argument

conda and the shell

`conda env list` / `conda activate NAME`	environments / switch to one
`conda list`, `conda env export`	what is installed (exact versions and builds)
`conda create -n NAME pkg=1.2`	`=1.2`: any 1.2.x; `==1.2`: exactly 1.2
`pip freeze`	Python packages as `name==version`
`set -euo pipefail`	bash strict mode: stop at the first error
`md5sum FILES > MD5SUMS`, `md5sum -c MD5SUMS`	record / check checksums
`chmod a-w FILE`	make read-only
`diff A B`	compare two text files line by line

Words

Pipeline / workflow	a set of steps executed in a defined order to process data
Workflow engine	the software that runs a pipeline (Snakemake, Nextflow, Galaxy)
DAG	directed acyclic graph: the jobs and the files that connect them
Re-entrancy	continuing where a pipeline stopped, redoing only what is needed
Provenance	the record of how a result was made: data, code, parameters, software
Environment	a set of installed programs and libraries with particular versions
Pinning	fixing the exact version of a package