Whole Exome Sequencing (WES): Variant Calling & VAF

Example environment: Ubuntu 24.04, Bash 5.1 (BWA-MEM 0.7.17, GATK 4.4.0)

Follow the data contract from paired reads to a filtered single-sample VCF, then calculate and inspect read-based allele fractions. The commands are a workflow template: paths, indexes, target intervals, software versions, and filtering decisions must match your experiment.

Format: Advanced template. For learners comfortable with Linux and sequencing files, this walkthrough traces the decisions and checks in a single-sample WES workflow. It supplies command templates, not sequencing data or reference resources; readers must obtain a mutually compatible FASTQ, reference, and capture-target set.

Learning objectives and prerequisites

By the end of this lesson you will be able to:

Trace an end-to-end germline Whole Exome Sequencing (WES) workflow from paired FASTQ through a filtered single-sample VCF, identifying the inputs and decisions needed to adapt the template.
Compute Variant Allele Frequency (VAF) from a VCF and interpret VAF distributions for germline WES data.
Inspect FASTQ quality, alignment summaries, duplicate metrics, and target-coverage metrics, and recognize when the tutorial lacks enough evidence to validate a callset.

Prerequisites: Linux command-line experience and familiarity with FASTQ, BAM, and VCF. This is an advanced template, not a bundled-data lab: the example names assume paired-end reads from one sample and a GRCh38 reference; the files are not included. The commands below require resources built for the exact FASTA and capture design, plus adaptation to installed software versions.

BWA-MEM (0.7.17), samtools (1.16+), FastQC, GATK (4.4+), Picard (via GATK or standalone), bcftools (1.15+), bedtools, liftover utilities (if needed), and R or Python for plotting.
Reference FASTA (GRCh38 primary assembly or GRCh37 as appropriate) and matching index files, and exome interval BED (e.g., Agilent SureSelect V6 intervals).

Why this matters: biological and computational context

WES enriches selected genomic intervals, usually including many protein-coding exons. Coverage is uneven: capture design, GC content, alignment, and library preparation affect which bases are callable. HaplotypeCaller performs local haplotype assembly in active regions; marking duplicates and, where appropriate, BQSR are distinct preparation steps. None of these steps removes the need to inspect coverage and call evidence.

Variant Allele Frequency (VAF) is the fraction of sequencing reads supporting the alternate allele at a locus. For a diploid germline heterozygous single-nucleotide variant (SNV) with no bias, VAF ≈ 0.5; homozygous alternate VAF ≈ 1.0; low VAFs may indicate mosaicism, contamination, subclonal variation, or sequencing error. Interpreting VAF requires considering depth, strand balance, read position, and mapping quality.

WES interrogates a capture-design-defined subset of the genome, not every coding base: some targets may have insufficient coverage for reliable calls, and most noncoding sequence is not systematically assessed. Detection of structural or other complex variants also depends on the assay and pipeline. Interpret a negative result within the validated target and callable regions, not as evidence that all relevant variation is absent.

Required input files and directory layout

Before beginning, create a working directory and ensure you have these files (examples paths below):

Raw paired FASTQ: raw/NA12878_R1.fastq.gz and raw/NA12878_R2.fastq.gz
Reference FASTA: refs/GRCh38.fa (with fasta index .fai and dictionary .dict); the BWA FM-index files (.amb, .ann, .bwt, .pac, .sa) are checked and prepared in the alignment step. The .fai and dictionary are not BWA indexes.
Exome interval file: refs/SureSelect_V6_intervals.bed (BED coordinate system matching the reference)
Optional: known-sites for BQSR (e.g., dbSNP, Mills indels) in GRCh38 coordinates

To source public teaching data, begin at the IGSR NA12878 sample page, which links data collections, and the NIST Genome in a Bottle resources. Use the IGSR portal’s technology and data-type filters to investigate candidate files. Select paired-end exome FASTQs (not WGS reads), record the run/sample and reference build, and obtain the matching capture BED from the assay provider or study. These pages are discovery points, not a preassembled compatible bundle: verify file manifests/checksums, read pairing, FASTQ sample identity, build, and BED coordinate convention before using any files here. If a suitable exome FASTQ pair or its matching capture BED cannot be verified, stop rather than substituting another assay or target file. Do not substitute a GIAB truth-region BED for the experiment's capture-target BED.

Reproducible workflow (sequential commands)

From the project root (the directory containing raw/ and refs/), make the destinations and fail early if the required input paths are absent:

