Variant Calling from FASTQ to VCF with GATK4 Best Practices

Turning raw sequencing reads into a trustworthy list of genetic variants takes more than one command. The GATK4 Best Practices workflow from the Broad Institute chains together alignment, duplicate marking, base quality recalibration, and statistical variant calling into a pipeline that is reproducible and defensible in a paper or a clinic. This tutorial walks through the full germline short-variant pipeline end to end, from paired-end FASTQ files to a hard-filtered VCF, with the exact command for every stage.

The Pipeline at a Glance

Every stage below reads the output of the previous one, so a mistake early on (a missing read group, an unindexed reference) propagates all the way to the final VCF. It helps to keep the whole shape in view before diving into individual commands.

Prerequisites: Reference, Tools, and Known Sites

You need three aligned tools working before touching a single FASTQ: bwa for alignment, samtools for BAM handling, and gatk4 itself. If you have not set these up yet, follow Install BWA and Bowtie2 with Conda and Install GATK4 with Conda first, then come back here.

You also need a reference genome and, for the BQSR step, a "known sites" VCF of already-catalogued variants, typically dbSNP plus the Mills and 1000 Genomes gold-standard indels, both distributed in the public GATK resource bundle. Before anything else, index the reference three ways, since bwa, samtools, and gatk each expect their own index format:

bwa index reference.fasta
samtools faidx reference.fasta
gatk CreateSequenceDictionary -R reference.fasta -O reference.dict

That gives you a BWA index, a .fai FASTA index, and a .dict sequence dictionary that GATK requires for every downstream tool.

Step 1: Align Reads with BWA-MEM

bwa mem is the standard aligner for Illumina paired-end reads longer than about 70 bp, and it is what the GATK Best Practices workflow assumes. The critical detail most people miss on their first run is the -R read group string: GATK refuses to process a BAM without one, since it uses the SM (sample) and ID fields to track which reads came from which sample and sequencing run.

bwa mem -t 8 \
  -R "@RG\tID:sample1\tSM:sample1\tPL:ILLUMINA\tLB:lib1" \
  reference.fasta sample_R1.fastq.gz sample_R2.fastq.gz \
  | samtools sort -@ 8 -o sample.sorted.bam -
samtools index sample.sorted.bam

Piping straight into samtools sort avoids writing a huge intermediate unsorted BAM to disk. The result, sample.sorted.bam, is coordinate-sorted and indexed, which every subsequent GATK step expects.

Step 2: Mark Duplicates

PCR amplification during library prep produces multiple read pairs that are optical or PCR copies of the same original DNA fragment, not independent evidence. Left uncorrected, these duplicates inflate coverage and bias variant allele fractions. MarkDuplicates flags them (it does not delete them) so downstream tools can ignore them.

gatk MarkDuplicates \
  -I sample.sorted.bam \
  -O sample.dedup.bam \
  -M sample.dedup.metrics.txt

samtools index sample.dedup.bam

The sample.dedup.metrics.txt file reports the duplication rate. A well-prepared library typically shows single-digit to low double-digit percent duplication; a rate above roughly 30-40 percent usually points to low input DNA or over-amplification during library prep.

Step 3: Base Quality Score Recalibration (BQSR)

Sequencer-reported base quality scores are systematically biased in predictable, machine-specific ways (by cycle, by dinucleotide context, by lane). BQSR models that bias empirically by comparing observed mismatches against a set of known variant sites, on the assumption that mismatches at known-variant positions are real biology while mismatches everywhere else are sequencing error. It then rewrites the quality scores to be statistically accurate, which matters because HaplotypeCaller's genotype likelihoods depend directly on those scores.

gatk BaseRecalibrator \
  -I sample.dedup.bam \
  -R reference.fasta \
  --known-sites dbsnp.vcf.gz \
  --known-sites Mills_and_1000G_gold_standard.indels.vcf.gz \
  -O sample.recal.table

gatk ApplyBQSR \
  -I sample.dedup.bam \
  -R reference.fasta \
  --bqsr-recal-file sample.recal.table \
  -O sample.recal.bam

