What WES Data Is and How to Work With It

Whole exome sequencing (WES) is short-read sequencing restricted to the protein-coding regions of the genome, about 30-35 Mb of target across roughly 20,000 genes, or 1-2% of the 3.1 Gb human genome. A capture kit of biotinylated probes pulls down fragments overlapping those regions before sequencing, so you spend your reads where the probes landed. What you get back is the same format as any short-read experiment: FASTQ, then an aligned CRAM or BAM, then a VCF or gVCF of variants, plus a BED file describing the capture target. The thing that makes WES data distinct from WGS data is that coverage is deliberately non-uniform. The file types are not what sets it apart. Every analysis you run has to account for where the probes did and did not work.

What is in a WES delivery, file by file

A complete WES delivery should contain:

• FASTQ (_R1.fastq.gz, _R2.fastq.gz), typically 2x100 or 2x150 paired-end. At 100x mean target depth on a 35 Mb target with roughly 50-70% on-target rate, expect 6-12 GB compressed per sample.
• BAM or CRAM aligned to a named reference. Insist on knowing which: GRCh38 with ALT contigs, GRCh38 no-ALT, or hg19. Mixing coordinate systems silently is the most common way to waste a week.
• Index files: .bai/.crai, .tbi for the VCF.
• VCF or gVCF. A gVCF is more useful because it records reference confidence blocks, so you can tell “this position is confidently reference” from “this position had no reads.”
• The capture target BED, from the kit vendor. Examples are IDT xGen Exome Research Panel v2, Twist Comprehensive Exome and Agilent SureSelect Human All Exon V8. Vendors publish both a “target” BED (what they intended to capture) and a “probe/bait” BED (where the probes sit). You want both.
• QC metrics, ideally Picard CollectHsMetrics output.

If the lab will not tell you the kit and the reference build, treat the data as unusable until they do.

The QC numbers that matter

Run picard CollectHsMetrics with BAIT_INTERVALS and TARGET_INTERVALS pointing at the kit’s interval lists, and read these fields:

• PCT_TARGET_BASES_20X: fraction of the target covered at 20x or more. Below 90% on a clinical-grade exome, ask why. Good libraries hit 95%+.
• MEAN_TARGET_COVERAGE: 100x is the usual clinical spec. Research exomes often ship at 50x, which materially degrades heterozygous detection at the coverage tail.
• PCT_SELECTED_BASES: on-bait fraction. 60-80% is normal. Below 50% means you paid for whole-genome reads at exome prices.
• FOLD_80_BASE_PENALTY: how much you would have to over-sequence to bring 80% of targets to the mean. Under 2.0 is good. Above 3 means the capture was uneven and your low-coverage tail is large.
• Duplicate rate from MarkDuplicates. Above 20% on a PCR-based library means low input DNA complexity, and your effective depth is lower than the raw number suggests.

Also compute per-sample contamination with VerifyBamID2 and check the sex chromosomes and a quick relatedness matrix (somalier is fast and works on exomes) before doing anything else. Sample swaps happen.

A pipeline we would run

Nothing exotic. The point is to make choices you can defend later.

# Align. -Y for soft-clipped supplementary alignments; -K for reproducible output
bwa-mem2 mem -t 32 -K 100000000 -Y -R "@RG\tID:S1\tSM:S1\tPL:ILLUMINA\tLB:lib1" \
  GRCh38_full_analysis_set_plus_decoy_hla.fa S1_R1.fq.gz S1_R2.fq.gz \
  | samtools sort -@8 -o S1.bam -

gatk MarkDuplicates -I S1.bam -O S1.md.bam -M S1.dup.txt
gatk BaseRecalibrator -I S1.md.bam -R ref.fa --known-sites dbsnp.vcf.gz -O S1.bqsr
gatk ApplyBQSR -I S1.md.bam -bqsr S1.bqsr -O S1.final.bam

# Call only inside the target, padded 100 bp to catch splice regions
gatk HaplotypeCaller -R ref.fa -I S1.final.bam -L target.bed --interval-padding 100 \
  -ERC GVCF -O S1.g.vcf.gz

For a single sample, hard filters beat VQSR (which needs many samples to train). Reasonable starting points for SNVs: QD < 2.0, FS > 60.0, MQ < 40.0, SOR > 3.0 and ReadPosRankSum < -8.0. For indels: QD < 2.0, FS > 200.0, ReadPosRankSum < -20.0. Then normalize before annotating, because otherwise your indel representations will not match any database:

bcftools norm -m -any -f ref.fa -Oz -o S1.norm.vcf.gz S1.vcf.gz

