What is Dr. Zubair Khalid's research focus?

Dr. Zubair Khalid specializes in molecular virology, mRNA vaccine development, and computational biology, with a focus on avian pathogens like IBDV and Avian Reovirus.

Where is Dr. Zubair Khalid currently working?

Dr. Zubair Khalid is a Postdoctoral Research Associate at the University of Maryland (UMD), specifically within the Department of Animal and Avian Sciences.

Variant Calling in Genomics: Methods, Pitfalls, and Best Practices

By Dr. Zubair Khalid, DVM, MS, PhD · 2026-08-11

Introduction to Variant Calling in Genomics

What is Variant Calling?

Variant calling is the computational process of identifying positions where an individual's genome differs from a reference genome. These differences—collectively termed genetic variants—are the raw material of genetic diversity, disease susceptibility, and evolutionary adaptation. In a typical human whole-genome sequencing experiment, variant calling identifies approximately 3–4 million single-nucleotide variants (SNVs), several hundred thousand insertions/deletions (indels), and thousands of larger structural variants relative to the GRCh38 reference.

The fundamental challenge of variant calling is distinguishing genuine biological variation from sequencing artifacts. A typical Illumina sequencing run produces reads with an error rate of approximately 0.1–1% per base, meaning that in a 30× genome, every position is covered by roughly 30 reads, and each read has a small but non-negligible chance of carrying an error. Variant callers use statistical models to determine whether the observed pattern of bases at a position is better explained by a true variant or by sequencing error.

Variant calling operates at multiple scales. At the nucleotide level, it detects SNVs and small indels (typically <50 bp). At the chromosomal level, it detects structural variants (SVs)—deletions, duplications, inversions, and translocations larger than 50 bp. While SNV calling is well-established and highly accurate, SV calling remains more challenging and often requires specialized algorithms that leverage read-pair orientation, split-read alignment, and read-depth signals.

Types of Genetic Variants

Single-nucleotide variants (SNVs) are substitutions of one base for another. They are the most common type of variant and the easiest to detect reliably. SNVs can be further classified by their functional impact: synonymous (no amino acid change), missense (amino acid change), or nonsense (premature stop codon). A well-known example is the sickle cell mutation in HBB (c.20A>T, p.Glu6Val), a single base substitution that alters the β-globin protein's hydrophobicity and causes hemoglobin polymerization.

Insertions and deletions (indels) involve the addition or removal of 1–50 nucleotides. Indels that are not multiples of three in coding regions cause frameshift mutations, which typically produce truncated, non-functional proteins. The CFTR gene provides a classic example: the deletion of phenylalanine at position 508 (c.1521_1523delCTT, p.Phe508del) is the most common cause of cystic fibrosis. Indel calling is more challenging than SNV calling because alignment algorithms must account for gaps, and repetitive regions often produce ambiguous alignments.

Structural variants (SVs) encompass deletions, duplications, inversions, and translocations larger than 50 bp. These variants can have profound phenotypic consequences—for example, the 22q11.2 deletion syndrome results from a ~3 Mb deletion that causes cardiac defects and palatal abnormalities. SV detection requires different computational strategies than SNV/indel calling, including analysis of discordant read pairs, split-read mapping, and read-depth ratios.

The Variant Calling Workflow: From Reads to Variants

Quality Control and Preprocessing

The variant calling pipeline begins before any variant is called. Raw sequencing data in FASTQ format must first pass quality control to identify problems that could compromise downstream analysis. The standard tool for this step is FastQC, which reports per-base quality scores, GC content, adapter contamination, and overrepresented sequences.

Quality scores, reported as Phred scores (Q), are logarithmically related to error probability: Q30 corresponds to a 1 in 1000 error rate, Q20 to 1 in 100. For variant calling, reads with mean quality below Q20 are typically flagged or removed. Adapter trimming is often necessary, particularly for reads longer than 150 bp where adapter sequences may be read through. Tools like Trimmomatic or cutadapt remove adapter sequences and low-quality bases from read ends.