mkdir -p raw refs aligned metrics vcf diagnostics
for f in raw/NA12878_R1.fastq.gz raw/NA12878_R2.fastq.gz refs/GRCh38.fa refs/GRCh38.fa.fai refs/GRCh38.dict refs/SureSelect_V6_intervals.bed; do
  test -s "$f" || { echo "Missing or empty required input: $f" >&2; exit 1; }
done

Before analysis, record the versions actually installed (the environment label above is illustrative, not a tested software lock):

date -Is
bwa 2>&1 | head -n 3
samtools --version | head -n 1
gatk --version
bcftools --version | head -n 1
R --version | head -n 1

Each command below states inputs and outputs. Adjust thread counts (-t, -@) for your environment. These commands are written to be safe (no destructive operations) and to document file provenance.

1. Alignment and on-the-fly sorting

# Input: raw/NA12878_R1.fastq.gz, raw/NA12878_R2.fastq.gz, refs/GRCh38.fa
# Output: aligned/NA12878.sorted.bam (coordinate-sorted BAM)
set -euo pipefail
# BWA-MEM needs its own FM-index for this exact FASTA; .fai/.dict are not substitutes.
# If refs/GRCh38.fa has changed since indexing, remove/rebuild its BWA index before running.
# GRCh38 is human-genome-sized; use BWA's bwtsw index builder rather than the <=2 GB IS algorithm.
bwa_index_missing=false
for suffix in amb ann bwt pac sa; do
  if [[ ! -s "refs/GRCh38.fa.${suffix}" ]]; then
    bwa_index_missing=true
  fi
done
if [[ "$bwa_index_missing" == true ]]; then
  echo "Preparing BWA FM-index from refs/GRCh38.fa (rebuild if this FASTA changes)." >&2
  bwa index -a bwtsw refs/GRCh38.fa
fi
for suffix in amb ann bwt pac sa; do
  test -s "refs/GRCh38.fa.${suffix}" || { echo "Missing BWA FM-index refs/GRCh38.fa.${suffix}; alignment stopped." >&2; exit 1; }
done
if ! bwa mem -t 16 -R '@RG\tID:NA12878\tSM:NA12878\tPL:ILLUMINA' refs/GRCh38.fa \
  raw/NA12878_R1.fastq.gz raw/NA12878_R2.fastq.gz \
  | samtools sort -@ 8 -o aligned/NA12878.sorted.bam -T aligned/tmp.NA12878; then
  echo "Alignment pipeline failed (BWA or samtools sort); stopping before BAM indexing." >&2
  exit 1
fi
# Index the sorted BAM
samtools index aligned/NA12878.sorted.bam

Notes: Run this alignment block as a Bash script: set -euo pipefail plus the explicit pipeline check propagates a BWA failure through samtools and stops before indexing. If the reference FASTA changes, rebuild the BWA FM-index from that exact FASTA; an existing .fai or dictionary does not replace the five BWA index files. Include a proper Read Group (-R) string: ID (read group ID), SM (sample name), PL (platform). GATK requires SM to be present in the BAM header.

2. Mark duplicates

# Input: aligned/NA12878.sorted.bam
# Output: aligned/NA12878.dedup.bam, metrics/NA12878_dup_metrics.txt
gatk MarkDuplicates \
  -I aligned/NA12878.sorted.bam \
  -O aligned/NA12878.dedup.bam \
  -M metrics/NA12878_dup_metrics.txt \
  --CREATE_INDEX true

What this produces: a coordinate-sorted, duplicate-marked BAM with an accompanying metrics file. The BAM will have the duplicate flag set (0x400) for reads inferred to be PCR/optical duplicates.

3. Base Quality Score Recalibration (BQSR): optional but recommended

# Input: aligned/NA12878.dedup.bam, known-sites VCFs (e.g., dbsnp, Mills indels)
# Output: aligned/NA12878.recal.table, aligned/NA12878.recal.bam
gatk BaseRecalibrator \
  -R refs/GRCh38.fa \
  -I aligned/NA12878.dedup.bam \
  --known-sites refs/dbsnp_146_GRCh38.vcf.gz \
  --known-sites refs/Mills_and_1000G_gold_standard.indels_GRCh38.vcf.gz \
  -O aligned/NA12878.recal.table \
  -L refs/SureSelect_V6_intervals.bed

gatk ApplyBQSR \
  -R refs/GRCh38.fa \
  -I aligned/NA12878.dedup.bam \
  --bqsr-recal-file aligned/NA12878.recal.table \
  -O aligned/NA12878.recal.bam \
  --create-output-bam-index true