Annotate with VEP (--everything --assembly GRCh38 --pick_allele_gene) or Ensembl’s offline cache plus gnomAD frequencies. DeepVariant is a reasonable alternative caller and tends to do better on indels in homopolymers. Use the WES model and do not use the WGS one. The models are trained on different coverage profiles and swapping them costs you accuracy.

Where WES quietly fails

The exome’s limits are structural. More reads will not fix them.

Regions the probes miss. GC-rich first exons, segmental duplications, and paralogous genes (PMS2, SMN1/SMN2, CYP2D6, HBA1/HBA2) capture badly or map ambiguously. A gene can be “on the panel” and still have exons at 5x.

Copy number. Exome CNV calling from depth ratios works, with tools such as ExomeDepth, CNVkit and GATK gCNV. It still needs a panel of normals processed with the identical kit and chemistry, and it resolves single-exon deletions poorly. Expect noisy calls and confirm anything important orthogonally.

Non-coding variation. Deep intronic splice variants, promoter and UTR variants, and regulatory elements are outside the target by construction. Structural variants with breakpoints in introns are invisible.

Repeat expansions and mitochondrial DNA. Short reads plus capture are the wrong instrument for STR expansions. Mitochondrial coverage from an exome is incidental and heteroplasmy estimates from it are unreliable.

This is the core tradeoff against WGS. WGS gives near-uniform coverage, better CNV and SV sensitivity, and the non-coding genome, at lower per-base depth for the same cost. Professional guidance has, for years, treated genome-scale sequencing as a clinical tool whose value depends on how the interpretation and reporting are handled rather than on the raw coverage number.¹ If your question is “what protein-coding variation do I carry,” WES answers it at high depth for less money. If your question is open-ended, WGS is the better instrument and the analysis burden is similar.

Interpretation is the hard part

Calling variants is a solved engineering problem. Deciding what any of them mean is not. A typical exome yields 20,000-25,000 coding variants per person, of which a few hundred are rare and predicted deleterious, and almost all of those are benign. Published best practices for clinical genome interpretation are explicit about the machinery required: a defined analytic validity threshold, ACMG/AMP-style evidence classification, phenotype-driven filtering, and confirmatory testing for anything reported.² Applying those rules to your own data usually turns a long candidate list into zero actionable findings, which is the correct outcome.

Two things follow. First, anything you might act on needs a clinician and an orthogonal confirmation in a CLIA-certified lab. A research-grade VCF is not a diagnosis, and we do not treat it as one. Second, sequencing surfaces findings you did not ask for, including carrier status and risk alleles unrelated to your original question. Sequencing also has implications for relatives who never consented to anything. That is a decision to make before you sequence. It should not wait until after you open the VCF.³

Questions people also ask

What is WES testing? A lab captures the protein-coding regions of your genome with hybridization probes, sequences them at high depth (typically 100x mean on target), aligns the reads to a reference, and calls variants in those regions. Clinically, it is used to look for coding variants that explain a specific phenotype. The output is a VCF plus a curated report.

How much data is a WES sample? Roughly 6-12 GB of compressed FASTQ per sample at 100x, 3-6 GB as CRAM, and a VCF in the tens of megabytes. Small enough to keep on a laptop, which is one of the real advantages over WGS.

Can I analyze WES data myself? Yes. A modern laptop can annotate and filter a VCF in minutes. Re-aligning from FASTQ wants 32 GB of RAM and a few hours of CPU, or a spot instance for a couple of dollars. The tooling is open source and well documented. The tools are bwa-mem2, GATK, DeepVariant, bcftools and VEP.

Does WES replace a genome? No. It covers 1-2% of the sequence and misses non-coding variation, most structural variation, repeat expansions, and reliable mitochondrial analysis. It is a deep look at a narrow slice.

Footnotes

1. Carla van El, Martina C. Cornel, Pascal Borry, et al. Whole-genome sequencing in health care. European Journal of Human Genetics, 2013. https://doi.org/10.1038/ejhg.2013.46 ↩

2. Christina Austin‐Tse, Vaidehi Jobanputra, Denise Perry, et al. Best practices for the interpretation and reporting of clinical whole genome sequencing. npj Genomic Medicine, 2022. https://doi.org/10.1038/s41525-022-00295-z ↩

3. Wim Pinxten, Heidi Howard. Ethical issues raised by whole genome sequencing. Best Practice & Research Clinical Gastroenterology, 2014. https://doi.org/10.1016/j.bpg.2014.02.004 ↩