A critical preprocessing step is read error correction or trimming of the 3′ ends of reads, where quality scores degrade. For Illumina data, a common practice is to trim bases with quality below Q20 from the 3′ end using a sliding window approach (e.g., a 4-base window with average quality below Q15 triggers trimming).

Read Alignment and BAM Files

After quality control, reads are aligned to a reference genome using an aligner such as BWA-MEM, Bowtie2, or minimap2. The choice of aligner matters: BWA-MEM is the standard for Illumina short reads, while minimap2 is preferred for long reads from PacBio or Oxford Nanopore platforms. The aligner produces a SAM (Sequence Alignment/Map) file, which is typically converted to its binary equivalent, BAM, using SAMtools.

Alignment is not a trivial mapping of reads to a linear reference. The aligner must handle:

Sequencing errors that create mismatches
Indels that require gap placement
Reads spanning exon junctions in RNA-seq data
Repetitive regions where reads map to multiple locations with equal confidence

For reads that map to multiple locations (multi-mappers), most pipelines discard them or assign them randomly. This is particularly important in regions of segmental duplication, where multi-mapping reads can create spurious variant calls.

The BAM file contains, for each read, its reference position, mapping quality (MAPQ), CIGAR string (which describes the alignment operations), and base qualities. The MAPQ score is critical: reads with MAPQ below 20 are typically excluded from variant calling because their alignment is uncertain.

Post-Alignment Processing (Mark Duplicates, BQSR)

Before variant calling, the BAM file requires additional processing to correct systematic artifacts.

Marking duplicates identifies reads that originate from the same DNA fragment during library preparation. During PCR amplification, a single fragment can produce many identical reads; if these are counted as independent observations, they artificially inflate confidence in the bases they carry. The Picard MarkDuplicates tool identifies read pairs with identical 5′ alignment positions and marks all but one as duplicates. For whole-genome sequencing, duplicate rates of 5–15% are typical; for amplicon-based targeted panels, rates can exceed 50%.

Base Quality Score Recalibration (BQSR) is a machine-learning step that adjusts base quality scores based on empirically observed error patterns. The Genome Analysis Toolkit (GATK) implements BQSR by building a model of errors as a function of:

Reported base quality
Position within the read (cycles 1–5 and the last few cycles have higher error rates)
Dinucleotide context (e.g., CpG sites have higher error rates)
Sequencing machine and flow cell

BQSR requires a set of known variant sites (e.g., from dbSNP) to exclude true variants from the error model. The recalibrated BAM file has adjusted quality scores that better reflect true error probabilities, which improves variant calling accuracy.

Core Algorithms for Variant Detection

Bayesian Genotyping

Most variant callers use a Bayesian framework to compute the probability of each possible genotype given the observed reads. The core equation is:

P(Genotype | Reads) ∝ P(Reads | Genotype) × P(Genotype)

where P(Genotype) is the prior probability (often based on population allele frequencies) and P(Reads | Genotype) is the likelihood of observing the read data given a particular genotype.

For a diploid organism, the possible genotypes at a position are homozygous reference (RR), heterozygous (RA), and homozygous alternate (AA). The likelihood calculation considers:

The base quality of each read at the position
The mapping quality of each read
The probability of sequencing error
The probability of a read being misaligned

The GATK HaplotypeCaller implements this Bayesian framework but adds a critical innovation: it performs local re-assembly of the region around each putative variant before genotyping. This haplotype-based approach improves indel calling accuracy because it can resolve variants that fall in repetitive or low-complexity regions where read alignment alone is ambiguous.

Haplotype-Based Calling

Haplotype-based callers, including GATK HaplotypeCaller and FreeBayes, work in several steps:

