🧬 Whole Genome Sequencing Variant Calling Pipeline

From Photons to Letters — Base-Calling and Raw Read Quality Control

Every whole genome sequencing pipeline begins with the sequencer's raw signal — fluorescence intensities captured over hundreds of imaging cycles on an Illumina NovaSeq X flow cell — being base-called into short nucleotide reads with per-base confidence scores. Before a single read touches the reference genome, quality control tools quantify sequencing error, adapter contamination, and systematic bias, because errors introduced here propagate silently through alignment and variant calling.

• 2×150bp: Read length (PE) (NovaSeq X Plus standard)
• ~90–120 GB: Raw data per genome (30× WGS, FASTQ.gz)
• >90%: Q30 target (bases with error prob <0.001)
• ~0.1%: Raw base error rate (Illumina sequencing-by-synthesis)

Phred quality scores and the FASTQ format

Each base call carries a Phred quality score Q = -10·log10(P_error), encoded as an ASCII character in the FASTQ fourth line:

Q10 → 1 in 10 chance of error (90% accuracy) Q20 → 1 in 100 chance of error (99% accuracy) Q30 → 1 in 1,000 chance of error (99.9% accuracy) — the de facto clinical threshold Q40 → 1 in 10,000 chance of error

A FASTQ record has four lines: read identifier (@instrument:run:flowcell:lane:tile:x:y), the base sequence, a "+" separator, and the quality string (same length as sequence). A typical 30× human WGS run produces ~800–900 million paired-end reads, generating 90–120 GB of compressed FASTQ per sample.

Quality decays toward the 3' end of each read as reagents degrade and phasing/pre-phasing error accumulates cycle over cycle — this is why adapter and quality trimming (fastp, Trimmomatic) removes the last few low-confidence bases before alignment.

Clinical labs typically require ≥85–90% of bases at Q30 or higher genome-wide, and >95% for genes on the ACMG secondary findings list, before releasing a run for downstream analysis.

FastQC / fastp metrics and common failure modes

Automated QC tools compute a standard panel of diagnostics per FASTQ file:

• Per-base sequence quality — box plot of Q score by read cycle; should stay above Q28 through 150bp • Per-sequence GC content — should match expected ~41% GC for the human genome; skew suggests contamination • Adapter content — Illumina TruSeq adapter read-through in short-insert fragments; trimmed by cutadapt/fastp • Overrepresented sequences — PCR duplicates or adapter dimers inflating specific k-mers • Sequence duplication levels — library complexity proxy; low complexity libraries waste sequencing depth on redundant reads • N content — fraction of no-call bases, should be near zero after cycle ~5

Failing QC triggers either re-trimming with stricter parameters, re-library-prep, or in clinical settings, re-draw and re-sequencing of the specimen before any variant is trusted.

BWA-MEM and minimap2 — Mapping Short Reads onto the Reference Genome

Alignment places each of the ~800 million reads onto its most likely genomic origin along the GRCh38 (or T2T-CHM13) reference, producing SAM/BAM records with a CIGAR string describing matches, mismatches, insertions, deletions, and soft-clips. This is the computational bottleneck of the pipeline and the foundation every downstream variant call depends on — a misplaced read is a variant call waiting to happen.

• >99%: Mapping rate (BWA-MEM) (reads mapped to GRCh38)
• >95%: Properly paired (concordant insert size)
• 5–15%: Duplicate rate (PCR + optical duplicates)
• ~3–5 h: Alignment runtime (30×) (BWA-MEM, 16-core node)

Seed-and-extend alignment and the CIGAR string

BWA-MEM (Li 2013) uses an FM-index built on the Burrows-Wheeler Transform of the reference to find maximal exact matches (MEMs) as seeds, then extends each seed with Smith-Waterman-Gotoh affine-gap local alignment to produce the final alignment and mapping quality (MAPQ). minimap2 uses a minimizer-based index instead, favored for long-read platforms (PacBio HiFi, Oxford Nanopore) and increasingly for short reads at scale due to speed.

Each alignment record encodes a CIGAR (Compact Idiosyncratic Gapped Alignment Report) string, e.g. "76M1I73M":

M — alignment match (may be a sequence match OR mismatch, base-level comparison needed) I — insertion relative to reference D — deletion relative to reference S — soft-clip (unaligned read portion, kept in record) H — hard-clip (unaligned portion, removed from record)

SAM/BAM records also carry MAPQ (0–60, Phred-scaled probability the mapping position is wrong), FLAG bits (paired, reverse strand, duplicate, supplementary), and the mate's position for insert-size calculations.

Duplicate marking and BQSR

PCR amplification during library prep and optical duplicates from adjacent flow-cell clusters create multiple identical read pairs that would otherwise inflate confidence at variant sites. Picard MarkDuplicates (or GATK MarkDuplicatesSpark) identifies read pairs sharing identical 5' mapping coordinates and marks all but one copy — typical duplicate rates run 5–15% for PCR-based libraries and <2% for PCR-free WGS preps.

