How to Analyze Your Own Exome Sequencing Data

Whole-exome sequencing (WES) captures the roughly 1–2% of the genome that codes for protein, about 30–35 Mb across ~180,000 exons, and sequences it deeply (typically 100x mean target coverage). If you have the raw data, the analysis is a well-defined pipeline. Align FASTQ to a reference with bwa-mem2, mark duplicates, and call variants with a haplotype-aware caller. Then annotate, and filter down from ~20,000–25,000 coding variants to a handful worth a second look. The hard part is knowing which of the remaining variants are real, which are capture artifacts, and which are real but meaningless. Running the tools is not the hard part.

What you get back, and what to check first

A lab should hand you paired-end FASTQ (or a CRAM). If all you get is a PDF report, ask for the underlying files. The minimum useful deliverable is FASTQ plus the capture kit’s target BED (IDT xGen, Agilent SureSelect, Twist, Illumina Nextera all have published targets and they differ substantially in what they include).

Before you call a single variant, check coverage against the BED. The number people quote is mean target coverage, which hides the distribution. What matters is the fraction of target bases at ≥20x, and how that fraction behaves in GC-rich first exons. Different capture platforms disagree meaningfully on which regions they cover well and on the variants they call, even from the same DNA¹. Run:

mosdepth --by targets.bed --thresholds 1,10,20,30 sample out.bam

and look at sample.thresholds.bed.gz. If ≥20x coverage of target is below ~90%, a large chunk of your exome is a blind spot, and “no variant found in gene X” means nothing there. Note also that exome capture is depleted at exon boundaries, so splice-site positions at ±1 to ±8 are often thinly covered even when the exon body looks fine.

The pipeline we would run

Reference: GRCh38 with the analysis set that includes ALT contigs and decoys (GCA_000001405.15_GRCh38_no_alt_analysis_set.fna if you want ALT-unaware alignment, which is simpler and what most callers assume). Avoid hg19 unless you are matching an old dataset.

bwa-mem2 mem -t 16 -R '@RG\tID:S1\tSM:S1\tPL:ILLUMINA\tLB:lib1' \
  GRCh38.fa R1.fq.gz R2.fq.gz \
| samtools sort -@4 -o S1.bam -

gatk MarkDuplicates -I S1.bam -O S1.md.bam -M dup.txt

gatk HaplotypeCaller -R GRCh38.fa -I S1.md.bam \
  -L targets.bed -ip 100 -ERC GVCF -O S1.g.vcf.gz

The -ip 100 interval padding matters. Capture probes pull in flanking sequence, and splice variants sit just outside the coding interval. Without padding you systematically miss them.

For a single sample, hard filters beat VQSR (variant quality score recalibration needs a large cohort to train on). A reasonable starting set, applied separately to SNVs and indels:

gatk VariantFiltration -V S1.vcf.gz \
  -filter "QD < 2.0" --filter-name QD2 \
  -filter "FS > 60.0" --filter-name FS60 \
  -filter "MQ < 40.0" --filter-name MQ40 \
  -filter "SOR > 3.0" --filter-name SOR3 \
  -filter "ReadPosRankSum < -8.0" --filter-name RPRS8 \
  -O S1.filt.vcf.gz

DeepVariant is the alternative and, in our experience, the better default for a single WES sample. It is a single container, needs no cohort for calibration, and handles indels in homopolymers better than HaplotypeCaller’s hard filters. Pass --model_type=WES and give it the same padded BED. The tradeoff: it is slower on CPU and gives you less interpretable per-variant annotations to reason about when something looks wrong.

Annotate with VEP or Ensembl’s offline cache:

vep -i S1.filt.vcf.gz --cache --offline --assembly GRCh38 \
  --everything --pick_allele_gene --vcf -o S1.vep.vcf

--everything pulls in gnomAD frequencies, CADD, SpliceAI (if you have the plugin data), ClinVar, and consequence terms. Use --pick_allele_gene rather than --pick, so you keep one consequence per gene rather than collapsing across genes.

Getting from 25,000 variants to something readable

A single human exome carries roughly 20,000–25,000 coding variants against the reference, of which a few hundred are rare and predicted damaging, and a handful are novel². The filtering cascade that survives contact with reality:

1. Frequency. Drop anything with gnomAD popmax AF above 0.01 for a recessive hypothesis, 0.0001 for a dominant one. Use popmax, not global AF, or you will misjudge variants common in one ancestry group and rare overall.
2. Consequence. Keep frameshift, stop-gained, canonical splice (±1,2), start-loss, and missense. Keep synonymous only if SpliceAI scores them above ~0.2.
3. Gene-level constraint. gnomAD’s pLI and LOEUF tell you whether loss of function in that gene is tolerated in the population. A stop-gained in a gene with LOEUF > 1.0 is usually noise.
4. Segregation, if you have relatives sequenced. A trio collapses the candidate list faster than any in-silico predictor, which is why trio design is standard in diagnostic work³⁴.

Expect to be left with 10–50 candidates from a singleton exome and single digits from a trio.

The false positives you will hit

