Calling Variants From Your Own Whole-Genome Sequencing Data

Oak
variant-calling wgs bioinformatics gatk deepvariant
A lab instrument compares two glowing DNA strands, flagging mismatched points with amber markers that fall into a tray.

By the end of this guide you will have a per-sample gVCF and a filtered VCF containing roughly 4 to 5 million small variants (single-nucleotide variants and short insertions and deletions) called against a human reference, plus a set of benchmark metrics telling you how much of that call set you should trust. You need your raw reads (paired-end FASTQ, or a CRAM/BAM you can revert), about 32 GB of RAM, 16 cores, and 1 TB of scratch disk for a 30x genome. Wall-clock time is 8 to 20 hours depending on the caller and your I/O. You also need a decision about the reference genome before you start, because everything downstream inherits it, and switching later means redoing the alignment.

1. Choose the reference and get the right auxiliary files

The reference choice is the single decision that will cause you the most pain if you get it wrong, because coordinates, annotation databases, and benchmark truth sets are all reference-specific. For a personal genome we use GRCh38 with the analysis set that includes the decoy contigs and excludes the alternate haplotype scaffolds from alignment scoring, specifically GCA_000001405.15_GRCh38_no_alt_analysis_set.fna. The no-alt analysis set matters because BWA will otherwise distribute reads across ALT contigs, drop mapping quality to zero, and silently erase variants in the HLA region and other polymorphic loci. T2T-CHM13 is a better reference in absolute terms and resolves regions GRCh38 cannot represent, but clinical annotation and most population frequency resources are still GRCh38-anchored, so we treat CHM13 as a second alignment you do later for specific questions rather than your primary.

# reference + index
samtools faidx GCA_000001405.15_GRCh38_no_alt_analysis_set.fna
bwa-mem2 index GCA_000001405.15_GRCh38_no_alt_analysis_set.fna
gatk CreateSequenceDictionary -R GCA_000001405.15_GRCh38_no_alt_analysis_set.fna

You will also want the GATK resource bundle files for base quality score recalibration (dbSNP 138, Mills and 1000G gold-standard indels, the 1000 Genomes phase 1 high-confidence indels) and, for benchmarking in step 7, the Genome in a Bottle truth sets. Download these once and keep them alongside the reference so that every run points at identical inputs. Reproducibility failures in personal pipelines are usually version drift in these auxiliary files rather than in the code.

2. Inspect and trim the reads only if you have to

Run fastqc or, faster on a whole genome, seqkit stats and fastp --disable_adapter_trimming --json to look at per-base quality, adapter content, duplication estimate, and GC distribution. For modern Illumina data on a well-prepared library, aggressive trimming is a mistake: the aligner soft-clips adapters and low-quality tails on its own, and trimming shortens reads, perturbs the insert-size distribution, and can bias indel representation. We trim only when adapter content exceeds a few percent, which happens with short inserts.

fastp -i R1.fq.gz -I R2.fq.gz -o R1.trim.fq.gz -O R2.trim.fq.gz \
  --detect_adapter_for_pe --length_required 50 \
  --json fastp.json --html fastp.html --thread 8

Read the fastp.json for duplication rate and insert-size peak. A duplication rate above roughly 20% at 30x usually means the library was over-amplified from low input, which caps the effective coverage you have and will show up later as unexpectedly low genotype quality in heterozygous calls.

3. Align with BWA-MEM2 and set read groups correctly

Alignment is where read-group metadata becomes non-negotiable, because GATK’s recalibration and duplicate marking both key on it. Give each sequencing run its own ID and PU, and keep SM identical across all runs from the same person.

bwa-mem2 mem -t 16 -K 100000000 -Y \
  -R '@RG\tID:L001\tSM:SUBJ01\tLB:lib1\tPL:ILLUMINA\tPU:HXXXX.1.ATCGAT' \
  GCA_000001405.15_GRCh38_no_alt_analysis_set.fna R1.fq.gz R2.fq.gz \
