DNA-Seq Variant Calling Workflow: Best Practices from FASTQ to Annotated VCF

Q: What is the difference between GATK HaplotypeCaller and DeepVariant?

GATK HaplotypeCaller uses classical probabilistic Bayesian modeling and local de Bruijn graph assembly. Google DeepVariant uses deep convolutional neural networks trained on truth datasets (Genome in a Bottle) to call variants as an image classification task, achieving high accuracy in difficult-to-map regions.

Q: What does the Transition/Transversion (Ti/Tv) ratio indicate in quality control?

Because methylcytosine deamination naturally elevates transition mutations (A G, C T) over transversions, the expected human genomic Ti/Tv ratio is ~2.1 (~3.0 in coding exomes). A significantly lower Ti/Tv indicates random sequencing error noise.

What Is DNA-Seq Variant Calling and Why Is It Foundational in Genomics?

Every individual organism possesses a unique genomic architecture containing millions of inherited or acquired sequence variations relative to species reference genomes. DNA sequencing (DNA-Seq) paired with computational variant calling identifies these variations, which include single nucleotide variants (SNVs), small insertions and deletions (indels), and large copy number and structural variations.

In medical genomics, variant calling identifies pathogenic mutations responsible for rare Mendelian disorders, characterizes pharmacogenomic drug metabolism profiles, and maps somatic driver mutations in cancer genomes. In agricultural and evolutionary biology, variant calling uncovers quantitative trait loci (QTLs) for crop improvement and illuminates population migration history.

How Do You Design DNA-Seq Experiments for Germline and Somatic Discovery?

Experimental design parameters differ fundamentally depending on whether the biological goal is germline variant discovery or somatic cancer mutation profiling:

Germline Variant Calling: Because germline variants exist in heterozygous (50%) or homozygous (100%) states in diploid genomes, Whole-Genome Sequencing (WGS) at 30x coverage or Whole-Exome Sequencing (WES) at 50x to 100x coverage provides ample statistical confidence for accurate genotyping.
Somatic Tumor Variant Calling: Tumors are heterogeneous mixtures of malignant clones, normal stromal cells, and immune infiltrates. Somatic mutations frequently occur at sub-clonal allele frequencies below 5% to 10%. Consequently, somatic profiling requires deep WES (150x to 300x) paired with a matched normal blood sample to subtract inherited germline variants.

Sequencing utilizes paired-end short reads (2x150 bp) to ensure confident alignment across repetitive elements and low-complexity regions.

How Are Reads Aligned, Duplicates Marked, and Base Scores Recalibrated?

Computational processing follows the globally established GATK Best Practices pipeline:

Alignment: Raw FASTQ reads are aligned to reference builds (such as GRCh38/hg38) using BWA-MEM2. Read group headers (@RG) with sample, library, and lane tags must be attached during alignment.
Duplicate Marking: PCR amplification during library construction produces duplicate fragments. Picard or GATK MarkDuplicates tags these duplicates so they are evaluated only once during variant calling, preventing PCR amplification bias from falsifying variant allele frequencies.
Base Quality Score Recalibration (BQSR): Sequencing machines produce systematic error patterns (e.g., quality score drops at specific dinucleotides or cycles). GATK BaseRecalibrator builds a covariate model against known polymorphic sites (dbSNP) and recalibrates base quality scores to reflect true physical error probabilities.

How Does GATK HaplotypeCaller Discover Germline and Somatic Variants?

Traditional variant callers evaluated single alignment columns independently. Modern callers like GATK HaplotypeCaller operate via local de novo assembly:

When HaplotypeCaller detects an active genomic region showing evidence of variation, it discards the existing alignment, builds a local de Bruijn graph from the reads, and assembles the most likely biological haplotypes. It subsequently aligns read sequences against each candidate haplotype using a Pair-HMM algorithm to calculate likelihoods for each possible genotype.

In multi-sample cohorts, HaplotypeCaller is executed in Genomic VCF (-ERC GVCF) mode for each individual. Joint genotyping across all GVCFs via GenotypeGVCFs distinguishes truly homozygous reference sites from sites lacking sequencing coverage, eliminating missing-data bias across cohorts.

How Do You Filter Raw Variants and Annotate Biological Consequences?

Raw VCF files contain false positives arising from sequencing errors, strand bias, and mapping artifacts. Filtering is executed using Variant Quality Score Recalibration (VQSR) in large cohorts, or rigorous hard filtering thresholds (e.g., QD < 2.0, FS > 60.0, MQ < 40.0, ReadPosRankSum < -8.0) for smaller datasets.

High-confidence filtered variants are annotated using Ensembl Variant Effect Predictor (VEP) or SnpEff. Annotation tools predict the transcript-level impact of each variant (synonymous, missense, nonsense, frameshift, splice site disruption) and append in silico pathogenicity scores (SIFT, PolyPhen-2, CADD, REVEL) alongside population allele frequencies from gnomAD, 1000 Genomes, and ClinVar clinical records.

What Are the Core Deliverables and Critical Pitfalls in DNA-Seq?

Standard deliverables include coordinate-sorted, deduplicated BAM files with index files, MultiQC alignment statistics reports, raw and filtered VCF files, comprehensive annotated variant spreadsheets, and transition/transversion (Ti/Tv) quality control ratios.

Critical pitfalls in DNA-Seq analysis include:

Mismatched Reference Builds: Mixing reference chromosome naming conventions (e.g., 'chr1' in UCSC vs '1' in Ensembl) causes alignment failures or silently drops chromosomes.
Calling Variants Without Read Groups: Failing to specify read group identifiers prevents GATK from executing duplicate marking and recalibration.
Abnormal Ti/Tv Ratios: In human whole-genome sequencing, the ratio of transitions to transversions should approximate 2.0 to 2.1 (and ~2.8 to 3.0 in exomes). A low Ti/Tv ratio indicates substantial false-positive error contamination.

Apply These Skills in Your Research

Master this computational workflow step by step with our hands-on mentoring, or collaborate with our data science team for production analysis:

Learn this in our 30-day workshop: DNA-Seq Variant Calling Workshop Need this analysis done for your data? WGS & WES Genomic Analysis Services

Frequently Asked Questions

What is the difference between GATK HaplotypeCaller and DeepVariant?

GATK HaplotypeCaller uses classical probabilistic Bayesian modeling and local de Bruijn graph assembly. Google DeepVariant uses deep convolutional neural networks trained on truth datasets (Genome in a Bottle) to call variants as an image classification task, achieving high accuracy in difficult-to-map regions.

What does the Transition/Transversion (Ti/Tv) ratio indicate in quality control?

Because methylcytosine deamination naturally elevates transition mutations (A<->G, C<->T) over transversions, the expected human genomic Ti/Tv ratio is ~2.1 (~3.0 in coding exomes). A significantly lower Ti/Tv indicates random sequencing error noise.

Why is joint genotyping preferred over single-sample variant calling in cohorts?

Joint genotyping pools sequencing evidence across all samples simultaneously. This empowers the algorithm to detect subtle variants in low-coverage individuals that would be missed in single-sample calling, and accurately separates non-variant homozygous sites from missing data.