Orchestrates the end-to-end germline short-variant pipeline from FASTQ to a filtered, normalized, benchmarked VCF, chaining QC/trim, BWA-MEM2 alignment, duplicate marking, optional BQSR, calling (bcftools/GATK HaplotypeCaller/DeepVariant/DRAGEN), normalization, site+genotype filtering, annotation, and hap.py/vcfeval benchmarking. Use when deciding the pipeline-wide reference-genome commitment (GRCh38 analysis set vs T2T, ALT/decoy handling), sequencing the steps in the defensible order (normalize BEFORE annotate, filter site- then genotype-level), choosing the calling engine and single-sample vs cohort joint-calling, picking a filtering strategy by cohort size, or benchmarking stratified within GIAB confident regions. Hands off mechanism to the variant-calling and read-alignment component skills; not a re-teach of any single step.
Orchestrates the end-to-end germline short-variant pipeline from FASTQ to a filtered, normalized, benchmarked VCF, chaining QC/trim, BWA-MEM2 alignment, duplicate marking, optional BQSR, calling (bcftools/GATK HaplotypeCaller/DeepVariant/DRAGEN), normalization, site+genotype filtering, annotation, and hap.py/vcfeval benchmarking. Use when deciding the pipeline-wide reference-genome commitment (GRCh38 analysis set vs T2T, ALT/decoy handling), sequencing the steps in the defensible order (normalize BEFORE annotate, filter site- then genotype-level), choosing the calling engine and single-sample vs cohort joint-calling, picking a filtering strategy by cohort size, or benchmarking stratified within GIAB confident regions. Hands off mechanism to the variant-calling and read-alignment component skills; not a re-teach of any single step.
[{"after_qc":"Q30 >85%, adapter content <1%"},{"after_alignment":"Mapping rate >95%, properly paired >90%"},{"after_dedup":"Duplication rate <30% for WGS, <50% for exome"},{"after_calling":"Ti/Tv ratio ~2.0-2.1 for WGS, ~3.0-3.3 for exome; dbSNP overlap >95% only after annotating the ID column"}]
Before using code patterns, verify installed versions match. If versions differ:
CLI: <tool> --version then <tool> --help to confirm flags
If code throws ImportError, AttributeError, or TypeError, introspect the installed
package and adapt the example to match the actual API rather than retrying.
Note: GenotypeGVCFs defaults (--max-alternate-alleles, --heterozygosity, --stand-call-conf), the GATK hard-filter thresholds, and DRAGEN speed/accuracy figures drift by version/vendor; confirm in-tool and against current GIAB benchmarks before quoting.
FASTQ to Variants Workflow
"Call variants from my whole-genome or exome FASTQ files" -> Chain QC/trim, alignment, duplicate marking, an engine-appropriate caller, normalization, filtering, annotation, and benchmarking into one filtered germline VCF.
This is a workflow skill: it owns the chaining decisions and hand-offs, not the internals of any one step. Every step below cross-references the component skill that teaches its mechanism.
The governing principle
A germline pipeline is a chain of commitments, and the two that decide whether the callset is trustworthy are made at the seams between steps, not inside them.
The reference is a pipeline-wide commitment made once and inherited by everything downstream. The build and analysis set chosen at alignment (step 2) fix the coordinates of every later comparison: the dbSNP/ClinVar/gnomAD records annotation matches against, the truth BED benchmarking scores within, and the cohort other samples are joint-called with must all be the SAME build. Changing it later means redoing alignment, calling, normalization, and annotation. Decide before aligning a single read.
The step order is not arbitrary; two orderings sink reviews. (a) Normalize BEFORE annotate. An indel that is not left-aligned to the database's canonical position is one base off, so the annotator silently misses the ClinVar/gnomAD record and reports a pathogenic variant as novel-absent -- a patient-safety failure that throws no error. (b) Never bcftools merge single-sample VCFs into a cohort. Absence of a record is then read as homozygous reference, fabricating genotypes; joint-genotype per-sample gVCFs instead so "confident hom-ref" is distinguished from "no data."
Filter site-level first, then genotype-level, then recompute cohort QC. Setting genotypes to no-call (./.) changes missingness, HWE, and allele frequencies, so those metrics must be computed on the genotype-filtered matrix, not before.
A single genome-wide F1 is nearly meaningless. Callers agree on easy SNPs (F1 > 0.999); they diverge in indels-in-repeats, segmental duplications, and MHC. Benchmark stratified, within the GIAB confident region, or do not claim accuracy (Krusche 2019 Nat Biotechnol 37:555-560).
Workflow overview
FASTQ
| [1] QC & trim -------------------> fastp (read-qc/fastp-workflow)
v
| [2] Align ----------------------> bwa-mem2 (read-alignment/bwa-alignment)
v ^-- reference commitment: GRCh38 analysis set / T2T, ALT/decoy handling
| [3] Mark duplicates ------------> samtools markdup (alignment-files/duplicate-handling)
v
| [4] (BQSR? optional on modern binned-quality instruments)
v
| [5] Call -----------------------> bcftools | GATK HaplotypeCaller | DeepVariant | DRAGEN
v ^-- single-sample OR per-sample gVCF -> joint-genotype (variant-calling/joint-calling)
| [6] Normalize (BEFORE annotate) -> bcftools norm -m-any -f ref (variant-calling/variant-normalization)
v
| [7] Filter: site-level THEN genotype-level (variant-calling/filtering-best-practices)
v
| [8] Recompute cohort QC on the genotype-filtered matrix (missingness, Ti/Tv, het/hom, excess-het HWE)
v
| [9] Annotate -------------------> VEP / SnpEff (variant-calling/variant-annotation)
v
| [10] Benchmark & QC ------------> hap.py/vcfeval, bcftools stats (variant-calling/vcf-statistics)
v
Filtered, normalized, benchmarked VCF
Reference genome: the pipeline-wide commitment
Decision made once, before alignment; everything downstream inherits it. Mechanism of building/indexing the reference lives in read-alignment/bwa-alignment; the deeper reasoning below is what a reviewer expects justified.
Choice
Commit to it when
Consequence inherited downstream
GRCh38 analysis set + decoys (hs38DH), ALT-aware (bwa-postalt)
Human germline, research or most clinical
Decoys soak up off-target reads; ALT-aware mapping recovers reads in MHC/segdup loci that ALT-unaware mapping force-fits to the primary, inflating false positives
Simpler and validated, but loses signal in ~5 Mb of ALT-bearing loci
Masked GRCh38 (false-duplication fix)
Calling in CBS, U2AF1, KCNE1B, KCNJ18 and other affected genes
Recovers reads whose mapQ collapsed across the phantom duplicate copy
T2T-CHM13
Research needing segdups/centromeres/dark genes; maximum accuracy
Reveals variants in newly resolved regions and removes GRCh38 false-duplication artifacts, but no lossless liftover to GRCh37/38, so the entire annotation/interpretation stack must be revalidated (Nurk 2022 Science 376:44-53; Aganezov 2022 Science 376:eabl3533)
The reason to fix this first: annotation databases, panel BEDs, benchmark truth sets, and any cohort a sample is joint-called with are all coordinate-specific. Mixing builds (e.g. normalizing to GRCh38 then annotating against a GRCh37 dbSNP) is a guaranteed silent miss. "We used GRCh38" is under-specified -- plain vs masked vs analysis-set-with-decoys materially changes results in named clinical genes.
The canonical order and why
Each step assumes the previous; the order is defensible under review (canonical preprocessing order, filtering/representation practice).
QC/trim -- remove adapters and low-quality tails before they corrupt alignment and duplicate detection.
Align -- to the committed reference, with read groups (SM/ID/PL/LB); read groups are a hard GATK requirement.
Mark duplicates -- PCR/optical duplicates are not independent evidence; marking (not removing) lets the caller down-weight them. Skip for amplicon/UMI data.
BQSR -- honestly optional on modern instruments. BQSR corrected context/cycle-dependent miscalibration on 2010-era continuous-quality Illumina. NovaSeq/NovaSeq X emit ~4 quality bins, leaving little to recalibrate; callsets are largely unchanged with vs without it. DeepVariant explicitly recommends NOT running BQSR (its CNN learned the raw-quality error model); DRAGEN handles quality internally (and uses DRAGSTR for STR/indel error modeling). Keep it for GATK-HaplotypeCaller consistency if a frozen pipeline demands it; otherwise the modern indel-accuracy lever is STR-aware error modeling (--dragen-mode), not BQSR.
Call -- per-sample VCF, or per-sample gVCF (-ERC GVCF) if a cohort will be joint-genotyped.
Normalize BEFORE annotate/compare -- bcftools norm -m-any -f ref.fa (split multiallelics, then left-align + parsimony), against the SAME reference used for annotation. Add -a (atomize) only when the downstream database is decomposed. This is the most common real ordering bug: annotate-then-normalize attaches consequences to a non-canonical representation that fails to match the database. The binding constraint is normalize-before-ANNOTATE/COMPARE, not normalize-before-filter: GATK convention runs VariantFiltration on the raw multiallelic records (its annotations are computed on that representation), then normalizes -- bwa_gatk_workflow.sh does exactly that, while bwa_bcftools_workflow.sh normalizes first. Both are correct; neither annotates before normalizing.
Filter site-level, then genotype-level -- site filters (VQSR/hard/ML) decide whether a site is real; genotype filters (GQ/DP/allele-balance) decide whether an individual genotype is trustworthy. SNPs and indels are filtered separately (different error processes and truth resources).
Recompute cohort QC on the genotype-filtered matrix -- missingness, Ti/Tv, het/hom, excess-het HWE. Doing HWE before genotype filtering lets low-GQ garbage drive spurious deviation.
Annotate on the normalized (and, where consequence matters, haplotype-resolved) representation.
Benchmark/validate after transforms (hap.py, bcftools stats).
Choosing the calling engine
Pipeline-level selection only; the mechanism and full decision table live in variant-calling/variant-calling. Hand off there to pick, then return here for chaining.
Situation
Lean toward
Hand off to
Auditable open-source, large cohort, joint calling
Maximum throughput on Illumina, hardware available
DRAGEN (or GATK --dragen-mode for the open equivalent)
variant-calling/variant-calling
Single-sample vs cohort is a chaining decision, not a caller feature. For a cohort, emit per-sample gVCFs and joint-genotype them so a variant seen in one sample is evaluated in all (cohort rescue of low-coverage hets, squared-off genotype matrix). This is what makes the pipeline forward-compatible with new samples (the N+1 problem). Full mechanism: variant-calling/joint-calling.
Primary path: BWA-MEM2 + bcftools
Fast, dependency-light, good for exploratory work and non-model organisms; weaker on indels in homopolymers than reassembly callers.
Read groups are mandatory; add -Y (soft-clip supplementary) if structural-variant calling is downstream, and -K 100000000 for thread-count-invariant output. Reference/analysis-set choice: read-alignment/bwa-alignment.
bwa-mem2 index reference.fa # once
bwa-mem2 mem -t 8 -K 100000000 \
-R "@RG\tID:sample\tSM:sample\tPL:ILLUMINA\tLB:lib1" \
reference.fa trimmed/sample_R1.fq.gz trimmed/sample_R2.fq.gz \
| samtools view -bS - > aligned/sample.bam
# collate groups mates by name; fixmate -m adds the ms/MC tags markdup needs; markdup needs coordinate order.# Do NOT coordinate-sort before fixmate, and do NOT markdup amplicon/PCR data (use UMIs there).
samtools collate -@ 8 -O -u aligned/sample.bam \
| samtools fixmate -m -@ 8 -u - - \
| samtools sort -@ 8 -u - \
| samtools markdup -@ 8 - aligned/sample.markdup.bam
samtools index aligned/sample.markdup.bam
Step 4: Call and normalize
# Single sample (mpileup passes MQ/BQ filters into the pileup)# -a FORMAT/DP,FORMAT/AD + call -f GQ emit the per-sample DP/GQ the Step 5 genotype filter needs.
bcftools mpileup -Ou -f reference.fa -a FORMAT/DP,FORMAT/AD --max-depth 250 --min-MQ 20 --min-BQ 20 \
aligned/sample.markdup.bam \
| bcftools call -mv -f GQ -Oz -o variants/sample.vcf.gz
# Normalize BEFORE any annotation or cross-callset comparison, against the SAME reference
bcftools norm -m-any -f reference.fa -Oz -o variants/sample.norm.vcf.gz variants/sample.vcf.gz
bcftools index variants/sample.norm.vcf.gz
For multi-sample cohorts, bcftools can call several BAMs jointly, but the GATK/DeepVariant gVCF path is preferred at scale (variant-calling/joint-calling).
Step 5: Filter (site then genotype)
# Site-level (bcftools flags rather than removes, so failures stay auditable)
bcftools filter -Oz -s LowQual \
-e 'QUAL<20 || INFO/DP<10 || MQ<30' \
-o variants/sample.siteflt.vcf.gz variants/sample.norm.vcf.gz
# Genotype-level: set low-confidence genotypes to no-call (NOT 0/0)
bcftools filter -Oz -S . \
-e 'FMT/GQ<20 | FMT/DP<8' \
-o variants/sample.filtered.vcf.gz variants/sample.siteflt.vcf.gz
bcftools index variants/sample.filtered.vcf.gz
Alternative path: BWA-MEM2 + GATK HaplotypeCaller
Local reassembly + PairHMM; the auditable reference implementation, strong on indels. Full mechanism: variant-calling/gatk-variant-calling.
The site-level filter is chosen by cohort size, platform, and organism; genotype-level filtering is always applied on top. Full mechanism and thresholds: variant-calling/filtering-best-practices.
Cohort / data
Site-level filter
Why
Large WGS cohort (~30+ jointly genotyped)
VQSR (or AS_VQSR for huge cohorts)
The Gaussian-mixture model needs tens of thousands of variants and truth-resource overlap to fit; unreliable below that
Single sample / small cohort
GATK hard filters or VETS/NVScoreVariants
VQSR is non-identifiable on few variants; a "converged" model on one exome is filtering on noise
Exome specifically
Hard filters (do NOT use DP as a VQSR annotation)
Capture-boundary coverage cliffs break the annotation manifold
Non-model organism
Hard filters or a bootstrapped truth set
No HapMap/Omni/Mills truth resources exist
DeepVariant / DRAGEN output
Use the caller's own calibration; do NOT re-apply GATK hard filters
Their error modes differ; classic annotations do not describe them
SNPs and indels are filtered separately (different error processes, truth resources, abundance). Hard-filter starting points (SNPs QD<2, FS>60, MQ<40, MQRankSum<-12.5, ReadPosRankSum<-8, SOR>3; indels loosen FS>200, tighten ReadPosRankSum<-20) are lenient heuristics to tune, not universal truth. RankSum annotations are only defined at het sites -- a hand-written filter must treat a missing annotation as PASS, or every hom-alt site vanishes.
Benchmarking the pipeline
Goal: a defensible accuracy statement, not a single number.
The only rigorous way to compare a callset to truth is haplotype-aware, stratified, and confined to the truth set's confident region. Two VCFs can encode the identical haplotype with different records, so a naive bcftools isec/line-diff overcounts errors; use hap.py wrapping the vcfeval engine, which replays variants onto the reference and matches at the haplotype level (Krusche 2019 Nat Biotechnol 37:555-560; GIAB truth, Zook 2019 Nat Biotechnol 37:561-566).
# Only meaningful when the sample IS a GIAB genome (HG001-HG007) with a truth VCF + confident BED.# -f = confident/callable region BED (TP/FP/FN counted ONLY inside it; calls outside are UNK, not FP)
hap.py truth.vcf.gz query.norm.vcf.gz \
-f HG002_confident.bed \
-r reference.fa \
-o bench/hg002 \
--engine=vcfeval \
--stratification stratification.tsv # GIAB region BEDs: low-complexity, segdup, MHC, GC-extreme
Discipline that separates a senior benchmark from a naive one:
Refuse the global F1. Report SNP and INDEL separately, and show the low-complexity/segmental-duplication/MHC rows explicitly -- hiding them behind an all-regions average is the most common soft cheat.
Do not benchmark an ML caller only on its training genome. DeepVariant and DRAGEN-ML train on GIAB coordinates, so scoring on HG002 alone partly measures memorization; score on a held-out or semi-blinded sample (HG003/HG004).
"We found variants GIAB missed" is almost always a category error -- calls outside the confident region, which GIAB declined to adjudicate, not accuracy. Real hard-region claims cite CMRG or an assembly-based benchmark.
When the sample is not a GIAB genome (the usual case), there is no truth VCF; fall back to proxy QC -- Ti/Tv (~2.0-2.1 WGS, ~3.0-3.3 exome), dbSNP overlap, het/hom by ancestry, and trio Mendelian concordance if a family is available. These are proxies for the absence of a benchmark, not a substitute for one (variant-calling/vcf-statistics).
Ti/Tv sliding toward 0.5 (random) signals false-positive inflation -- filters too loose. bcftools stats counts known sites from the ID column, and callers leave it as .: run bcftools annotate -c ID -a dbsnp.vcf.gz first or the dbSNP-overlap line reads 0% regardless of quality
Common Errors
Symptom
Cause
Fix
Annotation reports a known pathogenic variant as novel/absent
Annotated before normalizing; indel one base off the database coordinate
bcftools norm -m-any -f ref against the SAME build as the annotation DB, BEFORE annotation
Cohort has impossible all-hom-ref genotypes at variant sites
Built the cohort by bcftools merge of single-sample VCFs
Emit per-sample gVCFs and joint-genotype (variant-calling/joint-calling)
Every hom-alt site filtered out
Hand-written filter treats missing RankSum as failing
Treat missing annotation as PASS; RankSum is defined only at het sites
VQSR "converged" on one exome but the callset is garbage
VQSR needs tens of thousands of variants across ~30+ samples
Use hard filters or VETS/NVScoreVariants for single samples/exomes
Spurious variants in CBS/U2AF1/KCNE1B
GRCh38 false duplications collapse mapQ
Use a masked GRCh38 or T2T-CHM13; commit the reference before calling
GATK error "sample ... has no read group"
Read groups omitted at alignment
Re-run bwa-mem2 mem -R "@RG\t..." (SM/ID/PL/LB)
Different variant counts from vt vs bcftools on the same data
vt decomposes MNPs by default, bcftools does not
Standardize ONE normalization tool + flags across every cohort compared (variant-calling/variant-normalization)
Pipeline map (hand-offs)
read-qc/fastp-workflow -- QC/trim options and report interpretation
variant-calling/variant-calling - Engine selection and bcftools calling options
variant-calling/gatk-variant-calling - GATK HaplotypeCaller and DRAGEN mode
variant-calling/deepvariant - Deep-learning calling and GLnexus cohorts
variant-calling/joint-calling - Cohort joint genotyping and scaling
variant-calling/variant-normalization - Normalize before annotate/compare
variant-calling/filtering-best-practices - VQSR, hard filters, VETS
variant-calling/variant-annotation - Annotate variants with VEP
variant-calling/vcf-statistics - Ti/Tv, het/hom, and identity QC
References
Krusche P, Trigg L, Boutros PC, et al. (GA4GH Benchmarking Team). Best practices for benchmarking germline small-variant calls in human genomes. Nature Biotechnology 37:555-560 (2019). DOI 10.1038/s41587-019-0054-x. Stratified haplotype-aware benchmarking (hap.py/vcfeval).
Zook JM, McDaniel J, Olson ND, et al. An open resource for accurately benchmarking small variant and reference calls. Nature Biotechnology 37:561-566 (2019). DOI 10.1038/s41587-019-0074-6. GIAB truth set + confident regions.
DePristo MA, Banks E, Poplin R, et al. A framework for variation discovery and genotyping using next-generation DNA sequencing data. Nature Genetics 43:491-498 (2011). DOI 10.1038/ng.806. GATK framework.
Van der Auwera GA, Carneiro MO, Hartl C, et al. From FastQ Data to High-Confidence Variant Calls: The Genome Analysis Toolkit Best Practices Pipeline. Current Protocols in Bioinformatics 43:11.10.1-11.10.33 (2013). DOI 10.1002/0471250953.bi1110s43.
Poplin R, Ruano-Rubio V, DePristo MA, et al. Scaling accurate genetic variant discovery to tens of thousands of samples. bioRxiv 201178 (2018). DOI 10.1101/201178. Preprint only (never journal-published); the GVCF/joint-genotyping reference.
Poplin R, Chang P-C, Alexander D, et al. A universal SNP and small-indel variant caller using deep neural networks. Nature Biotechnology 36:983-987 (2018). DOI 10.1038/nbt.4235. DeepVariant.
Yun T, Li H, Chang P-C, et al. Accurate, scalable cohort variant calls using DeepVariant and GLnexus. Bioinformatics 36:5582-5589 (2020). DOI 10.1093/bioinformatics/btaa1081.
Danecek P, Bonfield JK, Liddle J, et al. Twelve years of SAMtools and BCFtools. GigaScience 10:giab008 (2021). DOI 10.1093/gigascience/giab008.
Nurk S, Koren S, Rhie A, et al. The complete sequence of a human genome. Science 376:44-53 (2022). DOI 10.1126/science.abj6987. T2T-CHM13.
Aganezov S, Yan SM, et al. A complete reference genome improves analysis of human genetic variation. Science 376:eabl3533 (2022). DOI 10.1126/science.abl3533. Reference-choice variant-calling payoff.