Notes: Known-sites must match the reference build. BQSR models systematic errors for the sequencer and adjusts base quality scores. BQSR requires known-sites resources matching the reference build and is not automatically beneficial for every data set. Compare the workflow with current project guidance; omitting it means use the duplicate-marked BAM consistently downstream.

4. Call variants with GATK HaplotypeCaller (per-sample GVCF or direct VCF)

# Input: aligned/NA12878.recal.bam (or aligned/NA12878.dedup.bam if skipping BQSR)
# Output (here): vcf/NA12878_raw.vcf.gz
gatk HaplotypeCaller \
  -R refs/GRCh38.fa \
  -I aligned/NA12878.recal.bam \
  -L refs/SureSelect_V6_intervals.bed \
  --native-pair-hmm-threads 8 \
  -O vcf/NA12878_raw.vcf.gz

This command produces a single-sample VCF. For cohort joint genotyping, emit one intermediate GVCF per sample with -ERC GVCF (also accepted as --emit-ref-confidence GVCF), then combine/genotype according to the GATK cohort workflow. A GVCF is not the final cohort VCF.

5. Hard-filtering (example filters; consider VQSR for larger sample sets)

# Input: vcf/NA12878_raw.vcf.gz
# Output: vcf/NA12878_filtered.vcf.gz
gatk VariantFiltration \
  -R refs/GRCh38.fa \
  -V vcf/NA12878_raw.vcf.gz \
  -O vcf/NA12878_filtered.vcf.gz \
  --filter-name "QD_lt_2" --filter-expression "QD < 2.0" \
  --filter-name "QUAL_lt_30" --filter-expression "QUAL < 30.0" \
  --filter-name "FS_gt_60" --filter-expression "FS > 60.0" \
  --filter-name "MQ_lt_40" --filter-expression "MQ < 40.0"

VariantFiltration annotates records; it does not remove failing records. A record with FILTER=. has not been evaluated by a filter and is not equivalent to PASS. Preserve the annotated file for review, then use the explicit PASS-only handoff below before describing downstream records as filtered. The expressions here are illustrative only: confirm annotations exist in the VCF header and tune against truth data or project-specific evidence. Hard filters act on sites and can discard true calls; VQSR also has data/resource requirements and is not automatically appropriate for every exome cohort.

5a. Create an explicit PASS-only VCF for downstream VAF

This handoff retains only records whose FILTER field is PASS; failing labels and unreviewed . records are excluded. The original annotated VCF remains the audit trail.

# Input: vcf/NA12878_filtered.vcf.gz (FILTER annotations retained)
# Output: vcf/NA12878_pass.vcf.gz (only records with FILTER=PASS)
set -euo pipefail
bcftools view --apply-filters PASS --output-type z \
  --output vcf/NA12878_pass.vcf.gz vcf/NA12878_filtered.vcf.gz

6. Optional: LiftOver variants between genome builds

The calling and filtering commands above use GRCh38, so their output does not need this conversion. If you instead start with a separate VCF already called against GRCh37, this optional Picard/GATK LiftoverVcf example converts that GRCh37 input to GRCh38 using an hg19-to-hg38 chain and the GRCh38 target reference. Always re-annotate and validate lifted variants.

# Separate input: this VCF was called against GRCh37; it is not produced above.
# Output: vcf/NA12878_GRCh37input_liftedToGRCh38.vcf.gz
gatk LiftoverVcf \
  -I vcf/NA12878_GRCh37_filtered.vcf.gz \
  -O vcf/NA12878_GRCh37input_liftedToGRCh38.vcf.gz \
  -CHAIN refs/hg19ToHg38.over.chain.gz \
  -R refs/GRCh38.fa \
  --REJECT vcf/NA12878_GRCh37_liftover_rejected.vcf.gz

7. Compute Variant Allele Frequency (VAF)

For a sample and a biallelic genotype, a common read-based estimate is alternate AD divided by the sum of reference and alternate AD. FORMAT/AD is allele-depth data; FORMAT/DP is a separate caller field and may not equal the sum of AD. This example reads the explicit PASS-only VCF and prints FILTER so the handoff stays inspectable. It is deliberately restricted to one sample and biallelic records: multiallelic records and missing values are reported as NA, not silently coerced into a biallelic estimate.

# Input: vcf/NA12878_pass.vcf.gz
# Output: vcf/NA12878_pass.withVAF.tsv (one sample; FILTER and FORMAT/DP shown)
set -euo pipefail
sample_count=$(bcftools query -l vcf/NA12878_pass.vcf.gz | wc -l)
if [[ "$sample_count" -ne 1 ]]; then
  echo "Expected exactly one VCF sample; select and review a sample before computing VAF." >&2
  exit 1