| samtools fixmate -m -u - - \
| samtools sort -@ 4 -m 2G -T /scratch/sort -o SUBJ01.sorted.bam -
samtools index -@ 8 SUBJ01.sorted.bam

The -K 100000000 flag fixes the batch size so that output is deterministic regardless of thread count, which you want if you ever intend to compare two runs. -Y uses soft clipping for supplementary alignments, which keeps the full read sequence available downstream. Piping straight into fixmate and sort avoids writing an intermediate SAM, which on a 30x genome saves several hundred gigabytes of writes.

Mark duplicates next. We use samtools markdup when speed matters and GATK MarkDuplicates when we want the metrics file and exact GATK-canonical behavior; on a single personal genome the difference in downstream call sets is small.

samtools markdup -@ 8 -s -f markdup.stats SUBJ01.sorted.bam SUBJ01.md.bam
samtools index -@ 8 SUBJ01.md.bam

4. Decide between a haplotype-based caller and a deep learning caller

This is the second consequential decision. The GATK Best Practices pipeline, which reassembles reads locally into candidate haplotypes and then genotypes them, is the reference implementation most of the field is built around and remains the approach we would use when we need gVCFs for joint calling or want to reason about every filtering knob 1. DeepVariant instead encodes pileups as images and classifies genotypes with a convolutional network, and it generally produces fewer false positives on Illumina data, particularly in indels and homopolymer regions, with far fewer hand-tuned parameters. Benchmarking work comparing callers consistently finds that the choice of caller interacts with coverage, variant type, and genomic region rather than producing one uniform winner 2. Our default for a single 30x Illumina personal genome is DeepVariant for the primary call set, with GATK HaplotypeCaller run in gVCF mode as a second call set for comparison.

The GATK path, including base quality score recalibration, looks like this:

gatk BaseRecalibrator -R $REF -I SUBJ01.md.bam \
  --known-sites Homo_sapiens_assembly38.dbsnp138.vcf \
  --known-sites Mills_and_1000G_gold_standard.indels.hg38.vcf.gz \
  --known-sites Homo_sapiens_assembly38.known_indels.vcf.gz \
  -O recal.table

gatk ApplyBQSR -R $REF -I SUBJ01.md.bam --bqsr-recal-file recal.table \
  -O SUBJ01.bqsr.bam

gatk --java-options "-Xmx8g" HaplotypeCaller \
  -R $REF -I SUBJ01.bqsr.bam -O SUBJ01.g.vcf.gz \
  -ERC GVCF -G StandardAnnotation -G AS_StandardAnnotation \
  --native-pair-hmm-threads 8

Run HaplotypeCaller per chromosome with -L chr1 and so on, in parallel, then merge the gVCFs with GatherVcfs. A single-threaded whole-genome run takes on the order of a day. Convert the gVCF to a genotyped VCF with GenotypeGVCFs.

The DeepVariant path is a single container invocation and needs no recalibration, since the model learns quality behavior from the pileup directly:

singularity run -B /data:/data docker://google/deepvariant:1.6.1 \
  /opt/deepvariant/bin/run_deepvariant \
  --model_type=WGS \
  --ref=/data/$REF \
  --reads=/data/SUBJ01.md.bam \
  --output_vcf=/data/SUBJ01.dv.vcf.gz \
  --output_gvcf=/data/SUBJ01.dv.g.vcf.gz \
  --num_shards=16 \
  --intermediate_results_dir=/scratch/dv_tmp

Set --model_type to match your data: WGS, WES, PACBIO, or ONT_R104. Running the WGS model on exome data inflates false positives at capture boundaries, and running it on long reads produces call sets that are close to useless.

5. Filter the call set

Raw VCFs contain a large tail of low-confidence calls, and how you filter determines whether your downstream analysis is interpretable. For GATK output on a single sample, Variant Quality Score Recalibration is poorly conditioned because it needs many samples to model the annotation distribution, so we use hard filters from the Best Practices recommendations instead, which were designed for exactly this case 1.

