Running a GATK Variant Calling Pipeline on Your Own Whole Genome

By the end of this guide you will have a single-sample germline VCF for your own genome. Expect roughly 4 to 5 million variants against GRCh38, of which about 3.5 million are SNVs (single-nucleotide variants) and 500,000 to 700,000 are small insertions and deletions. Each one is annotated with genomic context and population frequency. You also keep the intermediate CRAM alignment, which you can revisit for any region you care about. You need your raw reads as paired FASTQ or an unaligned or aligned CRAM. You also need a machine with at least 32 GB of RAM and 8 cores, plus about 500 GB of free disk for a 30x genome and its intermediates. Budget a few days of wall-clock patience if you run single-threaded HaplotypeCaller across the whole genome. Software: bwa-mem2 or bwa, samtools, and gatk 4.5 or later. You also need bcftools and either VEP or SnpEff for annotation. Everything below assumes germline short-read data from one person, which is the case GATK’s best practices were designed for and where they behave best.

A framing note before the commands. Variant calling produces a hypothesis about each position in your genome. It does not produce a fact. Pipelines that differ only in aligner or caller agree on well under 100% of calls: an early comparison of five exome and genome pipelines found concordance across all of them for only about 57% of SNVs called by any one pipeline.¹ Modern tools have narrowed that gap, but the lesson holds. Treat any single variant, especially one in a gene you recognize, as a lead that needs orthogonal confirmation. Nothing here is a clinical test, and interpreting a specific variant for your own health is a conversation with a clinician or a certified genetics lab. Such interpretation is not a bioinformatics exercise.

1. Get the reference and resource bundle right

Everything downstream depends on using one reference build consistently, and mixing builds is the single most common way to silently corrupt an analysis. We use GRCh38 with the analysis set that includes decoy and alt contigs, because reads from centromeric and other repeat-rich sequence that would otherwise pile up on real genes get absorbed by the decoy instead.

# GRCh38 analysis set + GATK resource bundle (Google Cloud public bucket)
gsutil -m cp \
  gs://genomics-public-data/resources/broad/hg38/v0/Homo_sapiens_assembly38.fasta* \
  gs://genomics-public-data/resources/broad/hg38/v0/Homo_sapiens_assembly38.dict \
  gs://genomics-public-data/resources/broad/hg38/v0/Homo_sapiens_assembly38.dbsnp138.vcf* \
  gs://genomics-public-data/resources/broad/hg38/v0/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz* \
  gs://genomics-public-data/resources/broad/hg38/v0/Homo_sapiens_assembly38.known_indels.vcf.gz* \
  ./ref/

Two decisions worth making deliberately here. First, use the Homo_sapiens_assembly38.fasta that ships with the GATK bundle rather than a UCSC or Ensembl FASTA, because the .dict sequence order and contig names must match your BAM header exactly or GATK will refuse to run. Second, if your data came back aligned to GRCh37, do not lift over the VCF. Realign from the CRAM. Liftover of coordinates is lossy at exactly the regions you would want to trust, and the reference-allele mismatches it introduces are hard to audit afterwards.

2. Align with BWA-MEM2 and mark duplicates

If you have FASTQ, align directly. The read group line is not optional: GATK’s base quality recalibration and HaplotypeCaller both read @RG fields, and a missing SM tag aborts the run several hours in.

bwa-mem2 index ref/Homo_sapiens_assembly38.fasta   # ~1 hr, ~80 GB RAM peak

bwa-mem2 mem -t 16 -K 100000000 -Y \
  -R '@RG\tID:L001\tSM:SAMPLE01\tLB:lib1\tPL:ILLUMINA\tPU:L001' \
  ref/Homo_sapiens_assembly38.fasta \
  reads_R1.fastq.gz reads_R2.fastq.gz \
  | samtools sort -@ 8 -m 2G -o aligned.bam -
samtools index -@ 8 aligned.bam

-K 100000000 fixes the chunk size so output is deterministic regardless of thread count, which matters when you want to reproduce a run. -Y uses soft clipping for supplementary alignments, which keeps the full read sequence available and is what GATK’s structural-variant-adjacent tools expect.

If you were given a CRAM aligned to GRCh38 by the sequencing provider, skip this and go straight to duplicate marking, but verify the header first with samtools view -H your.cram | grep '^@SQ' | head and confirm the M5 checksums match your reference.