Certain genes show up in nearly every exome as apparently loaded with rare damaging variants, and they are almost always artifacts: large multi-exon genes, genes with recent paralogs or pseudogenes, and regions where reads from a homologous locus mismap. Analyses of exome data found that a small set of genes accounts for a disproportionate share of spurious rare variant calls, and that filtering against an in-house set of previously observed calls removes most of them⁵. Practical version: maintain a blacklist. Start with the usual suspects below and add anything that recurs across unrelated samples you process.

• MUC family
• HLA region
• OR olfactory receptors
• PRAMEF
• NBPF
• TTN for sheer size

Second failure mode: indels in homopolymer and short tandem repeat tracts. Check the read pileup in IGV before believing any indel call. If the alternate allele fraction is 0.2 rather than near 0.5, and the supporting reads all share a strand or an end position, it is not a heterozygous germline indel.

Third: copy number. Standard SNV/indel calling is blind to exon-level deletions and duplications, which are a real fraction of coding disease variation. Tools like ExomeDepth, CNVkit, or GATK gCNV can recover some of this from read depth, but they need a panel of samples processed identically for normalization, and their resolution at single-exon scale is modest.

When the exome is the wrong instrument

Exome capture is non-uniform by construction. Coverage depends on probe design, GC content, and hybridization efficiency, so even at high mean depth you get regions with near-zero reads. A direct comparison found that whole-genome sequencing detects exonic variants that WES misses, with more uniform coverage across the coding regions themselves, and better sensitivity for indels and CNVs⁶. WGS also gives you the noncoding genome, mitochondrial variants at proper depth, and structural variants with breakpoint resolution.

The case for WES has always been economic and interpretive: fewer bases, cheaper, and the coding fraction is where the variants we know how to read live⁷⁸. That case has weakened as WGS costs fell. If you are choosing now and you want a dataset you will return to for years, we would sequence the genome. If you already have exome data, it is entirely workable, provided you keep the coverage map next to the variant list and treat uncovered regions as unknown rather than negative.

Anything that looks clinically consequential in your output needs confirmation by an accredited clinical laboratory and interpretation by a genetics professional. Research-grade calls are not diagnoses, and the counseling around an actionable result matters as much as the call itself.

Questions people also ask

What is an exome sequencing test? It sequences the protein-coding exons of roughly 20,000 genes, about 1–2% of the genome, at high depth. Clinically it is used to look for a single-gene cause of an unexplained phenotype, usually with parental samples for comparison³.

Does insurance cover exome sequencing? In the US, clinical exome sequencing is often covered when ordered by a physician for a specific diagnostic indication with prior authorization. Such indications include unexplained developmental delay, congenital anomalies and suspected Mendelian disease. Coverage for a healthy adult ordering it out of curiosity is uncommon, and the self-pay price for clinical WES is typically in the low thousands.

Exome or genome? Genome, if you are paying either way and want the data to age well. WGS covers exonic bases more evenly than capture does and finds exonic variants exomes miss⁶.

How much storage and compute? A 100x exome is roughly 6–12 GB as gzipped FASTQ, 4–8 GB as CRAM. Full alignment and calling runs in about 2–4 hours on 16 cores. The VCF is tens of megabytes. Keep the CRAM, not just the VCF, because callers improve and you will want to re-call.

Footnotes

1. Michael J Clark, Rui Chen, Hugo Y K Lam, et al. Performance comparison of exome DNA sequencing technologies. Nature Biotechnology, 2011. https://doi.org/10.1038/nbt.1975 ↩

2. Pauline C. Ng, Samuel Levy, Jiaqi Huang, et al. Genetic Variation in an Individual Human Exome. PLoS Genetics, 2008. https://doi.org/10.1371/journal.pgen.1000160 ↩

3. Eleanor G. Seaby, Reuben J. Pengelly, Sarah Ennis. Exome sequencing explained: a practical guide to its clinical application. Briefings in Functional Genomics, 2015. https://doi.org/10.1093/bfgp/elv054 ↩ ↩²

4. Gerald Goh, Murim Choi. Application of Whole Exome Sequencing to Identify Disease-Causing Variants in Inherited Human Diseases. Genomics & Informatics, 2012. https://doi.org/10.5808/gi.2012.10.4.214 ↩

5. Karin V. Fuentes Fajardo, David Adams, NISC Comparative Sequencing Program, et al. Detecting false-positive signals in exome sequencing. Human Mutation, 2012. https://doi.org/10.1002/humu.22033 ↩

6. Aziz Belkadi, Alexandre Bolze, Yuval Itan, et al. Whole-genome sequencing is more powerful than whole-exome sequencing for detecting exome variants. Proceedings of the National Academy of Sciences, 2015. https://doi.org/10.1073/pnas.1418631112 ↩ ↩²

7. Andrew B Singleton. Exome sequencing: a transformative technology. The Lancet Neurology, 2011. https://doi.org/10.1016/s1474-4422(11)70196-x ↩

8. Jacek Majewski, Jeremy Schwartzentruber, Emilie Lalonde, et al. What can exome sequencing do for you?. Journal of Medical Genetics, 2011. https://doi.org/10.1136/jmedgenet-2011-100223 ↩