gatk SelectVariants -V SUBJ01.hc.vcf.gz -select-type SNP -O snps.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 SelectVariants -V SUBJ01.hc.vcf.gz -select-type INDEL -O indels.vcf.gz
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

DeepVariant emits its own PASS/RefCall decisions and a calibrated GQ, so post-filtering is lighter. We apply a genotype-level depth and quality floor and drop anything in problem regions:

bcftools view -f PASS -i 'FMT/DP>=10 && FMT/GQ>=20' SUBJ01.dv.vcf.gz \
| bcftools view -T ^GRCh38_notinalldifficultregions.bed.gz -Oz -o SUBJ01.pass.vcf.gz
bcftools index SUBJ01.pass.vcf.gz

Sanity-check the filtered set with bcftools stats. On a 30x European-ancestry genome you should see roughly 4 to 5 million variants, a heterozygous-to-homozygous ratio near 1.5 to 2.0, a transition-to-transversion ratio around 2.0 to 2.1 genome-wide, and 15 to 25 thousand novel variants not in dbSNP. Deviations point at specific problems: a low Ti/Tv suggests false positives, an extreme het/hom ratio suggests contamination or reference mismatch, and a very high novel fraction suggests you called against the wrong reference build. Published quality-filtering guidelines for population-scale data formalize these checks and add per-sample metrics such as contamination, coverage uniformity, and sex-chromosome ploidy concordance 3.

6. Normalize and annotate

Before any annotation, normalize representation. The same indel can be written several ways, and mismatched representation is the most common reason a variant appears absent from a database that contains it.

bcftools norm -f $REF -m -both -c w -Oz -o SUBJ01.norm.vcf.gz SUBJ01.pass.vcf.gz
bcftools index SUBJ01.norm.vcf.gz

Then annotate with Ensembl VEP, picking a single consequence per variant only if you understand what you are discarding.

vep -i SUBJ01.norm.vcf.gz -o SUBJ01.vep.vcf.gz --vcf --compress_output bgzip \
  --cache --offline --assembly GRCh38 --fasta $REF \
  --everything --pick_allele_gene --fork 8 \
  --custom gnomad.genomes.v4.1.sites.vcf.gz,gnomAD,vcf,exact,0,AF,AF_popmax

The surrounding tooling ecosystem is large and changes quickly, and surveys of variant analysis tools are useful for orienting yourself when you need something specific like structural variant callers or copy-number tools that this pipeline does not cover 4. If your interest extends to somatic variants from tumor tissue, the analytical problem is different enough (subclonal allele fractions, matched normal requirements, tumor purity) that you should follow a dedicated somatic workflow rather than adapt this one 5.

7. Benchmark against a truth set so you know your error rate

Skipping this step is the most common mistake we see. Running the same pipeline you just built on a Genome in a Bottle sample, such as HG002, and comparing to the GIAB truth VCF inside the high-confidence BED tells you your own precision and recall rather than a number from a paper. Benchmarking practice has matured substantially, and current truth sets cover a much larger fraction of the genome than earlier versions, including many segmental duplications and other previously excluded regions 6.

hap.py HG002_GRCh38_1_22_v4.2.1_benchmark.vcf.gz SUBJ01_HG002.pass.vcf.gz \
  -f HG002_GRCh38_1_22_v4.2.1_benchmark_noinconsistent.bed \
  -r $REF -o happy_out --engine vcfeval --threads 16 \
  --stratification GRCh38_stratification.tsv

Use --engine vcfeval so that equivalent representations of the same indel are matched rather than counted as a false positive and a false negative. On 30x Illumina you should expect SNV F1 above 0.995 inside the high-confidence regions and indel F1 in the 0.99 range with DeepVariant, somewhat lower with hard-filtered GATK. The stratified output matters more than the headline number: performance in homopolymers, tandem repeats, and segmental duplications is dramatically worse than the genome-wide average, and knowing which stratum a given variant of interest falls into is what lets you decide how much weight to put on it 7.