fi
bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\t%FILTER[\t%GT\t%AD\t%DP]\n' vcf/NA12878_pass.vcf.gz \
  | awk 'BEGIN{OFS="\t"; print "CHROM","POS","REF","ALT","FILTER","GT","REF_AD","ALT_AD","FORMAT_DP","VAF"} \
    {split($7,ad,","); if ($5!="PASS" || $4 ~ /,/ || $6=="." || $6=="./." || $6==".|." || ad[1]=="" || ad[1]=="." || ad[2]=="" || ad[2]==".") {v="NA"} else {den=ad[1]+ad[2]; v=(den>0 ? ad[2]/den : "NA")} print $1,$2,$3,$4,$5,$6,ad[1],ad[2],$8,v}' \
  > vcf/NA12878_pass.withVAF.tsv

Notes: AD is emitted as one value per allele (REF, then each ALT); this simple output intentionally marks multiallelic records NA. FORMAT/DP is displayed for inspection but is not the VAF denominator. Confirm sample count, field declarations, and AD order in the VCF header.

Synthetic contract example (toy records, not biological data): A PASS biallelic row with GT:AD:DP=0/1:18,12:33 has AD sum 30 and VAF 12/(18+12)=0.4; DP remains 33. A QD_lt_2 record and a FILTER=. record are excluded by the PASS-only step. A PASS multiallelic row with AD=10,5,3 is retained for inspection but its VAF is NA, because the displayed two-allele calculation is not valid for it. The tiny VCF contract test uses these toy rows and checks this exact output:

CHROM  POS  REF  ALT  FILTER  GT   REF_AD  ALT_AD  FORMAT_DP  VAF
chrToy 2    C    T    PASS    0/1  18      12      33         0.4
chrToy 6    A    C,G  PASS    1/2  10      5       19         NA

Try the small contract before using research data

Download the invented four-record VCF. In a new disposable project, compress it with bcftools into vcf/NA12878_filtered.vcf.gz, then run only steps 5a and 7. Those example filenames do not make this toy file NA12878 research data. Do not run alignment or GATK on the toy fixture.

Check your answer: the table should contain only the two PASS rows shown above. The failed-filter and unreviewed-dot rows must be absent; the multiallelic VAF must be NA. Adding a second sample should stop the VAF block with an explicit error, not silently report the first sample. These checks demonstrate filter and table mechanics only.

8. Quick VAF distribution plot (example in R)

# Input: vcf/NA12878_pass.withVAF.tsv
# Output: diagnostics/VAF_distribution.png
Rscript -e '
d <- read.delim("vcf/NA12878_pass.withVAF.tsv", header = TRUE)
png("diagnostics/VAF_distribution.png", width = 800, height = 600)
hist(as.numeric(d$VAF[!is.na(d$VAF)]), breaks = 100, col = "steelblue",
     main = "VAF distribution", xlab = "VAF")
dev.off()
'

For a diploid, high-quality biallelic germline call, heterozygous and homozygous-alternate genotypes often have allele fractions near 0.5 and 1, respectively. A histogram is descriptive, not a genotype caller: depth, allele balance, ploidy, mapping/base quality, contamination, and capture bias affect the values. The included plot command assumes the table has a numeric VAF column and an existing diagnostics directory.

Output inspection and interpretation

When inspecting results, focus on both file-level metrics and per-variant signals. For quick read-level QC and BAM summaries, run from the project root:

# FASTQ quality summaries; inspect the FastQC HTML reports for per-base quality/adapters.
fastqc -o metrics raw/NA12878_R1.fastq.gz raw/NA12878_R2.fastq.gz
samtools flagstat aligned/NA12878.dedup.bam > metrics/NA12878_flagstat.txt
# Supply the target intervals in the format required by your Picard/GATK version.
gatk CollectHsMetrics \
  -I aligned/NA12878.dedup.bam \
  -O metrics/NA12878_hs_metrics.txt \
  -R refs/GRCh38.fa \
  -BI refs/SureSelect_V6_intervals.bed \
  -TI refs/SureSelect_V6_intervals.bed

Success means the reports exist and parse, not that the sample passes assay QC. Check that FastQC reports correspond to both mates, flagstat has mapped reads, and the HsMetrics report contains target/coverage fields rather than an error or empty result. Investigate missing outputs, zero mapped reads, incompatible interval input, or unexpected sample/read counts before continuing. BED is not accepted as the interval-list format by every Picard/GATK release; convert it with the matching Picard interval-list tooling and use that file if required. Interpret all metrics against assay-specific validated criteria.