Duplicate marking flags read pairs with identical 5’ positions, which arise from PCR amplification and optical artifacts and would otherwise be counted as independent evidence:

gatk MarkDuplicatesSpark \
  -I aligned.bam \
  -O markdup.bam \
  -M markdup.metrics.txt \
  --optical-duplicate-pixel-distance 2500 \
  --conf 'spark.executor.cores=8'

Use 2500 for patterned flowcells (NovaSeq, HiSeq X). The default of 100 is for unpatterned older instruments and will undercount optical duplicates. Expect a duplicate fraction of 2 to 10% for a well-prepared PCR-free 30x library. Above 20%, your library was over-amplified and effective coverage is lower than the nominal figure.

3. Base quality score recalibration

Base quality scores from the sequencer are systematically biased by machine cycle and sequence context, and callers weigh evidence by those scores.² BQSR builds an empirical error model from mismatches at non-polymorphic sites and rewrites the quality strings.

gatk BaseRecalibrator \
  -R ref/Homo_sapiens_assembly38.fasta \
  -I markdup.bam \
  --known-sites ref/Homo_sapiens_assembly38.dbsnp138.vcf \
  --known-sites ref/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz \
  --known-sites ref/Homo_sapiens_assembly38.known_indels.vcf.gz \
  -O recal.table

gatk ApplyBQSR \
  -R ref/Homo_sapiens_assembly38.fasta \
  -I markdup.bam \
  --bqsr-recal-file recal.table \
  -O analysis_ready.cram

Write CRAM, not BAM, at this point. With the reference on hand CRAM is roughly 40 to 60% the size of the equivalent BAM and every downstream GATK tool reads it natively. Inspect the before/after plots (gatk AnalyzeCovariates -bqsr recal.table -plots recal.pdf) once, so you know whether your instrument’s raw scores were well-calibrated to begin with. On recent NovaSeq data with only four binned quality values, BQSR changes very little, and skipping it costs you almost nothing. On older or unusual chemistry the correction is substantial.

4. Call variants with HaplotypeCaller

HaplotypeCaller finds active regions where the data disagree with the reference and assembles the local reads into a de Bruijn graph. It then realigns reads against each candidate haplotype with a pair hidden Markov model and genotypes them.² That local reassembly is why it handles indels and clustered substitutions far better than position-by-position pileup callers.

Run it in GVCF mode even for a single sample. The GVCF records confident reference blocks alongside variants, which lets you later distinguish “reference at this position” from “no data here”, and lets you re-genotype without re-assembling if you add samples.

# Parallelize by interval. Split the genome into ~50 shards.
gatk SplitIntervals \
  -R ref/Homo_sapiens_assembly38.fasta \
  -L ref/wgs_calling_regions.hg38.interval_list \
  --scatter-count 50 \
  -O intervals/