Common problems

Low mapping quality across an entire gene, most often in HLA, immunoglobulin loci, or the amylase cluster, usually means you aligned to a reference that includes ALT contigs without ALT-aware post-processing. Realign to the no-alt analysis set. A related symptom is a sudden block of missing calls in an otherwise well-covered region, which you should check in IGV before believing it.

Sample swaps and contamination are more common than people expect in personal sequencing. Run VerifyBamID2 or GATK CalculateContamination and check that estimated contamination is below about 1%, and confirm the genotype-inferred sex matches expectation by comparing chrX and chrY coverage. If you have any prior genotyping data, such as an array file, run bcftools gtcheck against it. A pipeline that produces a beautiful VCF from the wrong person’s DNA is a real failure mode.

Genotype quality that collapses in the last few megabases of each chromosome, or a Ti/Tv ratio drifting below 1.8, generally indicates coverage or base-quality problems upstream rather than caller error. Foundational treatments of SNP calling walk through how read depth, base quality, and mapping quality jointly determine confidence, and re-reading that logic is usually faster than tuning filters blindly 8.

Finally, a filtered, benchmarked VCF is a measurement, not an interpretation. Allele frequency, computational pathogenicity predictors, and even ClinVar assertions do not establish that a variant is affecting your health, and single-sample germline pipelines of the kind described here are not validated for clinical use. Clinical sequencing has its own requirements around confirmatory orthogonal testing, region-level coverage guarantees, and reporting standards that a research pipeline does not meet 9. If anything in your call set looks medically significant, take it to a clinical geneticist or genetic counselor for confirmatory testing and interpretation rather than acting on it yourself.

Oak builds longitudinal molecular profiles of individuals: whole-genome sequencing, RNA sequencing, proteomics, blood biomarkers, and continuous glucose data, integrated into one model of you. Build your profile.

Footnotes

  1. Geraldine A. Van der Auwera, Mauricio O. Carneiro, Christopher Hartl, et al. From FastQ Data to High‐Confidence Variant Calls: The Genome Analysis Toolkit Best Practices Pipeline. Current Protocols in Bioinformatics, 2013. https://doi.org/10.1002/0471250953.bi1110s43 ↩ ↩2

  2. Vera Pinto, Lisete Sousa, Carina Silva. Variant calling in genomics: A comparative performance analysis and decision guide. PLOS One, 2026. https://doi.org/10.1371/journal.pone.0339891 ↩

  3. Julia M. Sealock, Franjo Ivankovic, Calwing Liao, et al. Tutorial: guidelines for quality filtering of whole-exome and whole-genome sequencing data for population-scale association analyses. Nature Protocols, 2025. https://doi.org/10.1038/s41596-025-01169-1 ↩

  4. 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 ↩

  5. Riley J. Arseneau, Leah K. MacLean, Jeanette E. Boudreau, et al. Tutorial for variant interrogation in tumor samples. PLOS Computational Biology, 2026. https://doi.org/10.1371/journal.pcbi.1013924 ↩

  6. Nathan D. Olson, Justin Wagner, Nathan Dwarshuis, et al. Variant calling and benchmarking in an era of complete human genome sequences. Nature Reviews Genetics, 2023. https://doi.org/10.1038/s41576-023-00590-0 ↩

  7. Stepanka Zverinova, Victor Guryev. Variant calling: Considerations, practices, and developments. Human Mutation, 2021. https://doi.org/10.1002/humu.24311 ↩

  8. André Altmann, Peter Weber, Daniel Bader, et al. A beginners guide to SNP calling from high-throughput DNA-sequencing data. Human Genetics, 2012. https://doi.org/10.1007/s00439-012-1213-z ↩

  9. Daniel C. Koboldt. Best practices for variant calling in clinical sequencing. Genome Medicine, 2020. https://doi.org/10.1186/s13073-020-00791-w ↩