Base Quality Score Recalibration (BQSR) then models systematic, non-random errors in the reported Phred scores as a function of machine cycle, dinucleotide context, and reported quality bin, using known-variant sites (dbSNP) to distinguish real polymorphism from sequencing artifact — producing empirically corrected quality scores that materially improve downstream genotype likelihoods.

PCR-free library preparation is now standard for clinical WGS specifically because it removes GC-bias and duplicate-driven coverage gaps in GC-rich promoter and exon regions — regions that are disproportionately clinically relevant.

Local Reassembly and Genotyping — Calling SNPs and Indels from Pileups

Variant calling asks, at every position in the genome, whether the observed reads are better explained by the reference allele, an alternate allele, or both (heterozygous). GATK HaplotypeCaller solves this with local de novo assembly of active regions into a graph of candidate haplotypes; DeepVariant instead renders the pileup as a multi-channel image and classifies it with a convolutional neural network trained on GIAB truth sets.

• ~4.0–4.5M: SNPs called (30× WGS) (vs. GRCh38 reference)
• ~500–700k: Indels called (1–50bp insertions/deletions)
• 99.95%+: DeepVariant F1 (GIAB) (SNP accuracy, HG002 truth set)
• ≥20×: Recommended min. depth (for reliable heterozygote calls)

Active region detection and local De Bruijn graph assembly

HaplotypeCaller first scans the genome in windows, flagging "active regions" where pileup evidence (mismatches, soft-clips, indel signatures) departs from the reference above a noise threshold. Within each active region (typically 100–300bp), all overlapping reads are assembled de novo into a De Bruijn graph built from overlapping k-mers (default k=10 and k=25, multiple k-mer sizes to resolve repeats).

Paths through the graph that reconnect to the reference on both ends become candidate haplotypes. Each read is then re-aligned to every candidate haplotype using the PairHMM algorithm, producing per-read, per-haplotype likelihoods that feed a Bayesian genotyping model — this local reassembly step is precisely what allows accurate calling of indels sitting inside repetitive or low-complexity sequence, which naive pileup-counting methods get wrong.

Genotype likelihoods are computed for each candidate genotype (0/0, 0/1, 1/1, 1/2 …) using the standard diploid model, and the most probable genotype is emitted with a Phred-scaled quality (QUAL, GQ).

DeepVariant — pileups as images, genotyping as computer vision

DeepVariant (Poplin et al., Nature Biotechnology 2018) reframes variant calling as an image classification problem. For each candidate site, it renders a multi-channel "pileup image": rows are individual reads, columns are reference positions, and channel values encode base identity, base quality, mapping quality, and strand. An Inception-v3-derived CNN, trained on Genome in a Bottle (GIAB) truth sets (HG001–HG007), classifies each image into one of three genotypes.

Because it learns error patterns empirically rather than modeling them explicitly, DeepVariant is notably robust across sequencing platforms (Illumina, PacBio HiFi, Oxford Nanopore, Element) simply by retraining on platform-specific truth data, and has topped PrecisionFDA Truth Challenge leaderboards for SNP and indel accuracy since 2019.

At 30× coverage, a typical human genome yields ~4.8 million total variant calls relative to GRCh38 — roughly one variant every 620 bases — of which the overwhelming majority are common, benign population polymorphisms rather than disease-relevant findings.

Separating Signal from Noise — Variant Filtering and Functional Annotation

A raw call set of ~4.8 million variants contains real biology mixed with sequencing and alignment artifacts. Variant Quality Score Recalibration (VQSR) learns a statistical model of what a true variant looks like across several quality dimensions, while VEP or ANNOVAR annotate each surviving variant with its functional consequence, gene context, and population frequency in reference databases like gnomAD.

• 99.5% / 99.9%: VQSR sensitivity tranche (SNP / indel truth sensitivity)
• ~95–98%: Typical PASS rate (of raw HaplotypeCaller calls)
• 807,162: gnomAD v4 sample size (exomes + genomes, cross-ancestry)
• <0.1%: Rare disease AF cutoff (typical gnomAD max population AF)

VQSR — a Gaussian mixture model over annotation space

VQSR trains two Gaussian mixture models (positive: known true-positive sites from HapMap/1000 Genomes/Omni; negative: sites enriched for artifacts) over a multi-dimensional annotation space including:

QD (QualByDepth) — variant confidence normalized by depth, catches artifacts inflated by deep coverage FS (FisherStrand) — strand bias via Fisher's exact test; real variants shouldn't cluster on one strand MQ (RMSMappingQuality) — root-mean-square mapping quality of supporting reads MQRankSum, ReadPosRankSum — rank-sum tests comparing ref vs. alt supporting reads SOR (StrandOddsRatio) — an alternative strand-bias metric more robust at high depth