ls intervals/*.interval_list | parallel -j 8 '
  gatk --java-options "-Xmx6g" HaplotypeCaller \
    -R ref/Homo_sapiens_assembly38.fasta \
    -I analysis_ready.cram \
    -L {} \
    -ERC GVCF \
    -G StandardAnnotation -G AS_StandardAnnotation \
    --native-pair-hmm-threads 2 \
    -O gvcf/$(basename {} .interval_list).g.vcf.gz'

gatk GenomicsDBImport --genomicsdb-workspace-path gdb \
  $(for f in gvcf/*.g.vcf.gz; do echo -n "-V $f "; done) \
  -L ref/wgs_calling_regions.hg38.interval_list

gatk GenotypeGVCFs \
  -R ref/Homo_sapiens_assembly38.fasta \
  -V gendb://gdb \
  -O raw.vcf.gz

The wgs_calling_regions.hg38.interval_list from the bundle excludes long runs of N and the worst of the unplaceable sequence, which saves hours of compute on regions you would filter anyway. On 8 cores expect 12 to 20 hours for a 30x genome with this scatter. Setting --native-pair-hmm-threads 2 per shard is usually the throughput sweet spot, and higher values scale poorly.

5. Filter

Filtering is where a single-sample pipeline diverges most from the published best practices. VQSR is GATK’s Gaussian-mixture recalibration. It wants thousands of exomes or at least 30 genomes to fit a reliable model, and a single genome will not give it one. Two workable options remain.

The first is hard filtering on the annotations GATK emits, split by variant type because SNVs and indels have different error profiles:

gatk SelectVariants -V raw.vcf.gz -select-type SNP -O snps.vcf.gz
gatk SelectVariants -V raw.vcf.gz -select-type INDEL -select-type MIXED -O indels.vcf.gz

gatk VariantFiltration -V snps.vcf.gz -O snps.filt.vcf.gz \
  -filter "QD < 2.0"                --filter-name QD2 \
  -filter "QUAL < 30.0"             --filter-name QUAL30 \
  -filter "SOR > 3.0"               --filter-name SOR3 \
  -filter "FS > 60.0"               --filter-name FS60 \
  -filter "MQ < 40.0"               --filter-name MQ40 \
  -filter "MQRankSum < -12.5"       --filter-name MQRS \
  -filter "ReadPosRankSum < -8.0"   --filter-name RPRS

gatk VariantFiltration -V indels.vcf.gz -O indels.filt.vcf.gz \
  -filter "QD < 2.0"                --filter-name QD2 \
  -filter "QUAL < 30.0"             --filter-name QUAL30 \
  -filter "FS > 200.0"              --filter-name FS200 \
  -filter "ReadPosRankSum < -20.0"  --filter-name RPRS20

gatk MergeVcfs -I snps.filt.vcf.gz -I indels.filt.vcf.gz -O filtered.vcf.gz

Each of these thresholds targets a specific artifact. FS and SOR detect strand bias, where a putative variant appears on reads from only one strand. That pattern is characteristic of sequencing error rather than a real heterozygous site. MQ catches variants supported only by poorly mapped reads. ReadPosRankSum catches variants that always appear near read ends, the signature of alignment slippage around indels.

The second option, which we prefer when it is available, is GATK’s CNNScoreVariants followed by FilterVariantTranches. It scores each variant with a convolutional network over the read pileup and reference context and requires no cohort. Deep-learning callers and filters outperform hand-tuned heuristics on short-read data in head-to-head benchmarks. This holds particularly for indels, where DeepVariant has repeatedly achieved higher F1 than GATK on the same reads.³ If you want the best single-sample accuracy and are willing to step outside GATK entirely, run DeepVariant as a second caller. Take the intersection with your GATK calls as your high-confidence set and the symmetric difference as your review queue. Multi-caller consensus is the most reliable defense available against caller-specific systematic error.¹

Finally, add hard-won context. Mask your VCF against a low-complexity and segmental-duplication track, and treat anything inside as unverified:

bcftools annotate -a resources/GRCh38_notinalldifficultregions.bed.gz \
  -c CHROM,FROM,TO -m +EASY_REGION -Oz -o filtered.annot.vcf.gz filtered.vcf.gz

The GIAB stratification BEDs are the practical source for this. Variant discovery in highly divergent and repeat-rich regions requires workflow changes beyond default parameters, and accuracy there is materially worse than the genome-wide average.⁴

6. Annotate and make it queryable

Annotation attaches gene, consequence, and population frequency to each call. We use VEP with the offline cache because it is reproducible and does not depend on a web service being up.

vep --cache --offline --assembly GRCh38 --dir_cache ~/.vep \
  --fasta ref/Homo_sapiens_assembly38.fasta \
  -i filtered.annot.vcf.gz -o annotated.vcf.gz --vcf --compress_output bgzip \
  --everything --pick_allele_gene \
  --custom ~/.vep/gnomad.genomes.v4.1.sites.vcf.gz,gnomAD,vcf,exact,0,AF,AF_popmax \
  --plugin CADD,~/.vep/whole_genome_SNVs.tsv.gz \
  --fork 8

Then flatten to something you can query. A columnar table is far more useful than a VCF for exploration:

bcftools +split-vep annotated.vcf.gz -d -f \
  '%CHROM\t%POS\t%REF\t%ALT\t%FILTER\t%QUAL\t[%GT\t%DP\t%AD]\t%SYMBOL\t%Consequence\t%gnomAD_AF\t%CADD_PHRED\n' \
  | bgzip > variants.tsv.gz

Load that into DuckDB or Parquet and you can answer real questions in seconds: which protein-altering variants in your genome have a gnomAD allele frequency below 0.1%, how many homozygous loss-of-function calls you carry, and how your calls distribute across ClinVar review-status tiers. Integrated annotation pipelines that combine calling with consequence and frequency layers in one pass are becoming the standard packaging for exactly this reason.⁵ If any of those queries surfaces something in a gene with known clinical actionability, the next step is a clinician and a confirmatory assay in a certified lab. Further scripting is not the next step.

Common problems

Pseudogene cross-mapping is the error class that catches careful people. Reads from a processed pseudogene mismap onto the parent gene (or vice versa) and generate a cluster of high-quality-looking heterozygous calls with allele fractions near 0.5. A recent quantitative analysis of germline calling found these pseudogene-associated errors are systematic rather than random and concentrate in specific gene families.⁶ The tell is a tight cluster of variants within a few hundred bases, all matching the paralog’s sequence. Check whether your candidate gene has a known pseudogene before you believe anything in it.

Contig mismatch between reference and BAM produces the error Input files reference and reads have incompatible contigs. This almost always means a UCSC-style chr1 FASTA paired with an Ensembl-style 1 alignment, or the alt contigs present in one and absent from the other. Fix the reference, do not rename contigs in the BAM header.

Java heap exhaustion in GenotypeGVCFs with a GenomicsDB workspace shows up as an opaque OutOfMemoryError hours in. Pass --java-options "-Xmx16g" and --batch-size 50 --reader-threads 4 to GenomicsDBImport, and keep the workspace on local disk rather than a network mount.

Low-coverage regions producing confident-looking homozygous reference calls are a subtle trap. In GVCF mode, check the DP and GQ on reference blocks before concluding you lack a variant. A GQ of 3 over a region means the data were silent. It does not mean that the reference allele was confirmed. Any absence-of-variant conclusion needs a coverage check first.

Coverage itself sets a floor on what is possible. Below about 20x mean coverage, heterozygous sensitivity degrades quickly because one of the two alleles may be sampled only two or three times. At 30x PCR-free you should see mean coverage within 10% of the target and at least 95% of the callable genome at 10x or better. Compute this with gatk CollectWgsMetrics before you interpret anything.

Finally, pin your versions. GATK’s filtering thresholds and annotation names have changed across minor releases, and rerunning the same data with a newer version will shift your call set at the margins. Record the exact version of every tool, the reference checksum, and the resource bundle date alongside your VCF. Configuration-driven frameworks that capture the whole pipeline definition in a single versioned file exist for this purpose and are worth adopting if you plan to rerun as new tools appear.⁷ The broader tooling ecosystem has been fragmented and fast-moving since short-read sequencing became routine, and reproducibility discipline is what keeps your own results comparable to themselves over time.⁸

Footnotes

1. Jason O’Rawe, Tao Jiang, Guangqing Sun, et al. Low concordance of multiple variant-calling pipelines: practical implications for exome and genome sequencing. Genome Medicine, 2013. https://doi.org/10.1186/gm432 ↩ ↩²

2. Peter N. Robinson, Rosario M. Piro, Marten Jäger. Variant Calling and Quality-Based Filtering. Computational Exome and Genome Analysis, 2017. https://doi.org/10.1201/9781315154770-12 ↩ ↩²

3. Omar Abdelwahab, François Belzile, Davoud Torkamaneh. Performance analysis of conventional and AI-based variant callers using short and long reads. BMC Bioinformatics, 2023. https://doi.org/10.1186/s12859-023-05596-3 ↩

4. Shulan Tian, Huihuang Yan, Claudia Neuhauser, et al. An analytical workflow for accurate variant discovery in highly divergent regions. BMC Genomics, 2016. https://doi.org/10.1186/s12864-016-3045-z ↩

5. Karthikasri Karuppusamy, Shobana Sundar. A Comprehensive Integrated Pipeline for Detection and Annotation of Variants in Whole Exome Sequencing Data. Molecular Biotechnology, 2026. https://doi.org/10.1007/s12033-026-01599-6 ↩

6. Artem Podvalnyi, Arina Kopernik, Mariia Sayganova, et al. Quantitative Analysis of Pseudogene-Associated Errors During Germline Variant Calling. International Journal of Molecular Sciences, 2025. https://doi.org/10.3390/ijms26010363 ↩

7. Manuel Rueda, Dietmar Fernandez-Orth, Ivo G Gut. CBIcall: a configuration-driven framework for variant calling in large sequencing cohorts. Bioinformatics Advances, 2026. https://doi.org/10.1093/bioadv/vbag232 ↩

8. S. Pabinger, A. Dander, M. Fischer, et al. A survey of tools for variant analysis of next-generation genome sequencing data. Briefings in Bioinformatics, 2013. https://doi.org/10.1093/bib/bbs086 ↩