Skip to content

How to benchmark a variant-calling pipeline against GIAB

Short answer: run your pipeline on reads from a Genome in a Bottle (GIAB) reference sample such as HG002, compare the resulting VCF with the GIAB truth set using hap.py and the RTG vcfeval engine, restrict the comparison to the truth set's confident regions, and report recall, precision and F1 separately for SNVs and indels. Record every input by checksum so the same benchmark can be repeated exactly.

What is GIAB benchmarking?

Genome in a Bottle is a consortium that publishes thoroughly characterised human reference samples with "truth sets": high-confidence variant calls plus a BED file of the regions where those calls can be trusted. Because the truth is known, a pipeline's output for the same sample can be scored objectively. Laboratories use this to measure how well a germline variant-calling pipeline finds real variants (recall) and avoids false ones (precision).

GIAB benchmarking uses public reference materials only. No patient data is involved.

Which GIAB sample and truth set should I use?

GIAB publishes seven samples, HG001 to HG007 (see Concepts). HG002, the Ashkenazi son, has the most complete truth set and is the usual default. For GRCh38, Aver's default is the published HG002 v4.2.1 small-variant benchmark. GIAB's newer v5.0q set is still a draft, so Aver never selects it unless you ask for it explicitly.

Whichever release you choose, confirm it against the GIAB source and record the SHA-256 of every file you download. A benchmark is only repeatable if the truth set is identified exactly.

What are the steps of a GIAB benchmark?

  1. Get verified reference data. Download the truth VCF, the confident-regions BED and the matching reference genome build, and check each file against its published checksum.
  2. Get the reads. Use public GIAB sequencing data for the same sample. A single chromosome, such as chr20, keeps a first run laptop-sized.
  3. Run the pipeline exactly as validated. Use the same revision, parameters, containers and reference genome you use in production.
  4. Compare with hap.py using the vcfeval engine. Restrict the comparison to the confident regions, intersected with your panel or exome BED if you have one.
  5. Read the metrics per variant class. Report SNVs and indels separately, using the PASS rows as headline figures.
  6. Record everything. Keep the pipeline version, every container digest, the truth-set release, the comparator version and the outputs, so the benchmark can be re-run and audited.

Why use hap.py with the vcfeval engine?

The same variant can be written in more than one way in a VCF, particularly indels in repetitive sequence. A naive position-by-position comparison counts those as mismatches. hap.py with RTG vcfeval compares haplotypes, so equivalent representations match. It is the standard GA4GH benchmarking approach used with GIAB small-variant truth sets.

Run hap.py from a container pinned by digest rather than a mutable tag. The comparator's version is part of the result: a different version can change the numbers.

How do I read the results?

Metric Meaning
Recall TP / (TP + FN): the share of true variants the pipeline found
Precision TP / (TP + FP): the share of the pipeline's calls that are true
F1 The harmonic mean of recall and precision

For example, one real run of nf-core/sarek 3.10.0 on HG002 chr20 at roughly 30x coverage measured SNV recall of 99.4294% and indel recall of 97.8323%. The full record is in the worked example. Those figures describe one configuration. They are not a performance claim for sarek or GATK, and your laboratory sets its own acceptance criteria.

To see where a pipeline misses variants, add the GIAB stratifications. They split the metrics by genomic context such as low-mappability regions, segmental duplications and homopolymers. In the worked example, recall outside the difficult regions was above 99.9% for both classes, and most misses sat in low-mappability sequence and tandem repeats.

What are the common mistakes?

  • Scoring outside the confident regions. Calls outside them cannot be judged and distort both recall and precision.
  • Unverified reference files. A truth set or reference genome that differs from the published file silently changes every number. Check checksums before every use.
  • Unpinned comparator or containers. A rebuilt image under the same tag can change results without anyone noticing.
  • Ignoring filtering. If PASS and ALL rows are identical, your pipeline applied no variant filtering in that configuration. In the worked example, sarek with HaplotypeCaller filters only when --dbsnp or --known_indels is set.
  • Comparing against a different release or build from the one your baseline used, then reading the difference as a pipeline change.

How do I do this with Aver?

The open-source CLI automates the steps above, with checksum-verified reference data, a digest-pinned hap.py and a deterministic report:

pip install aver-bio
aver init --preset nf-core/sarek --clinical
aver resolve        # fingerprint every dependency
aver benchmark      # HG002 · GRCh38 · chr20 by default

If your pipeline runs elsewhere, compare an existing VCF directly:

aver compare --query out.vcf.gz --sample HG002 --stratification core

The VCF's sample name must be a registered reference material, so a clinical sample can never be compared by mistake. See the quickstart and the CLI reference.

Aver produces validation evidence. Decisions about whether a pipeline is fit for clinical use remain with the laboratory.