sample.recal.bam is now "analysis-ready" in GATK terminology: sorted, deduplicated, and recalibrated. This is the BAM every remaining step reads from.

Step 4: Call Variants with HaplotypeCaller

Unlike a naive pileup-based caller that looks at one position at a time, HaplotypeCaller locally reassembles the reads in each active region into candidate haplotypes using a de Bruijn-like graph, then computes genotype likelihoods for each haplotype against the data. That local reassembly is what makes it accurate around indels and closely spaced variants where simple pileup callers struggle.

For a single sample you can call genotypes directly. For anything beyond a single sample, the Best Practices recommend calling each sample in GVCF mode first, then joint-genotyping across the cohort, because it lets you add more samples later without recalling everyone from scratch.

gatk HaplotypeCaller \
  -R reference.fasta \
  -I sample.recal.bam \
  -O sample.g.vcf.gz \
  -ERC GVCF

gatk GenotypeGVCFs \
  -R reference.fasta \
  -V sample.g.vcf.gz \
  -O sample.vcf.gz

sample.vcf.gz now holds raw, unfiltered genotype calls for both SNPs and indels. "Raw" here means every candidate that passed HaplotypeCaller's internal calling threshold is present, including plenty of low-confidence artifacts that still need to be filtered out.

Step 5: Hard-Filter the Callset

The Broad recommends Variant Quality Score Recalibration (VQSR) when you have a large cohort, since it trains a Gaussian mixture model on true-positive resources like HapMap and Omni. For a single sample or a small project, VQSR does not have enough data to train reliably, so the documented fallback is hard-filtering: split SNPs and indels apart (their error profiles differ) and threshold on a handful of quality annotations.

gatk SelectVariants -R reference.fasta -V sample.vcf.gz \
  --select-type-to-include SNP -O sample.snps.vcf.gz

gatk SelectVariants -R reference.fasta -V sample.vcf.gz \
  --select-type-to-include INDEL -O sample.indels.vcf.gz

Filtering SNPs

gatk VariantFiltration \
  -R reference.fasta -V sample.snps.vcf.gz \
  --filter-expression "QD < 2.0" --filter-name "QD2" \
  --filter-expression "FS > 60.0" --filter-name "FS60" \
  --filter-expression "MQ < 40.0" --filter-name "MQ40" \
  --filter-expression "MQRankSum < -12.5" --filter-name "MQRankSum-12.5" \
  --filter-expression "ReadPosRankSum < -8.0" --filter-name "ReadPosRankSum-8" \
  --filter-expression "SOR > 3.0" --filter-name "SOR3" \
  -O sample.snps.filtered.vcf.gz

Filtering Indels

Indels get a looser strand-bias threshold because their error signature is noisier, and MQRankSum/MQ are dropped since mapping quality is less informative for gapped alignments.

gatk VariantFiltration \
  -R reference.fasta -V sample.indels.vcf.gz \
  --filter-expression "QD < 2.0" --filter-name "QD2" \
  --filter-expression "FS > 200.0" --filter-name "FS200" \
  --filter-expression "ReadPosRankSum < -20.0" --filter-name "ReadPosRankSum-20" \
  --filter-expression "SOR > 10.0" --filter-name "SOR10" \
  -O sample.indels.filtered.vcf.gz

Each command above tags failing records in the FILTER column with the matching filter name (e.g. QD2) rather than deleting them, so you always keep an audit trail. Merge the two filtered files back together with gatk MergeVcfs if you want a single output VCF, and in practice you would then subset to records where FILTER reads PASS for downstream analysis.

Common Pitfalls

Next Steps

This pipeline scales to many samples fastest when it is not run by hand for each one. A workflow manager such as Snakemake or Nextflow turns these commands into a reusable, parallelized pipeline, and pairing GATK4 with samtools/bcftools covers nearly everything you need for routine BAM and VCF manipulation. For cohorts beyond a handful of samples, look into GATK's joint-genotyping and VQSR workflows as the natural next step past the hard-filtering approach shown here.