BAM-level QC: use samtools flagstat for alignment summaries and a capture-aware metric such as Picard CollectHsMetrics for target coverage and on-target measures. Interpret coverage against the assay design and the lab/project’s validated callable-region criteria; there is no single depth threshold that guarantees reliable calling.
Variant counts and Ti/Tv: compare counts and transition/transversion ratios only with a like-for-like call set, genome build, target definition, variant class, and filtering protocol. These are broad QC signals, not universal thresholds or diagnoses of contamination.
Allele fractions: examine genotype, depth, allele-specific counts, mapping quality, and strand/read-position evidence together. Low fractions can arise from artifacts or biological mixtures; a histogram alone cannot distinguish them.
Per-variant evidence: inspect DP, AD, MQ, FS, and ReadPosRankSum when present, along with genotype and local read evidence. A fraction based on few reads is less stable than the same fraction supported by many reads; do not treat a fixed depth cutoff as a universal confidence rule.

What a plausible result means and does NOT prove:

Finding a heterozygous variant (VAF ~0.5) in a gene associated with disease is only an initial observation; it does not prove pathogenicity. It requires clinical interpretation, inheritance analysis, orthogonal validation (Sanger) for critical cases, and annotation (population frequency, predicted impact).
Low-VAF variants may be true somatic/mosaic events or artifacts. They do not prove somatic origin without further experimental or cohort context.

Troubleshooting (symptoms, likely causes, and fixes)

Symptom: HaplotypeCaller crashes with "No samples in header"

Likely cause: BAM lacks Read Group (RG) or SM tag. Fix: re-run BWA-MEM with -R '@RG\tID:... \tSM:... \tPL:...' or use Picard AddOrReplaceReadGroups to add RGs.

Symptom: Excessive false positives clustered near exon boundaries

Likely cause: Capture-edge mapping artifacts or incomplete target BED specification. Fix: confirm target/reference build and coordinate conventions, then compare edge calls with read-level evidence and the assay’s validated interval policy. Padding changes the callable region and should be justified by the design; it is not a generic repair. HaplotypeCaller already performs local assembly, so legacy local realignment is not a default fix.

Symptom: Low on-target rate or highly uneven coverage

Likely cause: Library preparation or capture failure, incorrect BED file (wrong build or coordinate system), or mislabeled sample. Fix: verify BED build matches reference, re-run capture QC, inspect insert-size distribution, and consult wet-lab team.

Symptom: Unexpected Ti/Tv ratio (very low)

Likely cause: excessive false positives from permissive filters, contamination, or misaligned reads. Fix: first confirm the variant class, target region, sample, and filtering method are comparable to the benchmark. Review annotation distributions and inspect representative calls before changing thresholds; tightening filters blindly can remove true variants.

Limitations, version and reference-data caveats, and responsible use

Reference-build mismatch: Known-sites, interval BEDs, and chain files must match the reference FASTA. Using GRCh37 resources on GRCh38 BAMs will produce incorrect results.
Version differences: GATK, BWA, and samtools evolve. Command flags or default behaviors may change. Always test on a small subset and consult the tool's versioned documentation. Pipeline reproducibility requires recording exact versions.
Hard-filter thresholds are dataset-dependent. Variant Quality Score Recalibration (VQSR) is preferred for large cohorts; hard filters may be adequate for single-sample WES but require tuning.
Ethical/responsible use: Do not make clinical claims based solely on this pipeline output. Clinical reporting requires validated pipelines, accredited labs, confirmatory testing, and genetic counseling.
Liftover caveats: Some variants do not map cleanly between builds and will be rejected by LiftoverVcf. Re-annotation and manual review are required for lifted variants in clinically relevant loci.

Use the results as a QC record

Run BAM-level QC (samtools flagstat, Picard CollectHsMetrics) and produce a MultiQC report.
Confirm RG tags and reference-build consistency for all inputs.
Run HaplotypeCaller in GVCF mode for cohort workflows, then joint genotype with GenotypeGVCFs.
Filter variants (VQSR for cohorts; tuned hard filters for single samples) and compute VAFs.
Annotate with VEP/ANNOVAR/Ensembl tools and interpret with population databases (gnomAD) and clinical resources (ClinVar).
For clinically actionable variants, plan orthogonal validation and integrate phenotype information before reporting.