Every variant receives a VQSLOD score (log-odds of being true vs. false under the two models), and variants are assigned to sensitivity tranches — eed 99.5% or 99.9% tranches mean the filter retains that fraction of known true variants from the training set, trading off some false-positive tolerance for near-complete sensitivity. Cohort-scale VQSR needs many samples (30+) to fit stable models; single-sample or small cohorts typically fall back to hard filters instead (QD<2.0, FS>60.0, MQ<40.0, ReadPosRankSum<-8.0 for SNPs).

Functional annotation with VEP and population frequency filtering with gnomAD

Ensembl's Variant Effect Predictor (VEP) and ANNOVAR both map each variant onto overlapping transcripts and predict its molecular consequence using Sequence Ontology terms:

missense_variant, synonymous_variant, stop_gained, frameshift_variant, splice_donor_variant, splice_acceptor_variant, 5_prime_UTR_variant, intron_variant, intergenic_variant …

Each annotated variant is cross-referenced against gnomAD (Genome Aggregation Database, v4: 807,162 exomes and genomes spanning diverse ancestries), which reports allele frequency per population (African/African-American, Latino/Admixed American, East Asian, European non-Finnish, European Finnish, South Asian, Ashkenazi Jewish, Middle Eastern). For rare disease diagnostics, variants with a maximum population AF above ~0.1–1% are filtered out under the assumption that a severe monogenic disease allele cannot be common in the general population — the exact cutoff depends on disease prevalence and inheritance model (dominant vs. recessive).

Additional in-silico predictors (CADD, REVEL, SpliceAI, AlphaMissense) are layered on to estimate deleteriousness of missense and splice-region variants before clinical review.

gnomAD population frequency filtering is the single most powerful noise-reduction step in the entire pipeline: it typically eliminates >99% of the ~20,000 coding variants in an individual genome, leaving only a few dozen to a few hundred rare, potentially disease-relevant candidates for manual review.

From Filtered VCF to Clinical Report — ACMG/AMP Variant Classification

The final stage translates a short list of rare, functionally plausible variants into a clinical verdict. The ACMG/AMP 2015 guidelines (Richards et al.) define 28 weighted evidence criteria — spanning population data, computational predictions, functional studies, segregation, and de novo status — that combine into a five-tier classification determining what, if anything, gets reported back to the ordering physician and patient.

• 28: ACMG/AMP criteria (PVS1, PS1–4, PM1–6, PP1–5, BA1, BS1–4, BP1–7)
• 5: Classification tiers (Pathogenic → Benign)
• 81 genes: ACMG SF v3.2 gene list (medically actionable secondary findings)
• ~30–40%: Diagnostic yield, rare disease WGS (trio-based clinical exome/genome)

The five-tier classification system

Every reviewed variant is sorted into one of five classes, each carrying different clinical action:

Pathogenic (P) — strong, multi-line evidence of disease causation; reported and acted on clinically Likely Pathogenic (LP) — substantial evidence, >90% confidence of pathogenicity; typically reported Variant of Uncertain Significance (VUS) — insufficient or conflicting evidence either way; NOT used for clinical decision-making, though often re-reviewed as evidence accumulates Likely Benign (LB) — evidence favors no clinical significance Benign (B) — common population variant or definitively demonstrated to be non-pathogenic (e.g., AF too high under BA1: >5% in any gnomAD population)

Evidence codes are combined with defined point-scoring rules: for example, PVS1 (null variant in a gene where loss-of-function is a known disease mechanism) is very strong evidence, while PP3 (multiple computational predictors agree on deleterious effect) is only supporting-level evidence. A variant needs a specific combination of strong + supporting evidence to cross into Pathogenic — a single weak criterion is never sufficient.

Secondary findings and the reportable range

Beyond the primary indication that prompted testing, ACMG maintains a curated Secondary Findings list (v3.2: 81 genes) — medically actionable conditions (hereditary cancer syndromes like BRCA1/2 and Lynch syndrome genes, cardiac conditions like hypertrophic cardiomyopathy and Long QT genes, familial hypercholesterolemia) that labs may optionally report back regardless of the original testing indication, if the patient has consented.

The final clinical VCF and accompanying report include: variant coordinates (genome build, chromosome, position, ref/alt), HGVS nomenclature (coding and protein), zygosity, ACMG classification with supporting evidence codes, and inheritance pattern — reviewed by a board-certified clinical molecular geneticist or genetic counselor before release, and typically Sanger-confirmed for any variant that will change patient management.

In clinically indicated rare disease cohorts, trio (proband + both parents) whole genome sequencing achieves a diagnostic yield of roughly 30–40% — meaningfully higher than exome sequencing alone, because WGS additionally captures structural variants, deep intronic splice variants, and regulatory regions inaccessible to exome capture.