Identify active regions: Scan the genome for regions where the read data deviate from the reference in ways that suggest variation.
Local assembly: For each active region, assemble the reads into haplotypes using a De Bruijn graph. This step determines the actual sequence of each chromosome copy in the region.
Read realignment: Align the original reads back to the assembled haplotypes, determining which haplotype each read supports.
Genotype likelihood calculation: For each variant site, compute the likelihood of each genotype given the read support for each haplotype.
Genotype assignment: Use Bayesian inference to assign the most probable genotype and compute a quality score (Phred-scaled probability that the genotype is wrong).

The key advantage of haplotype-based calling is that it treats indels and nearby SNVs as part of a single haplotype rather than independent events. This avoids the common error of double-counting a single variant as both an SNV and an indel in repetitive regions.

Germline vs. Somatic Calling

Germline variant calling identifies variants present in the germline genome—variants inherited from parents and present in all cells. The standard approach is to call variants from a single sample (or a family trio) using tools like GATK HaplotypeCaller or FreeBayes. Germline variants are expected to be present at allele fractions of 50% (heterozygous) or 100% (homozygous) in diploid samples.

Somatic variant calling identifies variants that arise during an individual's lifetime, typically in cancer. Somatic variants are present only in the tumor tissue and are often subclonal, meaning they are present in only a fraction of tumor cells. This creates two challenges:

Low allele fractions: A somatic variant present in 20% of tumor cells in a sample with 80% tumor purity will have an allele fraction of only 10%—well below the 50% expected for germline heterozygotes.
Contamination by normal cells: Tumor samples are always mixed with normal cells, diluting the variant signal.

Somatic callers such as Mutect2, Strelka2, and VarScan2 use a matched normal sample to distinguish somatic variants from germline variants and sequencing artifacts. The matched normal provides a baseline for the patient's germline genome; any variant present in the tumor but absent from the normal is a candidate somatic variant. Mutect2 uses a panel of normals (PoN) to model and filter common sequencing artifacts and germline variants that are missed in the matched normal.

Variant Filtering and Quality Control

Hard Filters vs. VQSR

Raw variant calls contain a substantial fraction of false positives. The GATK best practices recommend one of two filtering strategies:

Hard filters apply fixed thresholds to variant annotations. For example, a common hard filter set for SNVs requires:

Quality by depth (QD) > 2.0
Fisher strand bias (FS) < 60
Mapping quality (MQ) > 40
Strand odds ratio (SOR) < 3.0
Read position rank sum (ReadPosRankSum) > -2.5

Hard filters are simple, transparent, and reproducible, but they are not optimized for each dataset. A variant with QD of 1.9 might be a true positive in a low-coverage region, while a variant with QD of 5 might be a false positive in a repetitive region.

Variant Quality Score Recalibration (VQSR) is a machine-learning approach that builds a Gaussian mixture model from the distribution of variant annotations in the call set. It uses a set of known true variants (from HapMap, Omni, and 1000 Genomes) as training data and a set of known false positives (from the dbSNP "mismatch" set) as negative training data. VQSR produces a continuous quality score (VQSLOD) that can be thresholded to achieve a desired sensitivity/specificity trade-off.

VQSR requires a large call set (typically >30,000 variants) to build a reliable model, making it unsuitable for targeted panels or small exome datasets. For those applications, hard filters are the standard approach.

Common Filtering Criteria

The following annotations are commonly used for filtering:

Annotation	Description	Typical Threshold (SNVs)
QD	Quality by depth: variant quality divided by read depth	> 2.0
FS	Fisher strand bias: probability of strand bias	< 60
SOR	Strand odds ratio: strand bias measure	< 3.0
MQ	Mapping quality: root mean square of read mapping qualities	> 40.0
MQRankSum	Rank sum test for mapping quality of ref vs. alt reads	> -12.5
ReadPosRankSum	Rank sum test for read position of ref vs. alt reads	> -8.0
DP	Read depth at the variant site	> 10 (minimum)
GQ	Genotype quality: Phred-scaled confidence in genotype	> 20

Strand bias is a particularly important filter. If the alternate allele is observed predominantly on one strand, it is likely a sequencing artifact rather than a true variant. The FS and SOR annotations capture this signal.

Variant Annotation and Interpretation

Functional Effect Prediction

Once variants are called and filtered, they must be annotated with their functional consequences. Annotation tools such as ANNOVAR, SnpEff, and the Ensembl Variant Effect Predictor (VEP) map variants to genes and transcripts, then predict their effects.

For coding variants, the annotation includes:

Gene name and transcript ID (e.g., BRCA1, ENST00000357654)
Variant type (missense, nonsense, synonymous, frameshift, splice site)
Amino acid change (e.g., p.Val168Met)
Protein domain affected (e.g., kinase domain)
Conservation scores (e.g., PhyloP, GERP++)

For non-coding variants, annotation is more complex. Tools like VEP can annotate variants in promoters, enhancers, and other regulatory elements using databases such as ENCODE and Roadmap Epigenomics. The functional impact of non-coding variants is harder to predict, but tools like CADD (Combined Annotation-Dependent Depletion) integrate multiple annotations into a single deleteriousness score.

Population Databases (gnomAD, dbSNP)

Population frequency databases are essential for interpreting variant significance. The Genome Aggregation Database (gnomAD) contains variant frequencies from over 140,000 individuals across multiple populations. A variant present at high frequency in gnomAD (e.g., >1%) is unlikely to be a rare disease-causing variant, whereas a variant absent from gnomAD is more likely to be pathogenic.

dbSNP is a broader database that includes variants from many sources, including those of unknown clinical significance. It is important to note that dbSNP contains many common, benign variants; its presence in dbSNP does not imply pathogenicity.

For clinical interpretation, the American College of Medical Genetics and Genomics (ACMG) guidelines provide a framework for classifying variants as pathogenic, likely pathogenic, uncertain significance, likely benign, or benign. This classification integrates population frequency, functional prediction, segregation data, and functional assays.

Special Considerations for Different Sequencing Data

WGS vs. WES vs. Targeted Panels

Whole-genome sequencing (WGS) provides uniform coverage across the genome, including non-coding regions. The challenges for WGS variant calling include:

Computational cost: A 30× human genome produces ~90 GB of FASTQ data and requires substantial compute for alignment and calling.
Repetitive regions: Approximately 50% of the human genome is repetitive, and variant calling in these regions is unreliable.
Coverage uniformity: Despite improvements in library preparation, GC-rich and GC-poor regions are often undercovered. The sequencing coverage in these regions can drop below 10×, making variant calling unreliable.

Whole-exome sequencing (WES) targets the ~1–2% of the genome that is protein-coding. The advantages are lower cost and higher coverage in coding regions (typically 100–200×). The challenges include:

Capture bias: Exome capture probes have variable efficiency, leading to uneven coverage.
Off-target reads: A substantial fraction of reads (20–40%) map outside the targeted regions and are typically discarded.
Incomplete capture: Some exons are poorly captured or missed entirely, creating false-negative results.

Targeted panels focus on specific genes of interest (e.g., cancer gene panels). They offer the highest coverage (500–1000×) and lowest cost, but they cannot detect variants outside the targeted regions. The high coverage enables detection of low-allele-fraction somatic variants but also increases the chance of detecting PCR artifacts, making duplicate marking and error modeling critical.

RNA-seq Variant Calling

Variant calling from RNA-seq data is fundamentally different from DNA-based calling. The process of RNA-seq data introduces several challenges:

Splicing: Reads span exon junctions and cannot be aligned to the linear reference genome without a splice-aware aligner (STAR or HISAT2). Variant callers must account for the possibility that a read spanning a junction may not align contiguously.
Allele-specific expression: The expression level of each allele can differ, so the allele fraction at a variant site reflects both the genotype and the relative expression of the two alleles. A heterozygous variant may show 80:20 allele balance due to allele-specific expression, not because of a somatic mutation.
RNA editing: Post-transcriptional modifications, particularly A-to-I editing by ADAR enzymes, create apparent variants that are not present in the DNA. These are biological artifacts that must be distinguished from true variants.
Coverage bias: Genes with low expression have low coverage, making variant calling unreliable. Only variants in moderately to highly expressed genes can be called with confidence.

For RNA-seq variant calling, tools like GATK HaplotypeCaller in RNA-seq mode use a splice-aware approach and apply additional filtering to remove false positives from misaligned junction reads.

Common Pitfalls and How to Avoid Them

Reference Genome Artifacts

The choice of reference genome is critical. Using an outdated or incorrect reference can introduce systematic errors. For example, the GRCh37 (hg19) reference contains several known errors, including an inversion on chromosome 8 that is corrected in GRCh38. If you use GRCh37, variants in this region will be called incorrectly.

Reference genome artifacts also arise from:

Alternative contigs: GRCh38 includes alt contigs that represent alternate haplotypes of complex regions. Reads mapping to these contigs can confuse variant callers.
Decoy sequences: The hs37d5 decoy sequence (used with GRCh37) captures reads from common structural variants that would otherwise map incorrectly.
Blacklisted regions: The ENCODE project has defined blacklisted regions—areas with anomalous, unstructured read mapping that produce artifactual variant calls. These regions should be excluded from analysis.

Batch Effects and Confounders

Batch effects are systematic technical variations that affect samples processed in different batches. In variant calling, batch effects can arise from:

Different sequencing platforms (Illumina vs. Ion Torrent)
Different library preparation kits (e.g., different exome capture kits)
Different sequencing centers or operators
Different reagent lots or flow cells

These effects can create spurious associations in case-control studies. For example, if all cases are sequenced on one machine and all controls on another, variants that differ between the two machines will appear to be associated with the disease. The difference between genomics and epigenetics is relevant here: batch effects can mimic epigenetic differences if not properly controlled.

To mitigate batch effects:

Randomize samples across batches whenever possible.
Include control samples in every batch.
Use principal component analysis to identify and correct for batch effects.
Run joint variant calling across all samples to ensure consistent genotyping.

Validation with Sanger Sequencing

Sanger sequencing remains the gold standard for validating variant calls, particularly for clinically significant variants. The Sanger sequencing protocol produces high-quality, single-molecule sequencing that can confirm or refute a variant call.

However, Sanger validation has limitations:

It is expensive and low-throughput, making it impractical for validating thousands of variants.
It cannot validate variants in regions with complex secondary structure or high GC content.
It may fail for variants with allele fractions below ~15–20%, as the Sanger trace may not show the minor allele.

For research applications, orthogonal validation using a different sequencing technology (e.g., PacBio or Oxford Nanopore) is increasingly used instead of Sanger sequencing. For clinical applications, Sanger validation remains standard practice for medically actionable variants.

Practical Summary: Best Practices for Reliable Variant Calling

Key Steps Checklist

Start with high-quality input data: Assess FASTQ quality with FastQC; trim adapters and low-quality bases; verify sample identity and sex match.
Use the correct reference genome: Download the latest reference (GRCh38 for human) and its associated index files; use the same reference for all samples in a study.
Align with a modern aligner: Use BWA-MEM for short reads, minimap2 for long reads; set appropriate parameters for read length and error rate.
Process BAM files consistently: Mark duplicates, recalibrate base qualities (BQSR), and use the same pipeline for all samples.
Call variants with a haplotype-based caller: Use GATK HaplotypeCaller for germline, Mutect2 for somatic; run joint calling across samples for cohort studies.
Filter variants appropriately: Use VQSR for large datasets (>30,000 variants), hard filters for small datasets; always examine the distribution of annotations.
Annotate variants thoroughly: Use VEP or ANNOVAR with multiple transcript databases; include population frequencies from gnomAD.
Validate clinically significant variants: Confirm with Sanger sequencing or an orthogonal technology before reporting.

Resources and Tools

Alignment: BWA-MEM, Bowtie2, minimap2, STAR (RNA-seq)
BAM processing: SAMtools, Picard, GATK
Variant calling: GATK HaplotypeCaller, FreeBayes, Mutect2, Strelka2, VarScan2
Filtering: GATK VariantFiltration, VQSR, bcftools
Annotation: ANNOVAR, SnpEff, VEP, CADD
Visualization: IGV, Tablet

Frequently Asked Questions

What is variant calling in genomics?

Variant calling is the computational process of identifying positions in a genome where an individual's DNA sequence differs from a reference genome. It detects single-nucleotide variants (SNVs), small insertions and deletions (indels), and larger structural variants, and it assigns a genotype (homozygous reference, heterozygous, or homozygous alternate) at each variant site with an associated confidence score.

What are the main steps in variant calling?

The main steps are: (1) quality control of raw sequencing reads, (2) alignment of reads to a reference genome, (3) post-alignment processing including duplicate marking and base quality recalibration, (4) variant detection using a caller such as GATK HaplotypeCaller, (5) variant filtering to remove false positives, and (6) annotation of variants with functional and population frequency information.

What tools are commonly used for variant calling?

Common germline callers include GATK HaplotypeCaller, FreeBayes, and Strelka2. Common somatic callers include Mutect2, VarScan2, and Strelka2. Alignment is typically performed with BWA-MEM for short reads or minimap2 for long reads. Annotation is performed with ANNOVAR, SnpEff, or VEP.

What is the difference between germline and somatic variant calling?

Germline variant calling identifies variants inherited from parents and present in all cells of the body; these variants are expected at allele fractions of 50% (heterozygous) or 100% (homozygous). Somatic variant calling identifies variants that arise during an individual's lifetime, typically in cancer; these variants are present only in the tumor and may have allele fractions as low as 1–10%. Somatic calling requires a matched normal sample to distinguish somatic variants from germline variants.

How do I reduce false positives in variant calling?

Reduce false positives by: (1) using high-quality input reads with proper quality trimming, (2) marking PCR duplicates to avoid overcounting, (3) applying strand bias filters to remove orientation-specific artifacts, (4) using VQSR or well-calibrated hard filters, (5) excluding known problematic regions (e.g., ENCODE blacklisted regions), and (6) validating a subset of calls with an orthogonal method.

What is VQSR and when should I use it?

Variant Quality Score Recalibration (VQSR) is a machine-learning method that builds a model of variant quality from the distribution of annotations (QD, FS, MQ, etc.) in your call set, using known true and false variant sets as training data. It produces a continuous quality score that can be thresholded for optimal sensitivity/specificity. Use VQSR when you have a large call set (typically >30,000 variants, as in whole-genome or whole-exome data). For small targeted panels, use hard filters instead.

Why is variant calling more challenging in RNA-seq data?

RNA-seq variant calling is challenging because: (1) reads span exon junctions and require splice-aware alignment, (2) allele-specific expression can distort allele fractions, (3) RNA editing creates apparent variants not present in DNA, (4) coverage is highly variable and depends on gene expression levels, and (5) mapping bias can arise from sequence differences between the reference and the individual's genome.

Key Takeaways

Variant calling is a multi-step process requiring careful quality control, alignment, post-alignment processing, and statistical genotyping; shortcuts at any step compromise the final result.
Haplotype-based callers (GATK HaplotypeCaller, FreeBayes) outperform simple pileup-based callers for indels and variants in complex regions because they perform local assembly before genotyping.
Somatic variant calling is fundamentally different from germline calling: it requires a matched normal sample and specialized algorithms that can detect low-allele-fraction variants in the presence of normal cell contamination.
Filtering is essential: raw variant calls contain many false positives, and both hard filters and VQSR have appropriate use cases depending on dataset size.
The choice of reference genome, aligner, and variant caller introduces systematic biases; use the same pipeline consistently across all samples in a study to minimize batch effects.
RNA-seq variant calling requires splice-aware alignment and additional filtering to account for allele-specific expression and RNA editing; it cannot reliably detect variants in lowly expressed genes.
Always validate clinically significant variants with an orthogonal method such as Sanger sequencing before reporting results to patients or clinicians.